for an Exchangeabl
Solute in a
Layered Soil
II0I
. 0 *0 -
0 -I. ~ Ill 0 S-I l
November 1989
Agronomy and Soils Departmental Series No. 139
Alabama Agricultural Experiment Station
Auburn University
Lowell T. Frobish, Director
Auburn University, Alabama

ANALYTICAL AND NUMERICAL SOLUTIONS OF THE TRANSPORT EQUATION
FOR AN EXCHANGEABLE SOLUTE IN A LAYERED SOIL
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
LIST OF FIGURES........................ 
v
ABSTRACT........................... 
ix
INTRODUCTION.......................... 
1
THEORY . . . . . . . . . .
. . . . . . . . . . . . . . . . .0 
3
Analytical Solution of the ADE 
for a Nonreactive Solute in
a Homogeneous Medium.................... 
3
Analytical Solution of the ADE 
for a Nonreactive Solute in
a Two-layer Medium..................... 
8
Numerical Solution of the ADE for 
an Exchangeable Solute in
a Homogeneous Medium...................11
Numerical Solution of the 
ADE for an Exchangeable Solute 
in
a Layered Medium........................................ 
15
Additional Remarks about the Numerical 
Solution of the ADE . 20
RESULTS AND DISCUSSION.................... 
27
Numerical and Analytical Results 
for Transport of a
Non-reactive Solute...................................... 
28
The Effect of Nonlinear Exchange 
on Solute Transport .
43
The Influence of Soil Layers 
with Different Transport
Properties on Solute Transport.............. 
52
SUMMARY AND CONCLUSIONS4.................... 
63
LITERATURE CITED............................................65
APPENDIX A. Analytical Solution 
of the ADE for the Second 
Layer
of a Porous Medium Using 
a First-type Condition.......... 
68
APPENDIX B. Analytical Solution 
of the ADE for the Second Layer
of a Porous Medium Using 
a Third-type Condition........... 
73
APPENDIX C. Analytical Solution 
of the ADE for a Two-layer
Medium Using a Third-type Condition 
at the Inlet and a
Combined First- and Third-type 
Condition at the Interface 
. 77
APPENDIX D. Numerical Integration 
with the Gauss Chebyshev
formula..........................................80
iii
APPENDIX E. Expressions for a., b., 
c., and f. to Solve the
1 1 1
ADE for a Non-layered Soil and a 
Layered Soil
APPENDIX F. List of Parameter Values 
for Simulations
APPENDIX G. Listing of Computer Programs
FIRST PRINTING, NOVEMBER 1989
Information contained herein is 
available
to all without regard to race, color, 
sex,
or national origin
iv
81
. . . . 83
86
LIST OF FIGURES
FIG.1. Schematic of 
the problem and grid 
system for the
numerical solution of 
transport in a homogeneous 
soil.
The solid circles refer 
to grid points with 
known
values and the open 
circles refer to grid 
points
with unknown values 
................. 
13
FIG.2. Schematic of 
the problem and grid 
system at the
interface for the numerical 
solution of transport
in a layered soil. The 
solid circles and 
triangles
refer to grid points 
with known values and 
the open
circles and triangles 
refer to grid points
with unknown values 
.................. 
16
FIG.3A. The influence 
of the space step, 
Ax, on the numerical
solution of the ADE. 
The solid lines are based 
on
the analytical solution 
............... 22
FIG.3B. The influence of 
the time step, At, on 
the numerical
solution of the ADE. The 
solid lines are based-on
the analytical solution 
.... ............... 
22
FIG.4. Numerical solution 
of the ADE with zero 
dispersion
using two different 
approximations for 
8X/8t as
presented in Eq.(22) 
and (32). The solid line
represents the theoretical 
step front 
........ 24
FIG.SA. Numerical. solution 
of the ADE for a homogeneous 
soil
using a first-type condition 
............ 
29
FIG.5B. Numerical solution 
of the ADE for a homogeneous 
soil
using a third-type condition 
... ............ 
. 29
FIG.6A. Validation of 
the analytical solution 
of the ADE for a
two-layer medium using 
a first-type condition 
. . . . 31
FIG.6B. Validation of 
the analytical solution 
of the ADE for a
two-layer medium using 
a third-type condition 
. . . . 31
FIG.7A. Validation of 
the numerical solution 
of the ADE for
a two-layer medium 
using a first-type condition 
. . . 32
FIG.7B. Validation of 
the numerical solution 
of the ADE for
a two-layer medium using 
a third-type condition 
. . . 32
FIG.7C. Validation of the numerical solution 
of the ADE for
a two-layer medium using 
a third-type condition 
at the
inlet and a second-type 
condition at the interface 
. .33
V
FIG.8A. Numerical solution of the ADE 
for a two-layer medium
(D
1
<D
2
, v <v ) using a first-type condition
FIG.8B. Numerical solution of the ADE for a 
two-layer medium
(D
1
<D
2
, v
1
<v ) using a third-type condition .....
FIG.8C. Numerical solution of the ADE for a 
two-layer medium
(D <D v <v ) using a third-type condition 
at the
inlet and a second-type condition at the interface
FIG.9A. Numerical solution of the ADE 
for a two-layer medium
(D
1
>D
2
, v >v ) using a first-type 
condition
FIG.9B. Numerical solution of the ADE for a two-layer 
medium
(D >D
2
, v >v ) using a third-type condition .....
FIG.9C. Numerical solution of the ADE for a two-layer 
medium
(D >D , v >v
2
) using a first-type condition at the
inlet
2
and a second-type condition at the interface
FIG.10A. Analytical solution of the ADE for a two-layer 
medium
(D
1
<D
2
, v
1
=v
2
) using a first-type condition ....
FIG.10B. Analytical solution of the ADE for a two-layer 
medium
(D <D
2
, v
1
=v
2
) using a third-type condition .....
1 2'1
FIG.11A. Analytical solution of
(D >D
2
, v =v
2
) using a
FIG.11B. Analytical solution of
(D
1
>D
2
, v =v ) using a
FIG.12A. Analytical solution of
(D
1
=D
2
, v <v
2 
) 
using a
FIG.12B. Analytical solution of
(D
I
=D
2
, v <v ) using a
FIG.13A. Analytical solution of
(D
1
=D , v >v ) using a
FIG.13B. Analytical solution of
W =D , v >v ) using a
D12 12
the ADE for a two-layer medium
first-type condition
the ADE for a two-layer medium
third-type condition ....
the ADE for a two-layer medium
first-type condition
the ADE for a two-layer medium
third-type condition .....
the ADE for a two-layer medium
first-type condition ... .
the ADE for a two-layer medium
third-type condition .....
FIG.14A. The effect of front retardation on transport.
solution of the ADE using R = 1 .
FIG.14B. The effect of front retardation on transport.
solution of the ADE using R = 38.5 . . ..
Analytical
Analytical
a . . . .
34
34
35
36
36
37
38
38
39
39
41
41
42
42
44
44
vi
FIG.15A. The effect 
of nonlinear exchange 
on transport. Numerical
solution of the ADE 
for unfavorable 
exchange .... 
46
FIG.15B. The effect 
of nonlinear exchange 
on transport. Numerical
solution of the ADE 
for linear exchange. 
....... 46
FIG.15C. The effect 
of nonlinear exchange 
on transport. Numerical
solution of the ADE 
for favorable exchange 
47
FIG. 16. The effect 
of front sharpening 
on transport. Numerical
solution of the ADE 
for favorable exchange......... 
47
FIG.17A. The effect 
of nonlinear exchange 
on dispersion:
concentration profiles 
for various exchange 
isotherms
using a moving coordinate 
system after 1 d. 
..... 49
FIG.17B. The effect 
of nonlinear exchange 
on dispersion:
concentration profiles 
for various exchange 
isotherms
using a moving coordinate 
system after 2 d. 
..... 49
FIG.17C. The effect 
of nonlinear exchange 
on dispersion:
concentration profiles 
for various exchange 
isotherms
using a moving coordinate 
system after 3 d......... 
50
FIG.17D. The effect 
of nonlinear exchange 
on dispersion:
concentration profiles 
for various exchange 
isotherms
using a moving coordinate 
system after 4 d ...... 
50
FIG.18A. Numerical 
solution of the ADE 
for a pulse input
during unfavorable 
exchange...................... 
.51
FIG.18B. Numerical solution 
of the ADE for a pulse 
input
during linear exchange............................ 
.51
FIG. 18C. Numerical solution 
of the ADE for a pulse 
input
during favorable 
exchange.............. 
52
FIG. 19. The effect 
of layers with different 
dispersion
coefficients on solute 
transport. Numerical 
solution
of the ADE for a two-layer 
me'dium2using a third-type
condition with: (a)2D 
= D2= 10 cm /d, and 
(b) D1=10
cm-/d and D =100 cm/a 
2.............. 
54
2
vi
FIG.20A. The effect of layers with different exchange capacities
on solute transport. Numerical solution of the ADE for
a two-layer medium using a third-type condition with
(ST) 1 
=
(S
T
) 2=0 cmol/kg . ....... ................ 56
FIG.20B. The effect of layers with different exchange capacities
on solute transport. Numerical solution of the ADE for
a two-layer medium using a third-type condition with
(ST) =0 and (S
T
) 2=2 cmol/kg ............. .56
FIG.20C. The effect of layers with different exchange capacities
on solute transport. Numerical solution of the ADE for
a two-layer medium using a third-type condition at the
inlet and a second-type condition at the interface with
(ST) =0 and (S
T ) 
2=2 cmol/kg ............. 57
FIG.21. The effect of layers with different exchange isotherms
on solute transport. Numerical solution of the ADE for
a two-layer medium with: (a) (dY/dX)=1 and (dY/dX) =3X
2
1-2/3 2
and (b) (dY/dX) = 1 and (dY/dX) = _X - . . . . . . 58
1 2 3
FIG.22. Concentration profiles in a hypothetical sand with a
third-type inlet condition: nonlayered profile (a),
and profile with a hypothetical clay layer using a
third-type interface condition (b), and second-type
interface condition (c) ...................... ... 60
FIG.23. Concentration profiles in a hypothetical clay with
a third-type inlet condition: nonlayered profile (a),
and profile with a hypothetical sand layer using a
third-type interface condition (b), and second-type
interface condition (c) ..... ............... . 62
viii
ABSTRACT
The study of solute transport in layered media is important to
investigate how solute movement can be slowed by incorporating a
retarding layer in the soil profile and to examine the influence 
of
soil heterogeneity on solute transport. This is conveniently done by
comparing transport through layered media with transport 
through
homogeneous media.
A relatively straightforward implicit finite difference 
method
was used to simulate reactive solute transport in layered and
homogeneous soil profiles. The effects of numerical dispersion and the
differences between the use of a concentration and a flux boundary
condition at the inlet were demonstrated for a homogeneous soil.
Analytical solutions were obtained for a medium consisting of two
layers. Using a flux-type condition at the inlet boundary and
interface resulted in discontinuities in solute concentrations at these
posi tions.
The importance of the exchange capacity and the exchange isotherm
for the simulation of solute transport was investigated for some
hypothetical cases. Both of these phenomena have a strong impact on
solute transport..-Several examples of a solute pulse traveling through
a two-layer medium were presented. Simulation of solute transport
through a soil profile, in which a soil layer ,with low permeability
characteristics (soil barrier) was embedded, showed that the adsorption
capacity of the barrier was not fully utilized to retard the solute
because of dispersive transport.
ix
INTRODUCTION
Transport of chemicals in soils is of great interest in a number
of areas in agriculture and engineering. This interest is aroused by
the potential of chemicals (fertilizers, pesticides, wastes etc.) to
move from the soil surface through the unsaturated zone of the soil
towards the groundwater. This process generally results in a
deterioration of the quality of soil and water resources. Once a
chemical or solute reaches the groundwater, spread of the chemical is
greatly facilitated by streamflow within the aquifer system. Therefore,
the description of solute movement through the unsaturated zone is an
important research topic.
The two objectives of this study were to (1) simulate 1-D
transport of inorganic salts through homogeneous and layered soils and
(2) investigate the effects of linear versus nonlinear exchange.
Layering of soil, with layers of different physical and chemical
properties, is of interest due to the influence artificial or natural
barriers might have in slowing or preventing the movement of certain
chemicals. The effects of layering are also important to the study of
transport in heterogeneous media. A layered medium, consisting of a
series of homogeneous soil layers, is conceptually equal 'to a
heterogeneous soil. The second objective is important to determine how
much effort should be made to determine exchange properties. The
influence of the cation exchange capacity and the exchange isotherms
are investigated for some hypothetical cases with conditions similar to
our later experimental work (19).
Although solute transport in porous media is most often described
with the advection dispersion equation (ADE) a number of 
different
approaches have been proposed to model transport because 
the
application of the ADE might not be justified for many instances in the
field (15). However, the use of the ADE offers distinct advantages
because of its deterministic nature. For a limited number of
situations, analytical solutions are available. In most cases, however,
numerical methods need to be employed, because the specific conditions
for which analytical solutions exist are not encountered in the field.
Unfortunately, these methods do not permit as much insight into the
effects of physical and chemical processes on solute transport as
analytical solutions do, and they are susceptible to computational
errors. When long-term predictions of the behavior of chemicals in the
soil need to be made, careful consideration should be given to the
choice of the numerical model and the processes to be modeled.
A vast amount of work has been published on the numerical
solution of the ADE simulating chemical transport in soils. Finite
difference methods are quite popular for this purpose, because of their
relatively straightforward formulation of the problem (27). However,
the somewhat more Complicated finite element method has been used as
well (e.g., 26). The latter method has distinct advantages when systems
with irregular boundaries are considered. Another method is that of
characteristics, which is less vulnerable to the numerical problems
which plague finite difference solutions.
3
Because of the simple geometry 
of the media and the nature of
this study, viz. the examination of 
the importance of various physical
and chemical properties on solute 
transport, the numerical solution of
the ADE was carried out with a finite difference 
method. Verification
of the solution was accomplished with 
comparison to analytical
solutions and mass balance calculations.
THEORY
-The problem of transport of a nonlinearly exchanging solute 
is
formulated for a homogeneous porous medium and for a porous 
medium
consisting of distinct layers. Analytical solutions are given 
for
transport in a (homogeneous) one-layer and a two-layer medium subject
to a flux- and a concentration-type boundary condition at the inlet 
and
the interface. Numerical solutions of the ADE are also given 
for these
cases, but in addition solutions are presented for a 
flux-type
condition at the inlet of the soil and a combined flux- 
and
concentration-type at the interface of soil layers.
Analytical Solution of the ADE for a Nonreactive
Solute in a Homogeneous Medium
For a number of cases, involving simple boundary and initial
conditions and linearly exchanging solutes, analytical solutions 
of the
ADE are available. These solutions can be utilized to predict
transport, to determine transport parameters, or to validate numerical
solutions obtained for the same conditions. For 1-D steady-state flow,
transport of a nonreactive solute in a homogeneous soil column with
length L can be described as follows:
ac a
2
c aC
a-D-2 v 
O<x<L t>O (1)
at 2 8x
8x
where C is the solute concentration 
expressed in mass per volume of
solution [ML-3], t is time [T], 
D is the dispersion coefficient
2-1 
-1
[L2T 
- 1
, v is the mean pore 
water velocity [LT ], 
and x is the
coordinate in the direction of flow 
[L]. The dispersion coefficient
quantifies two mechanisms which contribute 
to solute spreading, viz.,
molecular diffusion and mechanical dispersion. 
A common expression for
the dispersion coefficient is:
D = D /A + alv 
(2)
where D is the coefficient of molecular 
diffusion in a free solution
0
[L2 T-1, X is the tortuosity factor, a is the dispersivity [LI, and Ivi
is the absolute value of v.
Although we assumed nonreactive solute transport, 
i.e., R=1, it
should be noted that transport of a linearly exchanging 
solute can be
described by making the substitutions v=v and 
D=D , where v =v/R and
*
D =D/R.
The choice of the initial and boundary conditions 
deserves
careful attention. If they do not represent the 
physical system, the
mathematical solution based on these conditions 
is of little practical
use. It is assumed that the initial concentration 
in the column is
constant:
C(x, o0) = C
(3)O<x<L
5
A great deal of attention has been paid to the proper 
formulation
of the boundary conditions for the ADE (11, 16, 34). We considered 
the
following boundary conditions at the inlet:
C1xNO = C
0  
(first-type) t>O (4-a)
(-D a+ vC)I V
0 
= VC (third-type) t>O (4-b)
where xO, sometimes denoted as x-O+, implies that x=O is approached
from inside the column. The first-type condition (Eq.(4-a)) suggests
that the concentration at the interface of influent solution (eluent)
is continuous, whereas the third-type condition (Eq.(4-b)) implies that
the flux is continuous at the interface. For solute displacement
experiments, the latter condition is the more realistic one.
In contrast with the first- or concentration-type condition, the
third- or flux-type condition leads to conservation of mass. However,
this condition has a few drawbacks. First, one might question the use
of a dispersion coefficient close to the interface (11). For practical
reasons we will accept its validity. Second, the third-type condition
leads to a discontinuity in solute concentration at the inlet.
Physically, this will not occur at the microscopic level because of
molecular diffusion.
To put these conditions in perspective, it is helpful to
formulate the theoretical boundary condition at the inlet as:
.= Cx,
0  
t>O (5-a)
00o 8-x 
+
ev C) = C-e D IC
C (x, t) = C ?xs--L t>O (5-c)
0 0
6
where e is the volumetric water content [L3L-3], 
L is the arbitrary
0
length of a boundary layer outside the soil 
column [L], and the
subscripts o and 1 refer to the 
eluent and the soil column,
respectively. All concentrations are resident concentrations 
(21). The
conditions expressed in Eq. (5) are of little practical 
use because the
location of L is generally not known. Therefore, Eq. 
(5-c) is omitted
0
and it. is assumed that the applied solution is well 
mixed with eluent
concentration C . Note that the concentration of the eluent is 
now a
0
flux-averaged concentration. Eq. (4-a) and (4-b) then follow from
Eq.(5-a) and (5-b). We prefer the use of Eq.(4-b) because of mass
preservation, accepting the discontinuity in concentration due 
to the
use of two different types of concentration, i.e., a flux-type for the
eluent and a resident-type for the soil. Solutions of the ADE subject
to Eq. (4-a) can also be useful (34). Therefore both Eq. (4-a) and (4-b)
will be used.
For the outlet, the boundary condition can be formulated as:
aC = 0 
t>O (6-a)
ax- xtL
where it is assumed that the concentration is macroscopically
continuous over the outlet boundary. This condition is invoked for
diffusive transport in finite media, e.g., mass transport by molecular
diffusion or heat transport with thermal insulation at x=L. In these
cases, no transport will occur across the surface at x=L. Physically,
it is difficult to imagine how dispersive transport will vanish during
transient advective transport across the outlet boundary. Furthermore,
van Genuchten and Parker (3.4) argue that in analogy to the macroscopic
7
discontinuity at the inlet, using Eq.(4-b), a similar discontinuity
might occur at the outlet. This invalidates Eq.(6-a). The problem of
formulating a boundary condition at x=L can be avoided by assuming that
the medium is in effect semi-infinite, i.e.:
aC
a 0 
t>0 (6-b)
ax x-)w
Based on the considerations outlined by van Genuchten and Parker (34),
we decided to use this outlet condition for our analytical solutions.
It should be noted that the differences between solutions for finite
and semi-infinite media are usually rather small (33).
The analytical solution of Eq.(1) subject to Eq.(3), (4-a) and
(6-b), i.e., for a semi-infinite medium with a first-type
condition at the inlet, is given by (18):
C(x,t) = C + (C -C -(+c IE + cca-)(7)
whereas for the third-type condition at the inlet (i.e., Eq.(1) subject
to Eq.(3), (4-b) and (6-b)), the solution is given by:
C(x,t) =Ci +(Co-Ci) + 4v t emp- 4Dt +
o 2V/-2 LxD
2 vx vx+vtl(
- (i + vx + t) xL 2J (8)
The evaluation of enc with standard tables is rather tedious. The
function is available as a library function for a number of computer
languages. However, it can also be approximated with a series expansion
(1). Van Genuchten and Alves (33) provided approximations for the
function eap(A)e4z:(B) which can be programmed easily.
8
Analytical Solution of the ADE for a Nonreactive
Solute in a Two-layer Medium
Analytical solutions can also be developed for transport in media
consisting of homogeneous layers. A variety of solutions does already
exist for heat conduction in composite media (21). Solutions for
transport have been reported by Shamir and Harleman (28), and Al-Niami
and Rushton (3). Shamir and Harleman (28) used a first-type inlet
condition and an outlet condition at infinity for each layer, whereas
Al-Niami and Rushton (3) used a first-type inlet condition and a finite
outlet condition for each layer. Neither solution leads to conservation
of mass. We solved this problem assuming that each layer is in effect
semi-infinite, using a first or a third-type condition at the inlet of
each layer. The solution procedure becomes more tedious with the
increase of the number of layers. Therefore, we only considered a
medium consisting of two homogeneous layers, with an interface located
at x=L . The flow is perpendicular to the interface.
1
The solution of the ADE for the first layer was already discussed
in the section on analytical solutions for a homogeneous medium.
Mathematically, the presence of a second layer does not influence the
solution for the first layer because the outlet condition of the first
layer was chosen at infinity. Physically, this seems a reasonable
assumption for steady flow if no extreme differences in dispersion
coefficients between the layers exist, i.e., upstream dispersion is
negligible.
The ADE for the second layer is given by:
9
aC D 
2
C_ 8C-c - -C v -
L <x<L 
t>O 
(9)
at 2 2 2ax 1 2
ax
subject to:
C(xO) = C L,<x<L (10-a)
i1 2
C xL
1 
CIxL t>O (first-type) (10-b)
1 1 ax 1 1IxT"L 
2 2 x +2v2CIx",L 
t (hr-tp)(Oc
1
c = 0 t>O (10-d)
a x x- w
where the subscripts 1 and 2 refer to the first and second layer. We
assumed that the initial concentration, C is equal for both layers.
For a first-type condition, the problem is solved in conjunction with
Eq.(7), whereas Eq.(8) was used to solve the problem for a third-type
condition. From now on it is understood that the term condition refers
to the boundary condition at the inlet of a homogeneous layer or a
nonlayered soil. Using Eq.(11) from Parker and van Genuchten (22), it
can be shown that Eq.(10-b) implies that the volume averaged or
resident concentration is continuous at the interface, while Eq.(10-c)
implies that the flux-averaged concentration is continuous.
The analytical solution of Eq. (9), subject to (10-a), (10-d), and
(10-b) or (10-c), was obtained with Laplace transforms. Appendices A
and B contain the details of the respective solution procedures. For a
first- and third-type condition, the results are, respectively:
10
(C -C ) (x-L ) tL-vT vL 1L +vT
C(x,t) = 
o 
1 + eapD1
l/6tD [4DT 
1 1
2 0 1 1
1 
(x-L-v 
(t-T))
2
x3/2 4D(-) dr + C 
(11)
(t-T) 
3/2 
4D2 
(
t-T)
v (C -C ) L -vT 
vL L+vT
C(x,t) = 2 o 
4 D I+X
4DD 04D 4D T
2 1 1
22
1 
ep 
(x-L 
-v 
(t-T))2
/(t-T) 
4D2 
( t - T )
v v (x-L) x-v(t)
eeA 
-L+v(t-T) 
-dT + C 
(12)
S4D 
2  
4D (t-T)
2 2
The integrals appearing in 
these solutions were evaluated 
with
help of the Gauss-Chebyshev formula, 
discussed in Appendix D.
If we deal with resident-type 
concentrations, both concentration
and flux are continuous across 
the interface and the solutions 
given by
Eq.(11) and (12) will not yield 
the exact concentration profile. 
The
transport problem where both the 
flux and concentration are continuous
resembles the one of heat flow 
in composite media with no contact
resistance at the interface (8). 
To solve this problem analytically,
the (simplifying) assumption 
that the first layer is, in 
effect,
semi-infinite cannot be used. 
The exit and inlet condition 
of
respectively the first and second 
layer consist of the two interface
conditions (Eq. (10-b) and (10-c)). 
Appendix C contains the mathematical
statement of the problem as well 
as the solution for the concentration
of both layers in the Laplace domain.
11
Numerical Solution of 
the ADE for an Exchangeable
Solute in a Homogeneous 
Medium
To describe transport 
of a reactive solute 
we rewrite Eq. (1) as
follows:
8C 8s 8
2
C ac (13)
e -f+ p T= 
eD - ev 
ax
8x
-3
where p is the dry bulk 
density of the soil 
[ML ] and S is the
adsorbed concentration of 
the solute expressed as 
mass of solute per
-1
mass of dry soil [MM ].
Equation (13) needs to be rewritten 
to solve for the dependent
variable, C. We assumed that 
a unique relationship, the 
exchange
isotherm, exists between the 
adsorbed ions and those present 
in the
liquid phase. It should be noted 
that the relationship between S 
and C
generally also depends on the 
total ionic strength and other 
chemical
species present (23). Such 
dependency can be included implicitly 
by
experimentally determining the 
exchange isotherm under conditions
similar to those for which the 
ADE needs to be solved. The adsorbed
concentration can be expressed in 
terms of the liquid concentration, 
C,
using the chainrule:
8S dS8C 
(14)
at dC at
It is convenient to introduce 
the following dimensionless
concentrations to describe transport 
of nonlinearly exchanging solutes
(17):
X=C/CT 
( 15-a)
y=S/S
T
..... (15-b)
12
where X is the dimensionless concentration of the 
solute in the liquid
phase, CT is the total amount of solute in the liquid 
phase per volume
-3
of solution [ML -, Y is the dimensionless 
concentration of the solute
in the adsorbed phase, and ST is the total 
amount of solute in the
adsorbed phase per mass of soil1 [MM -1 
Substitution of these
dimensionless quantities in Eq.(13) and the use 
of Eq.(14) yields:
R x D 
a 
ax
at ax2 ax
where the dimensionless retardation factor, R, is defined 
as follows:
R=1 _PS
5
T dY 
(
8)C TdX
The familiar form of the ADE is obtained by dividing 
Eq. (16) by R:
6)
.7)
ax D x 
v* a(18)
at 2 ax
8x
where v =v/R and D =D/R.
Eq.(16) will be so1ved for-the following initial 
and
boundary conditions:
X(x,O0) =X (x
(-D av X) vX
7x x
49
0
ax0
8x IxtL=0
The numerical
Crank-Nicolson method.
temporal spacing At,
j=0, 1, .. the central
0<x<L
(first-type)
(third-type)
t>0
t>0
t>0
(19-a)
(19-b)
v
(19-C)
( 19-d)
solution is achieved with the implicit
Using a grid system with spatial spacing Ax and
where x.=(i-1)Ax and t.=jAt (i=1,2,..,n and
1 j
difference approximation of Eq.(16) is*.
13
R D2j (xj
1
X+X) - D+j+1 
jXj ) +
Et i i 2 1+ 1 i 
i-i l 
(X1+1X 2X 
ii-
(X , j+x;+X -X )
1=2, 3p,4p.,*9.n-1
A schematic of the problem, 
along with the grid system 
used for the
numerical approximation, is shown 
in figure 1.
The approximation of the initial 
condition is:
x.= x
1
2-4i --n
(1-a
The inlet condition (x=0) was approximated 
by:
xi = X
1 0
D j+ lj+l+j 
.J+
-- x (X -XO +X-X X~ x
4Ax 2 0 2 0 2 1 1
vc
0
2
3
n-1
L ___
vce
(first-type)
SvX 0(third-type)
J-1 j J+1
FIG. 1. Schematic of the problem and grid
numerical solution of -transport in a homogeneous
circles refer to grid points with known values 
and
refer to grid points 'with unknown values.
system f or the
soil. The solid
the open circles
(20)
j>0
j>-0
(21-b)
(2 1-c)
Y1 Ax
At
14
where i=O is a fictituous node located at 
x=-Ax. The central difference
approximation for the outlet condition 
was written as:
X+1
-
X + +X -X = 0 jO0 (21-d)
n+1 n-1 n+1 n-1
where the column exit is located at (n-1)Ax=L 
and i=n+l represents a
fictituous node located at L+Ax. The value of X 
depends on the
condition at the inlet of the soil. For the first-type 
the value is
given by (21-a), whereas for the third-type condition 
X was
approximated according to Barry et al. (4).
To solve the above set of equations, all unknown variables 
at
time t need to be expressed in terms of known variables at the
previous time t
J . 
The variables at the fictituous nodes i=0 and n+1
were expressed in variables at column nodes using the central
difference approximation of the ADE at x=0 and x=L. It is convenient to
define G=D/(2Ax
2 ) 
and H=v/(4Ax) to rewrite Eq.(20) as:
r R'
(-G-H) Xi + i +2G X. + (-G+H) Xj+
1
11 At i+1
(22)
R .
(G+H) X
j  
+ [ 2G] xi 
+ (G-H)X
1-1 At 1 l+1
This equation needs to be solved for i=2,...,n-1 along with the two
equations for the inlet and outlet condition. The resulting system of
equations can be written in matrix form:
15
b X 
f
a
2  
b c . X
1
f
2b2 2 22
a
3  
b3 c3 X3 f
3 3 
3 3
an-i bn c" Xnf
.- n-1 n-1 n-1 n-1
0o............a b X f
n n n n
j+1
The n by n matrix is a tridiagonal band matrix for 
which the unknowns
X
j + l  
are conveniently solved with the Thomas algorithm (25).
Expressions for the coefficients a
i
, b
i
, c.i and f. are provided in
Appendix E.-
Numerical Solution of the ADE for an Exchangeable
Solute in a Layered Medium
In most soils, the physical and chemical properties will not be
constant with depth, e.g., ST might decrease with depth whereas 
p
generally increases with depth. This section is concerned with abrupt
changes in properties as they occur in naturally layered 
soils or
artificially constructed media. Each layer is assumed 
to be
homogeneous, with constant parameter values, and the layers 
are
separated by a sharp interface. A.schematic of the problem is provided
in figure 2A.
The ADE for transport within the k-th layer is 
written as:
R -D -v -(24)
k 8t k X2 k 8x
The same initial, inlet, and outlet conditions 
are used as for the
homogeneous medium in the previous section.,
16
0
INTERFACE CONDITION
first-type third-type
i -2
i-i
F'
- i -H---
i+1L
i +2
j+1 -
j j+i j-1
k
k+1
j j+1
fictituous node
b
C
d
FIG. 2. Schematic 
of the problem and 
grid system at -the 
interface for
the numerical solution 
of transport iri 
a layered soil. The 
solid
circles and triangles 
refer to grid points 
with known values 
and the
open circles and triangles 
refer to grid points 
with unknown values.
I
O vC
110
O vC
K K
a
L 
k
L k+l
LK-i1
L 
K
x
17
The interface condition can be expressed in 
a similar way as in
the section on the analytical solution of solute 
transport in a
two-layer medium. For the interface between layer 
k and k+1, where
spatial node i is located at x=Lk) these conditions are, respectively:
XlxtLk= X x4L 
(first-type) (25-a)
(-(D) -+(xv) X) =(-(eD) -+(Gv) X) (third-type) (25-b)
ik 8X (V)k x"L k+15X k+1 x4L
k k
In addition to a first- and third-type condition, a second-type
condition was introduced by combining Eq.(25-a) and (25-b):
(GD) x = (D) X (second-type) (25-c)
ODjk jx"Lk 
k+1 ax 
xLk
k 8xxk
The second-type condition, in conjunction with the ADE on both sides
of the interface, stipulates that both concentration and flux are
continuous at the interface.
The solution of Eq.(24) subject to Eq.(19-a), (19-d), and (1)
Eq.(19-b) and (25-a), (2) Eq.(19-c) and (25-b), or (3) Eq.(19-c) and
(25-c) was accomplished with the Crank-Nicolson method. The space-time
grid system is set up in such a way that a node is located at each
interface. The numerical approximation of the ADE for the nodes inside
each (homogeneous) layer is analogous to that discussed in the previous
section.
For a first-type (interface) condition, Eq.(24) was solved with a
central difference scheme where the concentration at the interface of
layers k and k+1 was solved by using weighted coefficients. The grid
system is shown in figure 2B. For example, for the coefficient a.
1
18
(Eq. (23)) we get ai=wkal+wk+l ai+ where the weights wk and wk1
depend on the solute residence time in the space increments 
on either
side of the interface. The solute flux at the interface 
in layers k,
J k and k+1, Jk1' at spatial node i and temporal 
node j, can be
approximated with the following simple explicit expressions:
(eD)
j ev) - k(XJ-x ) 
(26-a)
k k 1 Ax 1 i-1
(OD)k1 0
S= (ev) x _ D) (X -X
j
) (26-b)
k+1 k+1 1 Ax i+1 i
The weight factors for layers k and k+1 (w and w ), were then
k k+1
obtained from the solute fluxes at the interface (x=L), where the
k
concentration at the previous time equals X.:
1
W 1 (27-a)
k 1+(J /J )
k k+1
wk =1-w (27-b)
k+1 k
The solution procedure is similar as outlined for the numerical
solution of solute transport in a homogeneous medium. Regular and
weighted coefficients a, b, c, and f are listed in Appendix E.
For the third-type (interface) condition, two different
approximations for Eq. (25-b) were used to obtain the concentration at
lim Xj+1l=x j+1
the interface of layers k and k+1, Viz. xtLXiit (figure 2C) and
k
.
+1 j 1j+1 
X0+1
lim X.6
1
=X.J
4  
(figure 2D). A notation similar to that for X, and X.
x Lk
was used for the coefficients and concentrations involved in obtaining
these variables. The two approximations of Eq. (25-b) are formulated in
such a way that Eq. (23) could be used to formhlate all equations:
19
AX k i+1 -1 i+1 i 1 i
2 D X - X + Xj+1 
X (28-a)
-
( 
)k " -X_1, 
+ X +1 + X 
(28-b)
S-) xj+~ j+ + x - X + Xj+ + X
AX V k+1 L i+1 1-1 1
4
, 11J
It should be noted that X and X are concentrations at
fictituous nodes located at Lk+Ax and Lk-Ax, respectively. They are
obtained by using the central difference approximation of Eq. (24) for
layers k and k+1. A system of equations as given by Eq.(23) can be
obtained, except that n+K-1 equations need to be solved, where K is the
total number of layers.
Finally, the second-type interface condition was approximated by:
(eD) - X
1
+ X -X9x = (eD) [+- X +xXi - X (29)
k1- i i-1 k+1 I+I 1 +i+I 1
where the same grid system was used as for the first-type condition
(figure 2B).
The coefficients of the matrix system are listed in Appendix E.
The subscript i denotes the row number of the matrix. The coefficients
a, b, c, and f depend on location and time, because D, v, e, p, S
T
, and
dY/dX vary with each layer and dY/dX varies with time.
20
Additional Remarks about the Numerical Solution of the ADE
It is well known that the numerical solution of the ADE is
considerably more difficult than the solution of the heat 
equation
(14). In fact, the numerical solution of the ADE might 
not be accurate
due to numerical dispersion or because of numerical oscillations. These
numerical oscillations occur especially around the solute front,
whereas numerical dispersion occurs at low values of D when the
advection term is dominant. By adjusting the increments of the grid
system, some improvement can be achieved. Unfortunately, efforts to
reduce numerical dispersion might increase numerical oscillations and
vice versa. Huyakorn and Pinder (14) suggested keeping the Peclet
number (Pe=vAx/D) less than 2 and the Courant number (Cr=vAt/Ax) less
than 1. The latter condition is also referred to as the
Courant-Friedrichs-Lewy stability criterion (24).
Excessive rounding errors, which lead to unstable numerical
solutions, are avoided if the tridiagonal matrix, given by Eq.(23), is
diagonally dominant (7), i.e., if for elements on the diagonal of the
matrix the following holds:
n
Piil Pij1 
i=1,2. ,n 
(30)
j=1
j 
i
where p.. is an element at row i and column j of the matrix. In our
case, using a Crank-Nicolson scheme, this leads to the requirement that
b.lI IlI+c~I. ro the expressions provided in Appendix E, it i
obvious that no excessiye rounding error will occur in a homogeneous
soil, regardless of the size of Ax and At. However, the values of Ax
21
and At are important to obtain an accurate numerical approximation of
the differential equation.
As an illustration, we consider 1-D transport of a nonreactive
solute in a homogeneous soil with a dimensionless initial concentration
X =0.05. The ADE is solved numerically for a first-type or Dirichlet
i
boundary condition at the inlet for different values of Ax and At. The
soil properties, relevant for solute transport, used in this example
and for all later examples, are listed in Appendix F. These properties
were chosen somewhat arbitrarily. Although the units of the parameters
are listed, any consistent set of units would do. In fact, one might
prefer to state and solve the problem in dimensionless form.
Some results are presented in figure 3A, which shows numerical
solutions employing three different space increments (0.05, 0.5, and
2.5 cm) and a constant time step (0.005 d), as well as the analytical
solution (solid line). The three solutions are very similar for the
selected values of D and v. As expected, some deviations from the
analytical solution began to occur for Ax=2.5 when Pe=2.5 exceeded the
recommended maximum of 2. Figure 3B contains results for various values
of At and a constant Ax of 0.5 cm. The Courant condition is satisfied
if At<0.01 d. However, the results indicate that a smaller At is
desirable; the numerical and the analytical solutions are almost
identical when At:-0.001. Some deviations between analytical and
numerical solutions occur at the outlet (x=12 cm), because of the use
of a semi-infinite versus a finite length (Eq. (6-b) and Eq. (19-d),
respectively). It should be noted that the Crank-Nicolson method is a
22
0 2 4 6 8 0 
12 14
x [cm]
FIG.3A. The influence of the 
space step, Ax, on the numerical 
solution
of the ADE. The solid lines are 
based on the analytical solution.
x [cm]
FIG.3B. The influence of the time 
step, At, on the numerical solution
of the ADE. The solid lines are 
based on the analytical solution.
23
second-order method, which approximates the derivative 
of a function f
according to:
f(x
0
+h)-f(x
0
-h) 
h2 (3)
f(x0 2h f (3)() e [xo-hx0+h] 
(31)
f(Xo) 2h 6 6
(3)
where f is the third derivative of f. The 
last term of Eq.(31)
represents the error term associated with this approximation. 
A higher
order expansion leads to a reduction of 
the error.
Application of a higher order expansion for the 
time derivative
might partly correct the problem of numerical dispersion 
(32). A simple
example of reducing numerical dispersion by using a weighted expression
for the time derivative was provided by Stone and Brian 
(29):
xj+
1
-xj j+
1
-xj j+
1 
j
8X _ 1 i+1 i+1 2 1 1 i-1 i-1(32)
at 6 t 3 t 6 At
The effect of numerical dispersion is demonstrated by assuming
zero dispersion (figure 4). The results were obtained by using either
ax
Eq. (22) or a variation of Eq.(22) by substituting Eq. (32) for thea
term. The exact solution of the ADE resembles a step front 
(D=0 and
R=1). Both numerical solutions compare poorly with the actual 
solute
front, although the use of Eq.(32) leads to a somewhat steeper front.
Numerous other techniques have been used to deal with this 
phenomenon
(cf. 2).
24
+
0
+ Eq. (22)
o Eq. (32)
-?-- -6 - -
+ o
t = 0.1 d
+
0
0 2 4 6 8 
10 12 14
x [cm]
FIG. 4. Numerical solution of the 
ADE with zero dispersion using two
different approximations for 
ax/at as presented in Eq.(22) and (32).
The solid line represents the theoretical 
step front.
Besides comparison with analytical solutions, 
the accuracy of the
numerical solution can also be evaluated 
with the help of a mass
balance. We will do this by comparing 
the net influx of solute, denoted
as FLUX [MI, with the amount of solute 
which has accumulated in the
column according to the numerical solution, 
denoted as ACC [MI.
The net amount of solute which enters 
a hypothetical soil column
during steady flow in a time interval 
At=tj+ -t
j , 
is given by the
following flux term:
25
tj+
1
FLUX
1 
= veAC X(xt) 
X(xt)xL dt 
(33)
J
where A is the area of the column perpendicular 
to the direction of
flow [L2]. The net amount of solute which has 
accumulated in the column
during this period of time is given by:
.L
ACCj+1=A f CT X(xtj+l)-X(x
,t
j) +pST Y(x,tj+l)-Y(x
,t
J) dx (34)
0
where L is the length of the column in the direction of 
flow [L], and
X(x,t.) and Y(x,t.) denote the dimensionless solute concentrations 
in
J J
liquid and adsorbed phase, respectively. X(x,t.) can 
be determined
J
mathematically or experimentally. Because of the continuity principle,
the net accumulation should be equal to the net influx. This 
provided a
useful tool to check the accuracy of the numerically 
determined
X(x.,t j+).
1 j+1
Denoting the solute concentration of the eluent solution as Xin,
Eq.(33) was approximated as:
FLUX. AvC (XI - (X+XJ+1)] At (35-a)
j+1 T in 2 n n
X 0 < t < t
X. = o o (35-b)
in 
0  
elsewhere
where t is finite for pulse displacement and infinite 
for step
0
displacement. The integration to obtain the accumulation term 
is
performed by using the trapezoidal rule, which is convenient to use 
for
transport in layered media.
26
Two mass balance equations were used, viz. a relative 
mass
balance (RMB) over the time interval [t ,tj I and an overall 
relative
j' J+1
mass balance (ORMB) over the interval [0, t j+]. The error criterionj+1
based on the first balance is defined as:
=
4
FLX{j1I/c .+ " 100% (6
RMBj+1 FLUX j+1-ACC j+/ACCj+ 
100% 
(36)
To determine the mass balance over the whole time period, t j+1, the
amount of solute present at t=0 is subtracted from the amount present
at t j+ to yield the overall term SUMACC +. The total net supply of
j+1 ~J+1
solute to the column, SUMFLU+. 
,
and the overall mass balance are
j+1
respectively defined as:
j+1
SUMFLUj+ = E FLUX (37-a)
j+1 J
J=1
and:
ORMBj+. = SUMFLUj+ -SUMACC j+I/SUMACCj+ *100% (37-b)
j+1VSMF j+1 j+1 Ij+1f
As stated before, the accuracy of the numerical solution depends
on the size of the space and time increments, but also on the number of
calculations (rounding error). The question then arises, concerning
what the appropriate sizes of Ax and At should be to obtain an accurate
solution using a minimum number of calculations.
The value for At during the calculations was optimized with the
help of the overall relative mass balance criterion. Two bounds on ORMB
were used, namely a lower bound of 0.5% and an upper bound of 3%. The
value of At for the next time t was increased in case ORMB<0.5% and
j+1
27
decreased when ORMB>3%. A numerical scheme which adjusts both Ax 
and At
to solve the 1-D flow problem was reported by Dane and Mathis (12).
At this point we refer to Appendix G, which contains listings of
various WATFIV programs of numerical solutions used in this study. 
The
programs HOMSOIL and 
VARTIM calculate solute 
concentrations in
homogeneous media with a fixed and a variable time step, respectively.
Transport in layered media was solved with the programs LAY1 and LAY3
for a first- and third-type condition, respectively, with a fixed time
step. The program LAY1 was modified to solve the ADE for a second-type
interface condition (i.e., LAY2). The various parameters for each
problem are specified in the declaration section of the program. The
value of the retardation coefficient was determined with the subroutine
EXCHAN and the amount of adsorbed solute with the routine ADSORB.
Subroutine TRIDIA contains the Thomas algorithm to solve the
dimensionless liquid concentrations. For most problems, a fixed time
step was used, and values for ORMB were generally less than 1%,
depending on the type of problem.
RESULTS AND DISCUSSION
The first section of this chapter compares some of the analytical
and numerical solutions for layered with solutions for nonlayered
media. Furthermore, concentration profiles for layered media were
investigated with analytical solutions using a first- or a third-type
condition at the inlet of each layer. In the remaining sections we will
use only a third-type condition at the inlet of the soil and a second-
or third-type condition at the interface of soil layers (Appendix E).
28
Some basic effects of nonlinear exchange on solute transport will be
demonstrated for pulse and step input in the second section. The third
section discusses the effect of layers, with different values for their
transport parameters, on the transport of a solute pulse, 
along with
two different scenarios of reducing contaminant movement in a soil
profile for a step input. Most numerical values of the parameters were
chosen for illustrative purposes. All relevant values are listed in
Appendix F.
Numerical and Analytical Results for Transport
of a Nonreactive Solute
We start by examining the influence of the inlet boundary
condition on the mass balance. Numerical solutions of the ADE for a
first- and a third-type condition, obtained with the program HOMSOIL,
are presented in figure 5A and 5B, respectively. The values for At and
Ax were 0.005 d and 0.5 cm, respectively. These values were also used
for all later calculations, unless specified otherwise. Included in
figure 5A and 5B are values of RMB and ORMB. For a third-type
condition, X(O,t) gradually approaches X (figure 58), whereas for the
0
first-type condition, X(O,t) is instantly equal to X (figure 5A). In
0
the latter case, too much solute is forced into the column, resulting
in a significant overall mass balance error, particularly at the early
stage. In contrast, the third-type condition preserved mass. The mass
balance criterion to determine the value of a variable At can therefore
not be used for a first-type condition. Using a variable time step for
a third-type condition resulted in a CPU time which was one fifth of
29
0 2 4 6 8 
10 12 14
x [cm]
FIG. SA. Numerical solution
using a first-type condition.
of the ADE for a homogeneous
0 2 4 6 8 10 12 14
x [cm]
FIG.5B. Numerical solution of the ADE for a homogeneous soil using a
third-type condition.
soil
30
the CPU time for a fixed time step (0.005 d). The value for ORMB varied
during the calculations, but was well below 1% for both schemes
(HOMSOIL and VARTIM), except during the very early stages of
simulation.
Next, we considered transport in a two-layer medium, with the
interface at x=6 cm. One way of validating the analytical and numerical
solutions for this problem is to use the same parameter values for both
layers and compare the solution to the one for a homogeneous medium.
The analytical solutions according to Eq.(11) and (12) are shown in
figure 6A (first-type condition) and 6B (third-type condition),
respectively, along with the solutions for a homogeneous medium (Eq. (7)
and (8), respectively). Figure 7A, 7B, and 7C show the results obtained
with LAY1, LAY3, and LAY2, respectively, as well as the results
obtained with HOMSOIL (solid lines). The results indicate close
agreement between the solutions for the "homogeneous" and the
"two-layer" media.
31
*solution according to Eq.(1
solution according to Eq. (7)
2 4
8 10 12 
14
FIG.6A. Validation of
two-layer medium using
x cm]
the analytical solution of
a first-type condition.
the ADE for a
t=0.1 d
layer I
V= 50 cm/d
D= 50 cm2/d
*solution according to Eq. ( 12)
solution according to Eq. (8)
0.2 d
layer 11
v= 50 cm/d
D= 50 cm2/d
10 12 14
x[cm]
FIG.6B. Validation of the analytical solution of the ADE for 
a
two-layer medium using a third-type condition.
1.0-
0.8-
0.6-
0.4-
W
32
0 2 4 6 8 10 
12 14
x [cm]
FIG.7A. Validation of the numerical solution of
layer medium using a first-type condition.
the ADE for a
14
x [cm]
FIG.7B. Validation of the numerical 
solution of the ADE for a 
two-
layer medium using a third-type 
condition.
two-
33
0 2 4 6 8 10 12 14
x [cm]
FIG.7C. Validation of the numerical solution of the ADE for a two-
layer medium using a third-type condition at the inlet and a
second-type condition at the interface.
Transport during steady flow in a two-layer medium was
2
investigated for a medium with v =25 cm/d, D =12.5 cm /d, and e =0.4
1 1 1
(i.e., the first layer parameter values), and v =100 cm/d, D =100
2 2
cm /d, and e =0.1. Figure 8 shows the numerically calculated
2
concentration profiles for a first-, third-, and second-type condition
at the interface (VWC denotes 0). Because the D value is smaller in the
first layer, the solute front in the first layer is steeper than in the
second layer. The concentration profiles obtained with a second- and a
third-type condition are very similar except at the interface. The use
of a third-type condition results in a discontinuity of the (resident)
concentration, whereas the profile is continuous for a second-type
34
f irst type
OAZ
t = 0.3 d
layer 1
0.2-- V= 25 
cm/d
D = 12.5 cm2/d
0.1 VWC= 0.4
layer 11
V= 100 cm/d
D= 100 cm2/d
VWC= 0.1
0 2 4 6 8 10 12 1
x[cm]
FIG.8A. Numerical solution of the ADE for a
v1 <v 2 using a first-type condition.
third type
two-layer medium (D <D,
0.
0 24 6 8 10 
12 14
x [cm]
FIG.8B. Numerical solution of the
v 1(v2 ) using a third-type condition.
ADE f or a two-layer medium (D <D
35
second type
x
0 2 4 6 8
x [cm]
FIG. 8C. Numerical solution of the ADE for a
v 
1
<v ) using a third-type condition at the
condition at the interface.
two-layer medium (D
1
<D
2
,
inlet and a second-type
condition. For a third-type condition, the concentration increased with
position across the interface. Note that in both cases the flux and the
flux-averaged concentration are continuous at the interface. Figure 9
contains the results if the layers are interchanged. The concentration
profiles are now steeper in the second layer, because of the small D-
value compared to the first layer. For a third-type condition, the
(resident) concentration decreased across the interface.
36
first type
0 2 4 6 8 10 12 14
x [cm]
FIG.9A. Numerical solution of the ADE 
for
v >v ) using a first-type condition.
1 2
third type
a two-layer medium (D >D,
. 1 2
x [cm].
FIG.9B. Numerical solution pf the ADE for 
a two-layer medium (D >D
2
,
v >v ) using 
a third-type 
condition.
1 2
37
second type
J*0 2 4 
6 8 10 12 
14.
X [cm]
FIG. 9C. Numerical solution of the ADE for a two-layer medium (D >D
v >v ) using a first-type condition at the inlet ana 
2
a
S2
second-type at the interface.
We investigated the influence of different parameter values
further by using analytical solutions. First we examined the influence
of dispersion, using somewhat artificially large differences in
D-values. In figure 10, the D-value in the first layer is smaller than
in the second layer, whereas in figure 11 the D-values are
interchanged. The error in the mass balance can be evaluated by
comparing the first-type solution (figure iQA and 11A) with the
third-type solution (figure lOB and 1iB). For a layered medium, the
repeated application of the first-type condition leads to an increase
in the mass balance error of the solution compared to a homogeneous
medium (cf. figure 5). For a third-type condition, a discontinuity in
38
x [cm]
FIG.10A. Analytical solution of the ADE for a
v =v ) using a first-type condition.
1 2
10 12 14
two-layer medium (D <D
12
0 2 4 6 8 10 12 14
x [cm]
FIG. lOB..Analytical solution of the ADE for a two-layer medium (D <D
12
v =v ) using a third-type condition.
1 2
39
0 2 4 6 8 10 
12 14
x [cm]
FIG. 11A. Analytical solution of the ADE 
for a
v =v ) using a first-type condition.
1 2
two-layer medium
x [cm]
FIG. 11B. Analytical solution of the ADE 
for a
v =v ) using a third-type condition..
1 2
two-layer medium (D >D
1 2
(D >D
1 2
d
14
40
concentration can be observed at the interface. This discontinuity is
similar to the one at the inlet where flux-averaged concentrations,
rather than resident concentrations, are continuous. During heat
transport in composite media, a similar discontinuity arises if there
is a contact resistance. Heat flow across the interface is then
described with a 'radiation' boundary condition using a heat transfer
coeffcient.
Next, we investigated the influence of differences in v-values on
the concentration profile for steady-state flow. In figure 12, the
pore-water velocity is smaller in the first than in the second layer.
Figure 13 shows the concentration profiles if the layers are
interchanged. The solute front is steeper in the low velocity layer.
For a third-type condition, the concentration is lower in the low
velocity layer than in the high velocity layer at the interface. The
values of the advective and dispersive fluxes are, respectively, lower
and higher in the low velocity layer than in the high velocity layer
because of the differences in e. We note again that the flux-averaged
concentration is continuous at the interface.
From these examples it can be concluded that the use of a
first-type condition leads to large mass balance errors in layered
media. For a third-type condition, the discontinuity of the resident
concentration at the interface depends both on the differences in v and
D. If these differences are not too large, the analytical solution
using a third-type condition gives a good approximation of the actual
concentration profile. The latter can be obtained by assuming that both
concent rat ion and fl1ux are cont inuous at the i nt erface ( cf.
figures 8B and 8C and figures 9B and 9C).
41
0 2 4 6 8 10 12 14
x [cm]
FIG.12A. Analytical solution of the ADE for a 
two-layer
v <v ) using a first-type condition.
1 2
medium (D =D,
12
0 2 4 6 8 10 12 14
x [cm]
FIG.12B. Analytical solution of the
v <v ) using a third-type condition.
1 2
ADE for a two-layer medium (D =D,
12
42
0 2 4 6 8 10 12 14
x [cm]
FIG. 13A. Analytical solution of the ADE for a
v >v ) using a first-type condition.
1 2
two-layer medium (D =D
12
0 2 4 6 8 10 12 14
x [cm]
FIG.13B. Analytical solution of the ADE for a two-layer medium 
(D=D
12
v >v ) using a third-type condition.
1 2
43
The Effect of Nonlinear Exchange on Solute Transport
Reactions between the porous medium and the solute can have a
strong impact on the transport of a solute. The only reaction
considered here was ion exchange in a binary system. However, other
processes or systems could be combined as well with expressions like
Eq.(13).
For the purpose of comparison, we considered first the case of
linearly exchanging solutes. The relevant parameter values and the type
of inlet conditions are given in Appendix F. An arbitrary soil (ST=10
3 3 3
cmol c/kg, p=1.5 
g/cm, 8=0.4 
cm /cm ) is saturated 
with a 0.005 
M CaBr
2
solution which is displaced with a 0.01 M KC1 solution. The K movement,
assuming linear equilibrium exchange between Ca and K, is described
using R=38.5. The Cl displacement can be described using R=1, if
neither anion reacts with the soil. The profiles for Cl and K were
obtained with an adapted version of Eq.(8) by substituting v for v and
D for D and are shown in figure 14A and 14B, respectively. Although
this situation is somewhat idealized, it is clear that the ability of
the soil to adsorb the incoming cation greatly increases its travel
time through the soil.
In many instances the exchange is not governed by a linear but by
a nonlinear exchange isotherm. Therefore, transport of nonlinearly
exchanging solutes will be examined for a binary system. Two cases can
be distinguished, viz. preferential adsorption of the incoming cation
(favorable exchange) or of the resident cation (unfavorable exchange).
The shapes of the respective dimensionless exchange curves are convex
44
CI-profiles
400
x [cm]
FIG. 14A. The effect of front
solution of the ADE using R=1.
0.
retardation on transport. Analytical
K-profiles
t = 30 d
0 25 50 75
x [cm].
FIG. 14B. The effect of
solution of the ADE using
front retardation on transport. 
Analytical
R= 38.5.
100
45
and concave. A series of simple examples is provided to illustrate how
ion exchange influences the transport of solutes. Although complete
exchange rarely occurs, usually more than two species participate in
the exchange process, and S
T 
and C
T 
are not necessarily constant; the
illustrations are representative of the exchange behavior expected.
A number of authors have discussed nonlinear exchange during
transport. Bolt (5) provided analytical solutions for nonlinear
exchange. Lai and Jurinak (17) obtained numerical solutions using an
explicit finite difference method. Concentration profiles for
unfavorable, linear and favorable exchange are shown in figure 15. The
arbitrary exchange isotherms and initial concentration profiles are
included as well. All parameter values, as well as the
Y(X)-relationships, are given in Appendix F. Because of the choice of
an initial value for X>O, a substantial amount of (adsorbed) solute was
already present at t=O. For the soil with unfavorable exchange, the
adsorbed concentration was less than the solution concentration (figure
15A), for linear exchange the two were equal (figure 15B), and for
favorable exchange the adsorbed concentration was greater than the
solute concentration. It should be noted that Bruggenwert and Kamphorst
(6) provided an extensive survey of exchange data which might serve as
a guideline for many soil systems where solute transport needs to be
modeled. Examining the profile at t=2 d, it can be concluded, as was
shown by Bolt (5) and Valocchi et al. (30), that favorable exchange
(figure 15C) tends to steepen the solute front (counteracts spreading)
and nonfavorable exchange flattens the solute front (enhances
46
Unfavorable exchange
0 4 8 12
x [cm]
FIG. 15A. The effect of nonlinear exchange 
on transport. Numerical
solution of the ADE for unfavorable exchange.
Linear exchange
0 4 8 12
x [cm]
FIG. 15B.
solution
The effect of nonlinear exchange on transport. 
Numerical
of the ADE for linear exchange.
47
Favorable exchange
> 0.5
0 0. 0
x [cm]
FIG. 15C. The effect of nonlinear exchange on transport.
solution of the ADE for favorable exchange.
Numerical
x [cm]
FIG. 16. The effect of front sharpening on transport. Numerical solution
of the ADE for favorable exchange.
48
spreading). The self sharpening 
effect of favorable exchange 
becomes
even more apparent for 
lower values of v and D (figure 
16).
Differences in solute spreading 
caused by the exchange 
process
are conveniently demonstrated 
by using a moving coordinate 
system.
Figure 17 contains profiles 
shown in figure 15 plotted 
as a function of
x-v t, where v is the 
constant (retarded) advection 
term for a
linearly exchanging solute. 
The initial concentration 
profile is given
as a dashed line. The differences 
in steepness of the front 
can clearly
be observed. One of the implications 
is that the complete displacement
of the competing cation 
is accomplished most easily 
if the incoming
cation is favorably adsorbed. 
For a speedy, incomplete displacement 
of
the resident cation, an unfavorably 
adsorbed cation is to be 
preferred.
The effect of nonlinear 
exchange on transport is 
further
illustrated using a -solute 
pulse. Consider solute movement 
through a
medium with a constant initial- 
solute concentration (X =0.05). 
For an
unfavorably exchanging 
solute (figure 18A), considerable 
spreading
occurs at the front end of 
the pulse, whereas steepening 
takes place at
the tail end of the pulse. The 
reverse can be obvserved for 
a favorably
exchanging solute (figure 
18C), whereas for a linearly 
exchanging
solute the pulse remains symmetric 
(figure 18B).
49
Time = 1 d
-10 0 10 20
30
x-v*t 
[cm]
FIG.17A. The effect of nonlinear exchange on dispersion: concentration
profiles for various exchange isotherms using a moving coordinate
system after 1 d.
Time = 2 d
* : unfavorable
0.8+ : linear
08o : favorable
0.6
0.4
0.2t=
-10 0 10 
20 30
x-v* t [cm]
FIG.17B. The effect of nonlinear exchange on dispersion: concentration
profiles for various exchange isotherms using a moving coordinate
system after 2 d.
50
Time = 3 d
-10 0 10 20 30
x-v*t [cm]
FIG.17C. The effect of nonlinear exchange 
on dispersion: concentration
profiles for various exchange isotherms 
using a moving coordinate
system after 3 d.
Time = 4 d
1.0 -
I
* unfavorable
+:"linear
0.8 o 
favorableI
0.6
0.4
0.2
t=0 \
-10 0 
10 20 
30
. x-v*t [cm]
FIG.17D. The effect of 
nonlinear exchange on dispersion: 
concentration
profiles for various 
exchange isotherms using 
a moving coordinate
system after 4 d.
51
unfavorable exchange
ld
t =1.5 d
0 10 20 30
FIG. 18A. Numerical solution
during unfavorable exchange.
x [cm]
of the ADE for a pulse input
t=2d
0 10 20 30
x [cm]
FIG. 18B. Numerical solution of the ADE for a
during linear exchange.
pulse input
52
favorable exchange
t=2d
.-- t=3d
0 10 20 30
x [cm]
FIG. 18C. Numerical solution of the ADE for a pulse input
during favorable exchange.
The Influence of Soil Layers with Different Transport
Properties on Solute Transport
Transport of solutes can be influenced by having soil layers with
different transport properties. Examination of Eq.(13) shows that the
following mechanisms can alter the transport behavior during steady
one-dimensional flow.
1) Advection
The continuity condition requires that for an incompressible
fluid V*(Ov)=O throughout the medium. At the interface of layers
k and k+1 this leads to (ev) = (ev) . In a medium not
homogeneous with respect to e, the solute will, therefore, travel
at different velocities through the various layers.
53
2) Dispersion
Different soils will likely exhibit different mixing behavior. It
was already shown how transport is influenced by soil layers
having different dispersion values. Although the pore-water
velocity varies with each layer, differences in geometry of the
pores and the size distribution 
and orientation of the soil
particles is likely to have an even larger impact on spreading
(13). This also has implications for transport perpendicular to
the direction of flow.
3) Ion exchange
The exchange properties can be of major importance for the
transport of reactive. solutes. Consequently, differences in
exchange properties will result in differences in travel time for
the solute.
To illustrate the effect of chemical and physical properties
changing from layer to layer, numerical solutions were obtained for
transport in layered media with a second- or third-type condition at
the interface and a third-type condition at the inlet of the column. We
studied transport of a solute pulse in soils with layers that have
different values of D and ST and different exchange isotherms. The
influence of different values of D for a step change was already
illustrated in figure 10 and 11.
54
D = 10 cm 
2
/d
1.0
0.8
0.6
0.4
0.2
0.0
1.0
0.8
0.6
x
0.4
0.2
0.0
0 3010 20
X (cm)
FIG. 19. The effect of layers 
with different dispersion 
coefficients on
solute transport. Numerical 
solution of the ADE for 
a two-layer medium
usng a third-type ondition 
with: (a) D= D2= 10 
cm /d, and (b) D =10
usg o t D 
21
cm /d and D =100 cm /d.
2
D= 100cm 
2
/d
t=0.5d
0.1 d
It=0.3d
mb
I II YI \ I -~ III I1
55
resulting concentration profiles for a homogeneous medium, i.e. , the
bottom layer was assigned the same properties as the top one. The
initial concentration,X is indicated with a dashed line. For a
two-layer medium, the first layer (0<x<10) has the same properties as
the homogeneous medium, while the value for D of the second layer
(10<x<30) was, somewhat arbitrarily, increased ten-fold as compared to
the homogeneous case. The difference in solute spreading between the
two layers (figure 19B) and the discontinuity at the interface as the
front passes are apparent. At t=0.5 d, the peak concentrations occurred
at the same positions in figure 19A and 19B. However, the value of the
peak concentration in figure 19A was almost twice the value in figure
19B. A layer with increased dispersion effectively reduces the maximum
solute concentration without affecting the average solute movement.
A second example illustrates transport through two consecutive
layers with different ST values. Figure 20A demonstrates how the pulse
travels through a homogeneous medium without adsorption capacity (i. e.,
S T=0). Figure 20B and 20C show results for the medium with the second
layer (10<x<30) having the adsorption capacity S =2 cmol1 /kg, for a
T c
third- and second-type condition at the interface, respectively. The
movement of the solute is retarded considerably in the second layer.
Also, notice the discontinuity at the interface for a third-type
56
T=O 
cmol/kg
x [cm]
FIG.20A. The effect of layers with different 
exchange capacities on
solute transport. Numerical solution of the ADE 
for a two-layer medium
using a third-type condition with (ST 
) = (ST 
2
=0 cmol/kg.
ST=0 cmol/kg
t=0.2
10
S =2 cmol/kg
d t=t5 d
20
x [cm]
FIG.20B. The effect of layers with 
different exchange capacities on
solute transport. Numerical solution 
of the ADE for a two-layer medium
using a third-type condition with 
(ST )= 0 and (S )=2 cmol/kg.
T1 T 2
30
0.
57
ST=O cmol/kg ST=2 cmol/kg
1.0-
t-0.2d
0.81-
0.6--
0.4--
t= 1.0 d
0.2- t=05d t 1.5 d
0.0 '
0 10 20 30
x [cm]
FIG.20C. The effect of layers with different exchange capacities on
solute transport. Numerical solution of the ADE for a two-layer medium
using a third-type condition at the inlet and a second-type condition
at the interface with (S = 0 and (S
T
) =2 cmol/kg.
TiT 2
In conclusion, a layer with increased adsorption capacity effectively
reduces not only the solute movement but also the maximum solute
concentration in the liquid phase. This is attractive to minimize the
effects of solutes which pose a threat to the quality of ground water.
A third example treats transport of a reactive solute in a medium
having layers with different exchange isotherms. Linear exchange
occurred in the first layer, whereas the exchange was either
unfavorable (figure 21A) or favorable (figure 21B) for the incoming
cation in the second layer. The solute distributions were identical in
the first layer. However, in the second layer the incoming front of the
58
1.0 Unear 
exchange
01.0
0.8
0.6
0.4
0.2
0.0
Unfavorable exchange
t=0.2d I
If\=0.
5
d
t= 1.0d t= 1.5d
1.0 
near 
exchange
F- I
2d
t = 0.5 
d
i
, t= 1.0 d
Favorable exchange
1.5 d
t= 2.0 d
10 20
X (cm)
FIG.21. The effect of layers 
with different
solute transport. Numerical 
solution o the ADE
with: (a) (dY/dX) = 1 and 
(dY/dX) =3X and (b)
1 X
-
2/31 
2
3
exchange isotherms on
for a two-layer medium
(dY/dX) =1 and (dY/dX
1
x
0.8
0.6
0.4
0.2
0.0 L.
0
30
59
pulse showed spreading (figure 21A) and steepening (figure 
21B),
respectively, for the unfavorable and favorable exchange.
The next two examples were meant to illustrate the effect 
of a
"barrier" to solute (contaminant) transport. First, consider 
a sandy
soil which has a low S
T 
and where saturated or nearly saturated flow
conditions prevail. Evidently, any solute will be transported 
at a
relatively fast pace in such a medium (figure 22A). Movement of 
the
solute can be slowed by incorporating a layer which reduces the solute
flux via a decrease in the flow term and via an increase in adsorption.
A clay soil generally possesses these characteristics; the hydraulic
conductivity is relatively low during saturated conditions and the
value of S
T 
can be significant. Somewhat arbitrary parameter values,
listed in Appendix F, were chosen based on these considerations. The
effects on the solute concentration profiles of embedding a 2-cm-thick
clay layer, in the otherwise homogeneous profile, are shown for 
a
third-type condition (figure 22B) and a second-type condition (figure
22C) at the interfaces. The presence of the clay layer was assumed 
to
result in a ten-fold reduction of v (Appendix F). The value of the
dispersion coefficient was calculated according to Eq.(2), using D /A=2' 0
cm 2/d, while a was assumed equal to 0.5 cm for the sand and 2 cm for
the clay.
It can be observed in figure 22B that the concentration in the
clay layer was considerably lower than in the sand at either interface.
This was caused by the use of two different expressions for the ADE
* *
at the interface (i.e., with different values for D and v )
60
1.0 Sand
0.8
0.6
x t 0.8 d
0.4 t=0.6d
1 0.4 d
0.2
t 0.2 d
0.0
Sand-clay-sand
1.0-12
0.8 1t6d 1 I
0.I~
0.4 t=2d\
0.2 1 1 t 6d
0.0
Sand-clay-sand
1.0
0.8I
I t=l10d
0.6
0.4
0.2 
f=6 d
61
Physically, the concentration in 
the liquid phase should be continuous
at the interface. The use of a second-type condition 
is, therefore,
preferred, although the profiles 
in figure 22B and 22C were 
very
similar in the clay layer. To 
investigate the effectiveness of the
added adsorption capacity (the 
clay layer), the locations of the 
front
were also determined for zero 
dispersion by using the appropriate
retarded velocities in each 
layer. The locations of this
"Izero-dispersion" front are symbolized with the numerals 1 to 5 Cat
c/c =0.5) for the various times listed 
in figure 22B. The solute moved
0o
slower through the sand when no 
dispersion occurred. Becaus e of
dispersion, the adsorption capacity 
of the clay layer was not fully
utilized. It is clear that the magnitude 
of the dispersion coefficient
needs to be known to accurately 
predict solute movement under these
conditions.
Under some circumstances, it 
can be anticipated that "wet"
conditions rarely exist. Solute 
transport can then be slowed 
by
including a coarse textured layer 
in the porous medium. Winograd (35)
discussed the dry site disposal of 
waste. Concentration profiles in 
a
homogeneous clay (figure 23A) 
were compared with those in 
a
clay-sand-clay. system for a third- 
and second-type condition (figure
23B and 23C, respectively) under 
unsaturated soil water conditions. The
inclusion of the sand layer was 
assumed to result in a ten-fold
62
1.0 Clay
0.8
0.6 
t=1od
X
0.4
t 6 d
0.2
t2d
0.0
1.0 Clay-sand-clay
t= d 100 d
0.8 I
0.6
X
0.4 t=60d\
0.2
= 20 
d 
now
0.0 - - II 
V hi'Bn-
1.0 Clay-sand-clay
0.8
0.6
0.2
63
of the water. Because of the lower value of the Peclet number for the
3-layer system, considerable time was needed 
for the inlet
concentration to become equal to the eluent concentration. 
Notice that
the discontinuity at the interface behaved in the 
opposite manner as
for the clay layer in the previous example.
SUMMARY AND CONCLUSIONS
Analytical solutions of the 1-D ADE, subject to a first- 
or a
third-type inlet condition, were obtained to describe transport 
in a
two-layer medium. The 1-D ADE was solved numerically 
to describe
transport of linearly and nonlinearly exchanging solutes 
in layered
and nonlayered soils with the finite difference method. The accuracy
of the numerical solutions was checked for a number of situations where
analytical solutions were available or could be derived. Both numerical
and analytical solutions were used to study transport in layered 
soils.
The solution of the ADE subject to a concentration- 
or first-type
condition at the interface leads to substantial mass balance 
error.
These errors are avoided by using a flux- or third-type condition 
or a
combined concentration/flux- or second-type condition. For a 
third-type
condition, the resident concentration is not continuous 
at the
interface. This is in contrast to the resident concentration 
obtained
with a second-type solution. However, the flux-averaged 
concentration
is continuous at the interface in both cases.
Ion exchange and stratification appeared to have major impacts
on the movement of solutes in soils. The choice of cations 
in the
eluent solution is important when the resident cation needs 
to be
64
displaced. For complete displacement of the resident solute, an eluent
cation which is favorably exchanged with the resident cation is to be
preferred. In other instances, one might be interested in a relatively
fast displacement of only part of the resident solute and an
unfavorably exchanging solute should be applied. Spreading phenomena
are then advantageous. The examples which were presented offer some
qualitative insight in what might be expected. It was shown that
nonlinear exchange can become an important factor for reactive solute
transport, which implies that the description of the exchange reaction,
over the total concentration range, deserves careful attention. It was
also found that the dispersive flux can become quite important at lower
pore water velocities. These findings stress the importance of
obtaining a good knowledge of the values of physical and chemical
parameters.
Although the flow and transport conditions were relatively
simple, the developed computer programs can easily be modified to
account for different exchange curves, varying C
T 
and S
T
, multi species
exchange, sink/source terms, etc. Nonsteady flow conditions can also be
included by solving the flow equation along with the ADE.
LITERATURE CITED
(1) ABRAMOWITZ, M. and I.A. STEGUN. 1964. Handbook 
of Mathematical
Functions. National Bureau of Standards, 
Applied Mathematics
Series 55, U.S. Government Printing Office, Washington, D.C.
(2) ABRIOLA, L.M. 1987. Modeling Contaminant 
Transport in the
Subsurface: an Interdisciplinary Challenge. Rev. 
of Geophys.
25:125-134.
(3) AL-NIAMI, A.N.S and K.R. RUSHTON. 1979. Dispersion 
in Stratified
Porous Media: Analytical Solutions. Water 
Resour. Res.
15:1044-1048.
(4) BARRY, D.A., J.-Y. PARLANGE, and J.L. STARR. 1987. 
Numerical
Analysis of the Precursor Effect. Soil Sci. 
143:309-317.
(5) BOLT, G.H. 1982. Movement of Solutes in Soils: 
Principles of
Adsorption/exchange Chromatography. In G.H. 
Bolt (Ed). Soil
Chemistry: B.Physico-chemical Models. 
Elsevier Scientific
Publishing Co. Amsterdam.
(6) BRUGGENWERT, M.G.M. and A. KAMPHORST. 1982. Survey 
of Experimental
Information on Cation Exchange in Soil Systems. In 
G.H. Bolt (Ed).
Soil Chemistry: B.Physico-chemical Models. 
Elsevier Scientific
Publishing Co. Amsterdam.
(7) BURDEN, R.L., J. DOUGLAS FAIRES, and A.C. 
REYNOLDS. 1981.
Numerical Analysis. Prindle Weber & Schmidt. 
Boston.
(8) CARSLAW, H.S. and J.C. JAEGER. 1959. Conduction 
of Heat in Solids.
Clarendon Press, Oxford.
(9) CRUMP, K.S. 1976. Numerical Inversion of Laplace 
Transforms Using
a Fourier Series Approximation. J. ACM 23:89-96.
(10) DAGAN, G. and E. BRESLER. 1985. Comment 
on "Flux-averaged and
Volume-Averaged Concentrations in Continuum Approaches 
to Solute
Transport", by J.C. Parker and M.Th. van Genuchten. 
Water Resour.
Res. 21: 1299-1300.
(11) DANCKWERTS, P.V. 1953. Continuous Flow Systems. 
Chem. Eng. Sci.
2:1-13.
(12) DANE, J.H. and F.H. MATHIS. 1981. An Adaptive 
Finite Difference
Scheme for the One-dimensional Water Flow Equation. Soil Sci. Soc.
Am. J. 45:1048-1054.
65
66
(13) FRIED, J.J. and M.A. COMBARNOUS.
Adv. Hydroci. 9:169-282.
1971. Dispersion in Porous Media.
(14) HUYAKORN, P.S. and G.F. PINDER. 1983. Computational Methods in
Subsurface Flow. Academic Press. New York.
(15) JURY, W. A. 1982. Simulation of Solute Transport Using a Transfer
Function Model. Water Resour. Res.18:363-368.
(16) KREFT, A. and A. ZUBER. 1978. On the Physical Meaning of
Dispersion Equation and its Solutions for Different Initial
Boundary Conditions. Chem. Eng. Sci. 33:1471-1480.
the
and
(17) LAI, S.-H. and J.J. JURINAK. 1972. Cation Adsorption in One-
Dimensional Flow through Soils: a Numerical Solution. Water
Resour. 
Res. 8: 99-107.
(18) LAPIDUS, L and N.R. AMUNDSON. 1952. Mathematics of Adsorption in
Beds. 6. The Effect of Longitudinal Diffusion in Ion Exchange and
Chromatographic Columns. J. Phys. Chem. 56:984-988.
(19) Leij, F.J., and J.H. Dane. 1989. Determination of Transport
Parameters from Step and Pulse Displacement of Cations and Anions.
Ala. Agr. Exp. Sta. Dept of Agronomy and Soils Series 135.
(20) OBERHETTINGER, F. and L. BADII. 1973.
Transforms. Springer Verlag.
Tables of
(21) OZI$IK, M.H. 1980. Heat Conduction. Wiley-Interscience. New York.
(22) PARKER, J.C. and M.Th. VAN GENUCHTEN. 1984. Flux-averaged and
Volume-averaged Concentrations in Continuum Approaches to Solute
Transport. Water Resour Res. 20:866-872.
(23) PERSAUD, N. and P.J. WIERENGA. 1982.
One-dimensional Cation Transport
Ion-exchange Media. Soil Sci. Soc. Am.
A Differential
in Discrete
J.46:482-490.
Model for
Homoionic
(24) PRESS, W.H., B.P. FLANNERY, S.A. TEUKOLSKY, and W.T. VETTERLING.
1986. Numerical Recipes. Cambridge Univ. Press. New York.
(25) REMSON, I., G.M. HORNBERGER, and F.J. MOLZ. 1971. Numerical
Methods in Subsurface Hydrology. Wiley-Interscience. New York.
(26) RUBIN, J. 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.
Laplace
67
(27) SELIM, H.M., R. SCHULIN, and H. FLUHLER. 1987. Transport and Ion
Exchange of Calcium and Magnesium in an Aggregated Soil. Soil Sci.
Soc. Am. J. 51:876-884.
(28) SHAMIR, U.Y. and D.R.F. HARLEMAN. 1967. Dispersion in Layered
Porous Media. J. Hydraul. Div. Proc. ASCE 93(HYS):237-260.
(29) STONE, H.L. and P.T. BRIAN. 1963. Numerical Solution of Convective
Transport Problems. AIChE J. 9:681-688.
(30) VALOCCHI, A.J., 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.
(31) VANDERGRAFT, J.S. 1983. Introduction to Numerical Computations.
Academic Press. New York.
(32) VAN GENUCHTEN, M.Th. 1978. Mass Transfer in Saturated-unsaturated
Media. I. One-dimensional Solutions. Res Report No.78-WR-11.
Princeton Univ., Princeton.
(33) VAN GENUCHTEN, M.Th. and W.J. ALVES. 1982. Analytical Solutions of
the One-dimensional Convective-dispersive Solute Transport
Equation. U.S.D.A. Technical Bulletin No. 1661.
(34) VAN GENUCHTEN, M.Th. and J.C. PARKER. 1984. Boundary Conditions
for Displacement Experiments through Short Laboratory Soil
Columns. Soil Sci. Soc. Am. J.48:703-708.
(35) WINOGRAD, I.J. 1974. Radioactive Waste Storage in the Arid Zone.
Trans. Am. Geophys. Union 55:884-894.
APPENDIX A: Analytical Solution of the ADE for the Second Layer 
of
a Porous Medium Using a First-type Condition
The solution of the ADE for the first layer was given by
Eq.(7). We will derive more general solutions by considering pulse
displacement, and different initial concentrations for both layers. For
the second layer, the problem can be formulated as follows:
ac a
2
c ac
R -= D - - v L <x<L t>0 (Al)
2 8t 2 x2 2 X 1 2
8x
C(x,0) = g L <x<L (A2a)
2 1 2
C(x,t)IxL = C(x,t)IxL t>0 (A2b)
1 1
BC
= 0 
t>O (A2c)
8x x w
where the subscripts 1 and 2 denote the first and second layer,
respectively, and g
2 
is the prescribed initial concentration for the
second layer. An illustration of the problem for K layers was given in
figure 2.
Solution of (Al), subject to (A2), is achieved with the help of
Laplace transforms. The Laplace transform of C(x,t) with respect to t,
is defined as:
2 [C(x,t)] = C(x,t) exp(-st) dt = C(x,s) (A3)
J
0
which eliminates the dependency of C on time. By applying the Laplace
transform to (Al) and (A2), we obtain:
2- v - sR gR
dC 2d d_ 2- 2 2=
_ + - 0 (A4)
2 D dx D D
dx 2 2 2
C(L ,s) = F-ec(-tos) - exp (iL ) + -- (A5a)
1 [ O S
0
1 1 S
68
69
dC = 0 
(A5b)
dx x-)w
1 v 2 sR
where X = 1 and 
C is the prescribed eluent
1 2D 2D + D o
1 1 1
concentration (cf. Eq.(5a)). A general solution of 
(A4) is:
+ g
2
C(x,s) = a exp(XA x)+ 9 exp(Xx) + -
(A6)
2 2 s
v v 2 sR
+ 2 I2 1 2
where = 2 --2 + 2 and a and are 
coefficients to be
2 2D 2D D
2 2 2
determined by the boundary conditions. By using the 
outlet condition
(A5b), it follows that a=0. For convenience, the 
superscript of X will
2
be dropped. Application of the inlet condition (ASa) 
results in:
S= [3i-ep(-tos)-!expl(XL -A
L 
) + 12]ex(-L) (A7)
Substitution of (A7) in (A6) leads to the following 
solution for
C(x,s):
-g 1 C
C(x,s) s = O 3p 1 2( L+1 2 S) ezxp( L +X 
2
-ts) +
1 
2 
21 1 2
(
2
) 1 + C2 
+ C3 + C4 
(A8)
where g=x-L . By inverting the various terms of Eq.(A8), 
we can obtain
the desired solution for C(x,t). It is convenient 
to adopt the
following notation:
1 v
1  
1 19a)
alv/4R D (A9c)
k = /RD (A9d)
2 2 2
70
This allows us to write the inverse of C as:
C
1 
(x,t) = C ep(a ak 
-k a+s + a
2
k -k
2  
+S ) =
0 1 1 2 2
(Co-gl) erp(al k
1 
+a
2
k
2
) g*h (Al0)
where C depends on x because of k
2. 
The convolution g*h, defined as
t t
g(t)*h(t) = g(t-T)h(T) dr = g(T)h(t-T) dr, 
(All)
0 0
is used to carry out the inverse Laplace transform 
according to:
f- [g(s)h(s)] = g(t)*h(t) 
(A12)
The inverse of g(s) and h(s) can be obtained with the shift 
property of
the Laplace transform and a table of Laplace transforms 
(20, 33):
g(t) = exp(-a
2
t) -1 12 exp(-k/)]= 
(A13)
1
k -2a t 
k +2a t.
1 1 )
2 12 
22
k 
k
2
+4a
2
t
h(x,t) = exp(-a
2
t) -1 exp(-k
2  
s) = 2 exp- 22(A14)
2 
2 
4t
By applying Eq.(All), the following expression for C is found:
k -2aT 
k+2a T
C 1 (x, t) 2 (C-g ) [ec+ emp(2a1k 
)
e/,c [2 x
2 (k-2a (t-T))
2
The inverse of C follows from C
1
, using the properties of the Laplace
transform for a step function:
71
0 O<tst
I 0o
C2(xt) = C(A16)
S-g C (x,t-t) t>t
Co g11 o o
0 1
The inverse of C
3 
is determined according to:
C3(x,t) = 1[91 2]e(a k -k +s) =
k-2a t k2+2a t
2 1 22
(gg
2
){e4 2 + ecrp(2a
2
k
2
) 
4 t
]} A7
whereas C
4
=g
2
. The solution of the problem is:
C(x,t) = C
l(x
,t) + C2(x,t) + C3(x,t) + C
4  
(A18)
For a step-type displacement, i.e., t o, and g=g
2
=C. we obtain
0 1 2 1
C(x,t) 
= o 
1 2 
1 1 
+ ecp(2alk 
)enc
r -3/2 (k -2a Ct-/)
2
X(t-T)-3 ej- 22t-T) dT + C. (A19)
L-L 4(t-r) 1
The integrals in Eqs.(A15) and (A19) can be evaluated via standard
numerical methods.
The approach outlined in this Appendix is similar to that
originally described by Shamir and Harleman (28), which we were 
not
aware of at the time we derived (A19). It should be noted that their
Eq.(18) contains a misprint: the term 7D
L 
(t-7) should be read as
2
3
rD (t-) . Shamir and Harleman (28) extended their analysis to a
2
medium of K layers to obtain the concentration at the end of the K-th
layer. However, it can be verified that G(x,s) for an arbitrary
72
position in this layer is given by:
C -C K-1iC.
C(x,s) = 0 1 explltexp(A )+ 1 
(A20)
S k=1
where t is the thickness of the k-th layer and 
=x-L . Equation
k K 
K-1
(A20) can be inverted using the convolution theorem 
while making use of
the results for layer K-i. It is apparent from 
(A20) that the order in
which the layers appear does not affect the concentration in the final
layer K, as was already concluded by Shamir and Harleman 
(28).
APPENDIX B: Analytical solution of the ADE for the 
Second Layer of a
Porous Medium Using a Third-type Condition
In this Appendix we seek a solution of the ADE for the second
layer, using the principle of continuity of flux rather 
than
concentration, as in Appendix A. The use of a flux- or third-type
condition is preferred over a first-type condition for volume-averaged
concentrations because it leads to conservation of mass (34). 
The
solution in the first layer was given by Eq.(8). The problem 
can be
stated as follows:
ac a 
2
c ac
8C 82C BC
R - D - v - L <x<L t>0 (Bl)
2 8t 2 2 2ax 1 2
8x
C(x,0) = g L <x<L (B2a)
2 1 2
8C 8C
(-0 D 8c +ev C) = (-e D +vC)t> (B2b)
1 18x 11 xL 2 2 8x 2 2 xL t> B2b)
1
8C
S0 
t>0 (B2c)
8x x-)
Again, the method of Laplace transforms is employed to obtain a
solution of (Bl) subject to (B2). Transformation of (Bl) and (B2) leads
to:
2- v - sR gR
2dC 2 + - 0 (B3)
2 D dx D D
dx 2 2 2
dC dC
(-0 D + v C) = (-eD dC + v)(B4a)
1 1 dx 11 xL 22 dx 22 x4L
dC
d = 0 
(B4b)
dx x-
The Laplace transform of the concentration in the first layer is:
v -g C
0  
g
C(x,s) = ep(-ts) JcrpXx) + -- (BS)
v-DA - - 1 s
1 1 1
where A was defined in Appendix A.
1
73
74
The general solution of (B3) subject to (B4) 
is:
+ - 2
C(x,s) = a exp(A x) + 3 exp(A x) + -- 
(B6)
2 2 S
where a, 13, 
2 
and X2 were defined before. Boundary condition (B4b)
requires that a equals zero. An expression for 1 
is found by using the
inlet condition and assuming steady flow 
(i.e., V.(ve)=O):
V 2 g-g 2 C]
S= v 
2  
+ - 1-eC(-tos)exp(AL
) 
exp(-A L) (B7)
v -D X s so0 1 1 2 
1
2 2 2W
Substitution of Eq.(B7) into Eq.(B6) yields the solution 
of the
ordinary differential equation (B3). Using the 
notation according to
Eq.(A9) we obtain:
2 2
2C" 9s2)-- -- --
x exp(a -k -k a +s) + 2 = C (x,t) + C2(x t) + C3(xt) 
+ C
4  
(BS)
2 2 2 S 1 2' 
3' 4
To obtain C(x,t) we proceed by taking the inverse in a similar
manner as in Appendix A.
1 2 1 1 2 2
sa + /a +s
The inverse part of C is determined by using the 
expression for the
Laplace transform of the convolution of two functions 
g(t) and h(t),
for completeness, we also include the expressions 
for g(s) and h(s):
gs) = 1esp4 (a+s) 
(B10)
k -2a t k +2a t
g(t) = e { (-alk
1
)ere 1 ]+ e1p ~a k
1
)e [ ' (Bi)
2 1 11t /
75
x, s)=
h(x, t)=
22
a+p a- 
/s+
2 2
222
(4a t 2+k)k+2a t
L 4t 
]- 
2
ecp 
2 2
)eq{
(B12)
(B13)
After applying Eq. (All) we obtain:
t -2a T 
k +2aT
C (x,t)=a (C- {e_9[' +eaP2a k)e >21
(B14)
P-12 (k 2-2a 2Ct-z))
2  
k
2 
+2a 
2
(t T-r)
X .6(r(t-') e" [ 2 2 i-a 2 e4p(2a k ) e/,&c[J1l
L 4(t-14 4(t--T)
Using the properties of the Laplace transform:
O<t~st
0
C-g 
2, 
o)
01
(B15)
t>t
c C x, t) =2a 
2(g1-g2 ),ehp(a 
2k2
-
2
ro/16a t
-(g -g ) -(
21 2 II2
[earpt-k
2  
a +s]=
s (a
2
+ a +s]
,e(k2-2a 
2t)
2
,+ 
-k2 2a2t-
24tJ +eck 2]+
(l-2a k+4a-k +2a t
2+2a +a
2 
t) eocp(2a k
2
) e 2 
2]
4't
The solution of the problem is:
C~x, t)= C C x, t) 
+ c2 (x, t) + c 3Cxt 
t) + c4
For a step-type displacement, i.e., t --), and g =g =C. we obtain
CB16)
(B17)
(B18)
76
-2ak +2aT
C ( x , t ) 
= a 2 ( C 
o- C .) 
x 
4 q 
"4 T + 
a P ( a k 
e & / 
4
x{&ot -1/2e (k42-2 2(t-)) 
1
2k+2aa
2
(t-)
X- 7rJa) xpx
2 
exrpa 
2 2 
LA ~T+
6 4(-T 
4 Ct -T)
+ C.
1
(B19)
For the K-th layer, the solution in the Laplace domain is given
by:
v
C(x, 
S) 
-
c -C. K- c.
0 x 1 PL xE j exp(X )+ 
5
(v K-D KX K) (kl e.
Again, it appears that the order of layers 
between k=2 and K-1 is
irrelevant for the concentration 
in the K-th layer (LK- <x<L )
(B17)
APPENDIX C. Analytical Solution of the ADE for a Two-layer 
Porous
Medium Using a Third-Type Condition at the Inlet and a 
Combined
First- and Third-type Condition at the Interface
Physically, it is realistic that both concentration and flux 
are
continuous at the interface. This problem 
is similar to that of heat
flow in a composite medium where the heat flux across 
the interface is
continuous and intimate contact between different 
layers at the
interface, i.e., no thermal resistance, ensures 
continuity in
temperature (8). The application of continuity in both solute 
flux and
concentration at the inlet is not feasible, therefore a third-type
condition will be used at the inlet.
In this appendix, we seek a solution 
for the ADE in a medium
consisting of a finite layer (-L<x<O) and a semi-infinite 
layer (x>O)
(cf. 8, p.319). In contrast with the two previously 
investigated
problems, the concentration in the first layer depends 
on the
properties of the second layer. The problem can be stated as follows:
aC a2C 8C
RD- v -L<x<O t>O (Cla)
i at ax
2  
iax
22
8C 82C 8c
R 2 D 2 2
R D - -v - x>O t>O 
(Cib)
2 3t 2 2 2 ax
C (xO) = g -L<x<O 
(C2a)
1
C (x,O) = g2 x>Q (C2b)
8C rv C O<t-<t
(vC -D i) - v f 10 oC c
i i i ax xg-L 'L [ Cc
i 0 t0
77
78
C C t>O (C2d)
i Xto 2 X,
9C 8C
1 2
evC -D -=evC-eD t>0 (C2e)
111-1 8x x 222 22 a X0
ac
2 -=0 t>O (C2f)
8x x~oo
The Laplace transforms of (C1) and (C2) are:
dC v dC sR gR
1 1 1 1 11
1 + gO(C3a)
2 D dx D 1 D
dx 1 1 1
dC v dC sR gR
2 2 2 2 22
C + - 0 (C3b)
2 D dx D 2 D
dx 2 2 2
dC C
1 O 1P(tS
(vC- D 1 ) -vif S 1-exp(-t s) t>0 (C4a)
1 11 dx x -L 1 s L o
C xt= C 0t>0 (C4b)
dC dC
1 2
8D -x - D -t>0 (C4c)
1 1 dx x 0 2 2 dx x,0
dC
2 =0 
t>0 (C4d)
dx x-)
Note that Eq.(C4c) follows from Eq.(C2e) assuming steady flow and
continuity in concentration. The general solution of Eq.(C3) can be
written as:
C (x,s) = a exp(Ax) +t3 eXp(-Xx) exp 2D S+ 
(C5a)
v 2 x 92
C (x,s) = 7 exp(A
2
x) + p exp(-X
2
x)] exp2D + - (C5b)
2
where a, f3, 7' and ji are coefficients determined by the boundary
dv 2 sR '
conditions and =
Application of Eq. (C4d) gives y=0, whereas 
the use of CC4b) shows
79
that M=a+3. During steady flow, the following relationship between a
and 1 can be established with Eq. (C4c):
av DA +v DA = vD -v D (Ca)
or using q i=(D iAi)/v
a f3 (C6b)
q1+q2
whereas application of Eq. (C4a) yields:
a[-q exp(-AL) + 13+q exp(AXL) = C(-ts) exp (C7)
Explicit expressions for a and f3 can be found by using both
Eq.(C6) and (C7). It is convenient to transform the position coordinate
such that the interface is located at x=O. The resulting equations for
C and C are:
1 2
gl 
Co 
-g
1
-e
xp(
-tos )
C (x,s) = - +
Ss -q)(q -q 2)exp(-AIL) + (-+q )(ql+q ) exp(X L)
V X vx]}
s (q
1
-q
2
)e Xp- + AX(x-L) + (q +q
2
)rpX. - X(x-L) 
(C8)
2q [C 9 -eXP(-t s) ]eXP
- g
2  
2q1
C (xs) =--+x
2 s 1-q )(q 
-q2)exp(-XL) + ( 1+q 
)(q +q )exp(AXL)
v (x-L) (x-L) (C9)
sexp 
2D 
2
Inversion as described in Appendix A and B or application of the
inversion theorem or the expansion procedure (8) did not seem feasible.
Concentration profiles were therefore obtained via numerical inversion
of the Laplace transform according to the algorithm by Crump (9). As an
approximate analytical solution, Eq. (B14) can be used.
APPENDIX D. Numerical Integration with the Gauss Chebyshev Formula
Numerical integration can be performed in many ways (7). We
attempted Gaussian quadrature and the Gauss-Chebyshev method. The
latter one gave more accurate results at the end points and was easier
to code. In particular, the flexibility in the number of points is
attractive.
In order to evaluate the integrals appearing in Eq.(11) to (13),
the integration variable T, varying between 0 and t, needs to be
transformed to a variable n varying between -1 and +1. This was
accomplished by the linear transformation = 12T-t)] For a function
f(T) the integration changes as follows:
t 
+1
Jf(z) dr =J2 f( 2 ~t)diq (D1)
0
The Gauss-Chebyshev formula yields (31):
f (=) d 
1 a
1
f n
1
) + a
2
f (
2
) + ...... + anf(mn) 
+ ET 
(D2)
where E
T 
is the truncation error, and a and 71 are the coefficients
T
and points, respectively, which are defined for n points as:
= 
Coa( 
(
2k-l ) i)
k 2n
k=1,2,... 
,n 
(D3)
a =a
2  
ak= .. =an =n
In our case, the transformation of variables is straightforward.
Because the denominator does not contain a term V-{ ), both
denominator and numerator of the transformed function are multiplied
with this factor (cf. Appendix G).
80
APPENDIX E. Expressions for a
i
, b., c
i
, and f. to Solve
the ADE for a Non-layered and a Layered Soil
Non-layered soil
inlet condition (x=O) j>O
c =0 b = 1, f = X first-type
1 1 1 o
S= -2AtG
b = R + 2At(G + 2H(1 + H third-type
1 1 G
f, 
= 
(2R-b ) X + 2AtG X + 8AtH(1 + ) X
1 1 1 2 G o
a. =-G - H
1 .
R
b. =2G +
1 t i=2,3,.. ,n-1
c. = -G + H jaO
1
f. = -a.X + (b.-4G) X9 -c X
1 1 i-1 1 1 i i+1
outlet condition (x=L) j>O
a = -1,
R
n
b 
= 1 
+
n 2AtG
f = 
X 
+ (b 
-2) Xi
n n-1 n n
initial condition (t=O)
X i=1
X = first-type
SX. i>1
1
24H
= {4H2
=  
third-type
Note that G=D/(2Ax
2
) and H=v/(4Ax)
81
81
APPENDIX E. Continued
Layered soil
inlet condition (x=O) j>O,k=l
a = c = O, b = i, f = X first-type
1 1 o
al = , = -2AtG
H
b = R + 2At(G + 2H (1 + I)) third-type
f 
= 
(
2
R -b ) +2AtG X + 8AtH (1 + 1) X
1 i1 1 1 1 2 1 G o
1
coefficients for nodes not loacated at the boundaries or interface
a. =-G - H
1 k k
Re j=0, 1,2,....
1k
b. = 2G +
1 k At
k=1,2,....,K
c. = -G + H
1 k k
f. = -a.X + (b .-4G ) X - c.
1 1 1-1 i k 1 1 1+1
interface condition (x=L) j>O, k=1,2,..,K-1
k
a. =-w G -w G -w H-w H
1 k k k+1 k+1 k k k+1 k+1
wR. +w R.
k 1k k+1 1,k+1 + 2wG +2w 
G
i At k k k+1 k+1 first-type
c. = -w G -w G +wH +w H
1 k k k+l k+1 k k k+l k+1
f. = -a.X _ + (b.+4w G +4w G ) X -c.X
1 1 1 1 k k k+1 k+1 1 1 
1+1
2xv 
2 (D
k =2Ax k k Ax 
(
Vk+1
ai 
= -2G1
R
b 1k + 
2G + r 
(H-G )
i At k k k k third-type
c. = (k(Gk-Hk)
= k 1-1 At k k k k i+1 k+k 1
+ [k l (Gk-Hk )]X
k+1l k k +1
82
APPENDIX E. Continued.
a. = -D
1 k k
b. =e D +e D
1 k k k+1 k+1 second-type
c. =- D
1 k+1 k+1
f. = -a.X. - b.X. -c.X.
1 1 i-1 1 1 1 1+1
a. = - (G +H )
1 k+1 k+1 k+1
R.
b. = 9,k+1 + 
2
G + (G +H )
1 At k+1 k+1 k+1 k+1
c. =-2G third-type
1 k+1
f. = [q( (G +H)XL + (1-q) (G +H )lXi +
1 k k+1 k+1lk+1j11 L k k+1 k+1 k+1 1
RJ
E1 k+1G
2  
- (G 
+H )]X 
+ 2G X
At k+1 k+1 k+1 k+1 j 1, k+1 1+1
outlet condition (x=L ) j>0
K
a =-1, c = 0
n n
RJ
b =1 
+ 
n,
n 2AtG
K
f = X
j  
+ (bn-2) X
j
n n-i n n
initial condition (t=0) j=0 k=1,2,....,K
0 X0
X= first-type
i X. 
i>1
1
24H
k
Xc xi=1
X.
0 
= { 
24
Hk+
25
G it 
y
KOk k third-type
1
S X. i>1
1
2
Note that G =D /(2Ax ) and H =v /(4Ax)
k k k
APPENDIX F. List of Parameter Values Used for 
Simulations
Fig. k v D 8 p S 
dY X. t
condition 
T dX 
1 0
cm/d cm
2
/d g/cm
3 
cmol /kg d
3
a  
1 50 50 0.0 0.05 
00
3 1 50 50 .0.0 0.05 c
4 1 50 0 
0.0 0.05 00
5
a  
1 50 50 0.0 0.05 
00
5
b  
3 50 50 0.0 
0.05 o
6
a  
1 1 
50 50 
0.0 
0.05 
0
2 50 50 0.0 
0.05 00
6b 3 1 50 50 0.4 
0.0 0.05 0o
2 50 50 0.4 0.0 
0.05 0o
7
a  
1 1 50 50 0.0 0.05 
00
2 50 50 0.0 0.05 
00
7
b  
3 1 50 50 0.4 0.0 
0.05 00
2 50 50 0.4 0.0 
0.05 00
c  
1 50 50 0.4 0.0 0.05 00
2 50 50 0.4 0.0 
0.05 00
8
a 1 1 25 12.5 0.4 0.0 
0.05 00
2 100 100 0.1 0.0 
0.05 0o
8
b  
3 1 25 12.5 0.4 0.0 0.05 
0o
2 100 100 0.1 0.0 
0.05 o
8
c  
2 1 
25 12.5 
0.4 
0.0 
0.05 0o
2 100 100 0.1 0.0 0.05 
0
a 1 100 100 0.1 0.0 
0.05 00
9 1 2 25 12.5 0.4 0.0 
0.05 0o
b 1 100 100 0.1 0.0 
0.05 0
2 25 12.5 0.4 0.0 0.05 
0
9
c  
2 1 100 100 0.1 0.0 
0.05 0
2 25 12.5 0.4 0.0 0.05 
00
a  
1 1 10 10 0.0 0.05 
C
2 10 5 0.0 
0.05 0
b 3 1 10 10 0.4 0.0 
0.05 0o
10 3 2 
10 5 0.4 
0.0 
0.05 o
a 1 10 40 0.0 
0.05 0
2 10 5 0.0 
0.05 0
b 1 10 40 0.4 0.0 
0.05 0
3 2 10 5 0.4 0.0 0.05
12 1 20 20 0. 5 0. 0 
0.05 ci
2 50 20 0. 2 0. 0 
0.05 o
b 1 20 20 0. 5 0. 0 
0.05 Co
12 ~ 2 50 20 0. 2 0. 0 
0. 05
83
84
APPENDIX F. Continued.
Fig. k v D e p s dY 
X. t
condition 
T dX 
1 o
cm/d cm /d g/cm 
cmol /kg
13 a1
2
3b 
31
13 2
14
a 3
14 3
a
15 3
15 3
15 
c 3
16 3
17a3
17b3
17 3
18a3
b
18 3
18 
c 3
19 
a 
3 
1
2
1b 
31
19 ~ 2
2
2b 
31
20 3 2
a 1
21 
3 
2
21 b31
2
50 20
20 20
50 20
20 20
50 50
50 50
50 50
50 50
50 50
0.2
0.5
0.2
0.5
0.0
0.0
0.0
0.0
0.0
0.4 1.5 10.0
0.4 1.4 10.0
0.4 1.4 10.0
0.4 1.4 10.0
25 10 0.4 1.4 10.0
50 20
50 20
50 20
50 20
50 20
50 20
50 10
50 10
50 10
50 100
50 20
50 20
50 20
50 20
50 20
50 20
50 20
50 20
0.4 1.4 10.0
0.4 1.4 10.0
0.4 1.4 10.0
0.4 1.4 1.0
0.4 1.4 1.0
0.4 1.4 1.0
0.4 0.0
0.4 0.0
0.4 0.0
0.4 0.0
0.4 0.0
0.4 0.0
0.4 1.4 0.0
0.4 1.4 2.0
0.4 1.4 1.0
0.4 1.4 1.0
0.4 1.4 1.0
0.4 1.4 1.0
0.05 00
0.05 00
0.05 00
0.05 00
0.05 00
1 0.05 C0
3X
2  
shown C0
1 00
-x s 0
3
1 x-2/3
3
3X2
1 w
1 -2/3
-x 
n
00
00
3X
2  
0.05 0.4
1 0.05 0.4
1 -2/3
X- o .o5 0.4
0.05 0.1
0.05 0.1
0.05 0.1
0.05 0.1
0.05 0.2
0.05 0.2
1 0.05 0.2
1 0.05 0.2
1 0.05 0.2
3X
2  
0.05 0.2
1 0.05 0.2
1 -2/3
-X 0.05 0.2
LV Y vv uv v ~ r r
85
APPENDIX F. Continued.
dY
Fig. k v D p ST dX X. t
condition
cm/d cm
2
/d g/cm
3
cmol /kg
1
22a 3 2
3
1
b
22 3 2
3
1
2 3
a 3 2
3
b 1
23 3 2
3
100 52 0.4 1.3 1.0
100 52 0.4 1.3 1.0
100 52 0.4 1.3 1.0
10 7 0.4 1.3 1.0
8 18 0.5 1.3 5.0
10 7 0.4 1.3 1.0
10 7 0.2 1.3 1.0
10 7 0.2 1.3 1.0
10 7 0.2 1.3 1.0
1 4 0.2 1.3 1.0
0.5 2.5 0.4 1.3 0.0
1 4 0.2 1.3 1.0
0.05 c
0.05 o
0.05 00
0.05 00
0.05 00
0.05 o00
0.05 00
0.05 00
0.05 00
0.05 00
0.05 00
0.05 00
d
// In 11 17 1
APPENDIX G. Listing of Computer Programs
This Appendix contains the various WATFIV programs for 
the
IBM-mainframe used in this study. The values of the input parameters
are entered as program parameters in the declaration section, 
which is
followed by a brief explanation of these variables. Initial
concentration profiles can be entered as input variables (cf. HOMSOIL).
Any set of consisitent units can be used. The output format is
specified in the program. The structure is similar for all programs and
should be self explanatory.
FISEC : Analytical solution of the ADE for the second layer of
a two-layer medium using a first-type condition
THSEC : Analytical solution of the ADE for the second layer of
a two-layer medium using a third-type condition
HOMSOIL : Numerical solution of the ADE for a homogeneous medium
using a first- or a third-type condition
VARTIM Numerical solution of the ADE for a homogeneous medium
using a first- or a third-type condition with a
variable time step
LAY1 : Numerical solution of the ADE for a layered medium using
a first-type condition
LAY3 Numerical solution of the ADE for a layered medium using
a third-type condition
86
Listing of FISEC
//FIG1OA JOB (AYL59FL,124),'FEIKE LEIJ',MSGCLASS=P,
1/ NOTIFY=AYL59FL,MSGLEVEL=(1,1),REGION1O24KTIME(1,59)
/*ROUTE PRINT RMT4
/*JOBPARM LINES=10
//STEP1 EXEC WATFIV
//WATFIV.SYSIN DD *
/JOB FEIKE,TIME=(5) ,PAGES=60
C
INTEGER NUM/7/,I,J,UP/100/
REAL R/1.O/,POS(7)4VELO1/1O./,VELO2/1O./,DIFO1/100./,DIFO2/1O./
REAL A1,A2,C(30) ,F,H,M,N,K1,K2(25) ,DELPOS/1./,DUM
REAL ROOT(100),COEF,SUM,CON(12),DIMCON(12),X(12)
REAL NO/O.O/,CIN/.OOO5/,CZERO/.O1/,FACT,TEL1(100),TEL2(100)
REAL NOEM2(100),NOEM3(100),LEN1/6./,REPOS(25)
REAL EXF,FACT4(100),B(100),D(100),E(100),G(100),NOEM1(100)
REAL BOUND/170./,LOEF,DELTIM/O.1/,TIMMAX/4.O/
C
VARIABLES TO
NUM
VEL01
VEL02
DIF01
DIF02
CON
DIMCON
x
DELPOS
CIN
CZERO
LENi
DELTIM
.TIMMAX
BE SPECIFIED
NUMBER OF POSITIONS IN SECOND LAYER
PORE WATER VELOCITY IN FIRST LAYER
PORE WATER VELOCITY IN SECOND LAYER
DISPERSION COEFFICIENT IN FIRST LAYER
DISPERSION COEFFICIENT IN SECOND LAYER
CONCENTRATION
(C CI)/(COCI)
C/CO
SPACE INCREMENT FOR ANALYTICAL SOLUTION
INITIAL CONCENTRATION
ELUENT CONCENTRATION
POSTION OF INTERFACE
TIME INCREMENT FOR ANALYTICAL SOLUTION
MAXIMUM TIME FOR ANALYTICAL SOLUTION
L/T
L/T
L2/T
L2/T
M/L3
L
M/L3
M/ L3
L
T
T
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
c
87
PRINT,'* F I S E C
PRINT, '* 
A*I
PRINT,'* ANALYTICAL SOLUTION OF THE ADVECTION-DISPERSION 
'
PRINTI* EQUATION FOR 1-D TRANSPORT IN THE 2-ND LAYER OF
PRINT,'* A SEMI-INFINITE POROUS MEDIUM *
PRINT,'* 1-ST TYPE BOUNDARY CONDITION
PRINT,'* *
INITIALIZE SOME BASIC VARIABLES
A1=VELO1/(2.O*SQRT(DIFO1))
A2=VELO2/ (2. O*SQRT (DIFo2))
K1=LEN1/(SQRT(DIFO1))
POS(1 )=LEN1
REPOS (1) =0 .0
K2(1)=0.O
DO 10 1=2,NUM
POS(I)=POS(I-1)+DELPOS
REPOS(I)=P05(I)-LEN1
K2(I)=REPOS(I)/(SQRT(DIFO2)).
10 CONTINUE
C
PRINT,'
PRINT,' POSITIONS'
PRINT 113,'TIME',' ',(POS(I),I1,NUM)
113 FORMAt(l' ',A7,.A15,7F7.3)
PRINT,'
C
C INITIALIZE ROOTS AND COEFFICIENTS FOR GAUSS-CHEBYSHEV
C
DO 23 J=1,UP
DUM=((2*J-1)*1.57O96327)/FLOAT(UP)
ROOT(J)=C05(DUM)
23 CONTINUE
C
COEF=3. 141592654/FLOAT(UP)
C
TIME=DELTIM
WHILE(TIME.LE.TIMMAX) DO
C
C INITIALIZE VARIABLES FOR COMPLEMENTARY ERROR AND EXPONENTIAL FUNCTION4S
.C
DO 20 J=1,UP
TEL1(J)=A1*TIME*(ROOT(J)+1 .0)
TEL2(J)=A2*TIME*(1 .O-ROOT(J))
NOEMi (J)=2. O*TIME*(ROOT(J)+1)
NOEM2(J)=0.25*NOEM1(J)*S.QRT(O.25*NOEM1(J))
NOEM3(J)=SQRT(2.*TIME*(1.-ROOT(J)))
FACT4(J)=SQRT(1 .O-ROOT(J)*ROOT(J))
20 CONTINUE
C
C INITIALIZE POSITION DEPENDENT VARIABLES
C
DO 99 1=1,NUM
C.(I)-2*A2*K2(I)
DO 30 J1,fUP
B(J)=-((K1-TEL1(J) )**2)/NOEM1(J)
D(J)=(K2(I) -TEL2(J) )/NOEM3(J)
E(J)1 .O/NOEM2(J)
G(J)=(K2(I)+TEL2(J) )/NOEM3(J)
88
DO 40 J=1,UP
IF (B(J).LE.-BOUND) THEN
LOEFO0.O
ELSE
LOEF=EXP(B(J))
ENDIF
SUM=SUM+COEF*FACT4(J)*E(J)*LOEF*(EXF(NO,,D(J))+EXF(C(I),,G(J)))
40 CONTINUE
C
CON ( I )SUM*FACT+CIN
DIMCON(I)=(CON(I)-CIN)/(CZE'RO-CIN)
X(I)=CON(I)/CZERO
99 CONTINUE
C
PRINT 114,TIME, 'C-CI/CO-CI' ,(DIMCON(I),11I,NUM)
114 FORMAT(' 'IF7.3,AI17F7.3)
C
PRINT 113,' ','C/CO',(X(I),I=1,NUM)
C
PRINT,'1
C
TIME=TIME+DELTIM
ENDWHILE
C
STOP
END
C
REAL FUNCTION EXF(Q,Z)
C
REAL Q,Z,T,STUP1,STUP2,XX,YY,QQ
C
EXFO0.O
STUP1=ABS (Q)
C
IF (STUP1.GT.-170 .AND. Z.LE.O.O) GO TO 40
IF(Z.NE.O.) GO TO 31
EXF=EXP (Q)
GO TO 40
31 QQ=Q-Z*Z
STUP2-ABS (QQ)
IF((STUP2.GT.170) .AND. (Z.GT..0O)) GO TO 40
IF(QQ.LT.-170) GO TO 34
XX=ABS (Z)
IF (XX .GT. 3.0) GO TO 32
89
C
34 IF(Z.LT.O.O)THEN
EXF2. *EXP(Q)-EXF
ENDIF
C
40 CONTINUE
C
RETURN
END
/GO
90
Listing of THSEC
//FIG1OB JOB (AYL59FL,124),'FEIKE LEIJ',MSGCLASS=P,
// NOTIFY=AYL59FL,MSGLEVEL=(1,1),REGION=1024K,TIME=(0,59)
/*ROUTE PRINT RMT4
/*JOBPARM LINES=10
//STEP1 EXEC WATFIV
//WATFIV.SYSIN DD *
/JOB FEIKE,TIME=(05),PAGES=60
C
INTEGER NUM/7/,I,J,UP/100/
REAL /1./,POS(7),VEL01/1O./,VELO2/1O./,DIFO1/100./,DIFO2/10./
REAL TEL2(100),A1,A2,C,G(100),F(25),H,M,N,K1,K2(25),DELPOS/1.O/
REAL ROOT(100),COEF,TERM1,TERM2,CON(25),DIMCON(25),X(25)
REAL DUMMY,CIN/.0005/,CZERO/.01/,FACT1,FACT2,TEL1(100)
REAL NO/0.0/,NOEM1(100),NOEM2(100),LEN1/6./,REPOS(25),TIME
REAL DUM,LOEF,CHECK,EXFFACT3(100),DUMMY2,B(100),D(100),E(100)
REAL NOEM3(100),NOEM4(100),BOUND/150./,DELTIM/.1/,TIMMAX/4.0/
REAL GG(100),VMC1/0.5/,VMC2/0.5/
C
VARIABLES TO
NUM'
VELO1 :
VELO2
DIFO1
DIFO2
CIN
CZERO
LEN1
DELTIM
DELPOS
TIMMAX
VMC1
VMC2
BE SPECIFIED
NUMBER OF POINTS IN SECOND LAYER
PORE WATER VELOCITY IN FIRST LAYER
PORE WATER VELOCITY IN SECOND LAYER
DISPERSION COEFFICIENT IN FIRST LAYER
DISPERSION COEFFICIENT IN SECOND LAYER
INITIAL CONCENTRATION
ELUENT CONCENTRATION
POSITION OFINTERFACE
TIME STEP FOR ANALYTICAL SOLUTION
SPACE STEP FOR ANALYTICAL SOLUTION
MAXIMUM TIME FOR ANALYTICAL SOLUTION
VOLUMETRIC WATER CONTENT IN FIRST LAYER
VOLUMETRIC WATER CONTENT IN SECOND LAYER
L/T
L/T
L2/T
L2/T
M/L3
M/L3
L
T
L
T
L3/L3
L3/L3
PRINT,'***************************************************
PRINT,'* T H S E C
PRINT,'*
PRINT,'* ANALYTICAL SOLUTION OF THE ADVECTION-DISPERSION
PRINT,'* EQUATION FOR 1-D TRANSPORT 
IN THE 2-ND LAYER OF
PRINT,'* A SEMI-INFINITE 
POROUS MEDIUM
PRINT,'* 3-RD TYPE BOUNDARY 
CONDITION
PRINT,'*
TIME=DELTIM
INITIALIZE SOME BASIC VARIABLES
Al=VEL01/(2.0*SQRT(DIFO1))
A2=VELO2/(2.0*SQRT(DIFO2))
Kl=LEN1/(SQRT(DIFO1))
POS(1)=LEN1
91
C
C
C
C
C
C
C
C
C
C
C
C
C
C
c
C
C
C
REPOS (1)=0 .0
K2(1)0 .0
DO 10 12,1NUM
Pos(I)=POS(I-1)+DELPOS
REPOS(I)=POS(I)-LEN1
K2(I)=REPOS(I)/(SQRT(DIFO2))
10 CONTINUE
C
PRINT,,'DEPTH OF FIRST 
LAYER ',LEN1
PRINT,'DEPTH OF SECOND 
LAYER ',REPOS(NUM)
C
PRINT,'
PRINT,' 
POSITIONS.'
PRINT 113,'TIME',' 
',(POS(I),I=1,NUM)
113 FORMAT(' ',A7,A15,7F7.3)
PRINT,'
C
C INITIALIZE ROOTS 
AND COEFFICIENTS GAUSS-CHEBYSHEV
C
DO 23 J1,jUP
DUM=((2*J-1)*1.57O96327)/FLOAT(UP)
ROOT(J)=C05(DUM)
23 CONTINUE
C
COEF=3. 141592654/FLOAT(UP)
C
WHILE (TIME .LE .TIMMAX) 
DO
C
C INITIALIZE VARIABLES 
FOR COMPLEMENTARY 
ERROR AND EXPONENTIAL 
FUNCTIONS
C
DO 20 J=1,UP
TELl (J)=A1*TIME*(ROOT(J)+1)
TEL2(J)=A2*TIME*(1 
.O-ROOT(J))
NOEM1(J)=SQRT(2.O*TIME*(ROOT(J)+1))
NOEM2(J)=2.*TIME*(1.-ROOT(J))
NOEM3(J)=SQRT(NOEM2(J))
DUMMYO .785398163*NOEM2(J)
NOEM4(J)=SQRT(DUMMY)
FACT3(J)=SQRT(1.O-ROOT(J)*ROOT(J))
20 CONTINUE
C
C INITIALIZE POSITION 
DEPENDENT VARIABLES
C
DO 99 1=1,NUM
92
G(J)=-( (K2(I)-TEL2(J) )**2)/NOEM2(J)
30 CONTINUE
C
FACT2=(CIN*(VELO1*(VMC1,/VMC2)-VELO2))/(SQRT(DIFO2))
FACT1=O.2S*VMC1*VELO1*(CZERO-CIN)*TIME*COEF/(-VMC2*SQRT(DIFO2))
C
H=-((K2(I)-2.O*A2*TIME)**2)/(4.O*TIME)
MO.,5*(K2(I)-2.O*A2*TIME)/(SQRT(TIME)).
N=0.5*(K2(I)+2.O*A2*TIME)/(SQRT(TIME))
C
C EVALUATE THE INTEGRAL WITH GAUSS-CHEBYSHEV
C
TERM1O0.O
DO 40 J=1,UP
IF(G(J).LT.-BOUND)THEN
LOEF=0.0
ELSE
LOEF=EXP(G(J))
ENDIF
C
TERM1=TERM1+FACT3 (J)* (EXF (NO ,B (J) )+EXF (C, D(J) 
)) *
$ (E(J)*LOEF-A2*EXF(F(I),GG(J)))
40 CONTINUE
C
TERM1=FACT1t*TERM1
C
C EVALUATE SECOND PART OF THE ANALYTICAL 
SOLUTION
C
DUMMY=TIME/3.141592654
D UMMY=S QRT (D UMMY)
IF(H.LE. -BOUND) THEN
LOEF=O.
ELSE
LOEF=EXP (H.)
ENDIF
C
TERM2=FACT2*(DUMMY*LOEF+( (O.25AEXF(NO,M) 
)/A2)
$ -((O.25*(1.+F(I)+4.*A2*A2*TIME)*EXF(F(I) 
,N))/A2))
C
CON (I )=TERM1+TERM2+CIN
DIMCON(I)=(CON(I)-CIN)/(CZERO-CIN)
X(I)=CON(I )/CZERO
99 CONTINUE
C
93
STOP
END
c
REAL FUNCTION EXF(QZ)
c
REAL LOEF,Q,Z,T,STUP1ISTUP2,XX,YY,QQ,BOUND2/120./
c
EXFO0.O
STUP1=ABS (Q)
IF((STUP1.GT.BOUND2) .AND. (Z.LE..)) GO TO 40
IF(Z.NE.O.O) GO TO 31
EXF=EXP(Q)
GO TO 40
31 QQ=Q-Z*Z
STUP2=ABS (QQ)
IF((STUP2.GT.BOUND2) .AND. (Z.GT.0.)) GO TO 40
IF(QQ.LT.-BOUND2) GO TO 34
XX=ABS (Z)
IF(XX.GT.3.) GO TO 32
T=101(1 .0+0.3275911*XX)
YY=T*( .2548296-T*( .2844967-T*(1 .421414-
$ T*(1.453152-1.0614O5*T))))
GO TO 33
32 YY=O.56418958/(XX+.5/(XX+.I1(XX+1.5/(XX+2./
$ (XX+2.5/(XX+1.))))))
33 EXF=YY*EXP(QQ)
C
34 IF(Z.LT.0.0) THEN
EXF=2.*EXP(Q)-EXF
ENDIF
c
40 CONTINUE
C
RETURN
END
/GO
94
Listing of HOMSOIL
//FIG22 JOB (AYLS9FL,124),'FEIKE LEIJ',MSGCLASS=P,
// NOTIFY=AYL59FL,MSGLEVEL=(1,1),REGION=1024K,TIME(1,59)
/*ROUTE PRINT RMT4
/*JOBPARM LINES=10
//STEP1 EXEC WATFIV
//WATFIV.SYSIN DD *
/JOB FEIKE,TIME=(05),PAGES=60
C
C VARIABLE DECLARATION
INTEGER I,J,K,SIZE/25/,DUMBOU,INLET/3/
REAL DELTIM/0.005/,DELPOS/0.5/,COLD(63),CNEW(63),CTOT/O.O1/
REAL TIMMAX/5.0/,TIME,POS(63),CIN/1./
REAL DBD/1.4/,CEC/0.0/,VMC/0.4/,VELO/50./,DIFO/50./
REAL A(63),B(63),C(63),F(63),CINIT/0.05/
REAL DBDVMC,G,H,RETCO(63),ADS(63)
REAL ACCNEW,ACCOLD,ACCIN,ORMB,RMB,SUMFLU,SUMACC,FLUX
C
C VARIABLES
C INLET : INLET CONDITION(FIRST=1,THIRD=3)
C SIZE : NUMBER OF SPATIAL NODES
C K PRINT COUNTER
C DELTIM TIME STEP T
C DELPOS SPACE STEP L
C COLD : CONCENTRATION AT OLD TIME LEVEL
C CNEW : CONCENTRATION AT NEW TIME LEVEL
C CTOT : TOTAL SOLUTE CONCENTRATION M/L3
C TIMMAX : MAXIMUM TIME FOR NUMERICAL SIMULATION T
C CIN : ELUENT CONCENTRATION
C DBD : DRY BULK DENSITY M/L3
C CEC : CATION EXCHANGE CAPACITY (MMOL/G) M/M
C VMC : VOLUMETRIC WATER CONTENT L3/L3
C VELO : PORE WATER VELOCITY L/T
C DIFO DISPERSION COEFFICIENT L2/T
C CINIT INITIAL CONCENTRATION
PRINT,'
PRINT,'
PRINT,' * H 0 MSOIL
PRINT,' **I
PRINT,' * SOLUTION OF ADVECTION DISPERSION 
EQUATION *
PRINT,' * IN DIMENSIONLESS FORM, 
INCLUDING NON-LINEAR *1
PRINT,' * EXCHANGE, USING THE CRANK-NICHOLSON 
METHOD '
IF (INLET.EQ.1) THEN
PRINT,' * FIRST TYPE IBC HOMOGENEOUS SYSTEM 
*
ELSE
PRINT,' * THIRD TYPE IBC HOMOGENEOUS 
SYSTEM *
ENDIF
PRINT,' * *
PRINT,' ***************************************************
PRINT,'
C
PRINT,'TIMESTEP =' ,DELTIM
95
PRINT, 'SPACE-INCREMENTh' ,DELPOS
PRINT,'
C
TIMEO0.O
K1l
DUMBOU=SIZE- 1
C
G=DIFO/ (2 .*DELPOS*DELPOS)
H=VELO/ (4. *DELPOS)
C
DBDVMC=DBD*CEC/ (vMc*cToT)
C
C DEFINE INITIAL CONDITIONS AND BOUNDARY CONDITIONS FOR THE PROBLEM
C POSED, AND POSITIONS
C
READ, (COLD(I) ,11,1O)
READ, (COLD(I) ,I=11 ,20)
C READ,(COLD(I),I=21,SIZE)
C
DO 11 1=21,SIZE
COLD(I)=CINIT
11 CONTINUE
C
ACCOLDO .0
DO 26 1=2,SIZE
CALL ADSORB(I ,ADS, COLD)
ACCOLD=ACCOLD+O . *DELPOS*(CTOT*VMC*(COLD(I-J.)+COLD( I) )?
26 CONTINUE
C
SUMFLUO .0
ACCIN=ACCOLD
C PRINT, 'ACCIN=' ,ACCIN
ORMBO0.O
RMBO0.O
CIN1l.
C
C IF(INLET.EQ.1) THEN
C COLD(1)=CIN
C ELSE
C COLD(1)=12.*VELO*DELPOS*CIN/(12.*VELO*DELPOS+25*DIFO)
C ENDIF
C
DO 10 1=1,SIZE
96
115 FORMAT(' 1,3A6,16F6.2)
PRINT,'
9
PRINT 114,TIME,ORMB,RMB,(COLD(I),1I1,SIZE,4)
PRINT 115,' 1, 1 ',' ,(ADS(I),I=1,SIZE,4)
c
C INITIALIZE DUMMY VARIABLES AND TRIDIAGONAL MATRIX
C
A(SIZE)=-1.0
C(SIZE)=O.O
C
WHILE (TIME.LE.TIMMAX) DO
c
CALL EXCHAN(SIZE, COLD, DBDVMC,RETCO)
C
IF (INLET.EQ.1)-THEN
A(1)0O.O
B(1)=1.0
C(1)0O.O
F(1)=CIN
ELSE
A(1)0O.O
B(1)=RETCO(1)?2.*DELTIM*(G+2.*H*(1+(H/G)))
C(1 )=-G*2 .*DELTIM
F(1)=(-B(1)+2.*RETCO(1) )*COLD(1)+
$ 2.*DELTIM*(G*COLD(2)+4.*H*(1+(H/G))*CIN)
ENDIF
C
DO 25 I=2,DUMBOU
A(I )=-G-H
C(I)=-G+H
B(I)2 .O*G+(RETCO(I)/DELTIM)
$ C(I)*COLD(I+1)
25 CONTINUE
C
B(SIZE)1I.+O.5*RETCO(SIZE)/(G*DELTIM)
F(SIZE)=COLD(DUMBOU-)+(B(SIZE)-2. )*COLD'(SIZE)
C
C SOLVE THE DIFFERENTIAL EQUATION WITH THE CRANK-NICOLSON METHOD
C BY USING THE THOMAS ALGORITHM
C
CALL TRIDIA(A,B,C,CNEW,F,SIZE, INLET)
C
C REQUIRED PRINTOUT.
97
$ DBD*CEC*(ADS(I-1)+ADS(I)))
29 CONTINUE
C
FLUX=0. 5*VELO*DELTIMACTOTkVMC*(2*CIN-CNEW(SIZE)-COLD(SIZE))
SUMACC=ACCNEW-ACCIN
SUMFLU=SUMFLU+FLUX
RMB=100.*ABS( (ACCNIEW-ACCOLD-FLUX)/FLUX)
ORMB=100.*ABS( (SUMACC-SUMFLU)/SUMACC)
ACCOLD=ACCNEW
IF (K.EQ.PRINT) THEN
PRINT 114,TIME,ORMB,RMB,(CNEW(I),I=1,SIZE,4)
PRINT 115,1 1.1 1.1',(ADS(I),I=1,SIZE,4)
K=0
ENDIF
K=K+1
C
DO 60 I=1,SIZE
COLD(I )=CNEW(I)
60 CONTINUE
C
ENDWHILE
C
STOP
END
C
C
SUBROUTINE TRIDIA(A,B,C,CNEW,F,SIZE,I'NLET)
C
C THIS SUBROUTINE CONTAINS THE THOMAS-ALGORITHM, THE SOLUTION IS
C CONTAINED IN THE CNEW-ARRAY
C
REAL A(63),B(63),C(63),CNEW(63),F(63)
INTEGER SIZE, INLET
C
REAL GAMMA(63),BETA(63)
INTEGER I ,J,DUMBOU
C
DUMBOU=SI ZE-1
C
BETA(1)=B(1)
GAMMA (1) =F (1) /B (1)
DO 200 I=2,SI,ZE
98
CNEW(SIZE)=GAMMA(SIZE)
WHILE (I.GE.1) DO
CNEW(I)=GAMMA(I)-(C(I)*CNEW(I+1)/BETA(I))
I=I-1
ENDWHILE
C
RETURN
END
C
C *********************************************************
C
SUBROUTINE EXCHAN(SIZE,COLD,DBDVMC,RETCO)
C
REAL COLD(63),DBDVMC,RETCO(63)
INTEGER I,SIZE
REAL SLOPE(63)
C
DO 300 I=1,SIZE
SLOPE(I)=1.
RETCO(I)=1.+DBDVMC*SLOPE(I)
300 CONTINUE
C
RETURN
END
C
C *********************************************************************
c
SUBROUTINE ADSORB(I,ADS,CNEW)
C
REAL ADS(63),CNEW(63)
INTEGER I
C
ADS(I)=CNEW(I)
ADS(I-1)=CNEW(I-1)
C
RETURN
END
C
/GO
1.000 0.990 0.980 0.960 0.940 0.900 0.810 0.700 0.525 0.350
0.240 0.150 0.110 0.090 0.070 0.060 0.050 0.050 0.050 0.050
0.050 0.050 0.050 0.050 0.050 0.050 0.050 0.050 0.050 0.050
99
Listing of LAY1
//FIG19B JOB (AYL59FL,124),'FEIKE LEIJ',MSGCLASSP,
// NOTIFY=AYL59FL,MSGLEVEL=(1,1),REGION=1024K,TIME=(1,59)
/*ROUTE PRINT RMT4
/*JOBPARM LINES=10
//STEP1 EXEC WATFIV
//WATFIV.SYSIN DD *
/JOB FEIKE,TIME=(05),PAGES=60
C
C VARIABLE DECLARATION
INTEGER I,J,K,NUMLAY/2/,LAY(2)/11,70/,SIZE/61/,DUMBOU,PRINT/10/,PC
REAL DELTIM/0.005/,DELPOS/0.5/,COLD(63),CNEW(63),CTOT/0.O1/
REAL TIMMAX/5.0/,TIME,POS(63),CIN,ADS(63)
REAL DBD(2)/2*1.4/,CEC(2)/0.01,0.01/,PULSE/0.5/
REAL VMC(2)/2*0.4/,VELO(2)/50.,50./,DIFO(2)/20.,20./
REAL A(63),B(63),C(63),F(63),CINIT/0.05/
REAL DBDVMC(2),G(2),H(2),RETCO(63),W(2),FLU(2)
REAL ACCNEW,ACCOLD,ACCIN,ORMB,RMB,SUMFLU,SUMACC,FLUX
C
VARIABLES
NUMLAY
LAY
SIZE
PRINT
DELTIM
DELPOS
TIMMAX
DBD
CEC
PULSE
VMC
VELO
DIFO
CINIT
CTOT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT,
NUMBER OF LAYERS
INTERFACE NODES
NUMBER OF SPATIAL NODES
PRINT COUNTER
TIME STEP
SPACE STEP
MAXIMUM TIME FOR NUMERICAL SIMULATION
DRY BULK DENSITY
CATION EXCHANGE CAPACITY (MMOL/KG)
TIME OF PULSE DURATION
VOLUMETRIC WATER CONTENT
PORE WATER VELOCITY
DISPERSION COEFFICIENT
INITIAL CONCENTRATION
TOTAL CONCENTRATION
T
L
T
M/L3
M/M
T
L3/L3
L/T
L2/T
M/ L3
80,1
85,' *.******k********k*****************A******A.AA-AAA
85,' * LAYONE 
*1
85,' *
85,' * SOLUTION OF ADVECTION 
DISPERSION EQUATION '
85,' * IN DIMENSIONLESS FORM, INCLUDING 
NON-LINEAR '
85,' * EXCHANGE, USING THE 
CRANK-NICHOLSON METHOD *'
85,' * FIRST TYPE IBC 
A'
85,' * LINEAR & FAVORABLE 
EXCHANGE '
85,' ********************************************
PRINT 
85,'--------
DO 15 K=1,NUMLAY
PRINT 95,' LAYER # ',K
PRINT,'
100
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
c
F~I
PRINT 90, 'DISPERSION COEFFICIENT ',DIFO(K), 'CM2/D'
PRINT 90, 'PORE WATER VELOCITY ' ,VELO(K), 'CM/D'
PRINT 90, 'CATION EXCHANGE CAPACITY' ,CEC(K), 'ME/G'
PRINT 90,'DRY BULK DENSITY 'DBD(K),'G/CM3'
PRINT,.'
PE (K)=VELO(K)*DELPOS/DIFO(K)
CR(K)=VELO (K) *DELTIM/DELPOS
PRINT 90,1'PECLET NUMBER' ,PE(K),'
PRINT 90,.'COURANT NUMBER' ,CR(K),'
PRINT 85,' ----
15 CONTINUE
C
PRINT,'
PRINT 90, 'TIMESTEP' ,DELTIM, 'D'
PRINT 90, 'SPACE-INCREMENT' ,DELPOS, 'CM'
PRINT,'1
80 FORMAT(1',A60)
85 FORMAT(A60)
90 FORMAT(A32,3X,F7.3-,A1O)
95 FORMAT(A40,3X,1I2)
C
DUMBOU=SIZE- 1
CIN1.,
C
DO 21 K=1,NUMLAY
G(K)=DIFO(.K)/ (2. .DELPOS*DELPOS)
H(K)=VELO(K)/ (4. *DELPOS)
DBDVMC(K)=DBD(K)*CEC(K)/ (vMC(K)*CToT)
W(K)=O.5
21 CONTINUE
C
C DEFINE INITIAL CONDITIONS AND BOUNDARY CONDITIONS FOR THE PROBLEM
C POSED, AND POSITIONS
C
DO 11 1=1,SIZE
COLD(I )=CINIT
11 CONTINUE
C
DO 10 1=1,SIZE
POS(I)=(I.-1)*DELPOS
10 CONTINUE
C
C DETERMINE INITIAL AMOUNT OF SOLUTE
C
101
IF (I.EQ.LAY(K)) K=K+1
27 CONTINUE
C
SUMFLUO .0
ACCIN=ACCOLD
C
COLD(1 )=CIN
C
PRINT,' POSITIONS'
PRINT 113,' TIME.',' ORMB ',' RMB ',(POS(I),I=1,SIZE,4)
113 FORMAT(' ',3A6,16F6.1)
114 FORMAT(' ',19F6.2)
.115 FORMAT(' '3A6,16F6.2)
PRINT,'
C
TIME=O.
ORMB=O.
RMB=O.
C
PRINT 114,TIME,ORM4B,RMB,(COLD(I),I=1,SIZE,4)
PRINT 115,' ', 1'' ,(ADS(I),1I1,SIZE,4)
PC=1
c
C INITIALIZE DUMMY VARIABLES AND TRIDIAGONAL MATRIX
C
A(SIZE)=-1.O
C(SIzE)=O .0
C
WHILE (TIME.LE.TIMMAX) DO
C
CALL EXCHAN(SIZE ,W,LAY,COLD ,DBDVMCRETCO)
C
IF(TIME.LT.PULSE) THEN
CIN1l.
ELSE
CIN=CINIT
ENDIF
C
K1l
C
A(1)0O.O
B(1)=1.0
C(1)0O.O
F(1)=CIN
102
W(K)=1 .1(1.+(FLU(K)/FLU(K+1)))
B(I)=(RETCO(I)/DELTIM)+2.*(W(K)k*G(K)+(1-W(K))*G(K+1))
$ COLD(I>-C(I)*COLD(I+1)'
K=K+1.
ELSE
A(I)=-G(K)-H(K)
c(I)=-G(K)sH(K)
B(I)=2.O*G(K)+(RETCO(I)/DELTIM)
$ C(I)*COLD(I+1)
ENDIF
25 CONTINUE
C
B(SIZE)1 .+O 5*RETco(sIzE)/ (G(K)*DELTIM)
F(SIZE)=COLD(DUMBOU)+(B(SIZE)-2. )*COLD(SIZE)
C
C SOLVE THE DIFFERENTIAL EQUATION WITH 
THE CRANK-NICOLSON METHOD
C BY USING THE THOMAS ALGORITHM
C
CALL TRIDIA(A,B,C,CNEW,F SIZE, INLET)
C
CALL EXCHAN(SIZE-,W,,LAY,,CNEWDBDVMC,RETCO)
C
C REQUIRED PRINTOUTANDMASS BALANCE CALCULATIONS
C
TIME=TIME+DELTIM
C
ACCNEWO 0
K(1
DO 29 1=2,SIZE
CALL ADSORB(I,K,ADS,CNEW)
ACCNEW=ACCNEW+O. 5*DELPOS*(CTOT*VMC(K)*(CNEW(I-1 
)+CNEW(I) )?
$ DBD(K)*CEC(K)*(ADS(I'-1)+ADS(I)))
IF (I.EQ.LAY(K)) K=K+1
29 CONTINUE
C
FLUX=DELTIM*CTOT* (VELO (1) AVMC (1) CIN-O. 
5*VELO (NUMLAY) *
$ VMC(NUMLAY)*(CNEW(SIZE)+COLD(SIZE)))
SUMACC=ACCNEW-ACCIN
SUMFLU=SUMFLU+FLUX
RMB=100.*AB ( (ACCNEW-ACCOLD-FLUX) /FLUX)
103
ENDIF
PC=PC+1
DO 60 1=1,SIZE
COLD(I)=CNEW(I)
60 CONTINUE
ENDWHILE
STOP
END
SUBROUTINE TRIDIA(A,B,C,CNEW,F,SIZE,INLET)
C
C THIS SUBROUTINE CONTAINS THE THOMAS-ALGORITHM,
C CONTAINED IN THE CNEW-ARRAY
THE SOLUTION IS
REAL A(63),B(63),C(63),CNEW(63),F(63)
INTEGER SIZE, INLET
REAL GAMMA(63),BETA(63)
INTEGER I ,J,DUMBOU
DUMBOUSI ZE -1
BETA (1) =B( 1)
GAMMA(1)=F(1)/B(1)
DO 200 1=2,SIZE
200 CONTINUE
BEGIN BACKWARD SUBSTITUTION FROM LAST ROW, INCLUDE BOUNDARY 
CONDITIONS
I=SIZE-1
CNEW(SIZE)=GAMMA(SIZE)
WHILE (I.GE.1) DO
CNEW(I)=GAMMA(I)-(C(I)*CNEW(I+l)/BETA(I))
I1I-1
ENDWHILE
RETURN
END
SUBROUTINE EXCHAN(SIZE,W,LAY,COLD,DBDVMC,RETCO)
REAL COLD(63).,DBDVMC(2),RETCO(63),W(2)
104
C
INTEGER I,K,LAY(2),SIZE
REAL SLOPE(63)
C
K1l
DO 300 11I,SIZE
IF(I.LT.LAY(1))-THEN
SLOPE(I)1l.
RETCO(I)1 .+DBDVMC(K)*SLOPE(I)
ELSE IF (I.EQ.LAY(1)) THEN
RETCO(I)=W(K)*(1 .+DBDVMC(K) )+(1-W(K) )*(1 .+DBDVMC(K+1)*
$ (0.33333333333/(COLD(I)**O.666666666666)))
K=2
ELSE
SLOPE(I)=0.333333333333/(COLD(I)**0.6666666666666)
RETCO( I )1 .+DBDVMC(K)*SLOPE(I)-
ENDIF
300 CONTINUE
C
RETURN
END
C
C**********************************k
C
SUBROUTINE ADSORB(I ,K,ADS,CNEW)
C
REAL ADS(63),CNEW(63)
INTEGER I,K
C
IF(K.GT.1) THEN
ADS(I-1)=CNEW(I-1 )**0 .33333333333
ADS (I )=CNEW(I )**Q 3333333333
ELSE
ADS(I-1)=CNEW(I-1)
ADS(I)=CNEW(I)
ENDIF
.C
RETURN
END
/GO
105
Listing of LAY3
//FIG21B JOB (AYLS9FL,124),'FEIKE LEIJ',MSGCLASS=P,
1/ NOTIFY=AYLS9FLMSGLEVEL=(1,1) ,REGION=1024KTIME=(1,59)
/*ROUTE PRINT RMT4
/*JOBPARM LINES=10
//STEP1 EXEC WATFIV
//WATFIV.SYSIN DD *
/JOB FEIKE,TIME=(05),PAGES=60
C
C VARIABLE DECLARATION
INTEGER I,J,K,NUMLAY/3/,LAY(3)/21,26,80/,SIZE/63/,DUMBOUPRINT/20/
INTEGER UPLAY(3) ,PC
REAL DELTIM/O.05/,DELPOS/O.5/,COLD(63),CNEW(63) ,CTOT/O.O1/
REAL TIMMAX/100./,TIME,POS(63),CIN/1./,ADSOLD(63),ADSNEW(63)
REAL DBD(3)/3*1.3/,CEC(3)./O.O1,O.O,O.O1/
REAL PE(3),CR(3),,VMC(3)/O.2,.4,.2/,DIFO(3)/4..,2.5,4./
REAL A(63),B(63),C(63),F(63),CINIT/O.O5/,VELO(3)/1.,.5,1./
REAL DBDVMC(3),G(3),H(3),DV(3),VD(3),RETCO(63)
REAL ACCNEW,ACCOLD ,ACCIN,ORMB ,RMB ,SUMFLU, SUMACC, FLUX
VARIABLES
NUMLAY
LAY
SIZE
PRINT
DELTIM
DELPOS
CTOT
TIMMAX
CIN
CINIT
DBD
CEC
VMC
VELO
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT,
NUMBER OF LAYERS
INTERFACE NODES
NUMBER OF SPACE NODES
PRINT FREQUENCY
TIME STEP
SPACE STEP
TOTAL CONCENTRATION
MAXIMUM TIME OF SIMULATION
ELUENT CONCENTRATION
INITIAL CONCENTRATION
DRY BULK DENSITY
CATION EXCHANGE CAPACITY (MMOL/KG)
VOLUMETRIC WATER CONTENT
PORE WATER VELOCITY.
80,
85,'
85,'
85,'1
85,'
85,'1
851,'
85,'1
85,'
8-5,'
YK*
,,*
T
L
M/L3
T
M/L3
M/M
L3/L3
L/T
L A Y T HREE
SOLUTION OF ADVECTION DISPERSION EQUATION
IN DIMENSIONLESS FORM, INCLUDING NON-LINEAR
EXCHANGE, USING THE CRANK-NICHOLSON METHOD
THIRD TYPE IBC
SA ND
PRINT 85,' ---
DO 15 K1I,NUMLAY
PRINT 95,' LAYER
PRINT,'1
C
*t
*1
*1
*1
# K
106
PRINT 90, 'DISPERSION COEFFICIENT ',DIFO(K), 'CM2/D'
PRINT 90, 'PORE WATER VELOCITY ',VELO(K), 'CM/D'
PRINT 90, 'CATION EXCHANGE CAPACITY' ,CEC(K), 'ME/G'
PRINT 90,'DRY BULK DENSITY I DBD(K),#G/CM3'
PRINT 90, 'VOLUMETRIC WATER CONTENT' ,VMC(K), 'CM3/,CM3'-
PRINT,'
PE(K)=VELO(K)*DELPOS/DIFO(K)
CR (K)=VELO (K) *DELTIM/DELPOS
PRINT 90,'PECLET NUMB ER' ,PE(K),'
PRINT 90,'COURANT NUMBER' ,CR(K),'
PRINT 85,' - - - - - -
- -- - - - - - - - - -- 
- - - - - - - -
15 CONTINUE
C
PRINT,'
PRINT 90, 'TIMESTEP' ,DELTIM, 'D'
PRINT 90, 'SPACE-INCREMIENT' ,DELPOS, 'CM'
PRINT,'
80 FORMAT(1' ,A60)
85 FORMAT(A60)
90 FORMAT(A32,3X,F7.3,A1O)
95 FORMAT(A40,3X,12)
C
DUMBOU=SI ZE -1
C
DO 21 K=1,NUMLAY
G(K)=DIFO(K)/ (2,*DELPOS*DELPOS)
H(K)=VELO(K)/(4.*DELPOS)
DBDVMC(K)=DBD(K)*CEC.(K)/ (VMC(K)*CTOT)
UPLAY (K )=LAY(K) +1
VD(K)=2. *DELPOS*VELO(K)/DIFO(K)
DV(K)=2.*DIFO(K)/ (VELO(K)*DELPoS)
21 CONTINUE
c
C DEFINE*INITIAL CONDITIONS AND BOUNDARY CONDITIONS FOR THE PROBLEM
C POSED, AND POSITIONS
C
DO 11 1=1,SIZE
COLD(I )=CINIT
11 CONTINUE
c
K1l
DO 10 1=1,SIZE
POS(I)=(I-K)*DELPOS
IF (I.EQ.LAY(K)) K=K+1
107
I=2
WHILE (I.LE.SIZE) DO
CALL ADSORB(I ,K,ADSOLD,COLD)
ACCOLD=ACCOLD+0.5*DELPOS*(CTOT*VMC(K)*(COLD(I-1)+COLD(I) )+
$ DBD(K)*CEC(K)*(ADSOLD(I-1)+ADSOLD(I)))
IF (I.EQ.LAY(K)) THEN
K=K+1
1=1+1
ENDIF
1=1+1
ENDWHILE
C
SUMFLU=0 .0
ACCIN=ACCOLD
C
COLD(1)=12.*VELO(1)*DELPOS*CIN/(12.*VELO(1)*DELPOS+25.*DIFO(1))
C.
PRINT,' POSITIONS'
PRINT 113,' TIME 1,1 ORMB 1,1 RMB l,(POS(I),I=1,21,4),
$ (P05(I) ,I=22,26,2) ,(P05(I) ,I=27,SIZE,4)
113 FORMAT(' ',3A5,19F5.1)
114 FORMAT(' 1,22F5.2)
PRINT!'
c
PC=1
TIMEO0.O
C
C INITIALIZE DUMMY VARIABLES AND TRIDIAGONAL MATRIX
C
A(SIZE)=-1.0
C(SIZE)=O .0
C
CALL EXCHAN(SIZELAYCOLD ,DBDVMC ,RETCO)
C
WHILE (TIME.LE.TIMMAX) DO
C
K1l
C
A(1)0O.O
B(1)=RETCO(1)+2.*DELTIM*(G(1)+2.*H(1)*(1.+(H(1)/G(1))))
C(1)=-G(1)*2.*DELTIM
F(1)=(-B(1)+2.*RETCO(1) )*COLD(1)+
$ 2.*DELTIM*(G(1)*COLD(2)+4.*H(1)*(1 .+(H(1)/G(1) ))*CIN)
C
108
$ COLD(I+1)+( (G(K)-H(K) )*Dv(K+1)*vD(K) )*coLD(I+2)
ELSE IF(I.EQ.UPLAY(K)) THEN
A(I)=-(G(K+1)+H(K+1) )*vD(K+1)
B(I)=(RETCO(I)/DELTIM)+2.*G(K+1)+(G(K+1)+H(K+1))AVD(K+1)
c(I)=-2.*G(K+1)
F(I)=((G(K+1)+H(K+1))*DV(K)*VD(K+1))*COLD(I-2)+
$ ((G(K+1)+H(K+1) )*(1 .-DV(K))*VD(K+1))*COLD(I-1)+
$ ((RETCO(I)/DELTIM)-2.*G(K+1)-(G(K+1)-+H(K+1))*VD(K+1))*
$ COLD(I)+2.*G(K+1)*COLD(I+1)
K=K+1
ELSE
A(I)=-G(K)-H(K)
C(I).=-G(K)+H(K)
B(I )=2 O*G(K)+(RETCO(I )/DELTIM)
$ C(I)*coLD(I+1)
ENDIF
25 CONTINUE
C
B(SIZE)1l.+O.5*RETCO(SIZE)/(G(NUMLAY)*DELTIM)
F(SIZE)=COLD(DUMBOU)+(B(SIZE)-2. )*COLD(SIZE)
C
C SOLVE THE DIFFERENTIAL EQUATION WITH THE CRANK-NICOLSON METHOD
C BY USING THE THOMAS ALGORITHM
C
CALL TRIDIA(A,B,C,CNEW,F,SIZE,INLET)
C
CALL EXCHAN(SIZE ,LAY,CNEW,DBDVMC,RETCO)
C
C REQUIRED PRINTOUTANDMASS BALANCE CALCULATIONS
C
TIME=TIME+DELTIM
C
ACCNEW=0.0
I=2
K1l
WHILE (I.LE.SIZE) DO
CALL ADSORB(I ,K,ADSNEW,CNEW)
ACCNEW=ACCNEW+O. 5*DELPOS*(CTOTAVMC(K)*(CN4EW(I-I)+CNEW(I) )+
$ DBD(K)*CEC(K)*(ADSOLD(I-1)+ADSOLD(I)))
IF (I.EQ.LAY(K).) THEN
1=1+1
K=K+1
ENDIF
109
RIIB100.*ABS (ACCN4EW-ACCOLD-FLUX) /FLUX
ORMB=10O.*ABS(SUMlACC-SUMIFLU)/SUMACC
ACCOLD=ACCNEW
C
IF (PC.EQ.PRINT) THEN
PRINT 114,TIME,ORM4B,RMB,(CNEW(I),I1.21,4),
$ (CNEW(I) ,I=22,26,2) ,(CNEW(I) 
,I=27,SIZE,4)
PC=O
ENDIF
PC=PC+1
C
DO 60 I=1,SIZE
COLD(I )=CNEW(I)
ADSOLD( I)=ADSNEW(I)
60 CONTINUE
C
ENDWHILE
C
STOP
END
C
C
SUBROUTINE TRIDIA(A,B,C,CNEW,F,SIZE,INLET)
C
C THIS SUBROUTINE CONTAINS THE THOMAS-ALGORITHM,. THE SOLUTION IS
C CONTAINED IN THE CNEW-ARRAY
C
REAL A(63),B(63),C(63),CNEW(63),F(63)
INTEGER SIZE, INLET
C
REAL GAMMA(63),BETA(63)
INTEGER IJ,DUMBOU
C
DUMBOU=SIZE- 1
C
BETA (1)B ( 1)
GAMMA(1)=F(1)/B(1)
DO 200 1=2,SIZE
200 CONTINUE
C
C BEGIN BACKWARD SUBSTITUTION FROM LAST ROW, INCLUDE BOUNDARY CONDITIONS
110
C
RETURN
END
C
C *********************************
C
SUBROUTINE EXCHAN(SIZE ,LAY, COLD,DBDVMC,RETCO)
C
REAL COLD(63),DBDVMC(3),RETCO(63)
INTEGER I,K,LAY(3),SIZE
REAL SLOPE(63)
C
K1l
DO 300 1=1,SIZE
IF (I.LE.LAY(1)) THEN
SLOPE(I)1l.
RETCO(I)1 .+DBDVMC(l)*SLOPE(I)
ELSE IF (I.GT.LAY(2)) THEN
SLOPE(I)1.,
RETCO(IX1 .+DBDVMC(3)*SLOPE(I)
ELSE
SLOPE(I)1l.
RETCO(I)=1 .+DBDVMC(2)*SLOPE(I)
ENDIF
300 CONTINUE
C
RETURN
END
C
C **********************************
C
SUBROUTINE ADSORB (I ,K,ADS ,CNEW)
C
REAL ADS(63),CNEW(63)
INTEGER I,K
C
ADS(I)=CNEW(I)
ADS(I-1)=CNEW(I-1)
C
RETURN
END
/GO
i11
Listing of VARTIM
//F1G22 JOB (AYL59FL,124),'FEIKE LEIJ',MSGCLASSP,
1/ NOTIFYAYL59FL,MSGLEVEL=(1,1),REGION1O24KTIME=(159)
/*ROUTE PRINT RMT4
/*JOBPARM LINES=10
//STEP1 EXEC WATFIV
//WATFIV.SYSIN DD *
/JOB FEIKE,TIME=(OS),PAGES=6O
C
C VARIABLE DECLARATION
INTEGER I,JPC,SIZE/25/,DUMBOU,INLET/3/
REAL DELTIM,DELPOS/O.5/,COLD(33),CNlEW(33),CTOT/0OO/
REAL TIMMAX/4.O/,TIME,POS(33),CIN/.,/
REAL DBD/1.45/,CEC/O.O/,VMC/O.4/,VELO/50./,DIFO/50./
REAL A(33),-B(33),C(33),F(33),CINIT/O.05/
REAL DBDVMC,G,-H,RETCO(33),ADS(33),PRSTEP/O.2/,PRTIMIE,SFT
REAL ACCNEWACCOLD,ACCIN,OR14B,RIIB,SUIFLU,SU14ACC, 
FLUX
LOGICAL CHANGE
VARIABLES
SIZE
INLET
DELPOS
DELTIM
CTOT
CINIT
PRSTEP
CIN
TIMMAX
DBD
CEC
VELO
DI FO
VMC
PRINT,'
PRINT,'
PRINT,'
PRINT,'
PRINT,'
PRINT,'
PRINT,'
NUMBER OF SPATIAL NODES
CONDITION AT INLET BOUNDARY (FIRST=1, 
THIRD=3)
SPACE STEP
(VARIABLE) TIME STEP
TOTAL CONCENTRATION
INITIAL CONCENTRATION
TIME STEP FOR PRINTING
ELUENT CONCENTRATION
MAXIMUM TIME FOR NUMERICAL 
SIMULATION
DRY BULK DENSITY
CATION EXCHANGE CAPACITY 
(MOL/KG)
PORE WATER VELOCITY
DISPERSION COEFFICIENT
VOLUMETRIC WATER CONTENT
L
T
~M/L3
N /L
T
IL 3
M/ M
L/T
L2/T
L3/L3
*V VA R T I ME
4r 
*
SOLUTION OF ADVECTION DISPERSION 
EQUATION
IN DIMENSIONLESS FORM, INCLUDING 
NON-LINEAR
EXCHANGE, USING THE CRANK-NICHOLSON 
METHOD
IF (INLET.EQ.1) THEN
PRINT,' * FIRST 
TYPE IBC
ELSE
PRINT,' * THIRD TYPE 
IBC
ENDIF
PRINT,' *
PRINT,' ************
PRINT,'
PRINT, 'SPACE-INCREMENT' ,DELPOS
HOMOGENEOUS SYSTEM
HOMOGENEOUS SYSTEM
*1
*1
*1
*1
112
C
C
C
c
C
C
C
C
C
C
C
C
C
C
C
C
c
PRINT,'
C
DELTIMO0.005
TIMEO0.O
PC=1
DUMIBOU=SIZE-1
CHANGE . TRUE.
C
GDIFO/ (2 .*DELPOS*DELPOS)
H=VELO/ (4. *DELPOS)
C
DBDVMC=DBD*CEC/ (vMlc*cToT)
C
C DEFINE INITIAL CONDITIONS AND BOUNDARY CONDITIONS FOR THE PROBLEM
C POSED, AND POSITIONS
C
DO 11 11I,SIZE
COLD(I )=CINIT
11 CONTINUE
C
ACCOLDO .0
DO 26 1=2,SIZE
CALL ADSORB (IADSCOL)
ACCOLD=ACCOLD+O.5ADELPOS*(CTOTkVNC*(COLD(I-1 )+COLD(I) )+
$ DBD*CEC*(ADS(I-l1)+ADS(I)))
26 CONTINUE
C
SUMFLUO .0
ACCIN=ACCOLD
C
ORMBO0.O
RMBO0.O
CIN~1.
C
IF(INLET.EQ.1) THEN
COLD(1 )=CIN
ELSE
COLD (1 )=12 .*VELO*DELPOS*CIN/ (12 .*.VELOADELPOS+2 
5 ADI FO)
ENDIF
C
DO 10 11I,SIZE
POS(I)=(I-1)*DELPOS
10 CONTINUE
C
113
PRINT 115,1'f' I ',(ADS(I),I=1,SIZE,2)
C
C INITIALIZE DUMM4Y VARIABLES AND TRIDIAGONAL MATRIX
c
A(SIZE)=-1.0
C(SIzE)=O .0
C
WHILE (TIME.LE.TIMMAX) DO
C
CALL EXCHAN (SIZE,COLDDBDVMC, RETCO)
C
126 CONTINUE
C
IF (INLET.EQ.1) THEN
A(1)0O.0
B(1)=1.0
C(1)0O.0
F(1)=CIN
ELSE
A(1)0O.0
B(1>=RETCO(1)+2.kDELTIMA(G+2.*Hk(1+(H/G)))
C( 1)=-G*2.*DELTIII
F()=-B1) +2.*RETCO (1) )*COLD (1)+
$ 2.*DELTIM*(G*COLD(2)+4.*H*(1+(H/G))*CIN)
ENDIF
C
DO 25 1=2,DUMBOU
A(I )=-G-H
C(I)=-G+H
B(I )=2 .O*G+(RETCO(I)/DELTIM)
$ C(I)*COLD(I+1)
25 CONTINUE
C
B(SIZE)1 .+O .5*RETCO(SIZE)/ (G*DELTIM)
F(SIZE)=COLD(DUMBOU)+(B(SIZE)-2. )*COLD(SIZE)
C
C SOLVE THE DIFFERENTIAL EQUATION WITH THE CRANK-NICOLSON METHOD
C BY USING THE THOMAS ALGORITHM
C
CALL TRIDIA(A,B ,C,CNEW, F,SIZE,INLET)
C
C MASS BALANCE CALCULATIONS
c
114
FLUXO0.5*VELO*DELTIMk*CTOTAVMIC*(2*CINi-CNEW(SIZE)-COLD(SIZE))
SUMACC=ACCNEW-ACCIN
S FT=S UMFLU+ FLUX
RMB=100.*ABS( (ACCNEW-ACCOLD-FLUX)/FLUX)
ORMB=100. *ABS (SUMIACC-SUIMFLU) /SFT
C
IF(TIME.GT.O.1) THEN
C
IF(DELTIM.LT.O.OOOOOOO1) THEN4
PRINT, MASS BALANCE ERROR AT T=',TIME
STOP
ELSE IF(ORMB.LE.O.5 .AND. CHANGE) THEN
DELTIM=2 .*DELTIM
CHANGE=. FALSE.
GO TO 126
ELSE IF (ORMB.GT.3,.AND. CHANGE) THEN
DELTIM0. SADELTIM
CHANGE= .FALSE.
GO TO 126
ENDIF
C
ENDIF
C
CHANGE=. TRUE.
SUMFLU=SUMFLU+ FLUX
ACCOLD=ACCNEW
C
C REQUIRED PRINT OUT
C
TIME=TIME+DELTIM-
PRTIM=PC*PRSTEP
IF (TIME.GE.PRTIM) THEN
PRINT 114,TIME,ORMB,RMB,(CNEW(I),I=1,SIZE,2)
PRINT 115,.1', ',' ,(ADS(I),I=1,SIZE,2)
PC=PC+1
ENDIF
C
DO 60 1=1,SIZE
COLD(I)=CNEW(I)
60 CONTINUE
C.
ENDWHILE
C
STOP
115
C CONTAINED IN THE CINEW-ARRAY
c
REAL A(33),B(33),C(33"),CN\EW(33),F(33)
INTEGER SIZE, INLET
C
REAL GAMMA(.33),BETA(33)
INTEGER I ,J, DUNBOU
c
DUMBOU=SIZE -1
c
BETA(1)=B(1)
GAMMA(1)=F(l)/B(l)
DO 200 1=2,SIZE
200 CONTINUE
C
C BEGIN BACKWARD SUBSTITUTION FROM LAST ROW, INCLUDE BOUNDARY CONDITIONS
C
I=SIZE-1
CNEW(SIZE )=GAMMA(SIZE)
WHILE (I.GE.1) DO
CNEW(I)=GAMMA(I)-(C(I)*CNEW(I+1)/BETA(I))
ENDWHILE
C
RETURN
END
C
C
SUBROUTINE EXCHAN(SIZE,COLD,DBDVMC,RETCO)
C
REAL COLD(33)- ,DBDVMC,RETCO(33)
INTEGER I,SIZE
REAL SLOPE(33)
C
DO 300 I=1,SIZE
SLOPE(I)1l.
RETCO (I )=1 ?DBDVMC*S LOPE (I)
300 CONTINUE
C
RETURN
END
RETURN
END
c
/GO
116