A REVIEW OF PHYSICAL AND CHEMICAL PROCESSES PERTAINING TO
SOLUTE TRANSPORT
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
page
I. INTRODUCTION........................ .
II. FORMULATION OF THE TRANSPORT EQUATION............3
Non-reactive Solutes....................3
Reactive Solutes......................6
Solutions of the Transport Equation. .......... 9
Concentrations.................... 10
Boundary and Initial Conditions.............1
Analytical Solutions................. 
. 13
Non-equilibrium Conditions ...... .......... 18
Physical Non-equilibrium........................ 18
Chemical Non-equilibrium.................19
Apparent Non-equilibrium................. 22
III. DISPERSION AND DIFFUSION... ......................... 26
Molecular Diffusion................... 27
Dispersion........................................ 34
The Concept of Dispersion in Porous Media. ....... 40
Physical Non-equilibrium and Dispersion. ........ 44
Experimental Determination of Dispersion Coefficients . 50
IV. ION EXCHANGE. ........................ 62
Equilibrium Chromatography................ 62
The Combined Effect of Ion Exchange and Dispersion 65
Numerical Solutions of Transport Involving Equilibrium
Exchange................. ................. ........ 69
V. TRANSPORT IN STRUCTURED AND LAYERED SOILS....... 78
Mobile and Stagnant Regions. ............... 79
Aggregated Media..................... 82
Fractured Media........................90
Layered Media.........................94
iii
VI. PREDICTION OF SOLUTE TRANSPORT 
UNDER FIELD CONDITIONS
Deterministic Models
Stochastic Models .........
Mechanistic Models ..........
Non-mechanistic Models
Concluding Remarks
VII. TRANSPORT DURING UNSATURATED 
FLOW CONDITIONS
Steady Flow
Transient Flow
Experimental 
Considerations 
. .
VIII. MULTI-PHASE TRANSPORT .
.... ..
Transport in Vapor and Liquid 
Phase . .
Transport in Immiscible Liquid 
Phases .
IX. EPILOGUE . .. . . .... .
NOTATION . . . . . .. . . . . . . . . .
BIBLIOGRAPHY .. 0...........
S. . . 103
.. . 105
. . 105
.. . . 110
.. . 116
.. . 117
.. . 118
. . . . 120
. . . 122
. .. . . . . 126
. . .. . . . . 126
. ... . . . 129
S .. . . . . 131
. . .. . . 134
. . .. . . 140
FIRST PRINTING, APRIL 1989.
Information contained herein 
is available to all without
regard to race, color, sex, 
or national origin.
iv
a . 102
1
I. INTRODUCTION
Solute transport in soils has received 
considerable attention
during the past decades. Agriculturalists are interested 
in the
behavior and effectiveness of applied 
fertilizers and pesticides.
Leaching of salts (solutes) through the 
soil profile, which depends on
exchange, exclusion, dissolution and precipitation, 
volatilization,
chemical and biological transformations, and 
other processes, affects
both the chemical and physical condition of 
a soil. Other areas of
interest are miscible displacement, such as 
the mixing of fresh ground
water with salt sea water and the secondary recovery 
of oil.
The main reason for the increased interest in 
the transport of
solutes is the concern about contamination 
of the soil, in particular
the movement of contaminants to the ground 
water. A wide variety of
substances is involved, including agricultural chemicals, 
industrial
compounds, radio-active materials, and domestic waste.
Solute transport through porous media has 
been studied by workers
from many disciplines: chemical, civil, and petroleum 
engineering;
hydrology, hydrogeology, and geochemistry; 
and soil science and
agronomy. Because of this diversity, research areas 
include
mathematical approaches to solve flow and 
transport equations,
experimental work concerning homogeneous media 
under laboratory
conditions, and field experiments with 
inherent heterogeneity.
Laboratory work has focused on exchange chromatography, the influence
2
Of biological, hydrodynamic, and geochemical processes, and the
transport during unsaturated flow conditions.
This literature review covers, in very general terms, what has
been done to date in the area of solute transport in soils. In
particular, it examines both chemical and physical processes and
explores transport during saturated as well as. unsaturated flow
conditions.
II. FORMULATION OF THE TRANSPORT EQUATION
Non-reactive Solutes
As a starting point we will consider the transport of a chemical
species in an isotropic porous medium, which is homogeneous with
respect to the relevant transport and flow parameters. The porous
medium consists of a rigid solid phase, of which the pores are filled
with liquid. It is assumed that the species does not react or interact
with the medium and that the solute is completely miscible with the
solvent. At the macroscopic level, the equation of conservation of
mass, assuming no source/sink terms, leads to:
a- -V* JD + = -V*J (2-1)
at D V3 
-
where 0 is the volumetric water content [L3L-3I, C is the concentration
of the solute expressed in mass of solute per volume solution phase
-3
[ML ], t is time [T], JD is the diffusive-dispersive mass flux
-2-1 3 -2-1
[ML T ], J is the volumetric flux of the carrier [LL T ], J is
V s
the total solute mass flux [ML-2T
-
1] and V* denotes the divergence
-1
[L I.
The autonomous flux J is generally expressed as a Fick-type
equation:
J = -OD (2-2)
D 8x
where D is the coefficient of hydrodynamic dispersion [L2T
- 1] 
and x is
distance [L]. D represents the random effects due to molecular
3
diffusion and mechanical dispersion, which occur during transport 
or
flow. Many references exist with regard to 
this topic (15, 53).
Although the mechanisms behind these two processes 
are quite different,
they behave similarly and Fick's law for diffusion 
can be used to
describe both effects. Therefore, the coefficient 
of mechanical
2 -1
dispersion D [L2T-], and the effective coefficient 
of molecular
dis
2-1
diffusion, D [L2T-], are added to give an effective 
coefficient of
e
dispersion:
D=Ddis + D 
(2-3)
d e
D depends, among other things, on water content 
and pore-water
velocity. Because of the intimate relation 
between microscopic
variation of the water velocity and solute spreading, 
D is also called
coefficient of hydrodynamic dispersion (15). One might 
argue that this
is somewhat ambiguous because only Ddi depends on the pore-water
velocity. The use of one coefficient is especially 
convenient when
solving the transport equation analytically. 
Sometimes a third
mechanism, which contributes to dispersion, is included; 
the diffusive
transfer of solute between mobile and immobile 
regions of the liquid
phase (135).
Substituting Eq.(2-2) into Eq.(2-1) and using the relation
between the Darcy flux and the pore-water velocity 
(Jv=vO), we obtain
for one-dimensional transport:
aOC_ at X Ov Oacn -- vOC] ( 24 )
-1.
where v is the average pore-water velocity [LT ]. Assuming v and 0 and
5
to be uniform with respect to x and assuming steady-state flow,
Eq.(2-4) leads to:
ac a
2
c ac
at - D 2 v (2-5)
8x
which is a special form of the Fokker-Planck equation 
for one variable,
namely with v and D independent of x and t (157). It is a linear
second-order partial differential equation of parabolic type, which
will be referred to as the advection-dispersion equation (ADE); it is
also known as convection-dispersion or convection-diffusion equation.
It has been widely used in miscible displacement studies (e.g., 121).
For multidimensional, anisotropic flow, Eq.(2-S) needs to be rewritten
(cf. 165) as:
av.C
aC _a D aC 
1(vc
= D 1(2-6)at 
ax. 
ij ax. 
ax.
1 1
where D.. is the dispersion tensor and i,j are indices denoting
ij
directions.
For laboratory conditions the 
ADE is very useful for modeling
solute movement. The equation 
can be adjusted to include 
phenomena such
as cation-exchange (107), adsorption-desorption of pesticides (198),
anion exclusion (29), precipitation-dissolution (159), nutrient uptake
by plants (47), microbial induced transformations (41) and radioactive
decay (154).
For field conditions, Eq.(2-5) is less useful to describe
transport, because most soils are macroscopically heterogeneous 
due to
factors such as structural development, heterogeneity of the 
soil
profile, swelling and shrinking, and biological activity.
Reactive Solutes
Solutes and porous media often interact with each other. Many
soils, due to the presence of negatively charged clay minerals,
variably charged organic matter, hydrous oxides and mineral edges,
exhibit adsorption and exclusion phenomena. These generally cause the
solute and the carrier to travel at a different speed.
It is assumed that the chemical species is present only in the
liquid and the sorbed phase, and that there is only one liquid phase.
The process of adsorption/desorption can be included into Eq.(2-1) to
yield:
a .C + = -V* JD (2-7)
where q is the concentration in the sorbed phase expressed in mass of
-3
solute per volume of soil [ML
-3 ]
. Sorbed refers to adsorption and
precipitation, although the latter process will be ignored. Actually q
can be regarded as a source/sink term. It should be noted that, under
unsaturated flow conditions, the vapor phase needs to be considered if
volatile solutes are present. For one-dimensional transport, Eq. (2-7)
leads to:
aa c 1 avec
Cat + q xDe x ax(2
If an exchange isotherm is known, which relates the concentration
of the species in the adsorbed phase to the concentration of all
species in the liquid phase, the dependent variable q can be expressed
in terms of C. However, this is only justified in case of a single
valued relationship between q and C, Implying that the exchange
reaction reaches instantaneous equilibrium and that 
no hysteresis
exists. The importance of hysteresis in the q(C) relationship 
was
illustrated by van Genuchten and Cleary (198). In case of 
a binary
system, S
k
) the adsorbed mass of species k per mass of soil [MM-], can
be expressed as:
S
k 
= f(CT C
k
) (2-9)
where C
T  
is the total electrolyte concentration and C
k  
is the
-3
concentration of species k in solution [ML
-
1; both C and Ckare];bt T an k e
expressed per volume of solution phase. To express the 
adsorbed
concentration on the basis of volume of soil, the 
following
relationship is used:
q = SkPb (2-10)
-3
where pb 
is the soil 
bulk density 
[ML .
For conditions similar to those stated for the derivation of
Eq. (2-5), Eq. (2-8) can be written as:
8C Pb80S 82C
- -D 
-v 
(2-11)
at o at 2ax
8x
where C and S refer to the concentrations of solute species k in a
binary system. Applying the chain rule and assuming a constant total
electrolyte level with respect to time, i.e., Sk=f(Ck) yields:
k k'
8S _dS8C 
(-2
where dS/dC is the slope of the exchange isotherm, which is constant in
case of linear adsorption or exchange and which is usually referred to
3 -1
as the distribution coefficient K
d 
[L3M-!. This quantity indicates how
the species is distributed between the liquid and 
adsorbed phase:
S = KC (2-13)
d
Using Eq. (2-12), Eq.(2-11) can be rewritten in terms 
of one dependent
variable:
ac a
2
C ac2-14)
Ra= D 2 V-x(2
ax
where the dimensionless retardation factor R is defined as 
follows:
Pbd
R = 1 + dS (2-15)
0wdC
Depending on the assumed adsorption mechanism, different
transport equations were given by Gupta and Greenkorn (75):
1. Linear adsorption isotherm: q= k C
1
k 2
1ac a C ac
1+ o a-=D 2-V j(2-16)
at - x2 a -x
2. Freundlich 
adsorption 
isotherm: 
qk
2
Cn
,-nk
2  
2
[1 n-1a ]ac a c ac
I1+0(1+bC) 2J at = ax__ -2 
(218
In the above three examples k, n, a and b are empirical constants.
Dividing Eq. (2-14) by the retardation factor leads again
to the ADE, as for the non- reactive case:
t D 2 -(2-19)8t 82 Ox
ax
with D =D/R and v =v/R. In case of non-linear exchange 
(R=f(C.)),
1
Eq. (2-19) is a non-linear partial differential 
equation for which
analytical solutions are difficult to obtain.
A more general equation should include sink/source 
terms to
account for phenomena such as radioactive decay, precipitation,
dissolution and chemical reactions. Parker and van Genuchten 
(130)
presented the following, more general expression:
P sAPb sPb
b as + ac Da2 C C S Ssb(2-20)
e at at O - O w w
+  
0
where j and W are rate constants for first-order decay 
in liquid and
w s
-1 -3-1 -
solid phase [T
-
1, respectively, and 7w [ML-T-] and [T are rate
constants for zero-order production in liquid and solid 
phase,
respectively.
Solutions of the Transport Equation
Before presenting some analytical solutions for Eq.(2-19), 
it is
worthwhile to mention the different types of concentrations that 
can be
considered and to pay attention to the appropriate boundary 
and initial
conditions (131, 200).
10
Concentrat ions
The time average of C is defined (60) as:
At
o 2
'(x yZt ) = At C(x,y,z,t) dt (2-21)
At
t --
o 2
A spatial or volume average of c, a microscopic variable in 
this case,
at a position (xoYozo) is defined as
T cdV
C (x yo ,z t) - lim AVet(2-22)
v o o o AV->V
fdV
fAVe
where the volume AV and AV represent a "large" chunk of the 
porous
medium and the liquid phase, respectively, and dV and 6V correspond 
to
a microscopic differential volume element and the Representative
3
Elementary Volume (16), respectively. All volumes are expressed 
in L
Furthermore, x, y and z are fixed positions (the center 
of a
0 0 0
physical or material point). It should be noted that the macroscopic
concentration, C, is generally taken as the volume-averaged
concentration, C
V.
At a certain position, the flux-averaged concentration, Cf, is
given as the ratio of the transport and flow term:
C (t) = Js/Jv (2-23)
f s V
Cf(t) represents the mass 
of solute per unit volume of 
fluid passing
through a given cross section 
during an arbitrary time interval 
(102).
As was stated by Parker 
and van Genuchten (131), 
solute flux
11
distributions are in many cases of more interest than pore fluid
concentrations. The relationship between flux-averaged concentration,
Cf, and volume-averaged concentration, C
V ,  
the two types of
concentrations pertinent to many displacement experiments, can be
derived from Eq.(2-23):
aC
Cf = C D V (2-24)
f V v8x
Parker and van Genuchten (131) discussed the importance of
distinguishing between both concentration modes in conjunction with the
way an experiment is performed and the results are analyzed (cf. 34,
102). Unless stated otherwise, we will assume, as was done previously,
that all concentrations 
are volume averaged, 
i.e., C=CV
V.
Boundary and Initial Conditions
In order to solve Eq.(2-S), or for that matter (2-19), specific
boundary and initial conditions are needed. Following van Genuchten and
Alves (197), the initial condition is:
C(x,t) = f(x) x>O , t=0 (2-25)
where f(x) is an arbitrary 
function.
For the boundary condition at x=O, the first- or concentration-
type boundary condition, the Dirichlet problem, is given by:
C(x,t) = g(t) x
4
O , t>O (2-26)
and the third- or flux-type boundary condition is given by:
BC
(-D -- + vC) = vg(t) t>O (2-27)
xx8O
12
where g(t) is an arbitrary function describing the concentration of the
incoming solution. The exit boundary condition can be defined for 
a
semi-infinite or for a finite soil column. For a semi-infinite 
column
the appropriate boundary condition is:
ac (x t) = 0 
x-) P, t>O 
(2-28)
ax
In case of a finite system the following continuity condition at the
exit boundary needs to be satisfied:
(-D a- + vC) =vC t>O (2-29)
xtL ex
where C is the exit concentration, which is assumed to be equal 
to
ex
CIxtL Hence:
ac(x, t) = 0 xtL, t>0 (2-30)
ax
Analytical solutions of Eq. (2-5) for conditions (2-28) and (2-30) are
approximately equal.
For the inlet boundary, a third-type condition is generally
preferred, because the mathematical solution for this condition leads
to conservation of mass. Eq. (2-26) implies equal concentrations in the
feed solution and in the porous medium at the boundary. However, in
reality it takes some time to attain the same concentration. First, the
input solution might not be well mixed resulting in a boundary layer
outside the porous medium (24). More important, because at 
the
microscopic level the concentration varies gradually, a transition zone
should be taken into account. Obviously a volume averaged concentration
of a representative elementary volume (REV) will not immediately be
equal to the concentration of the feed solution. In other words, since
13
a (tracer) solution can only be injected at a certain rate, a
prescribed concentration will not be established instantaneously. For
displacement experiments, involving a step change in concentration
(g(t<O)=O and g(t>O)=C ), a third-type condition should therefore be
0
used:
aC
(-D -+ vC) =vC (2-31)
8x o
x 0
where C is the concentration of the feed solution. The discontinuity
0
in the concentration across the inlet boundary increases with
increasing D/v. Van Genuchten and Parker (200) showed that solutions of
the ADE, subject to a first-type boundary condition, lead to flux
average concentrations, whereas solutions for a third-type boundary
condition lead to volume average concentrations. One is also referred
to the work by Parlange and Starr (133, 134), Parlange et al. (132),
and Pandey and Gupta (128).
Analytical Solutions
In case of an infinite system (-o<x<oo), a solution of Eq.(2-5)
can be obtained by making the following transformation of coordinates:
= x-vt
(2-32)
T t
which transforms Eq. (2-5) to:
C a2C
-= D a__C(2-33)
aT1T 2
The original boundary and initial conditions
1 (2-34)
C = C x<O t<O, x- t>O
0
14
need to be transformed as well. Note that C. is the initial
1
concentration in the soil. The solution of 
Eq. (2-33) can be found in
the areas of heat flow (37) and diffusion 
(45).
Using an alternative transformation:
x-vt=(2-35)
f4 Dt
results in an ordinary differential equation:
2
dC+ 2 dC=0 
(2-36)
d 
2  d=
With the transformed initial and 
boundary concentrations
described by:
c = c. -o
1 
(2-37)
C =C 0 -00
the solution of Eq. (2-36) is:
C = (C-C.)/(C -C ) = erfc(2-38)
1 0 1
where C is the dimensionless concentration 
and erfc is the
complementary error function. For other 
boundary conditions, involving
semi-infinite or finite media, a coordinate 
transformation cannot be
employed or is inconvenient. Carslaw and Jaeger 
(37, p. 388) showed how
Laplace transforms can be used to solve 
Eq.(2-5). Van Genuchten and
Alves (197) provided a compendium of available 
analytical solutions to
Eq. (2-14), many of them obtained with the help of Laplace transforms.
The solutions for four combinations of inlet and exit conditions 
are
listed in table 1.
15
Table 1. Analytical Solutions 
of Eq.(2-14) for Various
Boundary Conditions 
with Constant R, after van
Genuchten and Wierenga 
(204)
I nlet
boundary
condition
C(0. t) = C
1)0,ac- tC
0
ax X
0
= )C
0
Analytical solution
crfc 2(L-flt' Ir? KexpAt7I erfc
~ Lu;) t' r%)ep (~ut)'j
1 1 + ? + u't ep( erfc Re t
2 DR11? D J1 2(DRt)'J
ti x FIX Wt 1Dt1
213o, s1n (LCxp12 4R L'-RJ
Fri 02+(2L +v
;; [cs(lLx) + I;s(azx]ex ty x '
C - I~_T~ -U
02 + 21 D m 21
Om cot (03 +) + - 0
p, 4!)
References: A-i: (109), A-2: 
(113), A-3- (42), A-4: (27)
By using the following dimensionless 
variables
T =vt /L
Z =x/L
P =vL/D
C=(C-C. 
)/(C 0-C.)
Eq. (2-14) can be rewritten as
2-
R a8- 1aC 
ac,
= P 2 32
(2-39)
(2-40)
(2-41)
(2-42)
(2-43)
where T is the number of pore 
volumes, L is the column length 
[L], Z is
the dimensionless distance 
and P is usually referred 
to as the Peclet
number. Analytical expressions 
for C =C the dimensionless 
exit
Z=1 ve
concentration, are provided 
in table 2 for different 
boundary
Case
Al1
A-2
A-3
A-4
(3ctI3, L
0, 2o1),)4
I \
16
Table 2. Expressions for the Relative Effluent Concentration, 
C (T), in
e
Terms of the Column Peclet Number (P) and Pore Volume
(T) for the Four Analytical Solutions Listed in table 1, after
van Genuchten and Wierenga (204)
Ca e Relative effluent concentration
A-I ce(T) erfc 
)2 (R -T) + 
exp(P)erfc 
2((R T
A-2 Ce(T) = erfc [(l (R - T) + ( exp - R T) + P erfe (R T)
P P7T 02 T
S2/msin(3,) exp 
2
- - 4{ PR 
P
A-3 c,(T) = 1 - _ P1 4m cot(3m) + -= 0
02 + P+ -
0"22,sin(I) exp2 P2
2 t in/ 12ep -- H PR PO'm cot(fOm) -t3 + -- 0
A-4 c,(T) = 1 - S_2 4R Pm cot(m) 4 - 0
M=/3 
2 - - ? I)
4
conditions. The influence of the different boundary conditions on the
concentration profiles is graphically illustrated in figure 1. At small
times, the use of the first-type inlet condition results in a
considerable error.
A difference in the solutions for the semi-infinite and the
finite cases occurs once the solute concentration starts to increase at
the exit. Just as for the inlet, the concentration at the exit is not
likely to be continuous. Unless backmixing occurs from the effluent
solution to the soil, the semi-infinite solution is to be preferred.
These solutions are useful to predict the solute concentration in
homogeneous media and to determine transport parameters 
in conjunction
with an experimentally determined concentration distribution.
It should be noted that much more work on the analytical solution
of the ADE has been published (e.g., 65, 112, 127).
17
1.0
0.8
z
Z
H
z
Z
o
LUJ
z
0
u
0.6
0.4
0.2
REDUCED DISTANCE, Z
u
z
0
I.-
z
w
z
0
C
1.0
REDUCED DISTANCE, Z
FIG. 1. Calculated concentration distribution 
for R=1 and P-values of 5
and 20, respectively. The curves 
were obtained with the
analytical solutions listed in table 
2, after van Genuchten and
Alves (197).
18
Non-equilibrium Conditions
So far, it has been assumed that instantaneous local equilibrium
exists between the solute in the liquid and the adsorbed phase, the
distribution of ions in both phases being determined by 
the exchange
isotherm. However, instantaneous equilibrium 
might not always be
achieved. Two different kinds of non-equilibrium exist: physical 
and
chemical.
Physical Non-equilibrium
Physical non-equilibrium is of particular importance for
aggregated media, where flow occurs predominantly between aggregates,
and for unsaturated media, where the liquid phase in effect 
might be
discontinuous. During miscible displacement the solute in the feed
solution cannot reach all sorption sites immediately and no
instantaneous exchange equilibrium can be achieved. This situation is
generally classified as physical non-equilibrium. 
Conceptually,
advective flow occurs only in the so called mobile 
region of the fluid
phase. The concentration in the immobile region, 
both in the liquid and
the adsorbed phase, will therefore lag behind 
the concentration in the
mobile region. A large part of the sorption 
sites are usually only
accessible via the immobile or stagnant region 
of the liquid. Transfer
of solute from the mobile to the immobile region, 
and vice versa,
occurs via a diffusion controlled 
process. For physical
non-equilIibr i um, the transport equat ion was reformulated by van
Genuchten and Wierenga (202) following Coats and Smith (43):
19
ac ac. as as. a
2
c ac
mo im mo im mo mo
S+ +p. f 
- 
(1-f) = D---- v (2-44)
mo at im at b at +Pb t mo 2 mo mo ax
where the medium is supposed to be homogeneous with respect to 8 , 8.
mo im
and v and the flow is steady. The fraction of sorption sites in
mo
direct contact with the mobile region of the liquid is given by f.
Transfer of the species between mobile (mo) and immobile (im) regions
of the liquid phase was given by:
acm 
asm
im + ( (C -C.(2-45)
im at b( at mo im
-1
where a is the mass transfer coefficient [T-1
Anion exclusion (29, 103) can be viewed as a situation of
physical non-equilibrium, where the 
exclusion volume roughly
corresponds to the immobile region. Eq. (2-44) and (2-45) can be readily
adapted (207) to describe the transport of anions which are excluded
from certain regions of the liquid phase. A further discussion on the
concept of mobile and immobile regions of the liquid phase will be
given in the section on mobile and stagnant regions in Chapter V.
Chemical Non-equilibrium
Chemical non-equilibrium occurs when the adsorption/exchange
process requires some time to be completed. No instantaneous
equilibrium between the solute concentration in the liquid and the
adsorbed phase will occur. To account for chemical non-equilibrium
exchange, a kinetic approach can be taken by combining the transport
equation with the appropriate rate equation for adsorption 
of the
species by the medium (7, 8, 109).
20
Because of the variety of adsorption sites (clay minerals,
organic matter, oxydes) models with two kinds of adsorption sites have
been introduced (168). Adsorption occurs almost instantaneously for
"type- 1" sites, whereas for "type-2" sites adsorption is
time-dependent. Using somewhat arbitrary first-order kinetics, the
general sorption rates were given by Selim et al. (168) as:
aS
1 _
Stkl
C 
- k2S 
(2-46)
8t b 21
2
at - k3C - kS
2  
(247)
t 
pb3 
42
where S and S are the concentrations of solute sorbed to sites 1 and
1 2
2, respectively. Furthermore, k1 and k
3 
are forward and k and k are
backward reaction rate coefficients [T
-1
, respectively. It should be
noted that the same assumptions and restrictions apply as were made for
the derivation of Eqs.(2-5) and (2-11). At equilibrium, the sorbed
concentrations are given by the following linear isotherms:
k
_ 1
S = -- k2 C = KIC (2-48)
S - C == KC (2-49)
2 Pb k
4  
2
The sorbed concentration for all sites, at equilibrium, is given 
by:
S = S + S = (K + K2)C = KC (2-50)
A fraction F of the total sites belongs to "type-l" (i.e.,
F=S1/S=K /K). Since "type-i" sites are always at equilibrium, 
it
follows from Eq. (2-48) that:
21
as
at = FK a(2-51)
For "type-2" sites, the sorption rate may be given by a linear,
reversible first-order rate equation (cf. 109). The sorbed
concentration for these sites follows from the rate expression given
by:
at - (K
2
2C - $2 (2-52)
where a is a first-order rate coefficient [T-], as before, which is in
this case equal to k Eq. (2-51) and (2-52) can be substituted into
4.
Eq. (2-11) resulting in the following pair of equations to be solved:
Upb 
ac Pb 
as
2  a
2
C 
ac
(1+ b) -t+ &at = D -2-v 
(2-53)
e at2 8
ax
(2 4(1-F)KC-S
2  
(2-54)
Various authors (124, 168, 196) were able to describe the breakthrough
curves fairly well with a numerical solution of Eq. (2-53) in
conjunction with Eq. (2-54).
The "one-site" kinetic non-equilibrium model follows from the
previous model. In this case, F=O and Eq.(2-53) and (2-54) transform
to:
ac _PbasD a2 ac (2-55)
+ - -= D v (-x
T eat 2
ax
as = (KC - 5) 
(2-56)
at
where S still refers to the adsorbed concentration of a particular
species.
22
Apparent Non-equilibrium
Cameron and Klute (36), in effect, combined the physical and
chemical non-equilibrium models. The transition from a chemical
reaction at one site (the microscopic viewpoint) to the average
observed for a large number of pores (the macroscopic viewpoint) is
very complicated because macroscopically uniform soils are generally
not microscopically uniform. This involves both the physical 
and
chemical processes pertaining to sorption. In the section 
on physical
non-equilibrium the limited supply of the solute, due to immobile
water, was discussed. In fact, a whole range of supply 
rates exists.
The same holds for the actual sorption process; many rate 
equations
might be needed to describe sorption depending on 
the various soil
constituents and the solute. Furthermore, a distinction 
between
physical and chemical non-equilibrium is generally 
not possible.
Therefore, Cameron and Klute (36) proposed a "black box" 
approach to
describe the sorption process. They differentiated between two types 
of
sites, those which appear to react rapidly with the solute and those
which appear to react more slowly. Sorption at the latter type is
described with a kinetic type of reaction which is used to take into
account both chemical and physical non-equilibrium.
The sorption sites are divided, as in the section on chemical
non-equilibrium, into two fractions. S1 (equilibrium) and S (kinetic).
Exchange between ions in the adsorbed and liquid phase occurs via a
linear Freundlich (S
I )  
and , a kinetic (S
2 )  
type of process,
respect ively:
23
aS e 8C 0
- K -- k - C -kS (2-7
at -1Pbat 3% -4 2b Pb
where K
1  
(=k /k
2
), k
3 
and k are the equilibrium constant and the
adsorption and desorption rate, respectively. S=S +S is the adsorbed
1 2
concentration and C is the liquid concentration, expressed in mass per
volume of solution. The transport equation for combined kinetic and
equilibrium adsorption is:
Pbas
2
- ac 2c ac
- (1 +K)- .=D -V-5(2-58)
e.at i t 2 ax
ax
as
2
- k -C-k S ( 9
at 3 4S2 (2-59)? Pb
Cameron and Klute (36) expressed these equations in dimensionless form
and obtained a solution via Laplace transforms. Application of the
model to transport of atrazine, phosphorus, and silver was successful;
a purely equilibrium or kinetic model did not fit the data accurately.
Nkedi-Kizza et al. (124) fitted two non-equilibrium models, a
diffusion controlled and *a first-order reversible kinetic model,
through experimental breakthrough curves. How fast equilibrium is
attained in the ion exchange process, merely a redistribution of ions
rather than a typical chemical reaction, is determined by two
mechanisms: the supply of solute 
through the liquid phase to 
the
liquid/solid interface and the nature of the exchange reaction at that
interface. The actual ion exchange reaction is generally not the rate
limiting factor in most instances (84). Rather, the diffusion of ions
from solution to exchange sites, and vice versa, seems to be the rate
limiting step even if no immobile water is present and we deal with
24
chemical non-equilibrium. Presumably, both chemical and physical
non-equilibrium models might be used to represent non-equilibrium
exchange.
In summary, both models can be considered to have two types of
adsorption sites. The physical non-equilibrium model has "mobile" sites
with instantaneous equilibrium, while for the "immobile" sites exchange
is diffusion controlled. This model is described with Eq. (2-44) and
(2-45). The chemical non-equilibrium model is described with Eq. (2-53)
and (2-54). For "type-2" sites, instantaneous equilibrium is not
achieved because of the kinetic nature of the exchange process. When
expressed in dimensionless variables it can be shown that the transport
models for both models are identical and have equivalent breakthrough
curves. Both models can be used when describing ion exchange during
transport through aggregated sorbing media. Based on breakthrough
curves, obtained via curve fitting, Nkedi-Kizza et al. (124), however,
did not want to draw the conclusion that the two models were
conceptually similar.
Judgment whether 
local equilibrium 
exists can be 
based on
exchange data and the physical properties of the medium (cf. 24).
Valocchi (191) presented criteria for the validity of the local
equilibrium assumption. Presently, much work is being done in the area
of transport through structured soils. As far as physical
non-equilibrium is concerned, sometimes an effective dispersion
coefficient can be used for soils containing aggregates of a particular
size and shape. In this way, the sink/source term describing mass
25
transfer from mobile to stagnant regions of the liquid phase 
can be
omitted (cf. 57, 135, 150). On the other hand,
- 
a distinction can be
made between the inter-aggregate pore space (macropores), 
containing
the mobile liquid, and the intra-aggregate pore space 
(micropores),
containing the immobile liquid. In the 
intra-aggregate porespace the
predominant transport mechanism consists of 
diffusion. Generally it is
assumed that inside the aggregate instantaneous 
equilibrium is achieved
between the species in the adsorbed and liquid phase. 
A number of
analytical solutions have been published 
(cf. 153, 186, 201). Some of
these solutions will be discussed in Chapter V.
III. DISPERSION AND DIFFUSION
When a feed solution containing a given solute displaces a
resident solution containing another solute or the same solute at a
different concentration, a transition zone will develop in which a
variation of solute concentration will occur. During displacement
experiments, the solute concentration of the effluent is generally
monitored to obtain a breakthrough curve. The curve is indicative of
the amount of mixing between the two solutions. For a non-reactive
solute, the spreading is caused by dispersion and diffusion. In order
to describe solute transport, one generally needs to quantify both
processes. Solutes can also be useful in hydrological studies, where
they are used as tracers to study flow phenomena in porous media. Many
references exist about diffusion and dispersion (15, 23, 
24, 60, 61,
165).
Following Fried and Combarnous (61), two mechanisms of dispersion
will be distinguished:
1. Mechanical dispersion. When adopting a microscopic viewpoint,
nonuniform velocity profiles exist in a porous medium because 
of the
following reasons:
- The velocity of the soil solution is zero at the solid surfaces. 
Slip
flow does not occur because of the relatively high viscosity 
of water
and the small mean free path length (165).
26
27
- The pore dimensions vary, so different maximum 
velocities occur along
the axes of the pores.
- Streamlines fluctuate with respect to the mean direction 
of flow.
2. Physico-chemical dispersion: diffusion. 
This type of dispersion is
due to concentration gradients, 
actually gradients in
chemical potential. The following phenomena 
take place:
- If inside a streamtube a concentration 
gradient occurs, then
diffusion tries to annihilate the gradient.
- If concentration gradients exist between two 
adjacent streamtubes,
mass transfer between the streamtubes occurs by 
diffusion.
For both types of dispersion a longitudinal 
and transverse
component can be distinguished. Mechanical dispersion 
(dispersion) and
physico-chemical dispersion (diffusion) are closely 
related. Diffusion
occurs at the molecular level and dispersion at the pore 
level. Usually
it is not possible to distinguish between these 
levels. Because of
their similar nature, diffusion and 
dispersion are conveniently
described together by the coefficient of hydrodynamic 
dispersion.
Molecular Diffusion
The process of molecular diffusion is of interest 
for at least
two reasons. First, at low pore-water velocities 
transport of solutes
is dominated by the diffusion process, and 
second, an analogy exists
between molecular diffusion and mechanical dispersion. Because of this
analogy, knowledge of the diffusion process is helpful in understanding
dispersion. Let us first consider molecular diffusion in a (free)
solution, next diffusion in porous media, and conclude with a
discussion of mechanical dispersion.
28
The solute molecules in the solution possess random thermal
motion which causes an exchange of molecules between adjacent volume
elements. If isothermal and isobaric conditions exist, a net. transfer
of molecules of some species k occurs when the concentration of k in
the adjacent volume elements differs. More particles of k move from
elements of higher to elements of lower concentrations than vice versa.
Such a net transfer under the influence of a concentration gradient is
termed molecular diffusion. The process of molecular diffusion is
described by Fick's first law. In the case of one-dimensional diffusion
in a free liquid, it can be expressed as:
J =-D ac(3-1)
D o ax
where J is the solute mass flux due to diffusion [ML-2T -] and D is
D 0
the coefficient 
of molecular diffusion 
[L 2T_]1
Thermodynamically, it is the gradient in the chemical potential
that is the driving force for the diffusion process. The chemical
potential is defined as:
8G(n.)' P njk(3-2)
is -1
where p
1
is the chemical potential of species k [J mol Ink is the
number of moles of species k, G is the f ree enthalpy or Gibbs free
-1
energy [J mol I and 9'J and P are temperature and pressure,
29
+R9J1in C(3-3)
0
where pk is the chemical potential of species i in a chosen standard
state and the mole fraction Ck=Ck/CM the ratio of the molar
Mk Mk M
concentration of species k to the total molar concentration. Hence, for
an ideal thermodynamic solution, the gradient of the chemical potential
is proportional to the gradient of the natural logarithm of the
concentration:
ap a ln C
~
k  
= M(-k
ax]ideal =ax(3-4)
However, most solutions do not exhibit ideal mixing behavior. In
transport studies, the non-ideality of the solute needs to be addressed
for diffusion dominated transport if, for example, the total solute
concent rat ion is "high" (sea water) or for the thermodynamic
description of ion exchange. For non-ideal systems, Eq. (3-3) is
generally rewritten as:
=0 + RYJ ln ak(3-5)
where ak is the chemical activity of component k, which is related to
the concentration by the activity coefficient:
a = kCMk(3-6)
a
k  
k3-Mk
The value of the activity coefficient, 7 k, depends on the nature of the
microscopic interactions. Activity coefficients for ions can be
estimated with the Debye-Hiickel equation for solutions with ionic
strength less than 0.1 N (e.g., 125). The 'Davies' extension of the
Debye-H~ckel equation can be used to approximate k for systems
30
containing mainly small ions at ionic strengths less than 0.5 M. For
non-ideal solutions, the gradient of the chemical 
potential is:
3
1[ ] 
8
lnak
k 
= 
k(3-7)
a x Jnon-id eal ax
The diffusion flux for ideal solutions is given 
by Fick's first
law for a constant total molar concentration CM:
N =dC 
(3-8)
k CMD dx
where N
k 
is the molar flux density of species k due 
to diffusion
[ML T
-
1. The diffusion flux for non-ideal solutions 
is obtained
by inserting the ratio of the non-ideal 
to the ideal chemical
potential gradients in Eq. (3-8):
aIn a dC
N =-CDk Mk 
(9
Nk 
M-CMDo 
* 
dx
8 ln CMk
The contribution of molecular diffusion to the 
total molar flux
can be shown by substituting Eq.(3-9) into the 
continuity equation,
Eq.(2-1). The result is:
CMk _aalIn ak dCMk JvC(Mk3-10)
at jx oaIn Ck 1 0x (3-
Assuming an ideal system with a liquid velocity 
of zero and D
0
independent of position and concentration, the well known diffusion
equation, Fick's second law, is obtained:
ac a
2
c
-D 
(3-11)
3t 0 2
where the subscript notation for Ck 
has been dropped.
31
It can be shown by differentiation that (45):
C(x,t) A xexp (3-12)
, -4D t
to
is a solution of Eq.(3-11), where A is an arbitrary constant. The total
amount of substance, m
O, 
diffusing in an infinitely long cylinder with
unit cross section, is given by:
m= C dx (3-13)
-00
With an initial condition C(x,O) = mO6(x), 8(x) being the Dirac delta
function (i.e. a spike of solute in an otherwise solute free medium),
the concentration distribution becomes (60):
m 2
C(x,t) = exp 4D (3-14)
4D t Do
Substituting o = 2D t, yields a Gaussian or normal distribution:
0
m 2
C(x,t) = 
0  
exp[ x (3-15)
S2 22
where o- is the standard deviation of the solute concentration
distribution.
The process of molecular diffusion in porous media is similar 
to
that in a free solution, except that we must now define the mean flow
path for diffusion in terms of the structure of the 
medium. Let us
assume that diffusion in other phases than the liquid phase can be
ignored. If necessary, the additional contribution to molecular
diffusion from exchangeable ions can be accounted for (126). The
coefficient of molecular diffusion in soils can be expressed as (24):
32
D = D /A (3-16)
e o
where A is the tortuosity factor. In simple terms, tortuosity can be
viewed as the ratio of the "true" (i.e., circuitous) and straight line
flow path distances between two points in a porous medium.
Theoretically, its value should only depend on the geometry of the
medium. Tortuosity can also be characterized with the so called
formation factor, which is frequently used in petroleum engineering.
This factor is the ratio of the electrical conductivity of the pure
liquid phase to that of the porous medium (cf. 46). It depends on the
volume fraction of liquid and the tortuosity of the medium:
D /D = /(Mo) (3-17)
e o
where D is considered an effective coefficient of molecular diffusion
e
2 -1
[L T 
i
] and ? is the formation factor. It should be noted that the
validity and usefulness of this concept of tortuosity has been
questioned (165).
In media with a very low hydraulic conductivity, diffusion is the
main mechanism of solute transport. Therefore, the determination of
diffusion coefficients for low-permeability media has recently received
a considerable amount of interest (e.g., 173). For a concentration
profile following a Gaussian distribution, i.e., Fickian diffusion, the
coefficient of molecular diffusion of a stagnant liquid in a porous
medium can be determined from:
-_ 2D (3-18)
dt e
where D has been substituted for D . This relationship is frequently
e o
33
used to obtain the diffusion coefficient from the standard deviation of
the concentration distribution as a function of position. In porous
media, spreading is caused by various mechanisms. Integration of
Eq. (3-18) shows that the variance increases linearly with time for
constant diffusion coefficient. Thus, Eq.(3-17) can be used to evaluate
whether these mechanisms exhibit Fickian behavior.
An approach which is not restricted to profiles following a
Gaussian distribution is the method of moments. The moments of a
concentration distribution, resulting from a pulse type of solute
input, can be used to determine the mean and variance from the observed
or calculated distribution (9, 10). No assumption about the
distribution of the solute concentration needs to be made in order to
determine these moments. The p-th moment of a concentration
distribution with respect to x, m , is defined as:
m = x C(x,t) dx (3-19)
p 
I
--00
This operator m can be applied to each term of the
advection-dispersion equation. Aris (10) showed how these moments can
be determined without evaluating the integrals, which is convenient if
the distribution of C(x,t) cannot be described by a mathematical
relationship. Some properties which can be determined with these
2 2
moments are the mean, . or i , and the variance, o or , of the
x x
concentration distribution, where
= ml/mo (3-20)
2 2
a = (m /m ) - 4 (3-21)
34
Aris' moment method is also an Important tool In obtaining solutions
for the diffusion/dispersion equation in stratified media (80, 81,
115).
Dispersion
The description of dispersion as a diffusion type process has
been shown plausible in the classical paper by Taylor (189). Taylor's
analysis, concerning flow in a circular tube, is often discussed in the
literature (e.g., 60, 61, 73,- 121). Because dispersion in a free
solution provides a qualitative explanation of dispersion in a porous
medium, a brief treatment of Taylor's paper will be presented.
Taylor- (189) considered laminar flow in a circular tube with
radius a, having the following parabolic velocity profile:
u(r) =u r ](3-22)
a-1
wh Iere u0 is the maximum velocity at the axis [LT ] and r is the radial
distance from the axis [L]. It can be seen that the mean velocity over
a cross section of the tube, u, is equal to u /2.
m 0
First, transport of a solute by advection only will be
considered. When a solute is introduced into the tube, the
concentration profile develops similarly to the velocity profile, since
no transport in the radial direction occurs. Taylor (189) discussed
three case for whic the initil1coditins-ad2th1mea
35
V
Co
Al 
CO
0 a
L*- u --t A3
rA A2
[ ott 
>i
t+x
FIG. 2. Distribution of mean concentration in three cases in absence of
molecular diffusion, after Taylor (189).
In the first case, Al, the solute is originally only present in a
short segment of the tube, with width X and concentration C . The
0
solute will be distorted into a parabola with its shape depending on
u(r). The mean concentration over a cross-section is now given by:
C =CX/(ut) O<x<u t t>O
m 0 0 0(3-24)
C = 0 x>u t t>O
mo
In the second case, A2, a solute with concentration C enters the tube.
0
This case can be solved by assuming that the constant initial
concentration for x<O consists of a number of thin sections as
discussed before, leading to the following mean concentration:
36
C = C x<O
m o
C = C (l-x/(u t)) O<x<u t t>O (3-25)
m 0 010
C = 0 x>u t
m 0
The third case, A3, deals with a solute which is initially confined to
a greater distance X than for the first case (189).
In all three cases, the concentration profiles are determined
merely by advection, i.e., by the shape of the 
velocity profile.
However, at low velocities the concentration is determined by molecular
diffusion as well. Subsequently, the second step is to include
molecular diffusion in the transport equation.
Taylor assumed that diffusion was significant for radial
transport and negligible for longitudinal transport. Therefore, the
following respective restrictions were imposed:
2a
2
2D At-2a or At D -D(3-26-a)
00
2D At-L or At<L 
(3-26-b)
o 0t
2
D
0
where At is the time required for a particle to travel through a tube
with length L.
For a circular tube, solute transport via diffusion and advection
is given by:
+=D + a-2 
- u
?
1 - ac(3-27)
at Oar rr ax
where D is assumed to be independent of concentration. To derive an
0
expression for the dispersion coefficient, as a result 
of the combined
action of diffusion and advection, we follow the 
treatment by Fischer
37
et al. (60). Deviations from the mean velocity and concentration 
were
described by:
u'= u(r)-u 
(3-28-a)
m
C'= C(r)-C 
(3-28-b)
m
where u' and C' are deviations from the cross-sectional 
means u and
m
C, respectively. Substitution of Eq. (3-28) in 
(3-27) results in:
at (C+C') = Do[r (C +C') +-(C 
+C') + (C+C')J+
at m o[ ar 
2  
m r a-r m ax 2 m
a(c +c')
- (u +u') m (3-29)
m ax
This equation can be simplified. First, the 
contribution of molecular
diffusion to transport in the longitudinal direction can 
be assumed to
ac
be negligible. Second, note that m 0. Furthermore, 
it is convenient
ar
to apply the following transformation of coordinates:
=x- u t
m
T = t 
(3-30)
Z = r/a
Eq. (3-29) can now be rewritten as:
acm -ac, Dra
2
c, ac, U)ac + ac 
(3-31)
aT aT a2 az
2  
I az a- a-
If we apply the operator 1 [ ] dA, i.e., the mean over a cross
W 
A
section A, to Eq. (3-31) we obtain:
ac
mT +uac =0 
(3-32)
where the cross-sectional mean is denoted 
with an overbar. Subtraction
of Eq. (3-32) from (3-31) results in:
38
C' _o 0 C' 1 ac' u C cc' ,8c'
aT 
2  
c + - u + u (3-33)
2 2 z 
8z
a LO8
8
2 ~
ac' aC'
It is reasonable to assume that u' ~ u' hence Eq.(3-33) can be
rewritten as:
D 2 l
D ra c' + - u'a 
(3-34)
a 2 2 Z az
8 2  
a
a L82
Because u'=0, the net addition of solute is zero for the system, which
moves at velocity u . After sufficient time has elapsed, a steady
m
concentration profile will have developed: the advective and the
dispersive flux balance each other.
ac'
We now want to solve 
the steady concentration 
profile (i.e.,a
aT
= 0. Using Eq.(3-22) and u =u /2, we can rewrite Eq.(3-34) as:
m o
2
a2c) lac au oa2
Sa C 1 2 (3-35)
az2 Zaz D 8a 2
82 0
After multiplication with Z, integration with respect to 2 yields:
2
au 2 4
SZ o oc + ( (3-36)
az D 4 4
o
where a1 is an integration constant. Division by Z and a second
integration result in:
2
au rz2 Z
4
C, o 8c - + fn Z + a 
(3-37)
C D 8a 8 16 
1  
2
0
where a2 is a second integration constant. Because no solute transfer
aC'
occurs across the wall of the tube, we have =0 at Z=1 and hence
82
21=0. 
Note that 
a2=C'(0).
The dispersive flux, i.e., the mass transport relative to the
39
moving coordinate axis, is given by:
j Au'C' dA 
(3-38)
which can also be expressed according to Fick's law (Eq.(2-2)):
ac
JD =-D M (3-39)
and hence:
8C -1
MAD =-jA 
8
am 1 [u'C' dA (3-40)
Substitution of Eq. (3-37), using a2=0, in Eq. (3-40) and subsequent
integration results in:
2
a2u
D = o (3-41)
192D
0
which is generally referred to as the dispersion coefficient 
rather
than diffusivity coefficient since it quantifies the spreading 
of a
solute due to the combined effects of advection and (transverse)
diffusion. Note that the dispersion coefficient D is inversely
proportional to the coefficient of molecular diffusion, D . Using 
the
0
continuity equation along with Eq.(3-32) and (3-33), the equation
governing longitudinal dispersion becomes:
ac a2c
m D 
m 
(3-42)
aT 2
Solving this equation, for appropriate boundary and initial conditions,
results in concentration distributions which are assumed to be caused
bydispersion only. Aris (9) extended the theory, by using 
different
velocity profiles and a concentration dependent D . In many 
situations
0
in the field, it might take a long time to obtain steady 
state
conditions and initially no Fickian diffusion will occur.
40
The Concept of Dispersion in Porous Media
Unfortunately, the dispersion process in a porous medium does not
lend itself very well to a theoretical description like the one given
by Taylor. Since the geometrical boundaries in a porous medium are
generally not known, it is impossible to describe transport at the
microscopic level in an exact mathematical way. In most cases it is
assumed that Fick's law, originally used to describe diffusion in free
liquids, can be applied. This law provides merely an operational
definition of the diffusion/dispersion coefficient in soils (126).
Most of the dispersion in a porous medium is caused by the
presence of solids (meandering, changing pore sizes) with their
concomitant microscopic velocity profiles in the liquid. To study
dispersion phenomena at the microscopic level requires detailed
knowledge of the pore structure. Since such knowledge is only available
for simplified cases, porous media models or statistical mechanics must
be employed. A variety of geometrical models has been suggested.
Klinkenberg (98), using a straight, non-interconnected, capillary pore
model, showed that pore size distributions obtained via displacement
studies, were quite different for different displacement methods. The
model does not seem to be very satisfactory. Random capillary models,
consisting of a network of randomly oriented and distributed pores,
have been presented by de Josselin de Jong (54), Saffman (164) and Bear
and Bachmat (17).
De Josselin de Jong (54) considered a 3-D medium which consisted
of a network of interconnected straight equal channels of length L,
41
oriented at random. The overall flow occurred in one direction. The
movement of a solute particle in each channel depends on the deviation
of the channel from the overall direction of flow. At the junction of
various channels, the probability that a particle chooses a particular
channel is proportional to the ratio of the discharge to that channel
to the total discharge at the junction. De Josselin de Jong (54) then
determined the probability of a particle arriving at a given point for
a particular time after a large number of displacements. He came up
with a normal distribution in three dimensions to formulate expressions
for the longitudinal and transverse dispersion coefficient.
Saffman (164) used a very similar model, but he included
molecular diffusion. The movement of the particle by advection was
studied using the Lagrangian approach, i.e., the movement of the
particle is described with respect to its initial position, whereas the
movement by diffusion was described with Euler's method, i.e., the
movement of all particles was observed from a fixed position.
It can be argued that only a completely statistical approach will
predict the phenomenon of dispersion (165). One can attempt to predict
the path of a solute particle by using probability theory. This path is
determined by physical processes, but we are unable to describe these
processes at a microscopic level. The pathline for a solute particle is
considered as a vectorsum of elementary displacements. The elementary
displacement in a time interval At, for a certain position, 
depends on
the statistical model, for instance a random walk model 
or a
probability density function as a function of time and position. 
After
42
many intervals At, the total displacement of that particle obeys a
Gaussian distribution according to the Central 
Limit Theorem. This can
be extended to a large number of particles 
at that particular position
to predict a spatial distribution of 
the particles at a certain time
(i.e., a concentration profile).
Scheidegger (165) derived the following 
probability that a
specific molecule is at a certain position 
at time t:
(x,.t). = (4irDt)-3/2 exp ](3-43)
4Dt (-3
where q is the probability and g=x-x is 
the position with respect to
the time averaged position, x, of the molecule. 
Scheidegger developed
longitudinal and transverse dispersivity 
constants from statistical
concepts, ignoring molecular diffusion. 
He obtained a linear
relationship between pore-water velocity 
and longitudinal and
transverse dispersion coefficients. This 
serves as a justification for
the well known empirical expressions used 
for the longitudinal and
transverse dispersion coefficients (14), namely:
D =D + a vn (3-44)
L e 
L
n
D D +cTv 
(3-45)
T e T
where a and oT are the longitudinal 
and transverse dispersivity,
respectively, [L], D and D are the coefficient 
of longitudinal and
L T
2 -
transverse di.spersion, respectively, [L T I, v is the 
pore-water
velocity [LT
- I
] and n is an empirical coefficient, frequently 
taken
equal to 1. De, the part of dispersion which is 
independent of the
velocity, is equal to the effective coefficient 
of molecular diffusion,
43
viz. D 
/A [L2T
-1 ]
.
0
The significance of the space-time scale at which dispersion
occurs has been discussed theoretically by Bhattacharya and Gupta (19),
who distinguished three scales: a kinetic or molecular, a microscopic,
and a macroscopic or Darcy scale. Transitions from one scale to the
next higher scale can be made through application of the Central Limit
Theorem. Transformation to field scales could be attempted as well to
study the so called macrodispersivity.
Plumb and Whitaker (145) used two length scales to study the
effect of the heterogeneity of the porous medium on solute transport.
First, they used local volume averages, with the ADE developed by
combining "point" equations within a homogeneous layer and interfacial
boundary conditions. Second they used large-scale averaging, i.e., a
length scale was used that is large compared to the scale of the
heterogeneity. The resulting ADE for the latter approach contained
additional terms involving the time derivative. The theory was applied
to predict dispersion in stratified and 2-D spatially periodic media
(146). The dispersion coefficient determined with the large-scale model
for such media compared favorably with experimental observations. The
value of the dispersion coefficients was several orders of magnitude
larger than the value predicted for a medium which was assumed to be
homogeneous.
Finally, the concept of multiple length scales should be
mentioned. The resulting transport problem can then be solved with the
use of fractals (213, 215).
44
Physical Non-equilibrium and Dispersion
The existence of physical non-equilibrium, which enhances
spreading, is an important factor when describing solute transport.
This is particularly true for structured soils, which will be discussed
in Chapter V. In order to retain the 'simple' form of transport
equations such as Eq.(2-5), attempts have been made to include the
effects of physical non-equilibrium in the effective dispersion
coefficient. This can be done by lumping together the dispersion
processes already discussed, namely mechanical dispersion and molecular
diffusion, and non-equilibrium spreading (cf. 53, 135). Following Bolt
(24), a general approach will be presented here, which describes
molecular diffusion, mechanical dispersion and apparent (physical)
non-equilibrium in terms of their respective diffusion lengths. This
practice is quite common in the theory of chromatography (66).
Summation of these lengths allows the use of one apparent coefficient
of dispersion.
According to Bolt (24), the diffusion length, Ldiff [L], is given
by:
L =ODe/J (3-46)
diff e V
This length can be associated with the distance in a column to which a
solute front has advanced such that the concentration gradient at the
column entrance is zero. At larger times, the effect of spreading due
to diffusion can usually be ignored, in particular for small values of
Ldiff"
Next, a similar expression for the dispersion length needs to be
45
found. It is assumed that the flow is laminar, i.e. , streamlines
coincide with pathlines. For a typical porous medium, where the pore
space consists of an inter-capillary and an intra-capillary part, the
coefficient of mechanical dispersion should contain terms accounting
for transverse diffusion into the intra-capillary pore space, the
microscopic velocity distribution (i.e., advective dispersion), and
mixing in the inter-capillary pore space. The latter process was
referred to by Bear (15) as intensive mixing. It takes place at
positions where fluid streams from different channels are
interconnected. This provides an opportunity to mix via diffusion or by
the (partial) combination of fluid streams. Mixing counteracts
advective dispersion, which occurs because of differences in the mean
pore velocity for the various tubes. Differences in solute
concentration due to physical non-equilibrium, which in effect leads to
longitudinal dispersion, are reduced by transverse diffusion. Because
mixing and transverse diffusion limit dispersion, they will be examined
in more detail.
Assuming that liquid flow in a porous medium occurs in
semi-infinite, circular tubes with non-limited advective dispersion,
i.e., the velocity profile causes all dispersion, Bolt presented:
D = (o')2<v>2 At (3-47)
dis
where At is the mean residence time of a solute particle in a tube. The
relative spread of velocities, (a-,)2, is given by <Av>2/<v>
2
, which is
a constant for a particular medium with average velocity <v> and
variations about that velocity of Av =Iv,<v>I. This part of the
46
dispersion in porous media was first developed 
by Taylor (189). By
reversing the direction of flow, this type of dispersion 
can be
annihilated (cf. 167).
Obviously, dispersion in a porous medium is much 
more complex
because pore channels are interconnected and 
not straight. This
situation can be modeled by using so 
called mixing compartments, figure
3, where complete mixing occurs at points x=O, x=L, x=2L, 
etc. The
effective dispersion coefficient has a value of 
zero at x=O and reaches
a maximum at x=L. With the mean residence time of 
a solute particle in
a tube between two mixing points being 
At=L/<v>, an average dispersion
coefficient can be used:
2
D = (o-')<v>L/2 
(3-48)
dis
Figure 3 illustrates the effect of mixing, 
at regular intervals L, on
the coefficient of mechanical 
dispersion according to Eq. (3-48). 
As can
be seen, 'unlimited' spreading occurs only over a 
finite distance, L,
between mixing points. An effective, average dispersion coefficient 
for
such a distance is used for computational purposes.
If transport in separate channels 
is the main mechanism behind
dispersion, this dispersion will be limited by mixing 
at intervals with
length L. Then the mechanical dispersion coefficient is presumably of
the form:
D =L <v> (3-49)
whreLcan be approximated from Eq. (3-48), i.e. , from (a ')2L/2.
whr dis
Bolt (24) suggested that a value equal to the grain diameter be used.
47
--
0
0 L 2L 3L
x [L]
FIG. 3. The value for D determined by mixing and the averaged value
dis
for D according to Eq. (3-48) (dashed line).
dis
Usually, however, L has to be quantified by experimental 
procedures.
dis
Theoretically, if the flow were to be reversed along the 
same
steamlines as before the reversal, only the contribution of molecular
diffusion to mixing would not vanish. In reality, a larger part of this
mixing is not reversible.
In addition to mixing, advective spreading is limited by
transverse diffusion. This mechanism becomes important if flow occurs
in a single tube or if no mixing occurs because of the absence of
inter-capillary pore space, as inside the aggregate. It is still
assumed that the porous medium consists of simple tubes for which
Taylor's (189) results can be used. An expression for Ddi
s 
was obtained
by combining Eq. (3-47) and an appropriate expression for the diffusion
48
time needed to traverse the width R over which the velocity profile is
spread. Bolt proposed:
D (') 2 2v>2gR2/D (3-50)
dis e
where g is a geometry factor. This expression is similar to Taylor's
result for a circular tube (Eq.(3-41)). The dispersion length becomes:
L ( ') 2<v>gR /D (3-51)
dis e
In porous media, mixing and transverse diffusion (i.e., in a
single tube) occur simultaneously. The question then arises as to which
of the two prevails. This can be estimated by evaluating Eq. (3-48) and
(3-50). Since both processes counteract spreading, it can be assumed
that the process which results in the smallest value for Ddis
dominates.
Physical non-equilibrium occurs when flow in the intra-capillary
pore space is negligible. The stagnant phase of the liquid in the
intra-capillary pores does not participate in any (turbulent) mixing.
Rather, mixing between mobile and immobile phases of the liquid occurs
via transverse diffusion only. A description of dispersion in the
stagnant phase can be obtained by adapting Eq.(3-50) to account for the
limited accessibility of the liquid in the immobile phase. This was
done using a functional relationship between mobile and stagnant liquid
concentrations. This relationship can be presented, assuming a linear
increase of the surface concentration of the 
aggregate with time,
C =kt, as (45, 135):
mo
49
C -<C. > = kg(R )2/D 
(3-52)
mo im a 
e
where Cmo is the concentration 
in the mobile phase [ML 
3, <Cim > is themo 
' im
average concentration 
in the immobile or stagnant 
phase [ML-3], k is 
a
ratecoeficent[ M-3T
- 
1
rate coefficient 
[ML T ] and R is 
the radius of the aggregate 
[L].
a
The geometry factor 
of the aggregate, g, 
varies between 1/8 
for
infinitely extended 
cylinders to 1/15 for 
a sphere. By using 
the solute
flux over the entire 
liquid phase, Passioura 
(1971) obtained the
following dispersion 
coefficient, as caused 
by the stagnant phase
effect:
2 2
Ddis = go. R2 <v>2/OD 
(3-53)
dis im a e
with accompanying 
dispersion length:
L = ge. R
2
<v>/OD 
(3-54)
r Im a e
The effects of longitudinal 
diffusion, advective 
dispersion
(microscopic variation in 
the velocity profile in 
combination with
mixing and/or transverse 
diffusion) and the presence 
of a stagnant
phase can presumably 
be added since they act 
more or less
independently. It should 
be noted, however, that 
this last assumption
deserves more investigation 
(123). In the case of 
autonomous effects,
the diffusion/dispersion 
flux can be written as:
BC
D = -JVLD x (3-)
D diff dis r e 
V dis av 
e
50
It appears that, in many instances, mixing is the major
counteracting factor of advective dispersion. Ldis can therefore 
be
assumed to be constant ( 2R ) L is only important at small fluxes,
a dif
while L (dispersion length due to stagnant phase) becomes significant
r
for larger values of JV and increasing aggregate size R .a An increasing
aggregate size also causes the breakthrough curve to lose its sigmoidal
shape (22). Because of the relative importance of the larger
(inter-capillary) pores for transport, the tracer will appear sooner in
the effluent with larger aggregate sizes.
It should be noted that Ldiff is inversely proportional to JV
(Eq.(3-46)), and assuming that mixing prevails (Eq.(3-48)), Ldis 
is
approximately equal to the aggregate diameter, i.e., independent of the
flow, and Lr is proportional to JV (Eq. (3-54)). Despite the fact that a
number of simplifications and assumptions were made, the correspondence
between the experimental and theoretical behavior of L
D 
was fairly
good, figure 9.8 (24). In terms of the Pe-number, defined as Pe=vL/D
e
(L=2R), diffusion is the dominant term for Pe<1 and the stagnant phase
a
effect becomes dominant for Pe>120.
Experimental Determination of Dispersion Coefficients
To describe the dispersion process and to solve the ADE, a
numerical value for the longitudinal, and sometimes for the transverse,
dispersion coefficient needs to be determined. For one-dimensional flow
in a homogeneous medium with two-dimensional transport, the ADE can be
written as:
51
ac ac ac ac (357)
at - DL 2 x DT .a52
ax
2  
- ay2
D
L 
and D
T 
can be obtained by relating experimental results of the
concentration profile with explicit expressions for C. For this
purpose, displacement experiments are carried out, either in the
laboratory or in the field. For the determination of D these
experiments often involve a laboratory column filled with a packed
soil. The solute concentration of the resident solution in the column
is known, and a feed solution, with generally the same concentration of
a different solute, is leached through the column. In most cases the
effluent concentration is determined as a function of time, although
the concentration in the soil solution, at different locations along
the direction of flow, can be determined as well. A variety of methods,
used to obtain dispersion coefficients from these experimental dat.,
will be presented here. A review of some of these methods can also be
found in van Genuchten and Wierenga (203).
Because macroscopic variations in water content and pore water
velocity result in additional spreading, the column should be packed in
such a way that the soil is homogeneous with respect to the advection
term. Although this situation will never be reached totally, a constant
'effective' pore 
water velocity can 
be used for many 
practical
purposes. Deviations from this velocity during the experiment obviously
lead to incorrect results.
Mixing due to viscosity and density differences between resident
and feed solution are not to be included in the coefficient of
52
hydrodynamic dispersion. Therefore, many studies involve the so called
tracer case (e.g., 20). The low concentration of a non-reactive solute
(i.e., a tracer) does not affect the density and the viscosity of the
solvent. Displacing and displaced fluids have the same density and
viscosity. Differences in viscosity and density lead to instabilities
(fingering). Biggar and Nielsen (21) showed that unstable flow
-3 -3
dominates for a density difference of 3.4x10 g cm and a difference
in viscosity of 0.003 cP, obscuring the effect of molecular diffusion.
Rose and Passioura (161) demonstrated that, during horizontal
displacement experiments, small differences in density of the displaced
and displacing fluids lead to quite different dispersion coefficients.
The displacing liquids consisted of water, 0.1, 0.2, and 0.4 M NaCl
solutions with densities of 0.99707, 1.00116, 1.00523, and 1.0133 g
-3
cm , respectively. The amount of gravity segregation could be
characterized with a gravity segregation factor 0:
= gkhAp/(rqvL) 
(3-58)
-2
where g is the acceleration due to gravity [LT ], k is the
permeability [L2], h the height of the porous medium in the transverse
-3
direction of 
flow, [L], 
Ap is the difference 
in liquid density 
[ML ],
^t is the mean of the viscosities of the resident and feed solution
[ML T-1, v is the mean pore-water velocity [LT
-
] and L is the length
of the porous medium in the direction of flow, [L]. Since k decreases
rapidly during desaturation, gravity segregation is 
not much of a
factor during unsaturated flow. An increase in D by 
a factor 2 was
-4
reported by these authors for their experiments, with 
Ap=8.2x10 g
53
-3
cm compared to Ap=O. Gravity segregation might also be 
important when
contaminated salt water intrudes horizontally into a fresh water
aquifer. Because of the stratified flow and the subsequent mixing 
by
diffusion at the interface, contamination will extend further into the
fresh water than expected on the basis of the conventional hydrodynamic
dispersion process. Further discussions on viscosity and density
effects can be found in Krupp and Elrick (105), Bachmat (12), and
Scheidegger (165). Due to the temperature dependence of these
parameters, displacement experiments need to be carried out at constant
temperature.
One should also be aware of apparatus-induced dispersion (87),
which is especially important for short-column experiments and for
unsaturated media. These authors studied a conceptual porous medium
consisting of two different layers: the soil and the apparatus.
Apparatus-induced dispersion was assessed by carrying out displacement
experiments in the absence of the porous medium. 
The dispersion
coefficient of the porous medium was obtained by 
solving the
advection-dispersion equation for a two-layer system. 
The value of this
dispersion coefficient was as much as 40% less than 
obtained with the
one layer equation, where mixing was assumed to occur 
in the medium
only. The so called dead-volume, inside the apparatus 
but outside the
porous medium, should therefore always be kept to a minimum.
Rather than relying on effluent concentrations, in situ
determinations of the solute concentration can also be made. Harleman
and Rumer (82) measured electrical resistance with a conductivity probe
54
to determine NaCl-concentrations. Gupta et al. (76) used an Ag-AgCl
electrode to measure Cl concentrations in an unsaturated glass bead
medium. Electrical potential, in contrast with electrical resistance,
is not greatly affected by the water content of the medium. Kirda 
et
al. (95) used labeled Cl in the feed solution, of which the
concentration was measured in situ with a Geiger-Muller tube. Grismer
et al. (70) used a dual-source gamma-attenuation system to 
determine
salt and water content. In an error analysis, thse authors 
showed that
the accuracy of solute concentrations determined at low water 
contents
was limited. Grismer (69) used the same technique to determine SrCl
2
and Nal concentrations in a displacement study 
during horizontal
transient unsaturated flow. Initially, no salt was present 
and the feed
solutions had a molality of 0.205 and 0.1 m, respectively. 
No mention
was made of density or viscosity effects.
Following Fried and Combarnous (61), some 
methods to determine DL
and D will be discussed. The solute (tracer) concentration 
can be
determined as a function of position at a certain time or for various
times at a certain position, as with breakthrough experiments. For 
both
cases, analytical solutions exist. For a uniform medium and steady 
flow
(v and 8 are constant), the 1-D ADE was given as:
c D 2C -vc 
(3-59)
at L x2 8x
subject to the following initial and boundary conditions
C(x,t) = C x=0 t>O
0
C(x,t) = 0 x- t>0 (3-60)
C(x,t) = 0 x>0 t=0
55
Utilizing the definition of 
the complementary error function, 
the
solution of Eq. (3-59) may be written 
as (i.e., an approximation of the
solution by Lapidus and Amundson 
(109)):
1 x- vtl
C/C =2erfc [ 
(3-61)
/4D 
t
L
For a given time, the solution 
follows a normal distribution
1-N[(x-p)/r] with mean displacement 
=vt and standard deviation - =
-2D t. N[ J is the probability 
density function for a normal
L
distribution with values N[-1]=O.16 
and N[1]=0.84. The width of the
transition zone, 2o-, can be determined 
by plotting C/Co versus x, with:
2o = x0 -x0 = 8Dt 
(3-62)
0.16 084 L
where x and x denote positions 
for which C/C is equal to 0.16
0.16 0.84 
0
and 0.84, respectively. The value 
of the longitudinal dispersion
coefficient follows from:
D=(x -x 2/8t 
(3-63)
L 0.84 0.16
If the concentration is observed at 
a certain position as a function of
time, the following expression for D 
can be used:
L
x-vt x-vt
DL 0 i[/ .16 _ t0.84] 
(-4
t t
0.16 0.84
After the discussion of the determination 
of the longitudinal
dispersion coefficient, D, Fried and Combarnous 
(61) continue with the
determination of the transverse dispersion 
coefficient, D
T
. It should
be noted that relatively little in-depth 
work has been done in the area
56
of transverse dispersion, because of the increased experimental
difficulties and the much smaller effect of transverse dispersion in
comparison with longitudinal dispersion (e.g., 82).
For the determination of D in a semi-infinite 2-D system
T
consisting of a uniform soil, it is assumed that 1-D steady state flow
conditions exist and that steady 
solute transport is established
(C=C(x,y)), in which case the longitudinal dispersive flux is usually
negligible. Under these circumstances the transport equation becomes:
2
v D C (3-65)
8ax T 2
8y
where y is the distance in the direction transverse of the flow [L].
For the following boundary conditions:
a= 01x>O 
, y4?0
C(O,y) = C x=O , O<y<m 
(3-66)
0
C(O,y) = 0 x=O , -W<y<O
the analytical solution is:
C/C = 1+ erf[ y ]] 
(3-67)
0 2[ ((
4
xDT)/v)1
/ 2
The transverse dispersion coefficient determined from this solution is:
D v 2 
(3-68)
T 2x (
where the standard deviation at an arbitrary point x=x is given by:
0
= y(xC/C,=0.84) - y(x C/C =0.16)] 
(3-69)
2 O ' 0
Because relatively little work has been .devoted to transverse
dispersion, it is assumed in this review that we deal with longitudinal
dispersion, i.e., D=D
L
unless stated otherwise.
57
Rose and Passioura (160), making use of-the solution of Brenner
(27), determined D by plotting C =(C-C )/(C -C.) versus In T on
e i o 1
probability paper. Ce is the dimensionless exit concentration and T is
the number of pore volumes (T=vt/L). Almost straight lines for
particular values of the Brenner number were obtained. The Brenner
number is defined as:
B = vL/D (3-70)
where L is the column length [L]. Because of this linear relationship,
the following expression seems plausible:
inverfc (2C -1) = -a e T - ( 
(3-71)
e
where inverfc denotes the inverse of the complementary error function,
and a and ( are slope and intercept, respectively. For 16<B<640, the
following empirical relationship was proposed by Rose and 
Passioura
(160):
& B = 0. 1139 (og a)
3 
- 0.3504 (t a)2+
2.3623 to ga + 0.4732 
(3-72)
By plotting -tn T versus inverfc (2C -1) the slope a and the intercept
e
( can be obtained. The Brenner number follows then from Eq.(3-72) and D
can be calculated from Eq.(3-70). Van Genuchten and Wierenga (204)
presented various improvements of this method, viz. an inversion
formula for inverfc, a simplified expression of Eq.(3-72), and
approximations for various boundary conditions.
A popular technique to obtain dispersion coefficients is via
curve fitting. Parker and van Genuchten (130) provided a curve fitting
58
program, using a non-linear, least-squares inversion method, to
determine various transport parameters including 
dispersion
coefficients. The curve fitting program can be used 
for experiments
involving concentrations at a fixed location at different 
times or for
a certain time at different locations; 
other parameters to be
determined are retardation factors and degradation 
constants. The
program can also deal with physical and chemical 
non-equilibrium
transport models as discussed in the section on non-equilibrium
conditions, and a stochastic model of transport under 
log-normally
distributed flow conditions. An independent validation of 
the physical
model cannot be achieved by curve fitting (36, 149).
Amoozegar-Fard et al. (6) used the solution given by 
Eq. (3-61).
By taking the inverse complementary error function 
and applying an
error minimization technique, D and R were determined. 
Some other
methods, such as the point calculation, finding the variance, 
and the
slope method, are discussed in Levenspiel and Smith 
(114).
Smiles et al. (179) used Matano's method to 
study the
hydrodynamic dispersion during absorption of water by 
soil. This method
uses the Boltzmann transformation to transform a partial 
differential
equation into an ordinary differential equation (cf. 35, 
45). Because
of the relatively low velocity, the small value of 
the Peclet number
and the rapidly changing water content during infiltration, 
the
dispersion coefficient is assumed to be only dependent on water
content. By using the relations D (O)=SD(0), the product of water
5
content and the coefficient of hydrodynamic dispersion, and JvV, the
59
volumetric water flux, the one-dimensional transport equation
(Eq. (2-4)), can be rewritten as:
aeC 
a 
ac 
aJvC
- D() 
(3-73)
ax s - a x
with the following boundary and initial conditions:
C = C x>O t=O
i
C = C x=O t a-0 (3-74)
0
C = C. x-W t>0
1
where D (8) is an alternative expression for the dispersion
s
coefficient. It is assumed that the medium is semi-infinite; the
initial and exit boundary conditions transform to the same condition.
For a homogeneous soil with horizontal, non-hysteretic flow we
have:
Jv = -D ( ) e (3-75)
w ax
therefore Eq. (2-132) was rewritten as:
0
ac _ a [D () ac + D () ae ac (3-76)
t a Ds) a J w axax
where D is the soil water diffusivity [L2T-1. Using the Boltzmann
w
transformation (w =x/V/t), Eq. (3-76) becomes
d dC] dc
D (0) + -- 0 (3-77)
d- s d- + 2 dw
while the boundary and initial conditions, Eq. (3-73), become
C=C. -o
1 ( 3-78 )
The function g is defined as:
dO 8
g = Ow + 2D (o) -e = Iw O(-9
w d-w f dcC-9
60
The dispersion coefficient can then be obtained from
C
D (e) = d g(3-80)
s 2 - dC (-0
C.
1
in a way as described by Bruce and Klute (35) for D . D Ce) may be
w s
determined by carrying out experiments using segmented soil columns.
For each segment, average values for e, C, and g can be determined as a
function of w. Laryea et al. (110) determined dispersion coefficients
for both cations and anions during horizontal infiltration. Elrick et
al. (59) extended the above analysis to vertical flow. A power series
solution similar to that of Philip (139) was used to express the
solutions of flow and transport equations.
The dispersion process depends on viscosity, density, velocity,
water content, molecular diffusion, and permeability. In order to
assess the influence of velocity and particle size, experimental
results are commonly analyzed by plotting DL/Do versus the Peclet
number (Pe=vd/D) on a log-log scale graph. In this case, d is a
0
characteristic pore or particle size dimension. Bear (16) distinguished
the following dispersion regimes (61) based on values of the Peclet
number:
- Molecular diffusion is dominant for Pe<0.4.
- Molecular diffusion and mechanical dispersion are of the same order
for O.4<Pe<5. Both effects can be added up.
- Major mechanical dispersion with some molecular diffusion occurs in
the range 5<Pe<iO00. These effects interfere and cannot be added up.
61
- Dominant mechanical dispersion for 1000<Pe<1.5x1O
5 
with negligible
molecular diffusion.
- Mechanical dispersion for the flow regime being out of the domain of
5
Darcy's 
law for 
Pe>1.5x10 
.
Several authors have provided graphical relationships between velocity
or the Peclet number and the dispersion coefficient (e.g., 22, 24, 61,
138, 160).
Finally, the so-called numerical dispersion should be mentioned.
When the transport equation is solved by numerical methods, using
experimental values of the dispersion coefficient, the solute front
exhibits additional spreading. Especially at high values of the Peclet
number, 'smearing' of the numerical solution occurs around the front.
This is an artifact due to the numerical procedure (26). A number of
approaches can be taken to reduce this phenomenon (cf. 1, 29, 39, 195).
IV. ION EXCHANGE
Equilibrium transport of linearly exchanging solutes can be
described in a straightforward way with Eq.(2-19). However, for
transport of non-linearly exchanging solutes, the situation is somewhat
more complicated.
First, the influence of exchange on the solute front without
dispersion/diffusion is reviewed. Second, the effect of dispersion is
included. Several approximate analytical techniques can be used to
study transport under these conditions, which illustrate the effect of
non-linear exchange. In a majority of the studies, however, the
transport problem is solved numerically by using experimentally
determined exchange isotherms. Some of these studies will be discussed.
Equilibrium Chromatography
Following Reiniger and Bolt (156) and Bolt (24), the qualitative
influence of the exchange isotherm on the concentration front will be
studied for a relatively simple case. It is assumed that only two
different cation species are present during miscible displacement under
steady flow conditions. This provides a fairly representative picture
for many soil-water systems. The adsorbed concentration q of a solute
species depends on the liquid concentration of that species, C, and the
total concentration of all solute species, C It is convenient to
-3
express q in mol m , i.e., moles of charge per volume of porous
c
62
63
medium, and C in 
mol m , i.e., moles 
of charge per volume 
of
solution. The total 
concentration is constant 
during ion exchange. 
When
diffusion and dispersion 
are neglected, a step 
change in concentration
occurs at the "concentration 
front" located around 
Jvt/0, the
V
penetration depth. 
If C is assumed 
to be constant across 
this
T
concentration front, 
q depends only on C. 
For a particular 
species
Eq. (2-7) then becomes:
.(q' + ) acJ ac
at - v -a1
where q' (=dq/dC) is 
the differential capacity 
of the exchanger for 
the
exchanging species, 
JV is the Darcy flux 
which is assumed to 
be
positive, and x is 
the (positive) distance 
from the inlet. In order 
to
solve Eq.(4-1), and 
to justify neglecting 
the dispersive flux, 
aC/ax
needs to be finite. 
With the chain rule 
for partial derivatives,
ac ax at
(-) (c) (-x)c= -1, 
Eq.(4-1) can be rewritten 
as
(ax _ JV(4-2)
at-C -q' +-0
If the condition that 
ax be finite is violated, 
i.e. if the
Ox
concentration profile 
exhibits jumps, the 
conservation of mass
(Eq. (4-1)) needs to 
be expressed in an alternative 
way:
(Aq + 8AC) dx = AC dV 
(4-3)
where V is the input 
volume per unit area 
of the column. The rate 
of
propagation, given by Eq. 
(4-2), then needs to be 
rewritten as:
dx JV(4)
( d -) a = a q 
( 4 -4
AC
64
ac
The position of a 
particular solute concentration 
(for finite)
follows from integration 
of Eq. (4-2). For a 
soil with q and 0
homogeneous with respect 
to position and dV=Jvdt, 
the position for a
particular concentration 
is given by:
tj 
V-V 
(C)
Xc= J qi+ dt = o___
xC q+ t q')+0 
( 4-5 )
0
where the feed solution, 
with concentration C, 
enters a column having
an initial concentration 
C ; V (C) is an inverse 
feed function, which
0
is the volume of solution 
applied to the column at 
the moment that the
concentration at x=O reaches 
C (for step type displacement 
V =0). The
0
average depth of the 
concentration fronts in 
the adsorbed and liquid
phase is given by:
C C
.0 0
={XC (q'+ 0) dC 
/ (q'+ 0) dC 
(4-6)
C. C.
1 1
The rate of propagation 
of the solute in the liquid 
phase follows
from Eq.(4-2). Equilibrium 
chromatography enables 
us to determine the
propagation in the adsorbed 
phase as well by using 
the exchange
isotherm.
To study the shape 
of the solute front, 
Eq. (4-5) is
differentiated (V =0):
0
ax _ Vq'' 
(4.7)
OC V (q'+ 0)
At this point, a distinction 
needs to be made among favorable,
unfavorable, and linear 
exchange. If the incoming 
cation has a convex
65
2 2
isotherm, q' '=d q/dC2<0, we deal with so-called favorable 
exchange,
while the exchange is said to be unfavorable 
if q''>O and linear if
q' '=0.
In case of favorable exchange there 
is a dilemma. According to
Eq.(4-1), the slope of the solute front 
is negative (aC/ax<O), but
using Eq. (4-7) we find that aC/ax>O. The latter 
implies that the solute
travels faster at higher concentrations 
than at lower concentrations.
This is physically not possible 
for a step front and, as already
mentioned, the rate of propagation should 
then be determined based on
Eq.(4-4). For cases other than a step 
change in the concentration of
the eluent, a "self sharpening" effect 
will occur until a step front
has been established.
For linear exchange (q''=0), we can 
see from Eq. (4-7) that
ax/aC=O for any applied volume. The initial 
profile will therefore not
be altered during passage through the 
porous medium. If the initial
profile contains jumps, Eq. (4-3) needs to 
be used, otherwise Eq. (4-1)
suffices.
For unfavorable exchange we can conclude 
from Eq. (4-7) that for a
particular value of C, the front flattens 
with increasing V.
Unfavorable exchange induces solute spreading 
(i.e., decreases the
concentration gradient).
The Combined Effect of Ion Exchange and Dispersion
Bolit (24) examined the combined effect of exchange and
diffusion/dispersion using analytical techniques. Having expressed 
all
(physical) dispersion effects in one dispersion length, 
L
,
he
66
attempted to include exchange 
as well. The ADE in terms 
of the
dispersion length is obtained 
with Eq. (2-7) and (3-55):
raC aC a__(L 1C
(q'+e) -c +J a-c-( x -) 
1=0 
(4-8)
V[8x -X D axJ
Although the 'effective' value 
of L
D 
varies through the soil column,
Frissel et al. (62) showed that 
a column averaged value might be 
used.
For a constant LD Eq.(4-8) 
might be written as Eq.(2-14). 
If, in
addition, dq/dC=constant (linear 
exchange), the same analytical
solutions as for Eq. (2-19) can 
be used.
In case of favorable exchange, the 
solute front in the adsorbed
phase is generally located 
ahead of the solute front in the 
liquid
phase. In terms of the effective 
retarded velocity, v , the "high"
concentrations will travel at a velocity 
greater than v and the "low"
concentrations will travel at 
a velocity less than v . In contrast,
diffusion/dispersion causes just 
the opposite effect. Eventually, 
a
steady front will develop with 
respect to the moving coordinate 
x-v t
and all concentrations will travel 
at a velocity equal to v . For 
a
constant L Eq. (4-8) can be rewritten 
to obtain a general expression
for the propragation of the solute 
front:
( - qO _1-LD a(ac/ax) (4-9)
Once steady state has been established, 
we obtain:
* x J v
v = (8 - Aq 
(4-10)
AC
Based on these two expressions, Bolt 
(24) presented a method to
analytically determine adsorbed and liquid 
concentrations.
67
In case of unfavorable exchange, no steady front 
will develop
with respect to the moving coordinate, because 
dispersion/diffusion and
exchange both cause the "low" concentrations to travel faster 
than the
"high" concentrations. The ADE must now be solved numerically to get an
accurate picture of the solute front in the adsorbed and liquid 
phase
(Cf. 106, 107). For estimation purposes, one can add the two effects by
first solving the front position as determined by Eq.(4-5) (exchange
only) and then adding the effect of diffusion/dispersion. The latter
effect causes the solute front to spread around x =v t. With Eq.(2-38),
p
we can estimate this spreading according to:
x-x = 2 D t inverfc(2C) (4-11)
p
where C is the dimensionless concentration defined in Eq.(2-38). This
term can be treated as a perturbation term to be added to the front
position according to Eq.(4-3). In this way, a conservative estimate
was obtained for the maximum amount of spreading.
Finally, the occurrence of non-equilibrium during ion 
exchange
should be mentioned. In aggregated soils, the deviation 
from local
equilibrium is most likely to be caused by physical 
non-equilibrium,
i.e. the limited accessibility of exchange or adsorption 
sites. Bolt
(24) assumed that all adsorption takes place in the stagnant 
region of
the liquid. Examples of adsorption in both mobile and 
stagnant region
were treated by van Genuchten and Wierenga (202), Tang 
et al. (186) and
Selim et al. (170). Transport from 
the mobile to the stagnant region
occurs via a diffusion process (Eq. 
(2-44)), which was described by Bolt
(24) with:
68
a(q.m +m .mC
=im imCim k 0 (C -C.) 
(4-12)
at a mo mo im
where qim is the concentration in the adsorbed phase 
(immobile sites
only), k is a rate constant for diffusion 
inside the aggregate, and
a
instantaneous chemical equilibrium is assumed 
inside the aggregate. A
distribution ratio between mobile and immobile phase, 
defined as:
Kmi (Aq + 0. 0o(4-13)
D A AC im mo
can be used to rewrite Eq. (4-12) as
ac. k
im - a (C -C.) 
(4-14)
at mi mo im
Without longitudinal diffusion/dispersion 
the transport equation
becomes (Eq. (4-1)):
aC. aC aCSmi im mo mo
0 K + m-.0 =+ o 0 
(4-15)
moKD at mo at vax
Equations. (4-14) and (4-15) can be solved 
by introducing a position
dependent time and by using a transformation 
in order to obtain scaled
variables. The solution is expressed 
in the so called Goldstein
J-function (cf. 56, 196).
Instead of using the above approach, 
the stagnant phase effect
can be described by an equivalent length parameter 
L . For spherical
r
aggregates with radius Ra the expression 
of Crank (45) is adapted for
a situation where exchange occurs:
(q.m + 0. C.
im i m- iO. D (C -C. ) (4-16)
at R
2
im mo I
69
Combination of Eqs.(4-12) and (2-16) results in the following
expression for k?
a
0.
15D im (4-17)
k- 2 0to
R mo
a
With the use of Eq. (3-46), L can now be expressed in terms of k a Two
' r a"
different expressions were given by Bolt (24), depending on whether
adsorption/exchange takes place. The advantage of this approach is, 
as
was mentioned in the section on physical non-equilibrium 
and
dispersion, that all effects are accounted for by one effective value
of the dispersion length. Consequently, analytical solutions can be
used to solve the advection-dispersion equation.
Numerical Solutions of Transport Involving Equilibrium Exchange
The two preceding sections, discussed some of the fundamentals of
transport of exchanging solutes with and without dispersion. This 
was
done in a way so that explicit expressions for the position of the
solute front were obtained, and that (approximate) analytical solutions
were available. However, because of the non-linearity of the (measured)
exchange isotherms, numerical techniques are usually needed to solve
the transport equation. For that reason, this section focuses on 
the
(numerical) simulation of exchanging solutes.
One of the earlier works on transport of reactive solutes was
reported by Kay and Elrick (94). These authors determined adsorption
isotherms of lindane for various soils and soil fractions. Linear
adsorption was observed; the l indane was particularly strong adsorbed
by organic matter. A chromatographic model developed by Hashimoto et
70
al. (83) was used to predict the movement of lindane through the soil.
A reasonably good agreement between predicted and experimental
breakthrough curves was obtained. At higher pore water velocities, 
the
prediction was somewhat poorer, perhaps because of physical
non-equilibrium. Displacement of the lindane (elution curve) was
described poorly, possibly because of inadequate knowledge of the
desorpt ion curve.
Lai and Jurinak (106, 107) solved the transport equation for
non-linear adsorption with an explicit finite difference method. These
authors considered homovalent exchange in a binary system for
one-dimensional steady flow. The dimensionless variables Xk=Ck/C
T 
and
Y =S /S were introduced for the solute concentration in the liquid and
adsorbed phase, respectively. The following equation was obtained,
using a similar relation as Eq. (2-12):
2
aXk a2Xk aXk
at k D(X2k 
kV(Xk ax 
(4-18)
ax
where the dispersion coefficient and velocity were adapted according
to:
D(Xk) = D (4-19)
S1+ (PbST/CT) J (Xk)
v(X
k
) = v (4-20)
S1+ (PbST/OCT) f' (Xk)
with f' (Xk)=dYk/dXk
. 
The denominator of the latter two terms can be
recognized as the retardation factor; its value depends on X
k 
because
of non-linear exchange. Following Helfferich (84), various adsorption
functions Yk=f(Xk) were used. For a number of isotherms (linear,
71
concave and convex), numerical solutions were obtained for Eq. (4-18)
subject to specified boundary and initial conditions. Column studies
were carried out involving Ca and Mg, which exhibit a slightly
non-linear exchange isotherm. The comparison between experimental and
numerical results was good.
Rubin and James (163) presented a comprehensive analysis of
multi-component equilibrium exchange during one-dimensional steady flow
in layered and homogeneous profiles with a variabl e C T The equations
to describe exchange and transport were solved numerically with the
Galerkin method. These authors showed some interesting features, such
as multiple fronts and plateau zones during multi-species transport
with a varying C T. It was also demonstrated how ion exchange and
hydrodynamic dispersion influence solute transport..
Along the same lines, Valocchi et al. (192) presented an
analytical framework for transport of various species of ion exchanging
solutes, based upon the theory of chromatography. In contrast with some
previous work in this' area, the effect of hydrodynamic dispersion was
included and the total electrolyte level was not assumed to be
constant. The development of the transport equations is quite similar
to that presented by Rubin and James (163). In case of n exchanging
ions, the governing set of equations is:
ac k 
as k 
a 2 
a
0 atk OD k kJk=1) 2). . n(4-21)
nt bR 2 Va
72
C
T 
is the variable total solute concentration in.the liquid phase and
S is the constant, total solute concentration in the adsorbed phase or
T
the cation exchange capacity. In case an adsorbed species k, with
valence uk is exchanged for species j, with valence uj, the exchange
coefficient can be defined as:
Kjk = (Y/X.)k (Xk/Yk)k
j  
(4-24)
where Y and X again denote dimensionless concentrations. It should be
noted that a more accurate solution might be obtained by using chemical
activities. Furthermore, the assumption that Kjk is constant is
strictly not correct.
To solve the transport equations, we need to reduce the number of
dependent variables, by expressing S
k 
in terms of C
k . 
This can be done
by using the n-1 independent equilibrium expressions, and by using
Eq. (4-23). The multi-component exchange isotherm for a species k is
given as:
Sk = Fk(C
1 
C
2  
,. C ) k=1,2,....,n (4-25)
Eq. (4-25) can be substituted into Eq. (4-21), which leads to a set of n
transport equations of the following form:
8
Ck 
n 
ac. 
a2Ck 
8
Ck(4-
it bEjkD 
2tJDx(x4-26)
k=1 a 8 x
2
aFk 
ask
where fJk = aC. -C." This system has been solved for binary and
3 3
ternary exchange using the Galerkin finite element method (e.g., 163).
73
As was mentioned earlier, the value of 
Kjk is important as it
affects the shape of the front. For a concave 
isotherm, Kjk< 1 the
Xk-profile travels at a speed proportional to 
t instead of t ,
1/2 
as for
Fickian dispersion. For a convex isotherm, 
Kjk>1, the Xk-profile
becomes steeper when the front travels 
through the medium. Although the
situation with a varying CT is somewhat more 
complex, one can assume
that each particular concentration travels 
in a constant CT
environment.
A qualitatively good comparison of theoretically 
predicted
concentrations, using multi-component exchange, 
with experimental
results, involving the monitoring of Ca and 
Mg in a ground water
aquifer, was reported by Valocchi et al. (192). 
Finally, these authors
reported that dispersion induced exchange, 
a phenomenon reported by
Lake and Helfferich (108), which occurs 
as a result of dispersive
mixing across the solute front in case C
T 
varies and which leads to
changes in the adsorbed phase, was of minor importance.
Selim et al. (169) examined the concentration 
of 2,4-D during
infiltration and redistribution. This was 
one of the first studies
concerning reactive solute movement during 
transient, unsaturated flow.
Numerical solutions were obtained with an 
explicit-implicit finite
difference technique. The dispersion coefficient 
was obtained as a
function of pore-water velocity and the 
retardation factor was
described similarly as by Lai and Jurinak (106, 107). 
The concentration
profiles obtained from a field experiment, involving the application 
of
2 cm of an aqueous solution containing 50 ig/ml of 2,4-D, 
were in good
agreement with the simulated profiles.
74
Although the retardation factor, R, is not constant for
non-linear exchange, one may want to use an effective, constant value
for R which can be used over the whole concentration range. Valocchi
(190) introduced an effective distribution coefficient, Kd, as the
ratio of the step change in aqueous phase concentration and the step
change in sorbed phase concentration. The effective velocity of the
solute front can be given by (Eq. (4-2)):
v =v/R = dx /dt (4-27)
p
which provides a way to determine a (constant) value for R. The use of
a constant value for R is plausible for a number of cases. For "trace"
quantities of the species (i.e., C.<C), as for certain pollution
1 T'
problems, we can safely assume linear exchange. Valocchi (190) claimed
that for many situations involving non-linear exchange and
multi-species systems, a constant R can be used for predictive purposes
provided that dispersion is negligible. Valocchi (190) applied a mass
balance across the step change in concentrations to come up with:
-1 PbAS -1
v /v =R = (i +- -) (4-28)
0 AC
where C=C -C. and S=S -S., i.e. , the difference between final and
0 1 0 1
initial concentrations, can be used to define an effective K It
seems, however, that the dependence of K
d 
on the initial and final
concentration values limits its applicability to specific cases.
Bobbins et al. (159) used a model by Childs and Hanks (40) to
simulate water movement and solute transport. Transport was simulated
for a variety of cases: non-reactive solute transport, precipitation
75
and dissolution of lime and gypsum during transport, and multi-cation
exchange. The numerical model to simulate transport involving exchange
accounted for ion pair formation and it used activities rather than
concentrations. During concomitant lysimeter studies, soil solution
samples were obtained at various depths using porous cups attached to
sampling tubes. Solute concentration and electrical conductivity could
only be predicted satisfactorily by including precipitation/dissolution
and cation exchange. The adsorbed concentrations of Ca, Mg, Na, and K
followed from the CEC, the cation exchange capacity, and the
selectivity coefficients characterizing equilibrium exchange for
various pairs of cations. These relationships were (158):
CEC= SCa 
+ 
SMg + SNa + S
K  
(4-29)
1/2
C S C S C S
CCa Mg C
K 
Ca CNaMg
1 1/2 3 1/2 5 1/2
C S C S CS
Mg Ca Ca K Mg Na 
(4-30)
CNa Ca 
CK Mg 
CNa K
2 1/2 4 1/2 6
Ca Na 
Mg
S
K 
K Na
-3
where C is the concentration in the liquid phase [mol m- ] and S is the
-1
adsorbed concentration [cmol kg ]. By combining the above expressions
c
one can determine, for example, the amount of adsorbed Ca as:
C K C C
S CEC / Mg 1 Na + K + 1 (4-31)
Ca C 1/2 + C 
1
/2K C 
1
/2K
Ca Ca 2 Ca 3
In contrast with Eq.(4-25), a straightforward 
explicit expression is
now obtained.
76
The approach by Robbins et al. (159) to calculate complexation,
precipitation/dissolution, and 
cation exchange in separate subroutines
seems to be quite popular. 
The calculations are facilitated 
by the
availability of computer packages 
which can be used to simulate the
equilibrium chemistry involving 
many components and reactions. 
A
drawback is that convergence problems 
are reported to be common when
the transport equation and the 
model describing the equilibrium
chemistry (exchange, complexation, 
protonation etc.) are combined.
Jennings et al. (88) suggested, therefore, 
that the technique of direct
insertion into the main transport equation 
be used rather than solving
a number of equations simultaneously. 
This is basically the approach
discussed so far, i.e., the 
inclusion of a retardation factor 
in the
transport equation. The dependency 
of the solid phase concentration on
various quantities was expressed 
by Jennings et al. (88) as follows:
S = f(CST, t,x,aC/8t) 
(4-32)
A solution of the problem was accomplished 
by a Galerkin finite element
method (cf. 144). Kirkner et al. 
(97) elaborated on this technique.
It should be noted that the type of 
chemical reaction needs to be
considered when chosing the technique 
of numerical solution. Rubin
(162) distinguished six broad 
classes of chemical reactions during
transport, each with its own mathematical 
formulation. One is also
referred to the review article by 
Abriola (1), which discusses recent
work on the modeling of contaminant transport, 
including geochemical
and sorption models.
77
The exchange of K by Na, governed by the Gapon equation (63) was
studied by van Eykeren and Loch (194) for soil systems with 
a mobile
and an immobile liquid phase. Experimental results obtained 
with a
cation exchanger column compared reasonably well with 
numerical
solutions. Selim et al. (170) obtained breakthrough curves 
for
aggregated soils by also considering a mobile and immobile region 
of
the liquid phase. Experimentally determined exchange isotherms 
were
used to quantify the exchange reaction.
V. TRANSPORT IN STRUCTURED 
AND LAYERED SOILS
The structure of a 
soil, i.e., the spatial 
arrangement or
clustering of primary 
soil particles into 
secondary units called
aggregates, can have a 
great impact on the transport 
of solutes through
that soil. The REV, i.e., 
the smallest volume for 
which macroscopic
uniformity exists, is much 
larger for an aggregated 
soil than for a
non-aggregated soil. The 
dispersion in aggregated 
media is caused by
molecular diffusion, 
mechanical dispersion, 
and the previously
discussed "stagnant phase" 
effect. Spreading in a 
medium which consists
of homogeneous layers 
exhibits some additional 
characteristics in
comparison with non-layered 
soils where dispersion can 
be described
with a Fickian process. For 
aggregated and layered media, 
the medium
can conceptually be divided 
into homogeneous zones (mobile 
and immobile
regions, individual layers). 
In this way the transport 
problem lends
itself to an exact description.
First, the concept 
of mobile and immobile 
regions of the liquid
phase will be reviewed 
in more detail than 
previously, along with 
a
discussion of the 
velocity dependency 
of the dispersion 
coeffcient.
Next, some cases of 
transport through media 
with a particular structure
will be cited followed 
by some remarks about 
transport in layered
media. The discussion 
is quite general and 
merely concerns transport 
in
well defined systems. 
Much work is still needed 
to explain and predict
transport on a field 
scale, which will 
be discussed in the 
next
Chapter.
78
79
Mobile and Stagnant Regions
In order to describe the appearance of asymmetrical breakthrough
curves (tailing), van Genuchten and Wierenga (202) used the concept of
mobile and immobile regions of the liquid with diffusional transfer
between the two (see section on physical non-equilibrium). This concept
had been previously applied by Coats and Smith (43) and Villermaux and
van Swaay (206). A sensitivity analysis was carried out by van
Genuchten and Wierenga (206) with an analytical model, which indicated
that the tailing observed during unsaturated flow in an aggregated
medium was succesfully predicted. In two subsequent papers, the
experimental determination of breakthrough curves for aggregated,
unsaturated media was described (203, 205). Columns were packed with a
clay loam, aggregate size 2.0 or 6.3 mm, to study breakthrough curves.
The outlets of the columns were connected to a vacuum chamber
containing an automatic fraction collector. The solution was applied
with a syringe pump to maintain a constant flux. The solutes used in
these studies were H20 (203) and 2,4,5-T (205). It appeared that the
2
amount of immobile water increased with decreasing velocity and with
increasing aggregate size. A fairly good prediction of experimentally
determined breakthrough curves was obtained.
Since then, the concept of mobile and immobile regions of the
liquid phase has-widely been used to describe physical non-equilibrium
dispersion. This concept, however, applies not only to structured soils
according to de Smedt and Wierenga (55, 57). These authors showed that,
for unsaturated conditions in a non-structured medium (glass beads),
80
the effluent concentration could only be described by taking into
account an immobile water fraction 
as well. If no distinction 
was made
between those fractions, 
a very large value for the 
dispersion
coefficient was needed to 
fit the data with the advection-dispersion
equation.
On the other hand, it has been 
reported that no basis exists 
for
a subdivision into mobile and 
immobile regions of the liquid 
phase.
Philip (140, 141) performed a 
theoretical study of diffusion 
in a
semi-infinite porous medium 
containing dead-end pores filled 
with
immobile water. The cumulative 
diffusive fluxes in the mobile 
and
immobile region of the medium, 
J and J were defined as:
mo im'
Jmo = f {C
o 
(x t)-C m (xO))} dx, (5-1l-a)
0
00
Jim = { {Ci (xt)-C (xO)} dx(5-1-b)
0
where the concentrations are 
based on unit volume of porous 
medium.
These concentrations were solved 
analytically.
Philip (140) found that transient 
diffusion is only different for
mobile and immobile regions during 
the initial stage of the diffusion
process, when the relationship 
between time of infiltration 
and the
dimensionless cumulative diffusive 
flux obeys a (time)1
/2 
law, whereas
it obeys a (time)3
/ 2 
law in the stagnant region. Shortly 
after the
initiation of the diffusion process, however, the dimensionless
cumulative diffusive flux in both kind of pores obeys a (time)
I / 
law.
Phii ip ( 141 ) argued that the effect of dead-end-poros 
ity is
81
insignificant and the distinction of two regions of 
the liquid phase
seems unwarranted.
Philip's arguments were supported by Smiles and Philip (178) 
and
Smiles et al. (179) who observed a pistonlike 
displacement during
infiltration of a KCl solution in relatively 
dry soil. Their medium
contained 30% kaolinite in a matrix of fine sand. 
Work by Warrick et
al. (209) and Kirda et al. (95, 96) also seemed 
to support the
observation that, at least in non-aggregated media, no basis 
exists for
a subdivision into mobile and immobile regions. Awad 
(11), who did
similar displacement experiments as de Smedt and Wierenga 
(57) in a
sand, also concluded that no distinction between 
mobile and immobile
regions could be made.
Somewhat related to this issue is the dependence 
of the
dispersion coefficient, D, on the pore water velocity, 
v. Smiles et al.
(179) found no dependence of D on v during 
infiltration of a
KCl-solution into the dry kaolinite/sand medium. As was 
pointed out by
these authors, different results can be expected for 
coarse textured
soils. Smiles et al. (179) observed, that for the 
low Pe-numbers
encountered during infiltration into a dry soil, the water 
and solute
profiles preserved similarity in terms of W (=x/V/7) (see section on
experimental determination of dispersiuon coefficients). 
This implies
that D is velocity-independent. Saffman (164) and Pfannkuch 
(138) also
reported that D is independent of v, but Saffman (164) used Pe<1 as a
condition for D to be independent of v (Pe=vd/De, where d is the
average grain diameter). Values for v can be found by using Philip's
82
theory of infiltration, while d can be estimated with the Kozeny-Carman
relation. It appeared that for t>1.1 s, D was independent of velocity
under the conditions stated by Smiles and Philip (178). Groenevelt et
al. (72) studied dispersion in a dispersed clay paste, using a slit
model to represent the flow between clay platelets. The dispersion
coefficient was found to be virtually velocity independent and could be
equated to the coefficient of molecular diffusion. It should be noted
that in their case the Peclet number (Pe=wv/D where w is half the
width of the slit), was also very small (<1) because of the absence of
large pores.
Findings by Watson and Jones (212), who conducted similar
experiments as Smiles et al. (179), seem to indicate that hydrodynamic
dispersion is velocity dependent. In the literature, a linear
dependence of D on v has often been reported for (steady) flow at
higher values of the Pe-number. Yule and Gardner (219), for example,
reported a longitudinal dispersion coefficient for a Plainfield sand
which was linearly related to v. At large Pe-numbers (e.g. Pe>35) it is
plausible to assume a linear relationship between v and D (cf. 15, 82,
206).
Aggregated Media
The study of transport in aggregated media has received
considerable in recent years. Agronomists seek to optimally use
irrigation water and applied chemicals, which might be easily lost in
aggregated soils because of "bypass" flow. Chemical properties, such as
pe, p11, and pO
2
, might be quite different for inter- than for
83
intra-aggregate pores, leading to more complex sink/source terms. The
performance of ion exchange columns is the topic of many studies in the
chemical engineering literature (e.g., 155). Transport in these columns
is usually described in terms of a transport equation containing a
sink/source term which describes the exchange inside the solid
particles of the column. This section will examine how the transport
equation is formulated for aggregated media.
Passioura (135) distinguished between micropores, inside the soil
aggregates where solute movement occurs only via diffusion, and in
macropores, located between the aggregates where transport occurs via
viscous (advective) flow. The latter type of transport can be
completely dominant, and considerable work has been devoted towards a
better understanding of flow in macropores (e.g., 18). The liquid
present in the micropores is supposed to be immobile or stagnant,
whereas the liquid in the macropores is said to be mobile. Passioura
(135) used the following equation to describe transport in an
aggregated medium (cf. Eq. (2-44)):
o. aC.. a
2
C ac
im Im + mo =D mo mo (5-2)
o at at mo 2 mo ax
mo ax
where D is the dispersion coefficient of the mobile liquid phase and
mo
where the immobile phase, inside the aggregate, is treated as a
distributed sink. In general, C and C. are not known and one has to
mo 1im
resort to a simpler model (Eq. (2-5)) and use an effective value for D.
An expression for D should combine the advect ive dispersion in the
mobile region with the stagnant phase effect in the immobile region.
84
Therefore, Passioura (135) adapted the expression 
of Aris (9) for a
dispersion coefficient in a tube, to obtain 
the following expression
for "overall" dispersion in an aggregated porous 
medium:
D =D 0 /e + gv
2 
a2D 
(
mo mo T ma/D .(5-3)
where a is the radius of the aggregate [L], 
g is a coefficient
characterizing shape, 0
T  
is the total (liquid) volume fraction
available for the solute, 0 is the volume fraction 
of the mobile
mo
liquid phase, and Dei is the coefficient of 
molecular diffusion in
elm
the immobile 
liquid 
phase 
[L2T
-I
].
Following Taylor (189), one can introduce a moving coordinate
=x-vt. For large t, steady-state transport 
will occur so that C -C. =
mo im
constant. As advective dispersion can be neglected 
for a plane moving
at velocity v (g=constant), the solute flux 
across such a plane is:
0. 22 aC
im v a mo
dif 0
T  
15D 8(5-4ax
where Passioura (135) erroneously omitted the negative 
sign, and the
expression of Crank (45) was used for the 
concentration inside the
spherical aggregates. Because the steady-state 
assumption implies that
aC /ax=aC/ax, the following expression for the 
effective dispersion
mo
coefficient inside the aggregate was suggested by 
Passioura (135):
0. 2 2
D. - m va (5-5)
im 
T  15D eim
The effective dispersion coefficient for the liquid in the 
macropores,
which characterizes advective dispersion, was described with 
(cf.
Eq. (3-44)):
85
D =D + v a (5-6)
mo emo mo
where D is the coefficient of molecular diffusion in the mobile
emo
2 -1
phase [LT ] and v the velocity in the mobile phase (=vOT / m). This
mo T mo
leads to the following overall dispersion coefficient:
D = Dmo mo/T + Dim (-7)
Passioura and Rose (136) performed experiments to evaluate D
according to Eq.(S-7). Water retention data were used to estimate 
8
im,
which is approximately the volumetric water content at a suction of 75
cm of water. Pores which drained at suctions less than 75 cm were
considered to be macropores. D was obtained via the technique described
-6 -6
by Rose and Passioura (160). 
Values for Dei . between 4x10 
and 6x10
elm
2 -1
cm s were used, depending on the porous material. It is of interest
to examine the Brenner number (B) as a function of va. For low values
of va, the main contribution to B will be from advective dispersion,
whereas at higher values of va the stagnant phase effect becomes more
important. Therefore, we will have somewhere an absolute maximum for B.
1/2
This maximum will occur when va=(1D D 0 /8. ) . Displacement is
elm emo mo im
likely to be most efficient, i.e., closest to piston flow, for this
value of va.
The study of transport during flow of other fluids than water
might provide independent ways to characterize soil structure.
Millington and Shearer (118) discussed the effect of aggregation 
in
porous media on gas diffusion. These authors distinguished between 
a
solid phase and an intra- and inter-aggregate porespace filled 
with
86
waiter or gas. Various expressions for D eIDo were derived, where De and
D0 are diffusivities of the fluid in presence and absence of the porous
medium, respectively. Millington and Shearer (118) illustrated the
contrasting behavior of diffusion in the gas phase and ion diffusion in
the liquid phase. The ratio De /Do was found to be higher for
non-aggregated than for aggregated media as far as diffusion in the
liquid phase is concerned, while the opposite seemed to hold for gas
diffusion, except during very dry conditions.
Scotter (166), studying the preferential solute movement through
larger soil voids (cf. 174), considered two pore geometries:
cylindrical channels and planar voids. Theoretical breakthrough curves
for a non-reactive solute showed that a substantial amount of the
solute appeared in the effluent before one pore volume had leached
through, especially for large channel diameters. Van Genuchten et al.
(201) used the development by Scotter (166) to solve transport through
cylindrical macropores analytically. The cylindrical pore 'is supposed
to contain the mobile region of the liquid phase, while the surrounding
porous medium contains small pores in which the immobile region
resides. The adsorption sites were divided into a fraction in close
contact with the mobile liquid phase, and a fraction in contact with
the immobile region. A separate retardation factor was used for' both
fractions. A set of equations similar to Eq. (2-44) and (2-45) was
87
furthermore, dispersion in the macropore 
can usually be Ignored without
loss of accuracy.
The approach taken 
by van Genuchten 
et al. (201) seems 
to be
generally accepted 
to describe solute 
transport in aggregated 
porous
media. The overall 
transport equation 
contains both the concentration
in the mobile and 
immobile phase. 
The geometry of 
the aggregates
influences the diffusive 
transport between the 
two regions as well 
as
the average solute 
concentration inside 
the aggregates.
Valocchi (191) studied 
the validity of the 
assumption of local
physical and chemical 
equilibrium by comparing 
solutions of the
non-equilibrium and 
equilibrium type transport 
equations for aggregated
media. Two cases 
of physical non-equilibrium 
were considered, 
namely
diffusion in spherical 
aggregates and first 
order transport between
mobile and immobile 
regions of the liquid 
phase. The general 
equation
for reactive solute 
transport, assuming 
linear adsorption in 
both
mobile and immobile 
region, can be written 
as:
ac ac. 
a2c aceR 
mo iR m 
mo mo
0 R mot+0 mR 
=0im =0 D -0 v 
(5-8)
momoat imim at 
mo 2 mo 8x
8x
For spherical aggregates 
the physical non-equilibrium 
diffusion model
is:
a
Ci -- (x,r,t) 
r 
2
dr 
(5-9)
J0
Ca(x,r,t) r= C (x,t) 
(5-10)
AC AC
Rm a =D 1 (2 
a) 
(5-11)
im~t eim 2r Ar
r
88
where a is the aggregate radius and C is the local concentration
a
inside the aggregate. In case the transfer of 
solute between mobile and
immobile zones can be described as a first order 
process we get:
ac.
0. R(C - CRm)_(5-12)
im im at mo im
The transport problem is thus posed by Eq. (5-8) 
in combination with
Eq. (5-9) to (5-19) or in combination with Eq. (5-12). 
One case of
chemical non-equilibrium, involving first order kinetics, 
was studied
by Valocchi (191). The ADE was described as
0
ac as a
2
c ac (-30 - + P b at -OD 2 ev a 
(5-13)
ax
2
where. S is governed by the first-order linear kinetic 
expression
as - k
1
C - k
2
S (5-14)
with k and k as the forward and reverse 
rate coefficient,
1 2
respectively. This is referred to as the 
linear chemical
non-equilibrium model.
These models are well documented (150, 151). They were 
written in
dimensionless form and solved in the Laplace domain 
by Valocchi (191).
In order to quantify deviations from local equilibrium, 
a time moment
analysis was carried out. Aris (10) showed how absolute 
time moments,
m which are defined analagous to moments for position 
dependent
concentrations (Eq. (3-18)), can be determined with the solution 
in the
Lapl ace domain:
m= (-1)nl d- c(x s)] 
(5-15)
p s L ds
p
89
where X is the dimensionless 
position and C(X,s) is the 
Laplace
transform of C(X,t) with s as 
the transformation variable. The
advantage of this approach is that 
the solution in the Laplace domain
is fairly easy to obtain. 
Based on these moments, Valocchi
(191) obtained expressions for 
the first normalized moment and 
the
second and third central moment 
for the equilibrium and the three
non-equilibrium models for a 
Dirac-type input. In this way, 
mean
breakthrough time, spreading, and tailing 
could be characterized.
It was found that sorption rate 
limitations (non-equilibrium)
caused enhanced spreading and tailing. 
Non-equilibrium, however, does
not influence the mean breakthrough 
time. At large values for the
dimensionless rate parameter, F=k2L/v, 
where L is the column length,
2
the non-equilibrium and equilibrium 
results are similar. By comparing
the equilibrium solution with the 
more realistic non-equilibrium
solution, the error, which is made 
by using the local equilibrium
solution, can be established. This 
error increases linearly with
incresing Pe/F and decreases 
with increasing retardation factor.
Following Passioura (135), Rao et 
al. (151), and de Smedt and Wierenga
(57), Valocchi (191) obtained effective 
dispersion coefficients for the
non-equilibrium situations.
Before discussing transport in fractured 
media, which is very
similar to transport in aggregated 
media, a few other studies should be
mentioned. Rasmuson and Neret ni eks (153) obtained 
an anal yt ical
solution for the ADE in a porous medium 
consisting of spherical
particles. Rasmuson (152) extended this analytical 
solution to a case
90
where two-dimensional dispersion occurred. Other solutions 
involving
spherical particles were given by Skopp and 
Warrick (175) and Tang et
al. (186), while Sudicky and Frind (184) provided solutions 
for porous
media with rectangular voids. One of the first 
attempts to combine
analytical and experimental work was reported by 
Rao et al. (150), who
studied transport in a medium consisting of spherical aggregates.
Fractured Media
Because of interest in the disposal of hazardous 
wastes in
fractured bedrock, geologists and engineers have lately 
been studying
transport in fractured media. Fractures are important for the advective
transport of contaminants because of their relatively 
high water
conductance (bypass). Formulation of the transport problem is 
very
similar to that for aggregated media, i.e. , the two components are
advective transport in the fissure (the mobile region) and diffusion
into the rock (the immobile region). For a number of cases, analytical
solutions are available. In other instances, numerical methods need to
be used. The finite element method 
is particularly convenient to obtain,
numerical solutions for transport in fractured media (86).
Tang et al. (186) presented an analytical -solution for
contaminant transport in a single fracture. They considered a
radio-active contaminant with a finite migration distance because 
of
decay. Diffusion of the contaminant into the porous rock matrix slows
91
the fracture, (3) molecular diffusion within the fracture, in the
direction of the fracture axis, (4) molecular diffusion from the
fracture into the porous matrix, (5) adsorption onto the face of the
matrix, (6) adsorption within the matrix, and (7) radio-active decay.
For solute transport in the fracture, the following equation was used:
:C 82C OD . C. C
mo D mo + eim Im 
v mo
-t R - C + - (5-16)
at R 2 w mo bR 8x R 8z
mo az mo x=w mo
-1
where 2w is the fracture width [L], pw is a decay constant [T 1, R =1
w mo
+ - K is a face retardation factor, and 
K is the distribution
0 dmo dmo
coefficient in the fracture. For transport inside the matrix the
equation was:
ac. D . d2C.
im _ elm im 
C 
(6-17)
at R. 2 w im
Im ax
where D. is the effective diffusion coefficient in the matrix, R. = 1
elm im
+ K is the matrix retardation factor, and K is the
0 dim dim
distribution coefficient in the porous matrix. In both cases the
adsorption was assumed to follow a linear adsorption isotherm.
Equations (5-16) and (5-17) were solved analytically by Tang et
al. (186) for specified initial and boundary conditions with the help
of Laplace transforms. Some numerical examples were given as well. It
was found that the effect of dispersion cannot be ignored. If D=0
2 -1
m s , the solution lags substantially behind the general 
solution,
which used an expression for D as Eq. (3-36). This is in agreement 
with
the findings of Rasmuson and Neretnieks (154) who also 
investigated
the migration of decaying radionuclides in fissured rocks. 
They found
92
that advective dispersion in the fissures results 
in larger travel
distances for the radionuclides, therefore 
decreasing diffusion into
the matrix and subsequent decay inside 
the rock matrix. Agreement
between the analytical solution by 
Tang et al. (186) and a finite
element method solution by Grisak and 
Pickens (68) was generally good,
except for the transient case. The results 
demonstrated that matrix
diffusion can prevent severe contamination 
of underlying aquifers with
the (decaying) radioactive substances.
Neretnieks (119) treated dispersion in 
fissured rock, idealizing
the medium by assuming flow through 
parallel channels of different
size. The fissure width, w, is given by a 
distribution function f(w). A
pulse with concentration C
0
, introduced at the inlet, travels 
a
distance x in time t. in fissures with width 
w.. Therefore, w can be
11
expressed in terms of t. If the residence 
time for a fissure is less
than a certain time t, it delivers tracer. 
At the outlet, mixing among
solutions exiting from all fissures occurs, 
resulting in the following
dimensionless exit concentration:
T f(w)Jv(W) dw
C(t)/Co = wt - J(t)/J
v  
(5-18)
0 00 V V
T f(w)Jv(W) dw
V
0
The flow of water containing tracer, Jv(t), 
from fissures with widths
V
w(t)<w<o is diluted by the total flow of water, JV from 
all fissures.
An analytical expression for the mean residence time 
for a response
2
curve C(t) was Obtained. The spreading was quantified .with 
o- obtained
93
with the second moment of the time dependent concentration. The
dispersion coefficient depended strongly on the fissure width
distribution.
Neretnieks (119) referred to his approach as stratified flow,
i.e., the flow in each fissure depends on the size fissure just like
the flow in a stratified medium depends on the permeability of each
layer. In case of stratified flow, the width of the dispersion zone, 0-,
1/2
is proportional to the traveled distance x (117) instead of x as for
Fickian dispersion. Therefore, an increase in the observation distance
will yield a larger value for D if the ADE is used to describe the
transport process in stratified media (cf. 49). As stated by Neretnieks
(119), the implications of using the wrong mechanism (for instance
Fickian) when using a model for extrapolation to large distances may
have grave consequences in some applications.
Besides the mixing due to the advection term, mixing occurs due
to matrix diffusion. Diffusion of the tracer from the fissure into the
matrix of the rock is accompanied by sorption. Neretnieks (119) could
not obtain a mean and variance of the residence time with the moment
method because they appear to be unbounded (i.e., the curve exhibits
considerable tailing) even if the penetration depth in the rock is
small compared to the distance between fissures. The dispersion process
could be evaluated by comparing experimental data with theoretically
calculated response curves. By solving the analytical model, Neretnieks
(119) demonstrated that, for larger distances (larger residence times),
the effects of stratification and matrix diffusion become major" factors
in comparison with the effect of hydrodynamic disper-sion.
94
Layered Media
Solute transport in layered media has drawn a considerable amount
of attention (81, 116, 117, 142, 182, 185, 188). Almost all studies
concern stratified aquifers with flow parallel to the stratification.
Besides different soil-physical properties for each layer, differences
in transport parameters occur. Consequently, a different transport and
flow behavior might be expected, both perpendicular to and parallel to
the stratification, compared to a uniform medium.
Tanji (188) studied the leaching of boron in a .soil column
consisting of homogeneous layers, the stratification was perpendicular
to the direction of flow. Each layer had particular adsorption and
moisture characteristics. A scaled down field profile was used to
determine the amount of leaching experimentally. Calculated results for
the stratified columns were in good agreement with the measurements.
As mentioned earlier, layered media tend to enhance the
dispersion effect. Due to transverse diffusion from layers with a
significant advective transport to layers with minimal advective
transport, spreading of the solute front will occur in case of flow
parallel to the stratification (67). When using the conventional ADE, a
large value for the dispersion coefficient will be found. It was argued
by Gillham et al. (67) that under such heterogeneous conditions, no
Fickian dispersion will occur. The dispersion coefficient is time
dependent due to the transverse diffusion (64, 116). Gupta and
Bhattacharya (74) and Grven and Molz (80), on the other 
hand, pointed
out that Fickian dispersion will occur at large 
times for nonuniform
velocity fields.
95
Sudicky et al. (185) and Starr et 
al. (182) discussed the
migration of a reactive and a non-reactive 
solute in a layered porous
medium. A system consisting of a thin 
sand layer between two silt
layers was used to study the movement 
of a Cl solution applied to the
sand layer for horizontal flow 
parallel to the stratification. The
important processes were advection 
in the sand layer, transverse
diffusion from sand to silt layers, and 
molecular diffusion (transverse
to the flow direction) in the 
silt layers. The medium was water
saturated. For the sand, where complete 
mixing was supposed to occur,
the following transport equation applied:
C
2
C aC
a mo D 
ay - v 
ax 
(5-19)
at T ay2 ax
where the diffusive flux between sand 
and silt layers is represented as
transverse dispersion, C is the 
concentration in the sand layer, and
mo
v is the pore water velocity, in the 
sand layer, in the x direction
parallel to the stratification. Transport 
in the silt layers, which is
supposed to be purely diffusive, is 
described with:
ac. aDC.
im_ D 
im 
(5-20)
at e ay2
where C. is the concentration in 
the silt layers and D is the
im 
e
effective diffusion coefficient.
Solutions of these equations for appropriate 
boundary and initial
conditions, under the assumptions that the 
thickness of the sand layer
can be ignored, were given by Tang et 
al. (186). As the sand layer
becomes thicker, a transver-se concentration 
profile develops in the
96
sand layer. The analytical solution for the "thick" sand layer is more
complicated (185).
The exit concentration was determined as an 
average over the sand
layer with a fraction collector. The solute arrival 
time for the
system, with transverse diffusion from the 
sand into the silt layers,
was later than for pure plug flow. 
Strong dispersion seemed to occur
due to interlayer solute transfer. The experimental 
breakthrough curve
was poorly fitted when transverse transport 
was ignored.
Starr et al. (182) studied the same problem 
for a reactive
85
solute, Sr. Since the Sr-concentration was 
low, the exchange-isotherm
was assumed to be linear and non-hysteretic, 
hence a constant K
d 
could
be used. Application of the conventional 
ADE gave poor results for the
reactive case, especially at higher velocities. 
Neglecting longitudinal
dispersion, transport in the sand layer was 
described by:
8C 0. D ac. ac
mo _ im e im v mo C 
(5-21)
at 0 R w 8y R 8x mo
mo mo mo
y=w
where R is the retardation factor in 
the sand layer, which has a
mo
-1
width 2w, and W is the radio-active decay constant [T ]. Transport in
the silt layers was described by:
aC. D O2C
im 
_ 
e 
im(5-22
at R. 
2 -72 C im 
(522)
im ay
A diffusion coefficient for Sr in dilute solutions 
of 1.33xIO
2 -1
m s was used. Since a correct value for the retardation 
factor iS
essential for the prediction of the breakthrough curve, this 
factor was
determined in two independent ways for the silt material. 
The first
97
method determined K
d 
by the batch procedure, which resulted in values
of 1.65 mL/g for the sand and 12.5 mL/g for the silt. R could 
be
determined with Eq.(2-15), resulting in R =10 and R. =60. 
The second
mo im
method was via diffusion. Two diffusion half cells, one of them
containing a solution with 8sSr and the other free of 
85
Sr, were
brought in contact with each other. After 35 days the columns were
sliced in increments of 0.5 cm. The concentration in each segment was
determined and, via curve fitting of the analytical solution, an
effective diffusion coefficient was obtained. The ratio of the
effective diffusion coefficient and the diffusion coefficient in water
provides 1/R. A value of 70 was found for R. . The agreement between
the values obtained with the two methods was considered to be good. The
experimental breakthrough curve could, however, not be as well
predicted as for the non-reactive case. Better results were obtained if
the values of R were varied with v. However, no physical basis exists
for such an adjustment.
In the study of transport parallel to the stratification,
frequent use is made of the moment method of Aris 
(10), already
discussed previously. By using depth averaged moments, 
one can account
for the heterogeneity due to the stratification (115). 
Therefore, this
method, following Fischer et al. (60), will be discussed 
briefly. These
authors considered transport in a medium with depth 
h, bounded by
horizontal impermeable layers at y=0 and y=h. The magnitude 
of the
flow, in the x-direction, depends 
on y. For a moving coordinate system
the transport equation is:
98
-ac DrC + a 
-v' 
(5-23)
at 2 ay2 a
where v'=v-v and (=x-vt. Having the overbar 
notation indicate a depth
averaged value:
- 1h
v = d v dy 
(5-24)
Spatial moments can be defined according 
to (cf. Eq.(3-19)):
M P
m = C((,y,t) d( (5-25)
p
The operator defined by Eq. (5-25) can be applied 
to the individual
S 8C
terms of Eq.(5-23). Assuming that C=O anda = 0 for ( +m, we can find:
am a2m.
a -t D p(p-1) mp
2 
+ 2 +Vpm (5-26)
at [-2 y2 pp 
1
This equation can be averaged over depth, according 
to Eq.(5-24), and
the resulting equation is:
am
p = Dp(p-l)m 
2
+ pv'mp (5-27)
at p-2 p1
We can now solve m sequentially for p=0,1,2,... Using the 
theory of
P
statistical moments we can characterize the concentration distribution
by determining the mean, variance, skewness, kurtosis, etc. (81).
Fried and Combarnous (61) described the work by Marle et al.
(115) which showed how the moment method can be applied to a
multi-layered medium with flow in 
the horizontal x-direction and for
which the transport properties depend 
on y. The mathematical
formulation of the problem is:
99
CC C CC
e(y)a = (y)DL(Y) 2 +e DT(Y))v 
(5-28)
t 
2  
y Tay
ac 0 
Y = l and 
y = 2 
(5-29)
where it is assumed that v is independent of x. The boundaries 
of the
stratified medium are at y=y
1 
and y=y
2
. The problem can be solved in a
similar manner as described before, using depth averaged moments. Fried
and Combarnous (61) also described how an equivalent 
dispersion
coefficient can be used to represent dispersion in a layered medium.
Gtiven et al. (81) applied Aris' moment method 
to describe
dispersion of a tracer in a horizontally stratified 
aquifer of finite
thickness and infinite length and with vertical 
variations in the
hydraulic conductivity. Several idealized hydraulic 
conductivity
profiles were used. The experimentally and 
theoretically determined
macro-dispersivities were of the same order of 
magnitude. These authors
stressed the importance of detailed measurements 
of the hydraulic
conductivity profile at various positions along the aquifer. 
In a later
paper, Guiven and Molz (80) extended this work to an unbounded aquifer
with flow components parallel and perpendicular to the stratification.
The variation of the hydraulic conductivity, indicative of 
the
heterogeneity of the medium, strongly influences D
L 
via variations in v
(79, 183). Unfortunately, porous media generally do 
not consist of
homogeneous layers with distinct transport and flow parameters; v(y) is
not known in sufficient detail to predict the concentration profile
with the moment method. An even more complicated case arises if v
varies in both the longitudinal and the transverse direction.
100
A simplified approach for vertical transport might 
be to consider
only one-dimensional flow. Heterogeneity 
in the direction of flow is
most easily modeled by assuming 
that the soil consists of homogeneous
layers, i.e., horizontal layers (cf. 
171). Bresler and Dagan (30) used
the approach of representing a field 
as an ensemble of individual
columns, i.e., vertical layers, 
restricting the transport and 
flow
problem to one dimension.
Simmons (172) described dispersion during 
one-dimensional solute
transport in an arbitrary medium 
in terms of stochastic-advective
transport. He contributed the spreading 
in laboratory soils largely to
local variations in the hydraulic 
conductivity, which was incorrectly
interpreted as hydrodynamic dispersion. 
In the field, the heterogeneity
of the soil will be even larger, 
and the use of laboratory-measured
hydrodynamic dispersion coefficients 
is inappropriate. The stochastic
flow velocity depends on time and 
position, i.e., v(x,t). In contrast
with Bresler and Dagan (30), 
who used an ensemble of columns, 
the
description of the transport is truly 
one-dimensional. Simmons (172)
also questioned the validity of deterministic 
descriptions of transport
during transient (unsaturated) flow 
conditions because of the time
dependence of the dispersivity.
It has been suggested that for large distances 
or times, Fickian
behavior could be assumed (asymptotical 
approach of diffusive type of
transport). However, Matheron and de Marsily 
(116) showed that this
assumption does not always hold, even for a 
Gaussian covariance
function of the velocity. This was attributed to 
the lack of transverse
101
mixing. These authors studied horizontal transport in a medium 
with
horizontal layers. In case of a small constant vertical velocity, 
i.e.,
the flow vector is somewhat inclined with respect to the horizontal
stratification, transverse mixing is enhanced and diffusive behavior is
approached much faster.
Among other studies, dealing with stochastic modeling of
transport in heterogeneous media, the work by Gelhar (64) can be
mentioned. However, a further discussion of the. topic is beyond the
scope of this review.
VI. PREDICTION OF SOLUTE TRANSPORT UNDER FIELD CONDITIONS
Prediction of solute movement under field conditions 
is the
primary objective of many workers in the area of solute transport. 
The
application of the conventional advection-dispersion 
equation, ADE, has
proven to be quite successful for laboratory experiments (see 
section
on experimental determination of dispersion coefficients), 
where in
general macroscopically homogeneous conditions exist. Once 
the relevant
transport and flow parameters 
are known, the ADE is solved, 
generally
numerically, and a unique concentration distribution can be 
determined.
This is referred to as the deterministic approach. If necessary, 
the
structuredness of the soil can be taken into account by using various
non-equilibrium models (Chapter V).
Some of the earlier work 
involving the prediction of solute
transport under field conditions made use 
of the ADE, while the
corresponding transport parameters were measured 
in the laboratory.
This approach, also referred to as deterministic, 
leads to a uniquely
defined outcome for a given set of events. 
However, the description of
solute movement in the field, using this 
deterministic approach, was
often not very succesful. Three reasons 
were pointed out by van
Genuchten and Jury (199): (1) inaccurate 
or incomplete description of
the chemistry affecting transport of reactive 
solutes, (2) preferential
movement of water and solutes through macropores, 
and (3) spatial and
102
103
temporal variability of field scale flow and transport properties. In
order to deal with some of these problems, a variety 
of approaches has
been taken to model solute flow in the field. The variability 
on a
field scale has recently been investigated (48, 181).
Deterministic Models
An example of the use of experimentally determined 
data to
numerically predict solute profiles under different conditions 
is the
work by Bresler (28). Bresler used data from Warrick 
et al. (209), who
applied a pulse of 0.2 N CaCl during infiltration on a 
6.1 by 6.1 m
2
plot. The chloride concentrations in a 1 by 1 m square in the 
center of
the plot were monitored at depths of 30, 60, 90, 120, and 
150 cm, using
duplicate ceramic cups. The dispersion coefficient 
was obtained by
fitting the experimentally and analytically determined 
profiles.
Bresler used these data to validate a numerical model to 
predict solute
profiles. The model used an implicit method, which minimized 
the effect
of numerical dispersion. After validation, the model 
was used to
predict solute concentrations during infiltration, redistribution, 
and
evaporation. In a later paper, Bresler (29) discussed 
anion movement
under transient conditions. He included ionic diffusion, 
anion
exclusion, mechanical dispersion, and osmosis.
Selim et al. (169) determined the redistribution of 
2,4-D and
water both numerically and experimentally in a field 
study. The solute
2
concentration was monitored, at various depths, for an area of 1 m
using ceramic solution samplers connected to a vacuum pump. The water
retention curve, h(o), and the functional relationship between
104
volumetric water content and unsaturated hydraulic conductivity, K(O),
were determined in the laboratory on undisturbed soil cores. The
adsorption isotherm and the D(v) relationship were determined as well
in the laboratory. The water distribution was predicted fairly well,
whereas agreement between measured and predicted distributions of the
2,4-D was only fair. These authors were able to show that the
cumulative amount of irrigation water, and not the intensity of
application, determined the penetration depth of the herbicide.
The above examples concern rather small studies for which
homogeneous conditions might be expected. In larger field studies, this
is not likely to be true. The macroscopic heterogeneity ofthe advection
and dispersion terms, was illustrated in a large field study by Biggar
and Nielsen (22). They used an analytical solution of the
advection-dispersion equation to fit the experimental data. In their
field study, 20 plots were ponded with water containing C1 and NO
3
.
Solution samples were taken via extraction through porous ceramic cups.
They observed a large variation in v and D for the various plots with
both variables following a log-normal distribution.
Researchers who applied the conventional ADE under field
conditions found significantly.higher values for the dispersivity than
for column experiments in the laboratory. Gillham et al. (67) gave
-4 -2
values for 0L in the range of 10 to 10 m for laboratory
experiments, while for heterogeneous sand or gravel aquifers, a value
for the macrodispersivity of 1 to 100 m was found when the ADE was used
to fit the field data. Biggar and Nielsen (22) obtained a value of 2.93
105
cm f ro via regression of D versus v for 359 measurements 
at
different depths and plots in case of vertical transport.
.The concept of dispersion, as outlined in the introduction to
Chapter III, does not account for macroscopic non-uniformities 
in the
advect ion term, which is a major cause of mixing 
in the f ield. If one
wants to proceed with the deterministic approach, a detailed 
knowledge
of basic soil physical properties as K(O) and h(8) is needed. 
However,
in view of the field heterogeneity of these soil 
properties, a
considerable experimental effort is needed.
Stochastic Models
To circumvent the problem of a detailed quantification 
of the
flow term, several researchers have chosen for a different 
approach by
using stochastic models. These assume that the outcome of a 
series of
events is uncertain. Due to the uncertainty of the mechanisms 
and the
values of the parameters involved, a concentration distribution 
is
described in terms of probability. These stochastic models 
can be of a
mechanistic nature (49, 137) or a non-mechanistic nature (90). 
In the
first case, an attempt is made to describe the fundamental 
processes
contributing to transport, whereas in the latter case no assumption 
is
made for the mechanism behind the flow and transport processes.
Mechanistic models
106
result of a step input of a solute. The statistically homogeneous
field, i.e., stationarity exists as discussed by Journel and Huijbregts
(89), represents an ensemble of columns, each with uniform soil
properties. A generic notation of soil parameters, a., was used to
describe the soil properties in the field, e.g.,
6=O)z't' a , 2
"' '
.a )" In order to describe the distribution of
particular soil properties, a multivariate probability density
function, f=f(
1
,a 
2
,. .,an), was used. The joint probability that a<A1P
a 2<A 2.. a <A at position x is given by:
A
I
A A
P(A
1
A
2  
,Anx)= { J J (c a 2 .cn)daida
2
" da (6-1)
Because the medium is assumed to be statistically homogeneous, f
does not depend on x. Since dispersion is supposed to be an
insignificant factor for spreading compared to the heterogeneity of the
advective term, transport in the vertical direction is governed by
v(z,t,a.) and 0(z,t,a.). The independent variables z and t are of a
1 1
deterministic nature and a. represents random hydraulic parameters that
1
characterize h(e) and K(O). It should be noted that the condition of no
interaction between different columns might be too stringent. In order
to describe the hydraulic conductivity in terms of a., the scaling
approach of Warrick et al. (210) was followed.
Because of the random nature of v and 0, the solute concentration
is also a random variable. The dependency in the horizontal x,y-plane
is also given in terms of f(a.). The probability that C<zA, in a plane
1
with depth z and at time t, is:
107
A
P(z,t,A)= { f(z,t,C) dA (6-2)
The concentration over the whole field can be described in terms of the
frequency function, f(al,X
2
,..,an), as well. The probability that
A<C(x,y,z,t)<A+dA is determined by the frequency function f(z,t,C)
where C varies between 0 and 1. This function is generally not known,
but can be characterized by its moments:
<C(z,t)> = C f(z,t,C) dC (6-3)
JO
2C (z,t)= (C-<C>) f(z,t,C) dC (6-4)
C J
0
where <C> is the average concentration over the entire field in a given
x,y-plane. The randomness of the flow term at the input boundary
(rainfall, irrigation) was taken into account by introducing a random
leaching rate, which is linked to the random behavior of the hydraulic
conductivity to give a probability distribution, P(v), for the pore
water velocity.
Because dispersion/diffusion effects are neglected, the
concentration profile obeys a step function (C(z,t)=O or 1 for z>vt and
z<vt, respectively). The concentration over the field is the average of
all individual column values:
<C(zt)>= 1 - P(z/t) = 1 - P(v) (6-5)
where P(v) is expressed in terms of variables x.. The field scale
dispersion, caused by random changes in 0 and v across the field,
exceeds the mechanical dispersion at the pore scale. If desirable, the
108
pore scale dispersion could be taken into account by introducing an
'effective' dispersion coefficient to make up for vertical
heterogeneity. Dagan and Bresler (49) argued that the large transition
zone, (0< <C> <1), which is commonly encountered under field
conditions, can now be described correctly. When using the
advection-dispersion equation an unrealistically high value for the
dispersion coefficient would be needed.
Bresler and Dagan (30) compared the solute spreading effect,. due
to heterogeneity of soil properties, with Taylor's (189) theory of
dispersion during flow in a capillary tube. The analogy is that
(transverse) diffusion initially cannot annihilate the longitudinal
concentration gradient caused by the advective term. However, for
larger times, diffusion becomes the dominant process. For conditions
normally encountered during infiltration, they concluded that only the
advection term is needed to describe the transport process, but that
for very large travel distances, e.g. 100 m, dispersion should be
included. The model derived by Dagan and Bresler (49) was used by
Bresler and Dagan (30) to predict solute concentrations for various
rates of infiltration using data of Warrick et al. (210). For a travel
distance of 1 m, an effective dispersivity between 16.1 and 334.2 cm
was found. If one had used the ADE, assuming Fickian dispersion, the
transition zone (the width of the solute front)would grow proportional
with t 
/
. However, Bresler and Dagan (30) found this zone to grow
linearly with t. It should be noted that the work by these authors did
not include an independent field study to validate the model by
comparing measured and calculated solute distributions.
109
Bresler and Dagan (31) added pore scale dispersion, as measured
in- the laboratory, or the much larger average field 
value for
dispersion, which reflects mainly heterogeneity of the profile, to the
previously mentioned mechanisms of spreading (horizontal heterogeneity
of' soil properties and a variable infiltration rate). They concluded
that the use of the conventional ADE, with constant coefficients,
poorly predicted <C(z,t)> in a heterogeneous field. To model transport,
one should focus on the variability of the hydraulic conductivity, the
average rate of recharge, and the value for the (average) field
dispersion. The influence of laboratory-measured pore scale dispersion
is negligible, and the variability of recharge on the soil surface,
provided that it is of modest value, had little influence on the
predicted solute profiles.
In a subsequent series of papers, a more general approach was
given for the stochastic modeling of unsteady flow and transport. Dagan
and Bresler (50) derived a model for infiltration and redistribution in
a heterogeneous field, which was applied to two spatially variable
soils (32). Bresler and Dagan (33) derived a model to predict solute
transport using the results of their work on water flow modeling.
Expected solute concentrations were calculated by using either an
accurate numerical scheme or a simplified model, or by representing the
variable field by an equivalent uniform field. The simplified model
110
hydrodynamic dispersion process, do not seem to have an appreciable
improvement upon prediction of statistical moments in spatially
variable fields. Using an equivalent uniform medium resulted in a much
larger error. In the simplified model, the solute concentration was
predicted with statistical moments rather than deterministic values.
The advantage of this procedure is that, although it might not
accurately predict concentrations at specific positions, the averaging
procedure resulted in reasonably accurate concentrations for an entire
field. This is often of more interest than specific values at given
locations.
Non-mechanist ic Models
A recent approach to model solute transport in the field ignores
the actual mechanism of flow and transport and merely focuses on the
solute concentration as a result of these unknown processes. In this
"black box" approach, the transport and flow processes are represented
by a transfer function, which characterizes systems whose internal
mechanisms are unknown or unknowable (90). This approach has been used
for some time by hydrologists and engineers (51).
Raats (147) applied the concept of transfer functions to describe
advective solute transport during steady flow in soils. He described
the movement of a certain water parcel in the course of time by
focusing on the travel time, usually of primary interest in
contamination studies. The position, x, that a parcel, , occupies at
time t can be given by:
x = x (,t)
(6-6)
111
with the velocity for 
X:
v = 
(6-7)
Denoting s as the 
arc-length along the 
streamline [L], Raats 
(147)
characterized v by 
its magnitude, v, and unit 
tangent vector, T:
V =a 
(6-8-a)
4 = vi/v 
(6-8-b)
The travel time between 
two points on the streamline 
was given by:
t-t = i ds 
(6-9)
s
0
where s and s are an arbitrary 
and a reference position, 
respectively.
0
For steady flow without 
sinks and sources, Raats 
(147) obtained the
following alternative 
expression for the 
travel time, making 
use of
Eq. (6-8-b):
s s
t-t -1 expT 
ds ds 
(6-10)
0 (Ov) [ j
wher (O) 0is the flux at s 0 Raats (147) noted that 
the divergence of
the unit tangent vector' field, 
V-T,. is a characteristic of 
the flow
pattern and a measure 
of the divergence or convergence 
of infinitesimal
stream tubes. For such 
a stream tube, Raats (147) 
obtained:
A(s) = Aexp {V- ds 
(6-11)
00
derived:
112
1 s
t-t (ev)A0 Ads 
(6-12)
0
The integral on the right-hand side corresponds to the total amount of
water in the stream tube between s and s, whereas (Sv) A is the
0 00
steady input flux. Eq. (6-12) gives the ratio of these quantities, 
being
the travel time between s and s . Having established a link between the
0
exact analytical expression for the travel time in terms of r, and an
expression for approximate graphical analyes, viz. Eq. (6-12), Raats
(147) proceeded by considering more specific flow problems.
After evaluating the travel time for one parcel, the movement of
collections of particles forming a surface was studied. This is useful
to examine the effect of solute application as a pulse or step front.
In order to predict whether the solute passes a certain point, e.g.,
the outlet of a soil column, the time needed to pass the column
(transit time) and the time the solute resides in the column (residence
time) need to be known. Therefore, Raats (147) used expressions for the
transit time, T, for each streamline and the residence time for a
parcel in the system. In order to study the effect of a solute
application, one can follow a collection of particles which enter the
system simultaneously, forming an isochrone.
Obviously, the frequency distribution of transit times in the
medium, consisting of various tubes, is of importance. It is also
important to determine the contribution of each tube to the total
transport. This can be done by monitoring the output concentration as a
function of time. The cumulative transit time distribution function, q,
113
is defined as the fraction of the streamtubes with transit times
smaller than T. At any time, a fraction q of the output will be
"younger" than z and a fraction (l-q) will be "older" 
than T. The
transit time density distribution function is equal to dT/dq, 
i.e., the
frequency distribution of the volume of the streamtubes. 
It takes into
account the rate of solute movement in a particular 
tube and the
contribution of that tube to the total solute transport. 
Finally, the
transfer function is the product of the input and 
transit time density
distribution function. Raats (148) applied the theory 
to a number of
flow patterns.
Jury (90) determined the solute concentration at various 
depths
in the soil by means of an output function, which 
was the result of a
solute flux applied at the soil surface, i.e., an input 
function. These
two functions were related by using a transfer function, 
which could be
determined by using measurements of the solute concentration 
at a
particular depth. The concentration at other depths could 
be predicted
with the transfer function. A brief derivation of the transfer 
function
follows.
Instead of using time as the independent variable, Jury (90) 
used
the cumulative amount of surface-applied water, I, as the independent
random variable for solute transport at a particular location (214).
The probability that a tracer injected at the surface will reach 
a
depth L after an amount of water I has been applied, is
L fL()') dI' (6-13)
0
114
I' is a dummy integration variable and fL(I) is the probability density
function, denoting that the injected tracer will arrive at depth 
L upon
application of an amount of water between I and I+dI, i.e., 
a
distribution of travel times. The inlet concentration is 
given by
C. =C =8(I), a narrow pulse injected at I=O at the soil surface. 
The
in o
average concentration, at z=L, for arbitrary variations in C. 
and
in
spatially uniform water input is given by:
CO
C(I) = 
C (I-I
) f(I') d' (6-14)
where f (I') is the probability of the solute reaching z=L between
L
'time' I' and I'+dI', i.e., the cumulative infiltration, and 
C. (I-I')
in
is the concentration discharging at depth L if for a 'breakthrough
time' I'. Note that this concentration remains equal to the inlet
concentration, because front spreading is ignored.
The rate of water input is not uniformly distributed over the
field, therefore a second probability density function, g(i), is 
needed
to describe the variabilty of input rate for various locations. The
relation between time and cumulative infiltration follows from i=dI/dt.
The probability that a unit area of soil will receive a water input
flux between i and i+di is ig(i). The probability that the solute
reaches a depth z=L between t and t+dt depends on the soil properties
(via fL (I)) and on the input rate (via g(i)). This probability is given
by the travel time density function, hLCt), which is the joint
probability function:
115
hL t) = Jig(i)fL(it) di 
(6-15)
%0
The average concentration at a given depth L, for 
a spatially uniform
input concentration C. t), follows from:
in
Cout Ct) = {C
in
(t-t' )hLt') dt' (6-16)
where t' is a dummy integration variable. 
Jury (90) assumed that the
distribution of physical processes contributing to 
the density
function, fL(I), between z=0 and z=L is equal for all 
depths. In that
case, the probability that a tracer will reach a given 
depth z, after a
cumulative infiltration 
1=I, is equal to 
the probability 
of reaching a
reference depth z=L after a cumulative infiltration 
of I=I1L/z:
IL/z
Pz M = PL(IL/z) = { fLI') dl' (6-17)
0
and by combining Eq.(6-15) and (6-16), the concentration 
for a
spatially variable water application is given by:
00 
00
C(znt) { C.(t-t) {L ig(i)fL(it') di dt' 
(6-18)
0 n 0
for a water application I=it. To obtain f a calibration 
at one depth
is necessary. For this purpose, Jury et al. (93) 
obtained soil solution
samples by vacuum extraction at a depth of 30 cm 
at 14 locations.
Usually a lognormal distribution seemed to fit 
fL(I)
, 
but other
functional re lat ionships 
could be used as well. 
Once fL( I) was
obtained, the performance 
of the model was tested 
by comparing
predicted values of C (z 
I) with measured values 
at other depths than
out
for which the model was 
calibrated. Good agreement 
was found between
116
results obtained by using the transfer 
function and experimental
values. It is not clear how well the model works for 
stratified soils.
The model also used for unsaturated conditions 
(92).
Concluding Remarks
A review of some modeling approaches was presented by Addiscott
and Wagenet (4), who concluded that among the many 
models available,
few have been exhaustively tested under field conditions. 
Tests by
others than the original authors, under different circumstances, 
have
been rare.
Most analytical and numerical models are deterministic 
and
mechanistic by nature. Because of the spatial variability 
of soil
properties in the field, stochastic models seem to be more 
attractive.
One such stochastic model, involving transport of a reactive solute,
was presented by van der Zee and van Riemsdijk (193). For fundamental
laboratory studies, the mechanistic deterministic model 
is the most
appropriate. Combination of the deterministic 
ADE with a stochastic
flow model could be a useful approach in 
future field studies.
Of course, many more efforts have been 
made to model water flow
for field conditions. Some involved scaling 
(180), others incorporated
geostatistics or the theory of stochastic 
processes (99, 100, 114, 208,
216, 217, 218). Charbeneau (38) applied 
the kinematic theory, i.e., a
functional relationship between flux and 
water content. A similar
theory was developed for solute transport. The kinematic 
waves can be
distinguished in self sharpening, during infilitrat ion, 
and self
spreading, during drainage.

118
Although much of the material already discussed includes
transport during unsaturated flow, some specific work appearin in the
soil science literature is discussed in this section. In particular,
some of the laboratory experiments are reviewed. A distinction is made
between steady state and transient flow conditions. A few approaches to
control the flow and/or water content for steady state conditions will
be mentioned as well.
Steady Flow
One of the first results on miscible displacement during
unsaturated flow was reported by Nielsen and Biggar (120). They
observed that the magnitude of water not readily. displaced, the hold
back, increased with desaturation. The inclusion of dispersion
phenomena in the transport equation seems more appropriate for
unsaturated than for saturated flow conditions. Biggar and Nielsen (20)
emphasized the role of, diffusion, pointing out that, although
mechanical dispersion will decrease at lower velocities, the process of
molecular diffusion tends to enhance spreading.
Krupp and Elrick (104) explained tailing during unsaturated
steady flow with the 'stagnant liquid concept'. According to these
authors, the dispersion coefficient is not related to 0 in a simple
manner. Variation in the sequences of fully filled and part ially filled
pores and transport in the surface films contribute to spreading.
119
the dispersion coefficient accurately in this 
case. Also, at very low
water contents, the variance in pore-water velocity 
is likely to vanish
since there is no pore sequence- for liquid flow. Thus, 
the amount of
spreading is reduced (104).
Yule and Gardner (219) studied both longitudinal 
and transverse
dispersion under unsaturated conditions. Transverse 
dispersion, i.e.,
diffusion, gains importance during unsaturated flow due 
to the lower
pore water velocity. The value for the transverse 
dispersion
coefficient was found to be nearly independent of pore 
water velocity.
Van Genuchten and Wierenga (202) noted that more tailing 
is to be
expected in unsaturated sorbing media. Not only will 
the relative
amount of stagnant water increase, the fraction of sorption 
sites
located in the stagnant region will increase as well. 
Among others,
they used the concept of a mobile and 'immobile region 
of the liquid
phase to formulate the transport equations for unsaturated 
soils. De
Smedt and Wierenga (57) observed early breakthrough and prolonged
tailing during unsaturated condit ions. Dispersion coefficients obtained
from saturated displacement experiments, with the same pore-water
velocity, could be used to fit the results from unsaturated experiments
if the transport equation accounted for mobile and immobile water. 
When
this distinction was not made, dispersion coefficients many 
times
larger than for saturated conditions were needed to 
fit the
120
qualitatively explained the asymmetry of the experimental concentration
profile. By varying the velocity and water content independently, using
different gravitational heads, Gupta et al. (78) found that D increased
with v and g enerally decreased with increasing 8. However, no
functional relationship was derived.
Awad (11) examined dispersion phenomena i n a medium fine sand. He
found that the hydrodynamic dispersion coefficient appeared to increase
with decreasing water content, but found no basis for the presence of a
mobile and stagnant region of the liqui d phase. A theoretical
investigation of the relationship between the hydrodynamic dispersion
coefficient and water content was pursued,. but not established.
Transient Flow
This section briefly reviews transport during infiltration, the
development of leaching strategies and the exploration of the
dependence of D on v. Wierenga (214) showed that transient flow could
be represented by an effective steady state flow to describe solute
transport. This simplifies considerably the task of obtaining a correct
flow term to predict solute movement considerably (38).
Kirda et al. (95, 96) determined the displacement of chloride
during infiltration with chloride-free water in soil columns, both
experimentally and theoretically. For a. given amount of cumulative
infiltration, the Cl salt was leached more efficiently (i.e., with less
121
the salt than larger, less frequent water applications. Warrick et al.
(209) observed that the water content at the soil surface determines
the propagation of the solute present in the irrigation water, while
the influence of the initial water content in the remainder of the soil
profile was shown to be small. For equal quantities of water
infiltrated, the depth of the maximum solute concentration increases
with decreasing surface water content.
Smiles et al. (179) investigated the absorption of a KCl solution
by initially rather dry soils. They observed that for a low
Peclet-number, the dispersion coefficient, D, was virtually
s
velocity-independent and that piston displacement of the solute
occurred. Smiles and Philip (178), using low initial salt and water
contents, noticed that the C1 concentration profile was somewhat
asymmetrical due to changes in the advection-term across the front.
However, the prediction of the solute front was relatively insensitive
to changes in D . Quantitatively, diffusion and dispersion are
s
identical. Smiles et al. (177) transformed the transport equation to
describe (transient) absorption by using a constant surface flux v .
0
2
The solution, expressed in terms of position, v0x, and time, v t, was
unique, suggesting that D is velOcity-independent. Smiles and Gardiner
5
(176) also observed close to piston displacement of the solute in a
clay soil. Thin films of water on the clay surface were not accessible
to the invading solute. The Cl front was located ahead of the piston
front, which was attributed to anion exclusion. Bond (25) analytically
solved transport of a solute, applied as a pulse, during unsteady flow.
122
He used a velocity dependent D. Experimentally and analytically
determined curves matched fairly well for 
3
H 0 but not for Cl. No
2
explanation could be offered for the poor prediction of the Cl profile.
The initial volumetric water content was 0.176.
Finally, Grismer (69) studied the absorption of water containing
Nal or SrCl2 by a silt loam. The dispersion coefficients, determined
according to Smiles et al. (179), depended on 0 in a way similar to the
diffusion coefficient. C omplete displacement of the initial water was
observed, although for the initial water content of 0.124, a higher
flux was required than for an initial water content of 0.035 to achieve
complete displacement.
Apparently, solute transport during steady flow differs from
transport during transient flow, at least for infiltration processes.
In the latter case, the dispersion coefficient is reported to be
virtually independent of the velocity and a complete displacement of
the resident solution occurs (piston displacement). Since the solute
acts as a tracer of relative water. movement, it is assumed that the
same holds for the flow process in general.
Experimental Considerations
A. major task when obtaining- breakthrough curves for steady,
unsaturated flow conditions is the independent control of water content
and pore water velocity. Because both factors influence dispersion, it
123
The driving force for water flow is a hydraulic head gradient,
VH, with H given by:
H = h + z (7-1)
where h is the pressure head [L] and z the gravitational head [L]. The
osmotic head is ignored in this expression. The flow of water is given
by Darcy's law:
JV = -K(O) VH (7-2)
where K(e) is the hydraulic conductivity [LT I. It follows from
Eq.(7-1) and (7-2) that the water flow can be manipulated directly by
adjusting h or z, or indirectly via the water content which influences
K and h.
Nielsen and Biggar (120) used hanging water columns and a
negative air pressure at the outlet to obtain a constant average water
content and flowrate. The glass bead medium required a change in head
of only 3 cm of water to obtain identical flowrates at full and
approximately half of the saturated water content. Fritted glass
plates, which have a negligible exchange capacity, a large capillary
conductivity, and a narrow pore size distribution, were used to
maintain unsaturated conditions.
Elrick et al. (58) and Krupp and Elrick (104) regulated the water
content by placing a perforated sample holder in a pressure chamber,
while the desired flowrate was obtained with a constant volume pump and
a head tube. Both ends of the sample were in close contact with a
cellulose acetate membrane filter, supported by a porous screen and a
124
stainless steel end cap containing three plastic nipples. 
The sample
and pressure chamber were placed on a balance to allow continuous
monitoring of the sample weight and hence water content. A driving
head, equal to the sample length, was established by increasing the air
pressure in the chamber(unit gradient). The gradient in pressure head
was equal to zero; a constant water content was obtained throughout the
sample. This meant that only one flowrate was possible for each water
content.
Yule and Gardner (219) stated that variations in the hydraulic
gradient can only be achieved by variations in water content. This is
not quite true, as will be discussed shortly, but it accurately points
out the problem. Their experimental setup consisted of a rectangular
plexiglasss soil container with an inflow and outflow control section.
Each section consisted of 11 one-bar porous ceramic tubes. At the
inlet, the inflow in each tube was controlled by a Mariotte bottle and
at the outlet, the outflow from each tube was collected in a plastic
vial maintained under suction controlled with a bubbling tower. A unit
gradient (gravitational head) was established to obtain a uniform
moisture content throughout the column.
Gupta et al. (77) used tilted soil columns to study dispersion
phenomena in an unsaturated glass bead medium. Solutions with different
KCl concentrations were applied with a pump. The water flow could be
125
Awad (11) used a two phase steady-state flow system to obtain
independent control of the volumetric water 
content and pore water
velocity. A horizontal system was used. Awad 
(11) proposed to use
changes in the pneumatic head to vary the 
flowrate. The pressure head
is given by:
h =h + h 
(7-3)
m a
where h is the matric pressure head, which 
governs the volumetric
m
water content, and h is the pneumatic head 
which does not influence
a
the water content (101). Hence, by manipulating 
ha) the flow can be
regulated (with a pump) without affecting the water 
content. A flow
cell was constructed with fritted glass plates. 
The pneumatic head was
regulated by maintaining a gradient in 
air pressure in the medium
between the plates. A constant flowrate of water 
was established with
a pump. Water content and bulk density 
were measured with a gamma
attenuation system. The hydraulic head was 
obtained by using
tensiometers connected to pressure transducers. 
The method could not be
used at very low water contents, due to the 
occurence of a high
resistance boundary layer between plates and medium.
VIII. MULTI-PHASE TRANSPORT
A logical extension of the previously discussed unsaturated
transport is the study of transport of immiscible fluids. Multiphase
transport involves vapor and liquid phase or various liquid phases.
There is a need to study these processes in order to predict movement
of (volatile) substances such as pesticides, oil products, and other
(organic) compounds. First, pesticide movement in liquid and vapor
phase will be discussed based on the work by Jury et al. (91), followed
by a brief mention of research in the area of immiscible liquid
transport.
Transport in Vapor and Liquid Phase
Jury et al. (91) determined the loss of pesticide via
volatilization and diffusion in the vapor phase and via advection and
diffusion in the liquid phase. Their investigation involved the
movement of triallate in soil samples placed in a volatilization
chamber at relative humidities of 50 and 100%. The vaporized triallate
was collected on polyurethane plugs.
The total concentration of the solute in the soil was given by:
C
t
= PbS +%0 Ce +O09C 
(8-1)
g g
where C is the total concentration expresses as mass of solute per
t
-3
volume of soil [ML ], S is the adsorbed concentration, expressed as
126
127
mass of solute per mass of soil [MM_], 
C and C9 are liquid and gas
concentration, respectively, expressed 
as mass of solute per volume of
fluid [ML -], and 0 and 0 are 
the liquid' and gas content,
3 -3
respectively, expressed as volume of fluid 
per volume of soil [L L
The solute flux was expressed as
ac 9 
ac,
JS -Dg ax -etax JVct(82
where Js is the solute flux, D9 is the gas 
diffusion coefficient
[L 2T ],p D is the liquid diffusion coefficient 
[L 2T_ l1and J is the,
t V
3 -2 -1
volumetric water flux [L L T I. It should be 
noted that advective
transport in the vapor phase was considered 
to be negligible. Under
some circumstances, -such as non-isothermal flow, this 
might not be
justified. The continuity equation for one-dimensional 
transport in the
absence of a sink/source term is:
at+ a S=0 
(8-3)
at ax
Eq. (8-1) and (8-2) can be substituted in Eq. 
(8-3) to describe mass
transport. In order to solve the resulting equation, it 
should contain
only one dependent variable,, viz, a concentration 
in a particular
phase. The following relationships are useful for this purpose:
1) Gas and liquid concentration are related by Henry's law
128
2) The adsorption isotherm is linear, which seems reasonable 
for the-
trace amounts of solute encountered:
q = aCt + b (8-5)
Some additional assumptions were made by Jury et al. (91) to
rewrite Eq. (8-3) in terms of C :
g
ac a 2 cac
S g D _ g+ V g (8-6)
at e 2 e ax
ax
where c = pbaKH+0K +0, De=Dg+KHD and ve=JvK For appropriate
b H t H g e g H e V H
boundary and initial conditions, Eq. (8-6) can be solved analytically
for C . Of course, due to hysteresis and non-equilibrium, Eq. (8-4) and
g
(8-5) might not be representative. This would subsequently influence
the accuracy of the solutions of Eq. (8-6).
Jury et al. (91) considered two cases, viz. transport by
diffusion only (Jv=O) and transport by advection-diffusion (Jv O),
which obviously led to different solutions for C . An expression 
for
g
the solute flux at the surface, J (O,t), was found by using 
Eq. (8-2),
5
in conjunction with the solution for C and Eq.(2-4). The cumulative
g
loss was then obtained via integration of the flux over time. 
The
experimentally measured cumulative loss was used to determine the
effective diffusion coefficient, D . Knowledge of this D enables one
e e
to make theoretical predictions of the mass flow. It was shown that
advection (evaporation) caused the vapor loss to increase slightly.
Loss of' the pesticide was mainly due to depletion from the upper soil
layer. Because of the increase in the concentration gradient, diffusion
becomes more important close to the soil surface.
129
Transport in Immiscible Liquid Phases
Recently, an increasing amount of research has been initiated in
the area of non-aqueous phase liquids, NAPL (2, 129). The situation is
slightly different for this type of flow than for flow in systems
containing a single liquid phase. Advective transport plays a role in
all phases, and knowledge of the flow properties for each phase is
needed. Multiple phase flow is, among others, discussed by Scheidegger
(165). Once the flow problem has been solved, an attempt can be made to
describe transport. The distribution coefficients of the solute between
the various phases, as well as in some cases the exchange and
adsorption parameters, need to be known. In some instances, the
immiscible fluids may consist of multiple components with distinct
coefficients describing the distribution between the various fluids
(44). Baehr and Corapcioglu (13) studied transport of the various
components of gasoline for which they predicted quite different travel
times. Numerical solutions for multiphase flow and 
transport can be
found in Huyakorn and Pinder (86).
Abriola and Pinder (2, 3) treated the simultaneous transport 
of
organic compounds in three phases: the non-aqueous, aqueous, 
and vapor
phases. These authors derived a system of three non-linear 
partial
differential equations with five independent variables 
(i.e., two
capillary pressures and three mass fractions). The capillary 
pressures
occur as a result of interfacial tension. The mass balance equation for
a particular species within a certain phase contains a storage term, an
advective and autonomous flux, a sink or source term, and exchange
130
resulting from a phase change or interphase diffusion. The mass balance
needs to be derived for all species in each phase. Often a number of
terms can be neglected, like the advective transport in the gasphase 
as
mentioned before. These authors did not consider the (ad)sorption of
the solute to the solid phase. An implicit numerical model was used to
solve the resulting system of non-linear equations with Newton-Raphson
iteration. The study was restricted to one-dimensional tranport. The
primary mechanism by which the volatile component was transported
through the medium was via gas diffusion. Since it was assumed that
local equilibrium existed between gas and liquid phases at all times,
the concentration in the liquid phases changed accordingly.
IX. EPILOGUE
It is hoped that this review demonstrated that the problem 
of
solute transport is a very broad one with many interesting 
aspects.
Combined with the urgency of ground water pollution, it is 
therefore
not surprising that a great deal of work has been done in 
this area.
Obviously, it was impossible to cover all aspects of transport modeling
and to discuss fundamental processes in great detail in a review like
this. Instead, it was restricted to a discussion of some basic
processes which affect transport. These processes are of a chemical,
physical, and biological nature, which explains why 
often an
interdisciplinary approach is needed (Chapter I). It should be noted
that, because of thisinterdisciplinary nature, the terminology and the
focus of the research vary widely. Some useful references were provided
in the discussions on the basic transport equation 
(Chapter II and
III). The reader is encouraged to consult 
the more advanced and
comprehensive discussions in the literature.
The principles of transport, i.e. using 
the principle of
conservation of mass, can be applied to transport 
under "complicating"
conditions: with non-linear exchange (Chapter IV), 
in structured and
layered media (Chapter V), during unsaturated 
flow (Chapter VII), and
for multi-phase systems (Chapter VIII). Although 
the fundamental
processes are essentially the same, these (and 
other) conditions offer
challenges for additional research. As our understanding of transport
131
132
progresses, new questions will undoubtedly arise. Because this report's
objective was to formulate transport problems, or at least refer to 
it,
relatively little attention was paid to the mathematical, in particular
the numerical, solution of transport problems Many good references
exist on the mathematical tools which are available to solve transport
problems.
The ultimate goal of much of the reasearch is geared towards a
better understanding of solute transport in natural porous media. Some
of the basic research might be helpful to solve "real problems," some
might not. Chapter VI pointed out that many of the traditional
deterministic models are of limited value. It is clear that modeling
should be accompanied by measurements, which makes such studies very
challenging for researchers because of the amount of effort involved
(time, money, expertise). Furthermore, some of the research findings
are not applicable under different conditions.
To put research efforts in the area of solute 
transport in
perspective, it is noted that many social and economical 
aspects are
involved as well. Some pollutants can safely be dispersed 
in large
water bodies, whereas others need to be contained. For 
some materials
(e.g., radioactive waste), one might argue that nonproduction 
and
nonuse is the only sensible approach. The short term economic 
benefits
of, for instance, agricultural chemicals, do not warrant 
their use in
view of possible adverse effects in the long term. Transport modelers
can provide part of the information needed to make a rational decision.
Finally, in light of all the difficulties associated with transport
133
modeling, it seems only 
logical to improve management 
practices in
dealing with da ngerous 
chemicals and to educate 
potential developers,
producers, and users. This 
requires input from researchers 
in the area
of solute transport as 
well.
NOTATION
This section contains 
the description and units 
of most of the symbols
used in this review. Because the 
notation in the literature is
not uniform, occasionally 
some unusual symbols 
were used or symbols
were used twice.
Symbol 
Description 
Units
a (aggregate) radius 
L
ak chemical activity 
of species k
B Brenner number (column 
Peclet number)
C concentration of solute in liquid 
phase ML
3
C dimensionless 
concentration
C(x,s) concentration in Laplace domain 
ML
-3
Ca local concentration inside aggregate 
ML-3
C dimensionless 
exit concentration
e
C exit concentration 
ML
-3
ex
Cf flux-averaged concentration ML
-3
f-3
C concentration 
of solute in gas phase 
ML
-3
C. initial concentration 
ML
-3
C. inlet concentration 
ML
in
-3
C concentration of 
solute in liquid phase (=C) 
ML
C mean concentration 
ML
m
CimCmo C for immobile 
and mobile region 
ML
-3
mo3
C molar concentration 
of species 
ML
C mole fraction 
of species k (C k/C
M
M 
Mk M-3
C concentration of 
feed solution (eluent) 
ML
-
C outlet concentration ML
3
out
C
t
total mass of solute per bulk 
volume ML
134
135
Symbol Description Units
-3
C
t  
time average of C ML
C
T  
total concentration in liquid 
phase ML
-3
C
V  
volume averaged concentration ML
d depth or characteristic particle size L
D,D coefficient of (hydrodynamic) dispersion L2T
-
1
D effective (retarded) value for D (=D/R) LT 1
D coefficient of mechanical dispersion LT
-
dis
D (effective) coefficient of molecular diffusion LT
-
D D for immobile region LT
eim e
D D for mobile region LT
emo e 
2 -1
D gas diffusion coefficient LT
g 2-1
De liquid diffusion coefficient LT
2 -1
D
L  
coefficient of longitudinal dispersion L2T
D. ,D D for immobile and mobile region L2T
imt mo 
2-1
D coefficient of molecular diffusion in the L2T
free 
liquid
D (O) dispersion coefficient (=OD(e)) LT
-
D
T  
coefficient of transverse dispersion LT
-
D (0) soil water diffusivity 
LT
-
f probability density function
f (I) probability density function for
movement in a soil as
f fraction of sorption sites in direct 
contact
with the mobile region
f' slope of exchange 
isotherm (dS/dC or 
dY/dX)
F fraction of adsorption sites belonging to
"type-i" (=S
1
/S)
F dimensionless rate parameter (k2L/v)
g acceleration due to gravity 
LT
g geometry factor
g(i) probability density function for infiltration
G Gibbs free energy J mol-
136
Symbol Description Units
h(O) soil water retention function
h pressure head or height L
h pneumatic head L
hL(t) travel time density function
h matric head L
m
H hydraulic head L
0~- 1
i infiltration rate LT
I cumulative infiltration L
JD diffusive/dispersive mass flux ML T
D-2
J. ,J cumulative diffusive flux in immobile and ML
mobile region
J total solute mass flux ML-2T
-
1
S-1
JV volumetric water flux (Darcy flux) LT
-1
k rate coefficient T
k permeability L
2
-1
k rate constant in the aggregate T
a
K(8) hydraulic conductivity LT
K equilibrium constant
K
d  
distribution coefficient LM
Kdi
m  
K
d 
for immobile region 
LM
K K for mobile 
region 
LM
.dmo d
K distribution ratio of ions between mobile and
D immobile region
KH partition coefficient in Henry's law
K exchange coefficient
jk
L column or mixing length L
L
D  
effective diffusion/dispersion length L
Ldi
f  
diffusion length L
L d dispersion length L
L dispersion length due to stagnant phase effect L
r-2
m
0  
total mass of solute per area of column ML
m p-th moment of temporal or spatial
distribution
137
Symbolj Description 
Units
m
p
n
N[ ]
N.
1
P
P
P or Pe
q
q'
qim
r
R
R
R. ,R
im, mo
R
a
s
S
S. ,S
im 
mo
ST
t
17
T
u
u
o
v
v
V
V
v
mo
v
m
w
ML-2T-1
_1 _)
ML T6
depth-averaged moment
number of moles
normal distribution function
diffusive molar flux
probability
pressure
Peclet number (vd/D, vL/D)
concentration of solute in adsorbed phase
differential capacity of exchanger (=dq/dC)
q for immobile sites
radial distance
gas constant
retardation factor
R for immobile and mobile region
aggregate radius
position along a streamline
concentration of solute in adsorbed phase on
a mass basis
S for immobile and mobile region
total solute concentration in adsorbed phase
time
temperature
number of pore volumes (vt/L)
velocity in a free liquid
maximum value for u
pore-water velocity
effective (retarded) velocity (v/R)
volume
volume per unit cross sectional area
pore-water velocity in mobile region (veT /e )
mean
width
L-3
ML
3
-3
L
J mol K
L
L
MM
- 3
MM
- 3
MM
3
mm-3
T
K
LT-1
LT
1
-1
LT
-1
LT
-1
LT
L
n aL3L
-2
excha-1
regiLT
I , I
138
Symbol I Description Units
x Cartesian coordinate in direction of flow L
X width L
xC position at which a particular concentration L
resides
X
k  
dimensionless concentration in liquid phase
for species 
k (=Ck/C
T
x equivalent depth of penetration of a solute L
front
y space coordinate 
L
Yk dimensionless concentration in adsorbed phase
for species k (=Sk/S
z space coordinate L
z gravitational head L
Z dimensionless distance (r/a, x/L)
a mass transfer coefficient T
a. random hydraulic parameter
1
aL longitudinal 
dispersivity 
L
a
T  
transversal dispersivity 
L
3 gravity segregation 
factor
activity coefficient
w -s rate constant 
for production in 
liquid and 
MLT- 1
W adsorbed phase 
T
6 Dirac delta function
Sviscosity 
ML-
1T
-
1
0 volumetric water content or volume fraction LL
- 3
3 -3
0 
im
0
moe for immobile and mobile region 
LL
- 3
e ,9 volume fraction in liquid and gas phase L3L
- 3
g gS-3
0 total volume fraction LL
A tortuosity factor
-1
pk chemical potential of species k J mol
o -1
pk k 
in a chosen 
standard state 
J mol
mean of spatial (i ) or temporal (It) L or T
distribution
139
Symbol Description 
Units
Awl s rate 
constant 
for first 
order decay 
in liquid 
T-1
and adsorbed phase
Vk valence of species k
Cartesian coordinate in a moving system 
L
P density 
ML-3
Pb dry 
bulk density 
ML-3
T standard deviation of spatial ((Y 
) or temporal L or T
(0-t) distribution x
0' relative spread in velocities
T (transit) time 
T
-T unit vector field (v/v)
formation factor
probability function
Boltzmann variable (x/ff) LT
I
BIBLIOGRAPHY
(1) ABRIOLA, L.M. 1987. Modeling Contaminant Transport in the
Subsurface: An Interdisciplinary Challenge. Rev. of Geophys.
25:125-134.
(2) and G.F. PINDER. 1985. A Multiphase Approach to the
Modeling of Porous Media Contamination by Organic Compounds. 1.
Equation Development. Water Resour. Res. 21:11-18.
(3) . 1985. A Multiphase Approach to the
Modeling of Porous Media Contamination by Organic Compounds. 2.
Numerical Simulation. Water Resour. Res. 21:19-26.
(4) ADDISCOTT, T.M. and R.J. WAGENET. 1985. Concepts of Solute
Leaching in Soils: A Review of Modelling Approaches. J. of Soil
Sci. 36:411-424.
(5) AMOOZEGAR-FARD, A., D.R. NIELSEN, and A.W. WARRICK. 1982. Soil
Solute Concentration Distributions for Spatially Varying Pore
Water Velocities and Apparent Diffusion Coefficients. Soil Sci.
Soc. Amer. J. 46:3-9.
(6) , A.W. WARRICK, and W.H. FULLER. 1983. A
Simplified Model for Solute Movement through Soils. Soil Sci. Soc.
Amer. J. 47:1047-1049.
(7) AMUNDSON, N.R. 1948. A Note on the Mathematics of Adsorption in
Beds. J. Phys. Chem. 52:1153-1157.
(8) . 1950. Mathematics of Adsorption in Beds. II. J.
Phys. Chem. 54:812-820.
(9) ARIS, R. 1956. On the Dispersion of a Solute in a Fluid Flowing
through a Tube. Proc. R. Soc. London Ser. A. 235:67-77.
(10) . 1958. On the Dispersion of Linear Kinematic Waves. Proc.
R. Soc. London Ser. A. 245:268-277.
(11) AWAD, M.M. 1986. Effect of Solution Content and Seepage Velocity
on Hydrodynamic Dispersion. Ph.D. Thesis. Colorado State Univ.
Fort Collins.
(12) BACHMAT, Y. 1969. Hydrodynamic Dispersion in a Saturated
Homogeneous Porous Medium at Low Peclet Numbers and Nonhomogeneous
Solution. Water Resour. Res. 5:139-143.
140
141
(13) BAEHR, A.L. and M.Y. CORAPCIOGLU. 1987. A Compositional Multiphase
Model for Groundwater Contaminants by Petroleum Products. 2.
Numerical solution. Water Resour. Res. 23:201-213.
(14) BEAR, J. 1969. Hydrodynamic Dispersion in Flow Through Porous
Media. In Flow through Porous Media. J. deWiest (Ed.) Academic
Press. New York.
(15) . 1975. Dynamics of Fluids in Porous Media. Elsevier. New
York.
(16) . 1979. Hydraulics of Groundwater. McGraw-Hill, New York.
(17) and Y. BACHMAT. 1967. A Generalized Theory on Hydrodynamic
Dispersion in Porous Media. Proc. Int. Assoc. Sci. Hydrol. Symp.
Haifa. Publ. 72.
(18) BEVEN, K. and P. GERMANN. 1982. Macropores and Water Flow in
Soils. Water Resour. Res. 18:1311-1325.
(19) BHATTACHARYA, R.N. and V.K. GUPTA. 1983. A Theoretical Explanation
of Solute Dispersion in Saturated Porous Media at the Darcy Scale.
Water Resour. Res. 19:938-944.
(20) BIGGAR, J.W. and D.R. NIELSEN. 1962. Miscible Displacement. II.
Behavior of Tracers. Soil Sci. Soc. Amer. Proc. 26:125-128.
(21) . 1964. Chloride-36 Diffusion During
Stable and Unstable Flow Through Glass Beads. Soil Sci. Soc. Amer.
Proc. 28:591-594.
(22) . 1976. Spatial Variability of the
Leaching Characteristics of a Field Soil. Water Resour. Res.
12:78-84.
(23) BOAST, C.W. 1973. Modeling the Movement of Chemicals in Soils by
Water. Soil Sci. 115:224-230.
(24) BOLT, G.H. 1982. Movement of Solutes in Soils: Principles of
Adsorption/Exchange Chromatography. In Soil Chemistry. B.
Physico-Chemical Models. G.H. Bolt (Ed.). Elsevier. Amsterdam.
(25) BOND, W.J. 1987. Solute Transport During Unsteady Unsaturated Soil
Water Flow: The Pulse Input. Aust. J. Soil Res. 25: 233-241.
(26) BREDEHOEFT, J.D. 1971. Numerical Dispersion Comments. Water
Resour. Res. 7:755-756.
(27) BRENNER, H. 1962. The Diffusion Model of Longitudinal Mixing in
Beds of Finite Length. Numerical Values. Chem. Eng. Sci.
17: 229-243.
142
(28) BRESLER, E. 1973. Simultaneous Transport 
of Solutes and Water
Under Transient Unsaturated Flow 
Conditions. Water Resour. Res.
9:975-986.
(29) . 1973. Anion 
Exclusion and Coupling Effects 
in
Nonsteady Transport through Unsaturated 
Soil. I. Theory. Soil Sci.
Soc. Amer. Proc. 37:663-669.
(30) and G. DAGAN. 1979. 
Solute Dispersion in Unsaturated
Heterogeneous Soil at Field Scale. 
II. Application. Soil Sci. Soc.
Amer. J. 43: 467-472.
(31) . 1981. Convective 
and Pore Scale
Dispersive Solute Transport in Unsaturated 
Heterogeneous Fields.
Water Resour. Res. 17: 1683-1693.
(32) . 1983. 
Unsaturated Flow in Spatially
Variable Fields. 2. Application of 
Water Flow Models to Various
Fields. Water Resour. Res. 19:421-428.
(33) . 1983. Unsaturated 
Flow in Spatially
Variable Fields. 3. Solute Transport 
Models and Their Applications
to Two Fields. Water Resour. Res. 19:429-435.
(34) BRIGHAM, W.E. 1974. Mixing 
Equations in Short Laboratory Cores.
Soc. Pet. Eng. J. 14:91-99.
(35) BRUCE, R.R. and A. KLUTE. 
1956. The Measurement of Soil Moisture
Diffusivity. Soil Sci. Soc. Amer. 
Proc. 20:458-462.
(36) CAMERON, D.R. and A. KLUTE. 1977. 
Convective-Dispersive Solute
Transport with a Combined Equilibrium 
and Kinetic Adsorption
Model. Water Resour. Res. 13:183-188.
(37) CARSLAW, H.S. and J.C. JAEGER. 
1959. Conduction of Heat in
Solids. Clarendon Press, Oxford.
(38) CHARBENEAU, R.J. 1984. Kinematic 
Models for Soil Moisture and
Solute Transport. Water Resour. Res. 
20:699-706.
(39) CHAUDHARI, N.M. 1971. An Improved 
Numerical Technique for Solving
Multidimensional Miscible Displacement 
Equations. Soc. of Petrol.
Eng. J. 11:277-284.
(40) CHILDS, S.W. and R.J. HANKS. 1975. Model 
of Salinity Effects on
Crop Growth. Soil Sci. Soc. Amer. Proc. 39:617-622.
(41) CHO, C.M. 1971. Convective Tranport 
of Ammonium with Nitrification
in Soil. Can J. Soil Sci. 51:339-350.
143
(42) CLEARY, R.W. and D.D. ADRIAN. 1973. Analytical Solution of the
Convective-Dispersive Equation for Cation Adsorption in Soils.
Soil Sci. Soc. Amer. Proc. 37:197-199.
(43) COATS, K.H. and B.D. SMITH. 1964. Dead End Pore Volume 
and
Dispersion in Porous Media. Soc. Petrol. Eng. J. 4:73-84.
(44) CORAPCIOGLU, M.Y and A.L. BAEHR. 1987. A Compositional Multiphase
Model for Groundwater Contaminants by Petroleum Products. 
1.
Theoretical Considerations. Water Resour. Res. 
23:191-200.
(45) CRANK, J. 1975. The Mathemat'ics of Diffusion. Clarendon 
Press,
Oxford.
(46) CREMERS, A. and H. LAUDELOUT. 1965. On the 'Isoconductivity 
Value
of Clay Gels. Soil Sci. 100: 298-299.
(47) CUSHMAN, J.H. 1982. Nutrient Transport Inside and Outside 
the Root
Rhizosphere. Theory. Soil Sci. Soc. Amer. J. 46:704-709.
(48) DAGAN, G. 1984. Solute Transport in Heterogeneous 
Porous
Formations. J. Fluid Mech. 145:151-177.
(49) and E. BRESLER. 1979. Solute Dispersion 
in Unsaturated
Heterogeneous Soil at Field Scale. I. Theory. 
Soil Sci. Soc. Amer.
J. 43:461-467.
(50) . 1983. Unsaturated Flow in Spatially
Variable Fields. 1. Derivation of Models 
of Infiltration and
Redistribution. Water Resour. Res. 19:413-420.
(51) DANCKWERTS, P.V. 1953. Continuous Flow Systems. Chem. Eng. 
Sci.
2:1-13.
(52) DANE, J.H. and A. KLUTE. 1977. Salt Effects on 
the Hydraulic
Properties of a Swelling Soil. Soil Sci. Soc. 
Amer. J.
41:1043-1049.
(53) DAVIDSON, J.M., P.S.C. RAO, and P. NKEDI-KIZZA. 
1983. Physical
Processes Influencing Water and Solute Transport in Soils. 
In
Chemical Mobility and Reactivity in Soil Systems. Soil Sci. Soc.
Amer., Spec. Publ. 11. Madison, Wis.
(54) DE JOSSELIN DE JONG, G. 1958. Longitudinal 
and Transverse
Diffusion in Granular Deposits. Trans. Amer. Geophys. Union.
39:67-74.
(55) DE SMEDT, F. and P.J. WIERENGA. 1979. Mass Transfer in Porous
Media with Immobile Water. J. of Hydr. 41:59-67.
144
(56) . 1979. A Generalized Solution 
for
Solute Flow in Soils with Mobile and Immobile Water. Water Resour.
Res. 15: 1137-1141.
(57) . 1984. Solute Transfer 
through
Columns of Glass Beads. Water Resour. Res. 20: 225-232.
(58) ELRICK, D.E., K.T. ERH, and H.K. KRUPP. 1966. 
Applications of
Miscible Displacement Techniques to Soils. Water Resour. Res.
2:717-727.
(59) , K.B. LARYEA, and P.H. GROENEVELT. 1979. Hydrodynamic
Dispersion During Infiltration of Water into Soil. Soil Sci. Soc.
Amer. J. 43:856-865.
(60) FISCHER, H.B., E. LIST, R.C.Y. KOH, J. IMBERGER, and N.H. BROOKS.
1979. Mixing in Inland and Coastal Waters. Academic Press.
(61) FRIED, J.J. and M.A. COMBARNOUS. 1971. Dispersion in Porous
Media. In Adv. Hydrosci. 9:169-282. Academic, New York.
(62) FRISSEL, M.J., P. POELSTRA, and P. REINIGER. 1970. Chromatographic
Transport Through soils. III. A Simulation Model for the
Evaluation of the Apparent Diffusion Coefficient in Undisturbed
Soils with Titrated Water. Plant and Soil 33:161-176.
(63) CAPON, E.N. 1933. On the Theory of Exchange Adsorption in Soil. J.
Gen. Chem. U.S.S.R. 3:144-152.
(64) GELHAR, L.W., A.L. GUTJAHR, and R.L. NAFF. 1979. Stochastic
Analysis of Macrodispersion in a Stratified Aquifer. Water 
Resour.
Res. 15:1387-1396.
(65) GERSHON, N.D. and A. NIR. 1969. Effects of Boundary Conditions of
Models on Tracer Distributions in Flow through Porous Mediums.
Water Resour. Res. 5:830-839.
(66) GIDDINGS, J.C. 1965. Dynamics of Chromatography. Part I.
Principles and Theory. Marcel Dekker. New York.
(67) GILLHAM, R.W., E.A. SUDICKY, J.A. CHERRY, and E.O. FRIND. 1984. An
Advection-Diffusion Concept for Solute Transport in Heterogeneous
Unconsolidated Geological Deposits. Water Resour. Res. 20:369-378.
(68) GRISAK, G.E. and J.F. PICKENS. 1980. Solute Transport 
through
Fractured Media. 1. The Effect of Matrix Diffusion. Water 
Resour.
Res. 16:719-730.
(69) GRISMER, M.E. 1986. Nondestructive Observations 
of Solution
Displacement in Soils. Soil Sci. 141:185-189.
145
(70) , D.B. MCWHORTER, and A. KLUTE. 1986. Monitoring
Water and Salt Movement in Soils at Low Solution Contents. Soil
Sci. 141:163-171.
(71) GROENEVELT, P.H. and G.H. BOLT. 1969. Non-Equilibrium
Thermodynamics of the Soil-Water System. J. Hydrol. 7:358-388.
(72) , D.E. ELRICK, and K.B. LARYEA. 1980. Coupling
Phenomena in Saturated Homo-Ionic Montmorillonite. IV. The
Dispersion Coefficient. Soil Sci. Soc. Amer. J. 44:1168-1173.
(73) GUPTA, V.K. and R.N. BHATTACHARYA. 1983. A New Derivation in the
Taylor-Aris Theory of Solute Dispersion in a Capillary. Water
Resour. Res. 19:945-951.
(74) . 1986. Solute Dispersion in
Multidimensional Periodic Saturated Porous Media. Water Resour.
Res. 22:156-164.
(75) GUPTA, S.P. and R.A. GREENKORN. 1974, Determination of Dispersion
and Nonlinear Adsorption Parameters for Flow in Porous Media.
Water Resour. Res. 10:839-846.
(76) GUPTA, R.K., A. KLUTE, and R.J. MILLINGTON. 1972. An Electrode
System for Transient Dispersion Studies in Unsaturated Glass Bead
Mediums. J. Indian Soc. Soil Sci. 20: 297-304.
(77) , R.J. MILLINGTON, and A. KLUTE. 1973. Hydrodynamic
Dispersion in Unsaturated Porous Media. I. Concentration
Distribution during Dispersion. J. Indian Soc. Soil Sci. 21:1-7.
(78) . 1973. Hydrodynamic
Dispersion in Unsaturated Porous Media. II. The Stagnant Zone
Concept and the Dispersion Co-efficient. J. Indian Soc. Soil Sci.
21:121-125.
(79) COVEN, O., R.W. FALTA, F.J. MOLZ, and J.G. MELVILLE. 1985.
Analysis and Interpretation of Single Well Tracer Tests in
Stratified Aquifers. Water Resour. Res. 21:676-684.
(80) and F.J. Molz. 1986. Deterministic and Stochastic
Analyses of Dispersion in an Unbounded Stratified Porous Medium.
Water Resour. Res. 22:1566-1574.
(81) , and J.G. MELVILLE. 1984. An Analysis of
Dispersion in a Stratified Aquifer. Water Resour. Res.
20: 1337-1354.
(82) HARLEMAN, D.R.F. and R.R. RUMER. 1963. Longitudinal and Lateral
Dispersion in an Isotropic Porous Medium. J. Fluid Mech.
16: 385-394.
146
(83) HASHIMITO, I., DESHPANDE, K.B., and H.C. THOMAS. 1964. Peclet
Numbers and Retardation Factors for Ion Exchange Columns. Ind.
Eng. Chem. Fund. 3:213-218.
(84) HELFFERICH, F. 1962. Ion Exchange. McGraw-Hill. New York.
(85) HOPMANS, J.W. and J.H. DANE. 1985. Effect of Temperature-
dependent Hydraulic Properties on Soil Water Movement. Soil Sci.
Soc. Amer. J. 49:51-58.
(86) HUYAKORN, P.S. and G.F. PINDER. 1983. Computational Methods In
Subsurface Flow. Academic Press. New York.
(87) JAMES, R.V. and J. RUBIN. 1972. Accounting for Apparatus-induced
Dispersion in Analyses of Miscible Displacement Experiments. Water
Resour. Res. 8:717-721.
(88) JENNINGS, A.A., D.J. KIRKNER, and T.L. THEIS. 1982. Multicomponent
Equilibrium Chemistry in Groundwater Quality Models. Water Resour.
Res. 18:1089-1096.
(89) JOURNEL, A.G. and CH.J. HUIJBREGTS. 1978. Mining Geostatistics.
Academic Press. London.
(90) JURY, W.A. 1982. Simulation of Solute Transport Using a Transfer
Function Model. Water Resour. Res. 18:363-368.
(91) , R. GROVER, W.F. SPENCER, and W.J. FARMER. 1980.
Modeling Vapor Losses of Soil-incorporated Triallate. Soil Sci.
Soc. Amer. J. 44:445-450.
(92) and G. SPOSITO. 1985. Field Calibration and Validation
of Solute Transport Models for the Unsaturated Zone. Soil Sci.
Soc. Amer. J. 49:1331-1341.
(93) , L.H. Stolzy, and P. Shouse. 1982. A Field Test of the
Transfer Function Model for Predicting Solute Transport. Water
Resour. Res. 18:369-375.
(94) KAY, B.D. and D.E. ELRICK. 1967. Adsorption and Movement of
Lindane in Soils. Soil Sci. 104:314-322.
(95) KIRDA, C., D.R. NIELSEN, and J.W. BIGGAR. 1973. Simultaneous
Transport of Chloride and Water During Infiltration. Soil Sci.
Soc. Amer. Proc. 37:339-34S.
(96) . 1974. The Combined
Effects of Infiltration and Redistribution of Leaching. Soil 
Sci.
117: 323-329.
147
(97) KIRKNER, D.J., A.A. JENNINGS, and T.L. THEIS. 1985. Multlsolute
Mass Transport with Chemical Interaction Kinetics. J. Hydrol.
76:107-117.
(98) KLINKENBERG, L.J. 1957. Pore Size Distribution of Porous Media
and Displacement Experiments with Miscible Liquids. J. Petrol.
Techn. 210:63-66.
(99) KNIGHTON, R.E. and R.J. WAGENET. 1987. Simulation of Solute
Transport Using a Continuous Time Markow Process. 1. Theory and
Steady State Application. Water Resour.Res. 23:1911-1916,
(100) . Simulation of Solute
Transport Using a Continuous Time Markow Process. 2. Application
to Transient Field Conditions. Water Resour. Res. 23:1917-1925.
(101) KOOREVAAR, P., G. MENELIK, and C. DIRKSEN. 1983. Elements of Soil
Physics. Elsevier. Amsterdam.
(102) KREFT, A. and A. ZUBER. 1978. On the Physical Meaning of the
Dispersion Equations and its Solutions for Different Initial and
Boundary Conditions. Chem. Eng. Sci. 33:1471-1480.
(103) KRUPP, H.K., J.W. BIGGAR, and D.R. NIELSEN. 1972. Relative Flow
Rates of Salt and Water in Soil. Soil Sci. Soc. Amer. Proc.
36:412-417.
(104) and D.E. ELRICK. 1968. Miscible Displacement in an
Unsaturated Glass Bead Medium. Water Resour. Res. 4:809-815.
(105) . 1969. Density Effects in Miscible
Displacement Experiments. Soil Sci. 107:372-380.
(106) LAI, S.H. and J.J. JURINAK. 1971. Numerical Approximation of
Cation Exchange in Miscible Displacement through Soil Columns.
Soil Sci. Soc. Amer. Proc. 35:894-899.
(107) . 1972. Cation Adsorption in
One-dimensional Flow through Soils: A Numerical Solution. Water
Resour. Res. 8:99-107.
(108) LAKE, L.W. and F. HELFFERICH. 1978. Cation Exchange in Chemical
Flooding. 2. The Effects of Dispersion, Cation Exchange, and
Polymer/Surfactant Adsorption on Chemical Flood Environment. Soc.
Pet. Eng. J.:435-444.
(109) LAPIDUS, L. and N.R. AMUNDSON. 1952. Mathematics of Adsorption in
beds. VI. The Effect of Longitudinal Diffusion in Ion Exchange
and Chromatographic Columns. J. Phys. Chem. 56:984-988.
148
(110) LARYEA, K.B., D.E. 
ELRICK, and M.J.L. ROBIN. 
1982. Hydrodynamic
Dispersion Involving 
Cationic Adsorption during 
Unsaturated
Transient Water Flow in Soil. 
Soil Sci. Soc. Amer. Proc.
46:667-671.
(111) LEVENSPIEL, 0. and W.K. 
SMITH. 1957. Notes on the 
Diffusion-type
Model for the Longitudinal Mixing 
of Fluids in Flow. Chem. Eng.
Sci. 6:227-233.
(112) LINDSTROM, F.T. 
1976. Pulse Dispersion of 
Trace Chemical
Concentrations in a Saturated Sorbing 
Porous Medium. Water
Resour. Res. 12:229-238.
(113) , R. HAQUE, 
V.H. FREED, and L. BOERSMA. 
1967.
Theory on the Movement of Some 
Herbicides in Soils. Environ. 
Sci.
Techn. 1:561-565.
(114) LUXMOORE, R.J. and M.L. 
SHARMA. 1980. Runoff Response to 
Soil
Heterogeneity: Experimental 
and Simulation Comparisons for 
Two
Contrasting Watersheds. Water 
Resour. Res. 16:675-684.
(115) MARLE, C., P. SIMANDOUX, J. 
PACSIRSKY, and C. GAULIER. 1967.
Etude du D6placement de Fluides Miscibles 
en Milieu Poreux
Stratifie. Rev. Inst. Fr. Petrol. 22: 
272-294.
(116) MATHERON, G. and G. DE MARSILY. 
1980. Is Transport in Porous
Media Always Diffusive? A Counterexample. 
Water Resour. Res.
16:901-917.
(117) MERCADO, A. 1967. The Spreading 
Pattern of Injected Water in a
Permeability Stratified Aquifer. 
Int. Ass. Sci. Hydrol. Symp.
Haifa. Publ. 72.
(118) MILLINGTON, R.J. and 
R.C. SHEARER. 1971. Diffusion 
in Aggregated
Porous Media. Soil Sci. 111: 
372-378.
(119) NERETNIEKS, I. 1983. 
A Note on Fracture Flow 
Dispersion
Mechanisms in the Ground. Water 
Resour. Res. 19:364-370.
(120) NIELSEN, D.R. and J.W. 
BIGGAR. 1961. Miscible Displacement 
in
Soils. I. Experimental Information. 
Soil Sci. Soc. Amer. Proc.
25:1-5.
(121) 
. 1962. Miscible Displacement. 
III.
Theoretical Considerations. 
Soil Sci. Soc. Amer. 
Proc. 26:
216-221.
(122) - - , R.D. 
JACKSON, J.W. CARY, and 
D.D. EVANS (Eds.).
1972. Soil Water. Amer. 
Soc. of Agron. Madison, 
Wis.
149
(123) 
, M.Th. VAN GENUCHTEN, 
and J.W. BIGGAR. 
1986. Water
Flow and Solute Transport 
Processes in the 
Unsaturated Zone.
Water Resour. Res. 
22:89S-108S.
(124) NKEDI-KIZZA, P., J.W. 
BIGGAR, M.Th. VAN GENUCHTEN, 
P.J. WIERENGA,
H.M. SELIM, J.M. DAVIDSON, 
and D.R. NIELSEN. 
1983. Modeling
Tritium and Chloride 
36 Transport through 
an Aggregated Oxisol.
Water Resour. Res. 
19:691-700.
(125) NOVOZAMSKY, 
I., J. BEEK, and 
G.H. BOLT. 1978. 
Chemical
Equilibria. In Soil Chemistry. 
A. Basic Elements. G.H. 
Bolt and
M.G.M. Bruggenwert 
(Eds.). Elsevier. Amsterdam.
(126) NYE, P.H. 1979. 
Diffusion of Ions and Uncharged 
Solutes in Soils
and Soil Clays. Adv. 
Agron. 31: 225-272.
(127) OGATA, A. and R.B. 
BANKS. 1961. Solution 
of the Differential
Equation of Longitudinal 
Dispersion in Porous Media. 
U.S. Geol.
Surv. Prof. Paper 
411-A:1-7.
(128) PANDEY, R.N. and 
S.K. GUPTA. 1980. Comment 
on 'Dispersion in Soil
Columns: Effect of 
Boundary Conditions 
and Irreversible
Reactions'. Soil Sci. 
Soc. Amer. J. 44: 661.
(129) PARKER, J.C., R.J. 
LENHARD, and T. KUPPUSAMY. 
1987. A Parametric
Model for Constitutive Properties 
Governing Multiphase Flow 
in
Porous Media. Water. Resour. 
Res. 23:618-624.
(130) and M.TH. 
VAN GENUCHTEN. 1984. Determining 
Transport
Parameters from Laboratory 
and Field Tracer Experiments. 
Virg.
Agric. Exp. Stat. 
Bull. 84-3.
(131) 
. 1984. Flux-averaged
Concentrations in Continuum 
Approaches to Solute 
Transport. Water
Resour. Res. 20:866-872.
(132) PARLANGE, J.-Y., 
D.A. BARRY, and J.L. 
STARR. 1985. Comments 
on
"Boundary Conditions 
for Displacement Experiments 
through Short
Laboratory Columns". 
Soil Sci. Soc. Amer. 
J. 49:1325.
(133) 
and J.L. STARR. 1975. 
Linear Dispersion in 
Finite
Columns. Soil Sci. Soc. 
Amer. Proc. 39:817-819.
(134) 
. 1978. Dispersion in 
Soil Columns:
Effect of the Boundary 
Conditions and 
Irreversible Reactions.
Soil Sci. Soc. Amer. 
J. 42:15-18.
(135) PASSIOURA, 
J.B. 1971. Hydrodynamic 
Dispersion in Aggregated
Media. 1.Theory. Soil 
Sci. 111:339-344.
150
(136) and D.A. ROSE. 
1971. Hydrodynamic Dispersion in
Aggregated media. 2. Effects 
of Velocity and Aggregated Size.
Soil Sci. 111:345-351.
(137) PERSAUD, N., J.V. 
GIRALDEZ, and A.C. CHANG. 
1985. Monte-Carlo
Simulation of Noninteracting 
Solute Transport in a Spatially
Heterogeneous Soil. Soil Sci. 
Soc. Amer. J. 49:562-568.
(138) PFANNKUCH, H.O. 1963. 
Contribution a 1'Etude des D6placements 
de
Fluides dans un Milieu poreux. Rev. 
Inst. Fr. Petr. 18:215-270.
(139) PHILIP, J.R. 1957. Numerical 
Solution of Equations of the
Diffusion Type with Diffusivity 
Concentration-dependent. II.
Aust. J. Phys. 10: 29-42.
(140) . 1968. The Theory 
of Absorption in Aggregated Media.
Aust. J. Soil Res. 6:1-19.
(141) . 1968. Diffusion, Dead-end-pores 
and Linearized
Absorption in Aggregated Media. Aust. J. Soil 
Res. 6:21-30.
(142) PICKENS, J.F. and G.E. GRISAK. 
1981. Scale-dependent Dispersion
in a Stratified Granular Aquifer. 
Water Resour. Res.
17:1191-1211.
(143) .
1981. Modeling of Scale-dependent
Dispersion in Hydrogeologic Systems. Water 
Resour. Res.
17:1701-1711.
(144) PINDER, G.F. and W.G. GRAY. 1977. Finite 
Elements in Surface and
Subsurface Hydrology. Academic 
Press. New York.
(145) PLUMB, O.A. and S. WHITAKER. 
1988. Dispersion in Heterogeneous
Porous Media. 1. Local Volume 
Averaging and Large-scale
Averaging. Water Resour. Res. 
24:913-924.
(146) 
. 1988. Dispersion in Heterogeneous
Porous Media. 2. Predictions for Stratified 
and Two-dimensional
Spatially Periodic Systems. Water 
Resour. Res. 24:927-938.
(147) RAATS, P.A.C. 1978. Convective 
Transport of Solutes by Steady
Flows. 1.General Theory. Agr. Water 
Manage. 1:201-218.
(148) . 1978. Convective 
Transport of Solutes by Steady
Flows. 2. Specific Flow Problems. 
Agr. Water Manage. 1:219-232.
(149) RAO, P.S.C., J.M. DAVIDSON, 
R.E. JESSUP, and H.M. 
SELIM. 1979.
Evaluation of Conceptual 
Models for Describing 
Non-equilibrium
Adsorption-Desorption of 
Pesticides during Steady 
Flow in Soils.
Soil Sci. Soc. Amer. J. 
43:22-28.
151
(150) , R.E. JESSUP, D.E. ROLSTON, J.M. DAVIDSON, and D.P.
KILCREASE. 1980. Experimental and Mathematical Description of
Non-adsorbed Solute Transfer by Diffusion in Spherical
Aggregates. Soil Sci. Soc. Amer. J.44:684-688.
(151) , D.E. ROLSTON, R.E. JESSUP, and J.M. DAVIDSON. 1980.
Solute Transport in Aggregated Porous Media. Theoretical and
Experimental Evaluation. Soil Sci. Soc. Amer. J. 44:1139-1146.
(152) RASMUSON, A. 1981. Diffusion and Sorption in Particles and
Two-dimensional Dispersion in a Porous Medium. Water Resour. Res.
17:321-328.
(153) and I. NERETNIEKS. 1980. Exact Solution of a Model
for Diffusion in Particles and Longitudinal Dispersion in Packed
Beds. AIChE Journal 26:686-690.
(154) . 1981. Migration of Radionuclides
in Fissured Rock: The Influence of Micropore Diffusion and
Longitudinal Dispersion. J. of Geophys. Res. 86:3749-3758.
(155) RAZAVI, M.-S., B.J. MCCOY, and R.G. CARBONELL. 1978. Moment
Theory of Breakthrough Curves for Fixed-bed Adsorbers and
Reactors. Chem. Eng. J. 16:211-222.
(156) REINIGER, P. and G.H. BOLT. 1972. Theory of Chromatography and
its Application to Cation Exchange in Soils. Neth. J. Agric. Sci.
20:301-313.
(157) RISKEN, H. 1984. The Fokker-Planck Equation. Methods of Solution
and Application. Springer-Verlag. New York.
(158) ROBBINS, C.W., J.J. JURINAK, and R.J. WAGENET. 1980. Calculating
Cation Exchange in a Salt Transport Model. Soil Sci. Soc. Amer.
J. 44: 1195-1199.
(159) , R.J. WAGENET, and J.J. JURINAK. 1980. A Combined
Salt Transport-chemical Equilibrium Model for Calcareous and
Gypsiferous Soils. Soil Sci. Soc. Amer. J. 44:1191-1194.
(160) ROSE, D.A. and J.B. PASSIOURA. 1971. The Analysis of Experiments
on Hydrodynamic Dispersion. Soil Sci. 111:252-257.
(161) . 1971. Gravity Segregation during
Miscible Displacement. Soil Sci. 111:258-265.
(162) RUBIN, J. 1983. Transport of Reacting Solutes in Porous Media:
Relation between Mathematical Nature of Problem Formulation and
Chemical Nature of Reactions. Water Resour. Res. 19:1231-1252.
152
(163) and R.V. JAMES. 1973. 
Dispersion-affected Transport of
Reacting Solutes in Saturated Porous 
Media: Galerkin Method
Applied to Equilibrium-controlled Exchange in Unidirectional
Steady Water Flow. Water Resour. Res. 9:1332-1356.
(164) SAFFMAN, P.G. 1959. A Theory of Dispersion 
in a Porous Medium. J.
Fluid Mech. 6:321-349.
(165) SCHEIDEGGER, A.E. 1974. The Physics of 
Flow Through Porous Media.
University of Toronto Press.
(166) SCOTTER, D.R. 1978. Preferential Solute Movement 
through Larger
Soil Voids. I. Some Computations Using Simple Theory. Aust. J.
Soil Res. 16:257-267.
(167) and P.A.C. RAATS. 1968. Dispersion in 
Porous
Mediums Due to Oscillating Flow. Water Resour. Res. 6:1201-1206.
(168) SELIM, H.M., J.M. DAVIDSON, and R.S. MANSELL. 1976. Evaluation 
of
a Two-site Adsorption-Desorption Model for Describing Solute
Transport in Soil. Proc. Summer Conf., Washington D.C. pp
444-448.
(169) , R.S. MANSELL, and A. ELZEFTWAY. 1976 . Distribution
of 2,4-D and Water in Soil during Infiltration and
Redistribution. Soil Sci. 121:176-183.
(170) , R. SCHULIN, and H. FLUHLER. 1987. Transport 
and Ion
Exchange in an Aggregated Soil. Soil Sci. Soc. Amer. J.
51:876-884.
(171) SHAMIR, U.Y. and D.R.F. HARLEMAN. 1967. 
Dispersion in Layered
Porous Media. J. Hydraul. Div. ASCE. 93(HY5):237-260.
(172) SIMMONS, C.S. 1982. A Stochastic-Convective 
Transport
Representation of Dispersion in One-dimensional Porous 
Media
Systems. Water Resour. Res. 18:1193-1214.
(173) SKAGIUS, K. and I. NERETNIEKS. 1986. Porosities 
and Diffusivities
of Some Nonsorbing Species in Crystalline Rocks. 
Water Resour.
Res. 22: 389-398.
(174) SKOPP, J., W.R. GARDNER, and E.J. Tyler. 
1981. Solute Movement in
Structured Soils: Two-region Model with Small Interaction. 
Soil
Sci. Soc. Amer. J. 45:837-842.
(175) SKOPP, J. and A.W. WARRICK. 
1974. A Two-phase Model for the
Miscible Displacement of Reactive 
Solutes in Soil. Soil Sci. Soc.
Amer. J. 138:545-550.
153
(176) SMILES, D.E. and B.N. GARDINER. 
1982. Hydrodynamic Dispersion
during Unsteady, Unsaturated 
Water Flow in a Clay Soil. Soil 
Sci.
Soc. Amer. J. 46: 9-14.
(177) , K.M. PERROUX, 
S.J. ZEGELIN, and P.A.C. RAATS. 
1981.
Hydrodynamic Dispersion During 
Constant Rate Absorption of Water
by Soil. Soil Sci. Soc. Amer. J. 
42:229-234.
(178) and J.R. PHILIP. 
1978. Solute Transport during
Absorption of Water by Soil: Laboratory 
Studies and Their
Practical Implications. Soil Sci. Soc. 
Amer. J. 42:537-544.
(179) , J.H. KNIGHT, 
and D.E. ELRICK.
1978. Hydrodynamic Dispersion during Absorption 
of Water by Soil.
Soil Sci. Soc. Amer. J. 42:229-234.
(180) SPOSITO, G. and W.A. JURY. 1985. Inspectional 
Analysis in the
Theory of Water Flow through Unsaturated Soil. 
Soil Sci. Soc.
Amer. J. 49:791-798.
(181) , and V.K. GUPTA. 
1986. Fundamental
Problems in the Stochastic Convection-Dispersion Model of 
Solute
Transport in Aquifers and Field Soils. Water Resour. 
Res.
22:77-88.
(182) STARR, R.C., R.W. GILLHAM, and E.A. SUDICKY. 
1985. Experimental
Investigation of Solute Transport in Stratified Porous 
Media. 2.
The Reactive Case. Water Resour. Res. 21:1043-1050.
(183) SUDICKY, E.A. 1986. A Natural Gradient Experiment 
on Solute
Transport in a Sandy Aquifer: Spatial 
Variability of Hydraulic
Conductivity and its Role in the Dispersion 
Process. Water
Resour. Res. 22:2069-2082.
(184) and E.O. FRIND. 
1982. Contaminant Transport in
Fractured Porous Media: Analytical Solutions 
for a System of
Parallel Fractures. Water Resour. Res. 
18:1634-1642.
(185) , 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:1035-1041.
(186) TANG, D.H., E.O. FRIND, and E.A. 
SUDICKY. 1981. Contaminant
Transport in Fractured Porous Media: An Analytical 
Solution for a
Single Fracture. Water Resour. 
Res. 17:SSS-S64.
(187) , F.W. SCHWARTZ, 
and L. SMITH. 1982. Stochastic
Modeling of Mass Transport 
in a Random Velocity Field. 
Water
Resour. Res. 18:231-244.
154
(188) TANJI, K.K. 1970. A Computer Analysis on the Leaching of Boron
from Stratified Soil Columns. Soil Sci. 110:44-51.
(189) TAYLOR, G.I. 1953. Dispersion of Soluble Matter in Solvent
Flowing through a Tube. Proc. R. Soc. London Ser. A. 219:186-203.
(190) VALOCCHI, A.J. 1984. Describing the Transport of Ion-exchanging
Contaminants Using an Effective Kd Approach. 
Water Resour. Res.
20: 499-503.
(191) . 1985. Validity of the Local Equilibrium Assumption
for Modeling Sorbing Solute Transport through Homogeneous 
Soils.
Water Resour. Res. 21:808-820.
(192) , R.L. STREET, and P.V. ROBERTS. 1981. Transport of
Ion-exchanging Solutes in Groundwater: Chromatographic 
Theory and
Field Simulation. Water Resour. Res. 17:1517-1527.
(193) VAN DER ZEE, S.E.A.T.M. and W.H. VAN RIEMSDIJK. 1987. Transport
of Reactive Solute in Spatially Variable Soil Systems. Water
Resour. Res. 23:2059-2069.
(194) VAN EYKEREN, J.C.H. and J.P.G. LOCH. 1984. Transport of' Cationic
Solutes in Sorbing Porous Media. Water Resour. Res. 20:714-718.
(195) VAN GENUCHTEN, 
M.Th. 1978. Mass 
Transport in 
Saturated-
Unsaturated Media. I. One-dimensional Solutions. Res. Rep.
78-WR-11. Princeton Univ.
(196) . 1981. Non-equilibrium Transport 
Parameters
from Miscible Displacement Experiments. Research Report 
No.119,
U.S. Salinity Laboratory, Riverside, Calif.
(197) and W.J. ALVES. 1982. Analytical Solutions
of the One-dimensional Convective-Dispersive Solute Transport
Equation. USDA Tech. Bull. No. 1661.
(198) and R.W. CLEARY. 1982. Movement of 
Solutes
in Soil: Computer-Simulated and Laboratory Results. In Soil
Chemistry. B. Physico-chemical models. G.H. Bolt (Ed.). Elsevier.
Amsterdam.
(199) and W.A. JURY. 1987. Progress in Unsaturated
Flow and Transient Modeling. Rev. of Geophys. 25:101-105.
Amer. Geophys. Union.
(200) and J.C. PARKER. 1984. Boundary Conditions
for Displacement Experiments through Short Laboratory Columns.
Soil Sci. Soc. Amer. J. 48:703-708.
1ss
(201) , D.H. TANG, and R. GUENNELON. 1984. Some
Exact Solutions for Solute Transport through Soils Containing
Large Cylindrical Macropores. Water Resour. Res. 20:335-346.
(202) and P.J. WIERENGA. 1976. Mass Transfer
Studies in Sorbing Porous Media. I. Analytical Solutions. 
Soil
Sci. Soc. Amer. J. 40:473-480.
(203) . 1977. Mass Transfer
Studies in Sorbing Porous Media. II. Experimental Evaluation with
Tritium ( H 0). Soil Sci. Soc. Amer. J. 41:272-278.
2
(204) . 1986. Solute Dispersion
Coefficients and Retardation Factors. In Methods of Soil
Analysis. I. Physical and Mineralogical Methods. Soil Sci. Soc.
Madison, Wis.
(205) , and G.A. O'CONNOR. 1977.
Mass Transfer Studies in Sorbing Porous Media. III. Experimental
Evaluation with 2,4,5-T. Soil Sci. Soc. Amer. J. 41: 278-284.
(206) VILLERMAUX, J. and W.P.M. VAN SWAAY. 1969. Modele Representatif
de la Distribution des Temps de Sejour dans un Reacteur
Semi-infini a Dispersion Axiale avec Zones Stagnantes.
Application a 1'Ecoulement Ruisselant dans des Colonnes 
d'Anneaux
Raschig. Chem. Eng. Sci. 24:1097-1111
(207) WAGENET, R.J. 1983. Principles of Salt Movement in Soils. In
Chemical Mobility and Reactivity in Soil Systems. Soil 
Sci. Soc.
Amer. Publ. 11. Madison, Wis.
(208) WARRICK, A.W. and A. AMOOZEGAR-FARD. 1979. Infiltration and
Drainage Calculations Using Spatially Scaled Hydraulic
Properties. Water Resour. Res. 15:1116-1120.
(209) , J.W. BIGGAR, and D.R. NIELSEN. 1971. Simultaneous
Solute and Water Transfer for an Unsaturated Soil. Water Resour.
Res. 7:1216-1225.
(210) , G.J. MULLEN, and D.R. NIELSEN. 1977. Scaling
Field-measured Soil Hydraulic Properties Using a Similar Media
Concept. Water Resour. Res. 13:355-362.
(211) and D.R. NIELSEN. 1980. Spatial Variability of Soil
Physical Properties in the Field. In D. Hillel. Applications of
Soil Physics. Academic Press. New York.
(212) WATSON, K.K. and M.J. JONES. 1982. Hydrodynamic Dispersion during
Absorption in a Fine Sand. 2. The Constant Flux Case. Water
Resour. Res. 18:1435-1443.
156
(213) WHEATCRAFT, S.W. and S.W. TYLER. 1988. An Explanation of Scale
Dependent Dispersivity Using Concepts of Fractal Geometry. Water
Resour. Res. 24:566-578
(214) WIERENGA, P.J. 1977. Solute Distribution Profiles Computed with
Steady-state and Transient Water Movement. Soil Sci. Soc. Amer.
J. 41:1045-1051.
(215) WILLIAMS, J.K. and R.A. DAWE. 1986. Fractals-An Overview of
Potential Applications to Transport in Porous Media. Transport
Porous Media 1:201-210.
(216) YEH, T.C., L.W. GELHAR, and A.L. GUTJAHR. 1985. Stochastic
Analysis of Unsaturated Flow in Heterogeneous Soils. 1.
Statistically Isotropic Media. Water Resour. Res. 21:447-456.
(217) . 1985. Stochastic
Analysis of Unsaturated Flow in Heterogeneous Soils. 2.
Statistically Anisotropic Media with Variable a. Water Resour.
Res. 21:457-464.
(218) . 1985. Stochastic
Analysis of Unsaturated Flow in Heterogeneous Soils. 3.
Observations and Applications. Water Resour. Res. 21:465-471.
(219) YULE, D.F. and W.R. GARDNER. 1978. Longitudinal and Transverse
Dispersion Coefficients in Unsaturated Plainfield Sand. Water
Resour. Res. 14:582-S88.