for an Exchangeabl
Solute in a
Layered Soil
II0I
. 0 *0 
0 I. ~ Ill 0 SI 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 Twolayer 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
Nonreactive 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 Firsttype Condition..........
68
APPENDIX B. Analytical Solution
of the ADE for the Second Layer
of a Porous Medium Using
a Thirdtype Condition...........
73
APPENDIX C. Analytical Solution
of the ADE for a Twolayer
Medium Using a Thirdtype Condition
at the Inlet and a
Combined First and Thirdtype
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 Nonlayered 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 basedon
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 firsttype condition
............
29
FIG.5B. Numerical solution
of the ADE for a homogeneous
soil
using a thirdtype condition
... ............
. 29
FIG.6A. Validation of
the analytical solution
of the ADE for a
twolayer medium using
a firsttype condition
. . . . 31
FIG.6B. Validation of
the analytical solution
of the ADE for a
twolayer medium using
a thirdtype condition
. . . . 31
FIG.7A. Validation of
the numerical solution
of the ADE for
a twolayer medium
using a firsttype condition
. . . 32
FIG.7B. Validation of
the numerical solution
of the ADE for
a twolayer medium using
a thirdtype condition
. . . 32
FIG.7C. Validation of the numerical solution
of the ADE for
a twolayer medium using
a thirdtype condition
at the
inlet and a secondtype
condition at the interface
. .33
V
FIG.8A. Numerical solution of the ADE
for a twolayer medium
(D
1
D
2
, v >v ) using a firsttype
condition
FIG.9B. Numerical solution of the ADE for a twolayer
medium
(D >D
2
, v >v ) using a thirdtype condition .....
FIG.9C. Numerical solution of the ADE for a twolayer
medium
(D >D , v >v
2
) using a firsttype condition at the
inlet
2
and a secondtype condition at the interface
FIG.10A. Analytical solution of the ADE for a twolayer
medium
(D
1
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 ) using a
FIG.13B. Analytical solution of
W =D , v >v ) using a
D12 12
the ADE for a twolayer medium
firsttype condition
the ADE for a twolayer medium
thirdtype condition ....
the ADE for a twolayer medium
firsttype condition
the ADE for a twolayer medium
thirdtype condition .....
the ADE for a twolayer medium
firsttype condition ... .
the ADE for a twolayer medium
thirdtype 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 twolayer
me'dium2using a thirdtype
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 twolayer medium using a thirdtype 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 twolayer medium using a thirdtype 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 twolayer medium using a thirdtype condition at the
inlet and a secondtype 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 twolayer medium with: (a) (dY/dX)=1 and (dY/dX) =3X
2
12/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
thirdtype inlet condition: nonlayered profile (a),
and profile with a hypothetical clay layer using a
thirdtype interface condition (b), and secondtype
interface condition (c) ...................... ... 60
FIG.23. Concentration profiles in a hypothetical clay with
a thirdtype inlet condition: nonlayered profile (a),
and profile with a hypothetical sand layer using a
thirdtype interface condition (b), and secondtype
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 fluxtype 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 twolayer 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 1D
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 longterm 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) onelayer and a twolayer medium subject
to a flux and a concentrationtype 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
fluxtype
condition at the inlet of the soil and a combined flux
and
concentrationtype 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 1D steadystate flow,
transport of a nonreactive solute in a homogeneous soil column with
length L can be described as follows:
ac a
2
c aC
aD2 v
OO (1)
at 2 8x
8x
where C is the solute concentration
expressed in mass per volume of
solution [ML3], t is time [T],
D is the dispersion coefficient
21
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 T1, 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)OO (4a)
(D a+ vC)I V
0
= VC (thirdtype) t>O (4b)
where xO, sometimes denoted as xO+, implies that x=O is approached
from inside the column. The firsttype condition (Eq.(4a)) suggests
that the concentration at the interface of influent solution (eluent)
is continuous, whereas the thirdtype condition (Eq.(4b)) 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 concentrationtype condition, the
third or fluxtype 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 thirdtype 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 (5a)
00o 8x
+
ev C) = Ce D IC
C (x, t) = C ?xsL t>O (5c)
0 0
6
where e is the volumetric water content [L3L3],
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.
(5c) 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
fluxaveraged concentration. Eq. (4a) and (4b) then follow from
Eq.(5a) and (5b). We prefer the use of Eq.(4b) because of mass
preservation, accepting the discontinuity in concentration due
to the
use of two different types of concentration, i.e., a fluxtype for the
eluent and a residenttype for the soil. Solutions of the ADE subject
to Eq. (4a) can also be useful (34). Therefore both Eq. (4a) and (4b)
will be used.
For the outlet, the boundary condition can be formulated as:
aC = 0
t>O (6a)
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.(4b), a similar discontinuity
might occur at the outlet. This invalidates Eq.(6a). The problem of
formulating a boundary condition at x=L can be avoided by assuming that
the medium is in effect semiinfinite, i.e.:
aC
a 0
t>0 (6b)
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 semiinfinite media are usually rather small (33).
The analytical solution of Eq.(1) subject to Eq.(3), (4a) and
(6b), i.e., for a semiinfinite medium with a firsttype
condition at the inlet, is given by (18):
C(x,t) = C + (C C (+c IE + cca)(7)
whereas for the thirdtype condition at the inlet (i.e., Eq.(1) subject
to Eq.(3), (4b) and (6b)), the solution is given by:
C(x,t) =Ci +(CoCi) + 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 Twolayer 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 AlNiami
and Rushton (3). Shamir and Harleman (28) used a firsttype inlet
condition and an outlet condition at infinity for each layer, whereas
AlNiami and Rushton (3) used a firsttype 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
semiinfinite, using a first or a thirdtype 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_ 8Cc  C v 
L O
(9)
at 2 2 2ax 1 2
ax
subject to:
C(xO) = C L,O (firsttype) (10b)
1 1 ax 1 1IxT"L
2 2 x +2v2CIx",L
t (hrtp)(Oc
1
c = 0 t>O (10d)
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 firsttype condition, the problem is solved in conjunction with
Eq.(7), whereas Eq.(8) was used to solve the problem for a thirdtype
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.(10b) implies that the volume averaged or
resident concentration is continuous at the interface, while Eq.(10c)
implies that the fluxaveraged concentration is continuous.
The analytical solution of Eq. (9), subject to (10a), (10d), and
(10b) or (10c), was obtained with Laplace transforms. Appendices A
and B contain the details of the respective solution procedures. For a
first and thirdtype condition, the results are, respectively:
10
(C C ) (xL ) tLvT vL 1L +vT
C(x,t) =
o
1 + eapD1
l/6tD [4DT
1 1
2 0 1 1
1
(xLv
(tT))
2
x3/2 4D() dr + C
(11)
(tT)
3/2
4D2
(
tT)
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
(xL
v
(tT))2
/(tT)
4D2
( t  T )
v v (xL) xv(t)
eeA
L+v(tT)
dT + C
(12)
S4D
2
4D (tT)
2 2
The integrals appearing in
these solutions were evaluated
with
help of the GaussChebyshev formula,
discussed in Appendix D.
If we deal with residenttype
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,
semiinfinite cannot be used.
The exit and inlet condition
of
respectively the first and second
layer consist of the two interface
conditions (Eq. (10b) and (10c)).
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
( 15a)
y=S/S
T
..... (15b)
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 forthe following initial
and
boundary conditions:
X(x,O0) =X (x
(D av X) vX
7x x
49
0
ax0
8x IxtL=0
The numerical
CrankNicolson method.
temporal spacing At,
j=0, 1, .. the central
00
t>0
t>0
(19a)
(19b)
v
(19C)
( 19d)
solution is achieved with the implicit
Using a grid system with spatial spacing Ax and
where x.=(i1)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
ii l
(X1+1X 2X
ii
(X , j+x;+X X )
1=2, 3p,4p.,*9.n1
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
24i n
(1a
The inlet condition (x=0) was approximated
by:
xi = X
1 0
D j+ lj+l+j
.J+
 x (X XO +XX X~ x
4Ax 2 0 2 0 2 1 1
vc
0
2
3
n1
L ___
vce
(firsttype)
SvX 0(thirdtype)
J1 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
(21b)
(2 1c)
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 (21d)
n+1 n1 n+1 n1
where the column exit is located at (n1)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 firsttype
the value is
given by (21a), whereas for the thirdtype 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'
(GH) Xi + i +2G X. + (G+H) Xj+
1
11 At i+1
(22)
R .
(G+H) X
j
+ [ 2G] xi
+ (GH)X
11 At 1 l+1
This equation needs to be solved for i=2,...,n1 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
ani bn c" Xnf
. n1 n1 n1 n1
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 kth 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
firsttype thirdtype
i 2
ii
F'
 i H
i+1L
i +2
j+1 
j j+i j1
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
LKi1
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
twolayer 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
(firsttype) (25a)
((D) +(xv) X) =((eD) +(Gv) X) (thirdtype) (25b)
ik 8X (V)k x"L k+15X k+1 x4L
k k
In addition to a first and thirdtype condition, a secondtype
condition was introduced by combining Eq.(25a) and (25b):
(GD) x = (D) X (secondtype) (25c)
ODjk jx"Lk
k+1 ax
xLk
k 8xxk
The secondtype 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.(19a), (19d), and (1)
Eq.(19b) and (25a), (2) Eq.(19c) and (25b), or (3) Eq.(19c) and
(25c) was accomplished with the CrankNicolson method. The spacetime
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 firsttype (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(XJx )
(26a)
k k 1 Ax 1 i1
(OD)k1 0
S= (ev) x _ D) (X X
j
) (26b)
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 (27a)
k 1+(J /J )
k k+1
wk =1w (27b)
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 thirdtype (interface) condition, two different
approximations for Eq. (25b) 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. (25b) 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 (28a)

(
)k " X_1,
+ X +1 + X
(28b)
S) xj+~ j+ + x  X + Xj+ + X
AX V k+1 L i+1 11 1
4
, 11J
It should be noted that X and X are concentrations at
fictituous nodes located at Lk+Ax and LkAx, 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+K1 equations need to be solved, where K is the
total number of layers.
Finally, the secondtype interface condition was approximated by:
(eD)  X
1
+ X X9x = (eD) [+ X +xXi  X (29)
k1 i i1 k+1 I+I 1 +i+I 1
where the same grid system was used as for the firsttype 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
CourantFriedrichsLewy 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 CrankNicolson 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 1D transport of a nonreactive
solute in a homogeneous soil with a dimensionless initial concentration
X =0.05. The ADE is solved numerically for a firsttype 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 semiinfinite versus a finite length (Eq. (6b) and Eq. (19d),
respectively). It should be noted that the CrankNicolson 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
secondorder 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 [xohx0+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 i1 i1(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 (35a)
j+1 T in 2 n n
X 0 < t < t
X. = o o (35b)
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+1ACC 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 (37a)
j+1 J
J=1
and:
ORMBj+. = SUMFLUj+ SUMACC j+I/SUMACCj+ *100% (37b)
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 1D 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 thirdtype condition, respectively, with a fixed time
step. The program LAY1 was modified to solve the ADE for a secondtype
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 thirdtype
condition at the inlet of each layer. In the remaining sections we will
use only a thirdtype condition at the inlet of the soil and a second
or thirdtype 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 thirdtype 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 thirdtype
condition, X(O,t) gradually approaches X (figure 58), whereas for the
0
firsttype 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 thirdtype condition preserved mass. The mass
balance criterion to determine the value of a variable At can therefore
not be used for a firsttype condition. Using a variable time step for
a thirdtype 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 firsttype 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
thirdtype 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 twolayer 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 (firsttype condition) and 6B (thirdtype 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
"twolayer" media.
31
*solution according to Eq.(1
solution according to Eq. (7)
2 4
8 10 12
14
FIG.6A. Validation of
twolayer medium using
x cm]
the analytical solution of
a firsttype 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
twolayer medium using a thirdtype 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 firsttype 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 thirdtype
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 thirdtype condition at the inlet and a
secondtype condition at the interface.
Transport during steady flow in a twolayer 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 secondtype 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
thirdtype condition are very similar except at the interface. The use
of a thirdtype condition results in a discontinuity of the (resident)
concentration, whereas the profile is continuous for a secondtype
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 ) using a firsttype condition.
1 2
third type
a twolayer medium (D >D,
. 1 2
x [cm].
FIG.9B. Numerical solution pf the ADE for
a twolayer medium (D >D
2
,
v >v ) using
a thirdtype
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 twolayer medium (D >D
v >v ) using a firsttype condition at the inlet ana
2
a
S2
secondtype 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
Dvalues. In figure 10, the Dvalue in the first layer is smaller than
in the second layer, whereas in figure 11 the Dvalues are
interchanged. The error in the mass balance can be evaluated by
comparing the firsttype solution (figure iQA and 11A) with the
thirdtype solution (figure lOB and 1iB). For a layered medium, the
repeated application of the firsttype condition leads to an increase
in the mass balance error of the solution compared to a homogeneous
medium (cf. figure 5). For a thirdtype condition, a discontinuity in
38
x [cm]
FIG.10A. Analytical solution of the ADE for a
v =v ) using a firsttype condition.
1 2
10 12 14
twolayer 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 fluxaveraged 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 vvalues on
the concentration profile for steadystate flow. In figure 12, the
porewater 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 thirdtype 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 fluxaveraged
concentration is continuous at the interface.
From these examples it can be concluded that the use of a
firsttype condition leads to large mass balance errors in layered
media. For a thirdtype 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 thirdtype 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
twolayer
v v ) using a firsttype condition.
1 2
twolayer medium (D =D
12
0 2 4 6 8 10 12 14
x [cm]
FIG.13B. Analytical solution of the ADE for a twolayer medium
(D=D
12
v >v ) using a thirdtype 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
CIprofiles
400
x [cm]
FIG. 14A. The effect of front
solution of the ADE using R=1.
0.
retardation on transport. Analytical
Kprofiles
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
xv 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
xv*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
xv* 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
xv*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
. xv*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
onedimensional 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 porewater
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 thirdtype condition at
the interface and a thirdtype 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 twolayer medium
usng a thirdtype 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
twolayer medium, the first layer (00 (Al)
2 8t 2 x2 2 X 1 2
8x
C(x,0) = g L 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) = Fec(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= [3iep(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=xL . 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
(Cogl) 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(tT)h(T) dr = g(T)h(tT) 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 (Cg ) [ec+ emp(2a1k
)
e/,c [2 x
2 (k2a (tT))
2
The inverse of C follows from C
1
, using the properties of the Laplace
transform for a step function:
71
0 Ot
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) =
k2a 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 steptype 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(tT)3 ej 22tT) dT + C. (A19)
LL 4(tr) 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
(t7) 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 Kth
layer. However, it can be verified that G(x,s) for an arbitrary
72
position in this layer is given by:
C C K1iC.
C(x,s) = 0 1 explltexp(A )+ 1
(A20)
S k=1
where t is the thickness of the kth layer and
=xL . Equation
k K
K1
(A20) can be inverted using the convolution theorem
while making use of
the results for layer Ki. 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 Thirdtype 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 thirdtype
condition is preferred over a firsttype condition for volumeaveraged
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 0 (Bl)
2 8t 2 2 2ax 1 2
8x
C(x,0) = g L (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)
vDA   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 gg 2 C]
S= v
2
+  1eC(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)
P12 (k 22a 2Ctz))
2
k
2
+2a
2
(t Tr)
X .6(r(t') e" [ 2 2 ia 2 e4p(2a k ) e/,&c[J1l
L 4(t14 4(tT)
Using the properties of the Laplace transform:
Ot
c C x, t) =2a
2(g1g2 ),ehp(a
2k2

2
ro/16a t
(g g ) (
21 2 II2
[earptk
2
a +s]=
s (a
2
+ a +s]
,e(k22a
2t)
2
,+
k2 2a2t
24tJ +eck 2]+
(l2a k+4ak +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 steptype 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 (k422 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 Kth 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 KD KX K) (kl e.
Again, it appears that the order of layers
between k=2 and K1 is
irrelevant for the concentration
in the Kth layer (LK 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 LO (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 LQ (C2b)
8C rv C OO (C2d)
i Xto 2 X,
9C 8C
1 2
evC D =evCeD t>0 (C2e)
1111 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 1exp(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(xL) + (q +q
2
)rpX.  X(xL)
(C8)
2q [C 9 eXP(t s) ]eXP
 g
2
2q1
C (xs) =+x
2 s 1q )(q
q2)exp(XL) + ( 1+q
)(q +q )exp(AXL)
v (xL) (xL) (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 GaussChebyshev 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 = 12Tt)] For a function
f(T) the integration changes as follows:
t
+1
Jf(z) dr =J2 f( 2 ~t)diq (D1)
0
The GaussChebyshev 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(
(
2kl ) 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 Nonlayered and a Layered Soil
Nonlayered soil
inlet condition (x=O) j>O
c =0 b = 1, f = X firsttype
1 1 1 o
S= 2AtG
b = R + 2At(G + 2H(1 + H thirdtype
1 1 G
f,
=
(2Rb ) X + 2AtG X + 8AtH(1 + ) X
1 1 1 2 G o
a. =G  H
1 .
R
b. =2G +
1 t i=2,3,.. ,n1
c. = G + H jaO
1
f. = a.X + (b.4G) X9 c X
1 1 i1 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 n1 n n
initial condition (t=O)
X i=1
X = firsttype
SX. i>1
1
24H
= {4H2
=
thirdtype
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 firsttype
1 1 o
al = , = 2AtG
H
b = R + 2At(G + 2H (1 + I)) thirdtype
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 11 i k 1 1 1+1
interface condition (x=L) j>O, k=1,2,..,K1
k
a. =w G w G w Hw 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 firsttype
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
(HG )
i At k k k k thirdtype
c. = (k(GkHk)
= k 11 At k k k k i+1 k+k 1
+ [k l (GkHk )]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 secondtype
c. = D
1 k+1 k+1
f. = a.X.  b.X. c.X.
1 1 i1 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 thirdtype
1 k+1
f. = [q( (G +H)XL + (1q) (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
+ (bn2) X
j
n ni n n
initial condition (t=0) j=0 k=1,2,....,K
0 X0
X= firsttype
i X.
i>1
1
24H
k
Xc xi=1
X.
0
= {
24
Hk+
25
G it
y
KOk k thirdtype
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 x2/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
IBMmainframe 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 twolayer medium using a firsttype condition
THSEC : Analytical solution of the ADE for the second layer of
a twolayer medium using a thirdtype condition
HOMSOIL : Numerical solution of the ADE for a homogeneous medium
using a first or a thirdtype condition
VARTIM Numerical solution of the ADE for a homogeneous medium
using a first or a thirdtype condition with a
variable time step
LAY1 : Numerical solution of the ADE for a layered medium using
a firsttype condition
LAY3 Numerical solution of the ADE for a layered medium using
a thirdtype 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 ADVECTIONDISPERSION
'
PRINTI* EQUATION FOR 1D TRANSPORT IN THE 2ND LAYER OF
PRINT,'* A SEMIINFINITE POROUS MEDIUM *
PRINT,'* 1ST 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(I1)+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 GAUSSCHEBYSHEV
C
DO 23 J=1,UP
DUM=((2*J1)*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 .OROOT(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 .OROOT(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)=((K1TEL1(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'ROCIN)
X(I)=CON(I)/CZERO
99 CONTINUE
C
PRINT 114,TIME, 'CCI/COCI' ,(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=QZ*Z
STUP2ABS (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 ADVECTIONDISPERSION
PRINT,'* EQUATION FOR 1D TRANSPORT
IN THE 2ND LAYER OF
PRINT,'* A SEMIINFINITE
POROUS MEDIUM
PRINT,'* 3RD 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(I1)+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 GAUSSCHEBYSHEV
C
DO 23 J1,jUP
DUM=((2*J1)*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
.OROOT(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.OROOT(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*(CZEROCIN)*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 GAUSSCHEBYSHEV
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)*LOEFA2*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)/(CZEROCIN)
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=QZ*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*( .2548296T*( .2844967T*(1 .421414
$ T*(1.4531521.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 NONLINEAR *1
PRINT,' * EXCHANGE, USING THE CRANKNICHOLSON
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, 'SPACEINCREMENTh' ,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(IJ.)+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 )=GH
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 CRANKNICOLSON 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(I1)+ADS(I)))
29 CONTINUE
C
FLUX=0. 5*VELO*DELTIMACTOTkVMC*(2*CINCNEW(SIZE)COLD(SIZE))
SUMACC=ACCNEWACCIN
SUMFLU=SUMFLU+FLUX
RMB=100.*ABS( (ACCNIEWACCOLDFLUX)/FLUX)
ORMB=100.*ABS( (SUMACCSUMFLU)/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 THOMASALGORITHM, THE SOLUTION IS
C CONTAINED IN THE CNEWARRAY
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 ZE1
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=I1
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(I1)=CNEW(I1)
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.AAAAA
85,' * LAYONE
*1
85,' *
85,' * SOLUTION OF ADVECTION
DISPERSION EQUATION '
85,' * IN DIMENSIONLESS FORM, INCLUDING
NONLINEAR '
85,' * EXCHANGE, USING THE
CRANKNICHOLSON 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, 'SPACEINCREMENT' ,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)+(1W(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 CRANKNICOLSON 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(I1
)+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) CINO.
5*VELO (NUMLAY) *
$ VMC(NUMLAY)*(CNEW(SIZE)+COLD(SIZE)))
SUMACC=ACCNEWACCIN
SUMFLU=SUMFLU+FLUX
RMB=100.*AB ( (ACCNEWACCOLDFLUX) /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 THOMASALGORITHM,
C CONTAINED IN THE CNEWARRAY
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=SIZE1
CNEW(SIZE)=GAMMA(SIZE)
WHILE (I.GE.1) DO
CNEW(I)=GAMMA(I)(C(I)*CNEW(I+l)/BETA(I))
I1I1
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) )+(1W(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(I1)=CNEW(I1 )**0 .33333333333
ADS (I )=CNEW(I )**Q 3333333333
ELSE
ADS(I1)=CNEW(I1)
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,'
85,'
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 NONLINEAR
EXCHANGE, USING THE CRANKNICHOLSON 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, 'SPACEINCREMIENT' ,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)=(IK)*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(I1)+COLD(I) )+
$ DBD(K)*CEC(K)*(ADSOLD(I1)+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(I2)+
$ ((G(K+1)+H(K+1) )*(1 .DV(K))*VD(K+1))*COLD(I1)+
$ ((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 CRANKNICOLSON 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(II)+CNEW(I) )+
$ DBD(K)*CEC(K)*(ADSOLD(I1)+ADSOLD(I)))
IF (I.EQ.LAY(K).) THEN
1=1+1
K=K+1
ENDIF
109
RIIB100.*ABS (ACCN4EWACCOLDFLUX) /FLUX
ORMB=10O.*ABS(SUMlACCSUMIFLU)/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 THOMASALGORITHM,. THE SOLUTION IS
C CONTAINED IN THE CNEWARRAY
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(I1)=CNEW(I1)
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
NONLINEAR
EXCHANGE, USING THE CRANKNICHOLSON
METHOD
IF (INLET.EQ.1) THEN
PRINT,' * FIRST
TYPE IBC
ELSE
PRINT,' * THIRD TYPE
IBC
ENDIF
PRINT,' *
PRINT,' ************
PRINT,'
PRINT, 'SPACEINCREMENT' ,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=SIZE1
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(I1 )+COLD(I) )+
$ DBD*CEC*(ADS(Il1)+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)=(I1)*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 )=GH
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 CRANKNICOLSON 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*CINiCNEW(SIZE)COLD(SIZE))
SUMACC=ACCNEWACCIN
S FT=S UMFLU+ FLUX
RMB=100.*ABS( (ACCNEWACCOLDFLUX)/FLUX)
ORMB=100. *ABS (SUMIACCSUIMFLU) /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 CINEWARRAY
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=SIZE1
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