Vadose Zone Processes and Chemical Transport
Aging Effects on Cadmium Transport in Undisturbed Contaminated
Sandy Soil Columns
P. Seuntjens,* K. Tirez, J. Šimůnek, M.Th. van Genuchten, C. Cornelis, and P. Geuzens
ABSTRACT

reaction and multicomponent models to predict Cdtransport in soil columns. A competitive equilibrium
ion-exchange transport model predicted Cd breakthrough adequately, while a multireaction model withindependently measured rate constants largely underestimated the extent of Cd retention. Hinz and Selim(1994) predicted Cd and Zn transport in uniformly
packed soil columns under constant and variable ionicstrength conditions. They failed to predict Cd transport
under constant ionic strength, and from this concluded
that not all reactions were accounted for by the adsorption isotherms. An empirical Cd sorption model that
included pH and organic matter was proposed by van
der Zee and van Riemsdijk (1987) for acid sandy soils.
Temminghoff et al. (1995) extended this sorption model
to a three-species Freundlich equation that accounted
for Ca activity, pH, complexation, and ionic strength of
the soil solution. The modified sorption isotherms were
derived for a range of sandy soils differing in pH, organic
matter content, and soil solution composition. Modified
isotherms are especially advantageous in estimating sorption parameters from basic soil information, thus making them useful for large-scale (field and regional) applications.
Relatively little attention has thus far been paid to
the transport behavior of aged contaminants (i.e., forsorbing solutes that reside for relatively long time periods in the soil, from several years up to tens of years).
Such solutes can migrate into less accessible immobile
soil water regions and/or may diffuse into and become
strongly bound within the soil solid phase. Consequently, the time needed for a contaminant to dissolve
or be transferred from stagnant to mobile water may
increase. Discussions are still ongoing whether or not
irreversible Cd fixation in soils occurs. For example,
Christensen (1984b) did not observe irreversible Cdbinding to a loamy sand and a sandy loam soil in 0.01
M CaCl2 at pH ⫽ 6. He showed that 67 wk of aging of
Cd in the soil did not significantly change Cd desorption.
Ainsworth et al. (1994) studied desorption of 16-wk
aged Cd on hydrous ferric oxide and observed that 20%
did not desorb. Wilkens (1997) found that Cd desorption
from sandy soils deviated from adsorption and attributed the difference to sorption hysteresis. Smolders etal. (1999) measured the exchange between radio-labile109
Cd and indigenous Cd in 10 Belgian soils and foundthat 62 to 90% of the indigenous Cd was labile (i.e.,

Limited information is available on the effects of contaminant aging(i.e., the contact time of Cd with the soil) on Cd transport in soils.
We conducted displacement experiments in which indigenous Cd andfreshly applied Cd were leached simultaneously from undisturbed
samples of three Spodosol horizons. Sorption of Cd was describedusing Freundlich isotherms, whereas transport was described as aconvection–dispersion process. Parameter optimization analysis usinga mobile–immobile transport model applied to nonsorbing tracer displacement data showed that 16 to 22% of the water in the columns
was immobile. The low dimensionless mass transfer coefficients in themobile–immobile model were indicative of diffusion-limited transfer
between mobile and immobile water, and hence physical nonequilibrium. A two-site kinetic sorption model could be fitted closely to
breakthrough curves of the non-aged Cd for three soil horizons. No
conclusive evidence was found that contaminant aging in soil affects
cadmium transport. On the one hand, predictions of aged Cd leaching,using parameters estimated from displacement experiments with nonaged Cd, differed from those for the aged Cd in the E horizon. On
the other hand, no meaningful differences in transport behavior between aged and non-aged Cd were found for the humus Bh and Bh/C
horizons. The two-site kinetic rate coefficient ␣c was found to depend
on water flux, further indicating that mass transfer between sorption
sites and the liquid is limited by diffusion rather than by kinetic
sorption.

Heavy metals are being added to the environment
by atmospheric deposition and by fertilizer, sludge,
and sewage sludge application to the soil. Specific soil
conditions may favor heavy metal migration throughthe unsaturated zone to underlying aquifers. In the longterm, mobile heavy metals may adversely affect drinking
water supplies and the quality of surface waters connected to contaminated areas. The physical and chemical mechanisms affecting heavy metal transport andtransformations need to be accurately represented insolute transport models in order to correctly predict soil
and ground water quality, as well as to evaluate thefeasibility of soil cleanup programs.
Because of its mobility in the environment and its
harmful effects on humans, the behavior and transportof Cd in soils has been widely studied. For example,
Selim et al. (1992) compared the performance of multiP. Seuntjens, K. Tirez, C. Cornelis, and P. Geuzens, Vito, Flemish
Institute for Technological Research, Boeretang 200, B-2400 Mol,
Belgium. J. Šimůnek and M.Th. van Genuchten, USDA-ARS, George
E. Brown, Jr., Salinity Lab., 450 W. Big Springs Road, Riverside, CA
92507-4617. Received 14 June 2000. *Corresponding author (piet.
seuntjens@vito.be).

Abbreviations: 2SF, two-species Freundlich model; 3SF, three-species
Freundlich model; BTC, breakthrough curve; CDE, convection–
dispersion model; LEA, local equilibrium; TRM, two-region model.

Published in J. Environ. Qual. 30:1040–1050 (2001).

1040

1041

SEUNTJENS ET AL.: AGING EFFECTS ON CADMIUM TRANSPORT IN SOIL COLUMNS

isotopically exchangeable). A decrease in the soil pH
increased the labile Cd pool, thus indicating that changes
in soil solution composition may cause a release of previously sorbed Cd.
Most Cd transport studies have been carried out using
homogeneous saturated soil columns. However, field
and laboratory experiments with conservative tracers inundisturbed soils have shown that heterogeneous flow
and preferential or nonequilibrium transport conditions
often exist in unsaturated soils due to spatial heterogeneity of hydraulic properties, the presence of macropores, flow instabilities, and variable boundary conditions. Studies indicate that preferential solute transport
is common in many coarse-textured soils (Glass et al.,
1989; Baker and Hillel, 1990; Ritsema et al., 1993; Kung,
1993). Laboratory experiments with undisturbed soilcores or in situ field experiments are useful for studying
preferential flow and transport. Büchter et al. (1996)
conducted a migration experiment with Cd in an undisturbed stony soil column. They concluded that Cd
transport could be well predicted using a convective–
dispersive transport model with either a monocomponent isotherm or an ion-exchange isotherm. They observed sorption hysteresis, which they modeled as anirreversible sorption process. Unsaturated conditions
as such may also contribute to non-ideal tracer transport (e.g., Biggar and Nielsen, 1960).
Non-ideal reactive solute transport may be describedusing a bicontinuum concept by dividing the soil into
two pore water regions (physical nonequilibrium; van
Genuchten and Wierenga, 1976) or two sorption domains (chemical nonequilibrium; Selim et al., 1976).
Nonequilibrium is then either assumed to be caused byrate-limited diffusive mass transfer between mobile andimmobile water regions, or by kinetic sorption reactions.
Nkedi-Kizza et al. (1984) showed that the mathematicalformulation of physical and chemical nonequilibrium
processes is equivalent and that the two types of nonequilibrium cannot be easily distinguished on the basisof column outflow experiments. Although conceptuallymore elaborate approaches have been proposed to distinguish between different sources of nonequilibrium(e.g., Brusseau et al., 1989), they seem to be of limited
practical importance when evaluating macroscopic solute fluxes.
In this study, we evaluated the ability of two modified
Freundlich isotherms to describe sorption of Cd in different horizons of a layered sandy soil. We also compared the leaching behavior of Cd that resided for decades in the soil (i.e., indigenous Cd) with the simultaneous
breakthrough curve of a 111Cd isotope tracer added tothe same soil. The indigenous Cd originates from CdO
deposition. Under acidifying conditions, CdO has beenshown to be readily soluble (Smolders, unpublisheddata, 1999). This means that both the 111Cd isotope andthe indigenous Cd are added to the soil as a soluble
metal salt. In the displacement experiments, they are in
the same chemical form (i.e., in part as a divalent cationand in part as soluble chloride complex). The transport
behavior of indigenous and tracer Cd can be readily
compared.

MATERIALS AND METHODS
Cadmium Retention in Soils
Retention of metals in soils can be described by two dominant processes: (i) the metal is sorbed to the soil matrix as areadily exchangeable ion, and (ii) the metal is sorbed with ahigh affinity to specific sites on the soil matrix (Selim et al.,1992). The second reaction (i.e., specific sorption) causes themetal to be strongly bound to the soil solid phase. This typeof reaction includes the formation of inner- and outer-spherecomplexes, surface precipitation, and incorporation of heavy
metal ions into the mineral lattices of the soil. In acid sandy
soils the first process is generally assumed to be predominantsince specific sorption of Cd is unlikely to occur (Boekhold
and van der Zee, 1992).
Reversible Cd sorption can be modeled using either a macroscopic sorption isotherm describing reversible sorption reactions, or an ion-exchange isotherm describing competition of
metal ions for the same exchange sites on the solid phase (e.g.,
between Cd2⫹ and Ca2⫹ ). The former isotherm is often used
when predicting the retention of reactive solutes under conditions of constant ionic strength and a constant backgroundsolution, while the latter isotherm is more appropriate whensimulating conditions of variable ionic strength. The sorption
isotherm relates the sorbed concentration s (mg kg⫺1 ) to thesolution concentration c (mg L⫺1 ). Cadmium sorption is frequently described using a Freundlich sorption isotherm:s ⫽ kfc nCd
(1⫺n)

n

[1]

⫺1

where kf (mg
L kg ) and nCd are empirical constants.
Based on earlier work by Chardon (1984), van der Zee andvan Riemsdijk (1987) proposed an empirical isotherm thatconsiders the effects of pH and organic matter content on
sorption. This model is referred to here as the two-species
Freundlich (2SF) model because it considers two species inthe soil solution (H⫹ and Cd2⫹ ):s ⫽ k2SF · OC · (H⫹)m c nCd
(1⫺n)

[2]

where k2SF is the 2SF sorption constant (mgL kg ), (H⫹ )is the proton activity in the soil solution: 10⫺pH (mol L⫺3 ), OC
is the organic matter content (%), and m and nCd are empiricalconstants. The exponent m (⫽⫺0.5) was derived in otherstudies (Chardon, 1984) involving sorption experiments at different pHs. The modified Freundlich sorption constant k2SF
depends on the composition of the soil solution. van der Zeeand van Riemsdijk (1987) compiled isotherms from experiments of Chardon (1984) for a variety of soils. They derived
a modified sorption constant k2SF in Eq. [2] of 0.048 mg(1⫺n )Ln kg⫺1, and a Freundlich exponent n of 0.85.
From theoretical considerations of sorption onto heterogeneous soil surfaces leading to a non-ideal competitive adsorption (NICA) model (Koopal et al., 1994), Temminghoff et al.
(1995) proposed a simplified equation for multicomponent Ca
and Cd sorption, referred to as the three-species Freundlich(3SF) model:ns ⫽ k3SFOC · (H⫹)m · (Ca)nCa · (c)nCd

⫺1

[3]

where (c ) is the Cd activity in the soil solution; (Ca) is the
Ca activity in the soil solution (mol L⫺3 ); k3SF is the modified
Freundlich sorption constant (mg(1⫺n) Ln kg⫺1 ); and nCa is anempirically derived exponent. The Cd activity in the soil solution is determined by inorganic complexation (e.g., with Cl⫺
and NO⫺3) and by the ionic strength of the soil solution.
Cadmium and Ca activities were calculated from, respectively:(c) ⫽ cfCdF

[4]

1042

J. ENVIRON. QUAL., VOL. 30, MAY–JUNE 2001

(Ca) ⫽ [Ca] fCa[5]

port equation for physical nonequilibrium transport of a sorbing solute was formulated by van Genuchten and Wierenga(1976) as follows:

[6]

⳵c
⳵s
⳵s
⳵s
⳵2c
␪m m ⫹ ␳ fp m ⫹ ␪im im ⫹ ␳(1 ⫺ fp) im ⫽ ␪mD 2m
⳵t
⳵t
⳵t
⳵t
⳵x

with:

⫺log( fCd/Ca) ⫽ 2√I

冢1 ⫹ √I⫺ 0.3I冣⳵c⳵c⳵s⫺ ␪mvm m ␳(1 ⫺ fp) im ⫹ ␪im im ⫽ ␣p (cm ⫺ cim)
⳵x
⳵t
⳵t

and:

F⫽

2⫹

1
[Cd ]
⫽
[Cdtotal] 1 ⫹ 兺q KIq[X ]q[7]

where in Eq. [4] and [5], fCd and fCa are the activity coefficientsof Cd and Ca, respectively, and F is the complexation coefficient. The activity coefficient f is calculated from the Davies
extension of the Debye–Hückel equation (Bolt and Bruggenwert, 1978) represented by Eq. [6], in which I is the ionicstrength. The complexation coefficient F is calculated from
Eq. [7] in which KIq is the formation constant of the Cd complexes with X at ionic strength I and [X] is the Cl concentration.
Formation constants of Cl complexes were taken from Sillen
and Martell (1971) and Smith and Martell (1981). The sorption
isotherm (Eq. [3]) thus explicitly accounts for solution complexation with Cl⫺, competition with Ca, and pH. All remaining processes influencing Cd sorption are lumped into
the modified Freundlich constant k3SF. Temminghoff et al.
(1995) found the parameters nCd ⫽ 0.85, nCa ⫽ ⫺0.34, m ⫽⫺0.69, and k3SF ⫽ 10⫺3.81 mol(1⫺n ) Ln kg⫺1 for combined sorption
data obtained at various pH levels and ionic strengths.

Cadmium Transport
Equilibrium Transport
One-dimensional transport of a sorbing solute through ahomogeneous soil at constant water content is described here
using the convection–dispersion equation (e.g., Lapidus andAmundson, 1952):

⳵s
⳵c
⳵c
⳵2 c
␳ ⫹ ␪ ⫽ ␪D 2 ⫺ Jw
⳵t
⳵t
⳵x
⳵x

[8]

where t is time (d), x is distance (dm), ␳ is the soil bulk density(kg L⫺1 ), ␪ is the soil water content (dm3 dm⫺3 ), D is thedispersion coefficient (dm2 dm⫺1 ), and Jw is the volumetricwater flux density (dm2 dm⫺1 ). Cadmium sorption is introduced into Eq. [8] using a nonlinear Freundlich isotherm (Eq.
[1]). Equation [8] must be solved numerically unless the isotherm is linearized and relatively simple initial and boundary
conditions are invoked. The dispersion coefficient D in Eq.
[8] can be determined experimentally by performing a nonreactive tracer experiment at the same flow velocity and measuring
its breakthrough curve in the column effluent.

[9]
where the subscripts m and im refer to the mobile and immobile regions, respectively; ␣p is a first-order mass transfer coefficient (d⫺1 ), representing the rate of solute exchange between
the mobile and immobile regions; and fp is the fraction ofsorption sites in contact with mobile water. In the case of
linear sorption (i.e., nCd ⫽ 1 in Eq. [1]), Eq. [9] can be solved
analytically and used to optimize one or several of the parameters in the transport model (Toride et al., 1995; Šimůnek etal., 1999).
Another formulation for nonequilibrium transport of asorbing solute was proposed by Selim et al. (1976). Nonequilibrium in their case was assumed to be caused by the kinetics
of chemical reactions between the soil solution and the solid
phase of the soil. A first-order kinetic rate equation, applicableto heterogeneous ion exchange, was combined with Eq. [8]to obtain the chemical nonequilibrium model. The conceptualmodel of Selim et al. (1976) divides the exchange sites of thesolid soil phase into two domains: Type 1 sites that are in
equilibrium with the sorbing solute, and Type 2 sites that are
subject to kinetic sorption. For nonlinear Freundlich sorption(Eq. [1]), the nonequilibrium transport equation becomes (Šimůnek et al., 1998):

⳵2 c
⳵␪Rc
⳵c
⫽ ␪D 2 ⫺ Jw ⫺ ␣c␳[(1 ⫺ fc)kcnCd ⫺ s]
⳵t
⳵x
⳵x
⳵s
⫽ ␣c[(1 ⫺ fc)kc nCd ⫺ s]⳵t
R ⫽ 1 ⫹ fc

␳nCd kc nCd⫺1␪

[10]

where fc is the fraction of equilibrium sites and ␣c a first-orderkinetic rate coefficient (d⫺1 ). Sorption here is assumed to be
nonlinear on both Type 1 and Type 2 sites. We used in thisstudy the numerical solution of Eq. [10] coupled with nonlinear
optimization as provided by the HYDRUS-1D code (Šimůnek
et al., 1998) to optimize the sorption constant, the Freundlich
exponent (nCd ), the rate coefficient (␣c ), and the fraction ofequilibrium sites ( fc ).

Soil
Nonequilibrium Transport
In many displacement experiments high water fluxes areoften used in order to leach strongly sorbing chemicals. Thesehigh water fluxes may induce nonequilibrium transport. Twotypes of nonequilibrium are frequently assumed (i.e., physicaland chemical nonequilibrium). The conceptual model describing physical nonequilibrium divides the soil water phase into
mobile and immobile (or stagnant) regions. Convective–
dispersive transport is confined to the mobile water phase,
while transfer of solutes into and out of immobile water isassumed to be diffusion controlled. Aggregate size, pore watervelocity, and solution concentration affect the degree of nonequilibrium (Nkedi-Kizza et al., 1983). The governing trans-

For the batch and column experiments we used sandy soil
(sandy type Humod; Soil Survey Staff, 1998) samples takenfrom different horizons of a Spodosol profile located at the
Kattenbos experimental site (Lommel, Belgium). The profile,situated in the vicinity of non-ferrous industry, is contaminatedby atmospheric deposition with heavy metals (Cd, Pb, Zn,
and Cu). The (historical) Cd contamination consisted mainlyof CdO, which is shown to be readily soluble under acidifying
conditions (Smolders, unpublished data, 1999). The soil at theexperimental site is moderately contaminated as comparedwith the contamination at the industrial site itself (i.e., up to
30 mg kg⫺1 ). The Cd concentration is representative for the
diffuse pollution of the area. The studied Spodosols represent

1043

SEUNTJENS ET AL.: AGING EFFECTS ON CADMIUM TRANSPORT IN SOIL COLUMNS

Table 1. Properties of the five soil horizons used in the experiments (values represent means of three replicates).
Horizon

Cd

Zn

pH

CEC†

mg kg⫺1
A
E
Bh
Bh/C
C

1.70
1.07
1.93
1.03
0.80

OC‡

Fe2O3

Al2O3

3.66
0.78
3.03
0.93
0.20

0.20
0.09
1.19
0.68
0.57

0.64
0.55
2.32
1.76
1.53

kg⫺1

69.7
5.9
18.2
15.0
12.0

mmolc
26.6
12.3
96.7
15.3
8.0

3.7
3.9
3.7
4.4
4.6

Sand§

Silt§

Clay§

94.7
92.0
96.7
96.2
98.2

0.28
0.38
1.04
0.17
0.33

0.12
0.17
1.11
1.33
0.24

%

† Cation exchange capacity.
‡ Organic carbon content.
§ Sand (particle size ⱖ 50 ␮m); silt (2 ⱕ 50 ␮m); clay (particle size ⬍ 2 ␮m).

about 40% of the contaminated area. Table 1 summarizes themain characteristics of the five distinct soil horizons: a humus
A horizon, a eluvial E horizon, a Spodic Bh horizon, a transitional Bh/C horizon, and a C horizon (Seuntjens et al., 1999a).
Disturbed soil samples were taken from each horizon. The
samples were dried at 105⬚C and crushed with an agate mill
before determination of total metal concentration. Total concentrations of Cd, Zn, Al, and Fe were determined using an
inductively coupled plasma–atomic emission spectrometer(ICP–AES) after digestion with HNO3–HCl–H3PO4–
HBF4–HF in a microwave oven. The remaining analyses were
performed on air-dried soil. The soil pH was measured in a
1:5 (v/v) suspension of soil in water (pH H2O). Total organiccarbon was analyzed using a total organic carbon analyzer
with infrared spectrometer after combustion in a furnace at
500⬚C (total carbon) and acidification (total inorganic carbon).
The cation exchange capacity (CEC) was determined using
BaCl2 according to an experimental procedure adapted from
Gilman (1979).

Displacement Experiments
In the field, 0.12-m-long and 0.056-m-diam. acrylic sampling
cores with sharpened edges were driven into the soil to obtain
undisturbed soil samples from the E, Bh, and Bh/C horizons.
The soil cores were capped and transferred to the laboratory
where they were placed on a porous glass filter (pore size 16␮m). A 0.001 M CaCl2 solution was applied at the top boundary
at a constant flux using an adjustable-speed peristaltic pump.
A rotating lid ensured a homogeneous distribution of thesolute at the top of the column. The pH of the solution was
adjusted to the pH of the particular soil horizon by adding
small amounts of HNO3 (see Table 2). The resulting ionicstrength of the applied solution depended on the desired pH
and varied between 0.030 and 0.031 M. In order to establish
steady-state unsaturated flow in the soil column, a constantnegative pressure head was maintained below the glass filter
using a vacuum pump with adjustable pressure. General conditions of the displacement experiments are given in Table 2.
The experiments in the Bh/C column were conducted consecutively at two different flux densities to evaluate the effect ofwater flux on Cd transport.
A fraction collector allowed the collection of effluent samples from two columns simultaneously every 20 min. Indigenous Cd, Ca, and Cl concentrations and the pH of the effluent
were measured. When the pH reached a constant value, apulse of 111Cd stable isotope (dissolved in a 1% HNO3 solution,pH 1) in a 0.001 M CaCl2 background solution was added tothe soil column and leached with 0.001 M CaCl2. To adjust
the pH to the pH of the particular soil horizon, NaOH wasadded to the solution. This resulted in a 111Cd backgroundsolution of 0.05 M NaNO3 and 0.001 M CaCl2. Indigenous Cd
and 111Cd isotope concentrations were determined simultaneously in the collected sample fractions using an inductively
coupled plasma–mass spectrometer (ICP–MS). Details of thecalibration, the analytical analysis of the different Cd isotopes,
and the accompanying uncertainty analysis can be found in
Tirez et al. (1999).
A displacement experiment with a nonsorbing solute was
next conducted on the same soil columns to obtain physicaltransport parameters. A single pulse of 0.02 M CaCl2 (for the

Batch Experiments
Upon completing the transport experiments, soil samples
were taken from the E and Bh columns to determine thesorption isotherms. For the Bh/C horizon, four field soil samples were taken surrounding the location where the column
was sampled. After homogenization, 2 g of air-dry soil was
equilibrated with seven different Cd concentrations (5, 50,
100, 250, 500, 750, and 1000 ␮g L⫺1 ) prepared in a 1:10 (v/v)
suspension of soil in 0.001 M CaCl2. The suspension was
shaken head-over-end for 24 h to ensure equilibrium between
soil and solution. After separation of the solute using a 15min centrifuging at 10 000 ⫻ g, the pH was measured in the
supernatant. Concentrations of Cd, Ca, and Zn were measured
in the acidified (2% HNO3 ) supernatant using an inductively
coupled plasma–atomic emission spectrometer (ICP–AES).
Three replicate samples were obtained by repeating the sameprocedure. The amount of Cd sorbed was calculated from thedifference between the initial concentration of the solution
and the concentration after equilibration with the soil, divided
by the soil to solution ratio.
Table 2. Summary of the transport experiments.

Cl⫺/NO⫺3§
Horizon
E
Bh
Bh/C
Bh/C

L†

q

cm
9
8
10.5
10.5

cm d⫺1
32.9
31.7
100.7
11.4

␪
3

cm⫺3

cm
0.26
0.39
0.36
0.25

vo‡

␳

cm d⫺1g cm⫺3

127
81
280
46

1.59
1.30
1.43
1.43

pH

co
mol L⫺1

3.85
3.93
4.65
4.75

0.02
0.02
0.05
0.05

† Length of the soil column.
‡ vo ⫽ q/␪: the average pore water velocity.
§ co, To: input concentration and time pulse, respectively (Cl⫺/NO3⫺ ⫽ nonsorbing solute, Cd111 ⫽ sorbing solute).

Cd111§
To

co

To

d
0.017
0.017
0.00790.044

mol L⫺1

d
0.011
0.011
0.0039
0.026

4.46
4.46
4.95
4.95

1044

J. ENVIRON. QUAL., VOL. 30, MAY–JUNE 2001

Table 3. Nonlinear least-squares fit of the Freundlich parameters. The traditional Freundlich (Eq. [1]), two-species Freundlich (2SF,
Eq. [2]), and three-species Freundlich (3SF, Eq. [3]) isotherms were used in the analysis.
Freundlich (Eq. [1])Horizon
A
E
Bh
Bh/C
C
All horizons

2SF (Eq. [2])3SF (Eq. [3])

nCd

R

2

k2SF

nCd

R

mg(l⫺n ) Ln kg⫺1
16.7
0.31
7.2
0.75
26.1
0.86
5.8
0.77
1.5
0.53
no convergence

0.83
0.99
0.98
0.98
0.82

mg(l⫺n ) Ln kg⫺1
0.040
0.104
0.152
0.038
0.039
0.046

0.77
0.80
0.86
0.75
0.75
0.94

0.98
1.00
1.00
0.98
0.63
0.86

kf

E and Bh columns) or 0.05 M NaNO3 (for the Bh/C columns)
was applied at the top of the soil columns. Effluent conductivity was continuously monitored at the outlet of the soil columns using a flow-through conductivity electrode connected

2

k3SF

nCd

R2

mol(l⫺n ) Ln kg⫺1
2.33 ⫻ 10⫺43.76 ⫻ 10⫺4
7.58 ⫻ 10⫺43.00 ⫻ 10⫺5
8.53 ⫻ 10⫺6
1.68 ⫻ 10⫺4

0.88
0.83
0.87
0.74
0.65
0.86

0.98
1.00
1.00
0.98
0.94
0.82

to a conductivity meter and a personal computer for automatic
data collection.

RESULTS AND DISCUSSION
Sorption Isotherms

Fig. 1. Cadmium sorption isotherms: (a) Freundlich isotherm forthree different horizons, (b) field-scale averaged two-species
Freundlich isotherm, and (c) field-scale averaged three-species
Freundlich isotherm.

Sorption constants and exponents of the Freundlich
isotherm (Eq. [1]), the two-species Freundlich isotherm(2SF, Eq. [2]), and the three-species Freundlich equation (3SF, Eq. [3]) were optimized using a nonlinear
least-squares procedure taken from the STATISTICA
software package (Statsoft, 1996). Because the Bh/C
soil sample contained a considerable amount of indigenous Cd, the intercept of the Freundlich isotherm wasfitted as well. During optimization, the exponents m in
Eq. [2] and m and nCa in Eq. [3] were fixed at values
obtained for similar soils by van der Zee and van Riemsdijk (1987) (i.e., m ⫽ ⫺0.5) and Temminghoff et al.
(1995) (i.e., m ⫽ ⫺0.69 and nCa ⫽ ⫺0.34), respectively.
Results of the optimizations are given in Table 3. The
Freundlich isotherm (Eq. [1]) had kf values ranging from
1.5 mg(1⫺n ) Ln kg⫺1 for the C horizon to 26.1 mg(1⫺n ) Ln
kg⫺1 for the humus Bh horizon. The isotherms werestrongly nonlinear with exponents nCd between 0.31 and
0.86. The low nCd values for the A and C horizons werepartly due to a poor fit between measured and optimized
data. Figure 1a presents the Freundlich isotherms for
the E, Bh, and Bh/C horizons (i.e., those used in thedisplacement experiments). The negative intercept of
the Bh/C isotherm (i.e., ⫺0.45 mg kg⫺1 ) indicates desorption at low concentrations.
The optimized Freundlich exponents nCd of the 2SFmodel (Eq. [2]) ranged from 0.75 in the E horizon to 0.86
in the Bh horizon. Except for the humus Bh horizon, nCdvalues were considerably smaller than the value of 0.85
obtained by van der Zee and van Riemsdijk (1987). The
A horizon produced a significantly higher coefficient ofdetermination (R2 ) for the 2SF isotherm as comparedwith the Freundlich isotherm, while R2 for the C horizon
was much lower (Table 3). These results suggest thatdifferences in soil organic matter and pH explain thevariance between the different Freundlich isotherms for
the A horizon, but not for the C horizon. The low R2
for the C horizon indicates that the variance of thesorption constant cannot be adequately explained byvariations in soil organic matter and pH.
The optimized sorption parameter k3SF of the threespecies Freundlich isotherm varied from 8.53 ⫻ 10⫺6 for
the C horizon to 7.58 ⫻ 10⫺3 mol(1⫺n ) Ln kg⫺1 for the1045

SEUNTJENS ET AL.: AGING EFFECTS ON CADMIUM TRANSPORT IN SOIL COLUMNS

Table 4. Nonlinear least-squares fits of the parameters vm, Dm, ␪m, and ␣p in the two-region model (TRM) and the parameters v and D
in the convection-dispersion model (CDE) (values in brackets are standard errors of estimate).
CDE
Horizon
E
Bh
Bh/C
Bh/C

v

D

cm d⫺12

134(1.45)
99(0.46)
349(2.98)
31(0.20)

d⫺1

cm
23.3
(2.23)
19.3(0.67)
111(6.48)
11.6(0.45)

TRM
␭†

R2

cm
0.17

0.96

0.19

0.99

0.32

0.94

0.37

0.96

v m‡

D‡

cm d⫺12

126(1.49)
98(0.55)
308(2.42)
29(0.12)d⫺1cm
6.7(1.97)
1.9(9.5)
28(2.21)
3.7(0.21)

␭m§
cm
0.06
0.03
0.12
0.16

␪m/␪

␻

␣p§

R2

d⫺ 1
0.84
(0.019)0.65(0.23)
0.78(0.006)0.79(0.006)0.94(0.25)
5.6
(4.5)
0.78
(0.047)1.02(0.065)

2.90

0.99

2.54

0.99

7.78

0.99

0.64

0.99

† ␭ ⫽ D/v; ␭m ⫽ Dm/vm.
‡ vm ⫽ v␪/␪m, Dm ⫽ D␪/␪m.
§ ␣p ⫽ ␻␪v/L.

Bh horizon. This large variation in sorption constants
between soil layers indicates that Ca is presumably not
the only competing metal ion for the same sorption
sites. Based on the coefficients of determination, the3SF model describes sorption of Cd in the C horizon
significantly better than the 2SF model.
Isotherm data from all horizons (A, E, Bh, Bh/C, andC) were pooled in order to derive field-scale averaged
2SF and 3SF isotherms. The average 2SF and 3SF isotherms are shown in Fig. 1b and 1c, respectively. Resultsof the optimization suggest that 86% of the variance ofthe sorption constant k2SF between samples from different layers is explained by the organic matter content
and the pH. The average 3SF model, which additionally
includes Ca activity, ionic strength, and choride complexation, does not seem to predict sorbed concentrations better than the average 2SF model. The valueof 0.046 for the modified 2SF sorption constant k2SF
corresponds remarkably well with the modified Freundlich constant (k2SF ⫽ 0.048 mg(1⫺n ) Ln kg⫺1 ) derived byvan der Zee and van Riemsdijk (1987). Optimizing k3SF
and nCd of the modified 3SF equation to isotherm datafrom all horizons resulted in a modified k3SF of 10⫺3.77
mol(1⫺n ) Ln kg⫺1 and in a nCd value of 0.86. The modified
3SF constants are close to the values obtained by Temminghoff et al. (1995), assuming that the organic carbon
content in our soil is equal to 1%. Depth averages oforganic carbon content in our soil (i.e., OC ⫽ 1.01%)
indicate that this assumption is indeed valid. The results
show that sorption of cadmium in our layered sandy
soil can be reasonably well approximated by modified
isotherms valid for the entire soil profile.

Transport Experiments
Nonreactive Transport
The parameters v and D of the convection–dispersion
model (CDE) (Eq. [8] with ⳵s/⳵t ⫽ 0), as well as theparameters vm, D, and ␪m of the two-region model(TRM) (Eq. [9], with ⳵sm/⳵t ⫽ 0 and ⳵sim/⳵t ⫽ 0), were
optimized using the chloride breakthrough data for the
E and Bh horizons, and the nitrate breakthrough data
for the Bh/C horizon. The optimized CDE and TRM
parameters obtained using the CXTFIT curve-fitting
program (Toride et al., 1995) are given in Table 4. The
dispersivity (␭) calculated from the CDE parameters v

and D for the Bh/C horizon was somewhat larger forthe lower water content experiment. This is consistent
with observations by Biggar and Nielsen (1960), whosuggested that the flow path should become more tortuous when the soil water content decreases, thus leading
to higher dispersivities. Quadruplicate displacement experiments at three different fluxes for four different
horizons of the same soil confirmed this inverse relationship between the water content and dispersivity (results
not further shown here).
Figure 2 shows that the CDE model slightly underestimated the peak concentrations and poorly described
the tailing parts of the breakthrough curves for the E
and Bh/C horizons, while the TRM model describedboth the peak and the tailing well. Transport in the Bhhorizon was equally well described using the CDE andTRM models. Relatively large values of the dimensionless mass transfer coefficient (␻) (typically larger
than 10) are indicative of local equilibrium during transport, in which case Eq. [4] should be the appropriate
transport model. The optimized ␻ values shown in Table
4 are smaller than or equal to 1, except for the Bhhorizon. This indicates a significant degree of nonequilibrium in the E and Bh/C horizons for the imposedwater fluxes. The relative mobile water contents (␪m/␪)ranged from 0.65 to 0.84. These values are comparablewith those obtained in other studies on sandy soils(Gaudet et al., 1977; Vanclooster et al., 1993). Calculated mass transfer coefficients (␣p ) increased with increasing pore water velocity (v), which confirms findings
by Nkedi-Kizza et al. (1983), Reedy et al. (1996), andPang and Close (1999), among others. A higher porewater velocity or, equivalently, a shorter solute residence time, promotes incomplete mixing of solutes between mobile and stagnant pore water, and a larger
diffusive mass transfer rate coefficient ␣p.
In the following sections we will discuss Cd transport
for the same columns and the same boundary conditions
as the above tracer experiments. The v and D parameters, optimized using the (macroscopic) CDE model,
will be used to predict breakthrough of the sorbing
solute. Since physical nonequilibrium conditions prevailed in the E and Bh/C horizons, D optimized using
the CDE model is considered to be an effective model
parameter for these columns. The effective dispersion
coefficient lumps two processes (e.g., Valocchi, 1985):

Fig. 2. Observed chloride (E, Bh) and nitrate (Bh/C) breakthroughdata and fitted curves using the convection–dispersion (CDE) andtwo-region (TRM) models for (a) the E horizon, (b) the Bh horizon,(c) the Bh/C horizon with high flux, and (d ) the Bh/C horizon withlow flux.

Fig. 3. Non-aged 111Cd breakthrough data for (a) the E horizon, (b)the Bh horizon, (c) the Bh/C horizon with high flux, and (d ) the
Bh/C horizon with low flux. Curves represent predictions of 111Cd
concentrations using the convection–dispersion (CDE) model withindependently determined parameters (predicted local equilibrium[LEA]), using the optimized two-site model (TSM) with k and ndetermined from batch experiments (fitted two-site model [TSM],
batch k and n ), and using the optimized TSM model with k andn also fitted (fitted TSM, optimized k and n ).

SEUNTJENS ET AL.: AGING EFFECTS ON CADMIUM TRANSPORT IN SOIL COLUMNS

1047

dispersion owing to variations in the pore water velocity(vm ) and additional dispersion due to limited mass transfer between stagnant and mobile water phases.
Equilibrium Cadmium Transport
Measured breakthrough curves (BTCs) of 111Cd areshown in Fig. 3. The BTCs are asymmetrical, displaying
a relatively sharp breakthrough front and prolonged
tailing. Recovery rates calculated from observed 111Cd
concentrations in the column effluent were 96, 80, 66,
and 72% for the E and Bh horizons, and the Bh/C
horizons at the high and low fluxes, respectively. Given
the incomplete observed effluent curves for columns Bhand Bh/C, the recovery percentages indicate that the111
Cd isotope probably sorbs reversibly.
Independent predictions of the 111Cd breakthroughcurves (denoted as predicted local equilibrium [LEA]
in Fig. 3) were made using the dispersivity ␭ obtained
from the nonreactive tracer experiments and kf and nCdparameters derived from batch data using the Freundlich isotherm (Eq. [1]). The calculations were performed
using the HYDRUS-1D code (Šimůnek et al., 1998). As
explained before, physical nonequilibrium conditions
existed in the E and Bh/C horizons. Therefore, ourcalculations assume chemical equilibrium and physicalnonequilibrium in the E and Bh/C horizons, and chemical and physical equilibrium in the Bh horizon.
As shown in Fig. 3, the chemical equilibrium model
generally failed to predict the observed breakthroughdata of the 111Cd isotope. The mean breakthrough time
was substantially overestimated, while spreading of the
breakthrough curve was underestimated. The 111Cd isotope seemed to travel much faster than predicted withthe equilibrium model. The rather severe deviations
between the observed and predicted Cd breakthroughcurve are presumably caused by (i) overestimation ofretardation using the batch sorption parameters and/or(ii) inability of the chemical to reach equilibrium withthe solid phase during transport through the soil. Twomechanisms may explain this observed nonequilibrium(van Genuchten and Cleary, 1982). Chemical nonequilibrium can be due to kinetic sorption (i.e., the sorption
reaction is too slow compared with the transport velocity
to ensure chemical equilibrium between sorbate andsorbent). Second, physical constraints during transport
may prevent a chemical from reaching the exchangesites in the soil. The diffusion rate of the chemical fromthe mobile water fraction to the stagnant water surrounding the sorption sites is then slow compared withthe transport velocity in the mobile water fraction. Realchemical kinetic nonequilibrium may explain the fast
breakthrough in the Bh horizon, because no physicalnonequilibrium was demonstrated by the displacement
experiment using the nonsorbing tracer. In the othersoil columns, both types of nonequilibrium probably
were operative at the same time. The equilibrium model
predicted 111Cd breakthrough at the low flux in the
Bh/C horizon relatively well, although still relatively
fast breakthrough and tailing occurred.
Leaching of indigenous Cd is presented in Fig. 4.

Fig. 4. Indigenous aged Cd breakthrough data for (a) the E horizon,(b) the Bh horizon, and (c) the Bh/C horizon with high flux. Curves
represent predictions of Cd concentrations using the convection–
dispersion (CDE) model with independently determined parameters (predicted local equilibrium [LEA]), using the optimized twosite model (TSM) with k and n determined from batch experiments(predicted TSM, batch k and n ), and using the optimized TSM
model with k and n also fitted (predicted TSM, optimized k andn ). The TSM model parameters were optimized using the 111Cd
breakthrough curve.

The effluent curves initially show a sharp decrease inconcentration and a more gradual decrease toward theend of the experiments. The increase in Cd concentrations between 50 and 100 pore volumes is due to increased ionic strength resulting from the 111Cd tracer
pulse (0.05 M NaNO3, 0.001 M CaCl2 background). The
increase is followed by sharp decrease due to increased

1048

J. ENVIRON. QUAL., VOL. 30, MAY–JUNE 2001

Table 5. Two-parameter optimizations of the nonequilibrium parameters ␣c and ␤c using the two-site kinetic model and independently
measured parameters kf and nCd obtained from the batch experiments (values between parentheses are standard errors of estimate).
111

Fixed
HorizonE
Bh
Bh/C-high flux
Bh/C-low flux

Cd

kf

nCd

fc

␣c

␮g(l⫺n ) Ln kg⫺1
40.5
(6.7)
68.6
(7.5)
28.3
(9.5)
28.3
(9.5)

0.75
(0.03)
0.86(0.02)
0.77(0.05)
0.77(0.05)

0.60(0.010)0.43(0.013)0.41(0.007)0.63(0.005)d ⫺1
0.08(0.014)0.19(0.010)0.90(0.030)0.08(0.002)sorption of Cd and release of Na. This phenomenon is,
however, of short duration (except for column Bh) anddoes not seem to affect leaching of indigenous Cd. Recovery percentages, calculated from the ratio of the
cumulative Cd mass in the column effluent and theinitial total Cd mass in the soil columns, were 30, 38,
54, and 87% for the E, Bh, and high-flux and low-flux
Bh/C soil columns, respectively.
Predictions made using the equilibrium model (i.e.,
assuming local equilibrium [LEA]) strongly deviatedfrom the observed transport data. The chemical equilibrium model clearly failed to predict the sharp concentration decrease and overestimated the mass leachedwithin the calculated period. Apparently, desorption
and/or mass transfer of Cd to the transport flow regions
is too slow to assume equilibrium. We note here thatthe boundary conditions for the low-flux leaching caseof column Bh/C were different from those of the othercolumns. As shown in Fig. 4, the gradual decrease inconcentration between 150 and 300 pore volumes (lowflux leaching) is preceded by a sudden increase resulting
from a 15-d interrupted flow period between the highand the low-flux cases in the same column. Therefore,
the low-flux experiment cannot be considered as an
initial value problem similar to the other columns, andno predictions were made for this case.
Nonequilibrium Cadmium Transport
From the available observations we could not conclude conclusively whether irreversible sorption or nonequilibrium due to rate-limited diffusion or kinetic sorption was responsible for the low recovery percentages
and the overall poor prediction of the Cd breakthroughcurves using the equilibrium model. We pointed out
earlier that no clear evidence exists in the literaturethat irreversible Cd fixation occurs in acid sandy soils(Smolders et al., 1999). Therefore, in our calculations
below we will assume that nonequilibrium is responsiblefor deviations between predicted and observed Cd
breakthrough curves. Several researchers have demonstrated that actual Cd sorption is instantaneous in wellstirred and sieved soils (e.g., Christensen, 1984a). However, in systems with flowing water, and especially inaggregated or unsaturated soils, the global reaction rate
can be limited by the rate at which ions are transported
to the exchange sites, even when the actual sorption
reaction is instantaneous (Nkedi-Kizza et al., 1984). In
the next section we describe Cd transport using theIndigenous Cd
2

fc

␣c

0.81

0.34
(0.005)0.37(0.008)0.48(0.004)d ⫺1
0.28(0.014)0.15(0.014)1.21(0.040)

R

0.91
0.85

R2
0.99
0.94
0.96

0.98

chemical nonequilibrium transport model (Eq. [10]),which assumes that sorption is reversible, but that the
global reaction is not instantaneous at some fraction of
the sorption sites.
The fraction of sites fc for which sorption is instantaneous, and the first-order kinetic rate coefficient ␣c,were optimized using a numerical solution of Eq. [10](Šimůnek et al., 1998), while fixing the Freundlich parameters kf and nCd to the independently measured batchvalues. Results of the optimization are given in Table
5 for both the 111Cd isotope and the indigenous Cd.
Figure 3 shows that the nonequilibrium model adequately predicts the position of the 111Cd breakthroughcurve. The model, however, failed to correctly describe
the early breakthrough data and tailing for the E, Bh,and Bh/C (high flux) columns. Optimizing fc and ␣c forthe low-flux breakthrough curve for the Bh/C column
resulted in an almost perfect prediction of the adsorption front of the BTC, whereas tailing was somewhat
overestimated.
The description of the 111Cd breakthrough curve improved significantly when all four parameters (kf, nCd,
␣c, fc ) in Eq. [10] were optimized. Figure 3 illustrates
the close correspondence between model results andexperimental data. The optimized parameters are given
in Table 6. The optimized sorption constants kf fromthe column experiments were consistently smaller than
the kf values obtained from the batch experiments (Table 2), whereas the optimized exponents nCd were larger
than the batch values, except for the Bh horizon. Basedon the standard errors of estimate of measured andoptimized kf and nCd parameters, differences between
kf and nCd optimized using the breakthrough data andfitted to the batch data seem to be rather small. Only
for the E column did we find that the optimized nCdfrom the column breakthrough data was significantly
larger than the batch value. We believe that results ofthe four-parameter fit for this column are unreliable
and should be interpreted with care.
Independent predictions of indigenous Cd leaching
using model parameters obtained from the displacement
tests with 111Cd isotope are displayed in Fig. 4. Independent predictions were made using the batch kf and nCdvalues and fc and ␣c values estimated from the 111Cd
breakthrough curve (Table 5, depicted as predicted twosite model [TSM] batch k,n in Fig. 4), as well as using
all kf, nCd, ␣c, and fc parameters fitted to the 111Cd breakthrough curve (Table 6, depicted as predicted TSM opti-

SEUNTJENS ET AL.: AGING EFFECTS ON CADMIUM TRANSPORT IN SOIL COLUMNS

Table 6. Four-parameter optimizations of the nonequilibrium parameters ␣c and ␤c and sorption parameters kf and nCd using
the two-site kinetic model (values between parentheses are
standard errors of estimate).
Horizonkf(l⫺n )

␮gE
Bh
Bh/C-high flux
Bh/C-low flux

nCd
n

fc

kg⫺1

L
30.9
(0.83)
55.0
(2.5)
20.1
(0.17)
24.9(0.49)

␣c

R2

d⫺1
1.06
(0.04)
0.81(0.04)
0.90(0.001)
0.82(0.02)

0.69(0.02)
0.40(0.04)
0.44(0.004)0.62(0.01)

0.29(0.04)
0.31(0.04)
1.74(0.04)
0.10(0.008)

0.98
0.98
0.99
0.99

mized k,n in Fig. 4). A comparison between the observed
and as such predicted indigenous Cd leaching should
provide evidence whether or not the transport behaviorof 111Cd and indigenous Cd is similar. Fig. 4 shows thatleaching of indigenous Cd in the E horizon cannot be
described using parameters determined from the 111Cd
isotope. Much better predictions of Cd leaching wereobtained using the batch kf and nCd values (two-parameter fit) than when those two parameters were optimized
to the 111Cd breakthrough data (four-parameter fit). Thisagain shows that the four-parameter optimization may
yield unreliable results in this case. The results alsosuggest that transport of 111Cd in the E column is closer
to equilibrium than leaching of indigenous Cd. This isconfirmed by the larger fraction of equilibrium sites fc
for 111Cd than for indigenous Cd (Table 5). For the Bh
and Bh/C horizons, indigenous Cd leaching was described adequately using 111Cd parameters. No differences between the transport behavior of aged and nonaged Cd in these two layers could be inferred from theavailable observations.
The kinetic rate coefficient ␣c was smaller at low water
fluxes than at large water fluxes. Similar observationsof the dependence of the mass transfer rate coefficienton boundary conditions were reported by Gaber et al.
(1995) for atrazine. This indicates that diffusion, rather
than sorption kinetics, probably limits the transfer fromexchange sites to mobile soil water.

1049

Cd concentrations. Modifying the 3SF isotherm did not
improve the predictions as compared with the 2SF isotherm. Additional experiments at the scale of an entire
soil profile are needed to evaluate the usefulness ofmodified isotherms in predicting transport through layered soil profiles.
Unsaturated conditions during the displacement experiments for the two soil layers (i.e., the eluvial E
horizon and the transitional Bh/C horizon) produced
physical nonequilibrium. Mobile water fractions were
found to be about 80% of the total water content, indicating that 20% of the soil water did not contribute totransport. Physical equilibrium prevailed in the lowpermeability humus Bh horizon, which previously had
been shown to redistribute solute in the soil profile(Seuntjens et al., 1999b). The mass transfer coefficient␣p depended on water flux, thus showing that diffusion
limited the solute exchange between mobile and stagnant soil water.
The local equilibrium convection–dispersion model
failed to describe breakthrough curve of 111Cd isotope
and leaching of indigenous Cd. Parameters optimization
using the two-site kinetic nonequilibrium model showedthat retardation was somewhat overestimated by thebatch parameters, and that the global sorption rate was
not instantaneous at about 30 to 60% of the sorption
sites. Mass transfer coefficients were found to dependon the water flux, which suggests that diffusion rather
than kinetic sorption was responsible for the observed
nonequilibrium.
Predictions of Cd leaching using parameters representative for non-aged Cd compared favorably with observed leaching of aged Cd in the humus Bh horizon andthe transitional Bh/C horizon. The transport behavior of
non-aged and aged Cd hence was about the same inthese two layers. Aging seemed to affect the leaching
behavior of Cd in the E horizon as shown by a disparity
in optimized model parameters between aged and nonaged Cd.
ACKNOWLEDGMENTS

CONCLUSIONS

The authors wish to thank Dirk Mallants for revisions onan earlier version of this manuscript.

Sorption and transport of Cd in undisturbed unsaturated soil horizons of a multilayered sandy soil wasinvestigated. Sorption of Cd was nonlinear and could
be described well using a Freundlich isotherm. A twospecies Freundlich (2SF) isotherm, accounting for the(horizontal) variance of organic matter and pH, improved the prediction of sorption in the humus A horizon, but failed to predict sorption in the deeper C
horizon. A three-species Freundlich (3SF) isotherm adequately described sorption in all horizons. The 3SF isotherm predicted sorption in the C horizon better ascompared with the 2SF isotherm, thus showing that soil
water composition has to be taken into account whenestimating sorption of Cd in the deeper soil layers.
Modified 2SF isotherms that accounted for the variability of organic matter and pH in the entire soil profile
were found to explain 86% of the variance of sorbed

Ainsworth, C.C., J.L. Pilou, P.L. Gassman, and W.G. Van der Sluys.
1994. Cobalt, cadmium and lead sorption to hydrous iron oxide:
Residence time effect. Soil Sci. Soc. Am. J. 58:1615–1623.
Baker, R.S., and D. Hillel. 1990. Laboratory tests of a theory offingering during infiltration into layered soils. Soil Sci. Soc. Am.
J. 54:20–30.
Biggar, J.W., and D.R. Nielsen. 1960. Diffusion effects in saturated
and unsaturated porous materials. J. Geophys. Res. 65:2887–2891.
Boekhold, A.E., and S.E.A.T.M. van der Zee. 1992. A modified sorption model validated at the column scale to predict cadmium contents in a spatially variable field soil. Soil Sci. 154:105–112.
Bolt, G.H., and M.G.M. Bruggenwert. 1978. Soil chemsitry A: Basicelements. Elsevier Scientific Publ. Company, Amsterdam, theNetherlands.
Brusseau, M.L., R.E. Jessup, and P.S.C. Rao. 1989. Modeling thetransport of solutes influenced by multiprocess nonequilibrium.
Water Resour. Res. 25:1971–1988.
Büchter, B., C. Hinz, M. Gfeller, and H. Flühler. 1996. Cadmium

REFERENCES

1050

J. ENVIRON. QUAL., VOL. 30, MAY–JUNE 2001

transport in an unsaturated stony subsoil monolith. Soil Sci. Soc.
Am. J. 60:716–721.
Chardon, W.J. 1984. Mobiliteit van cadmium in de bodem. Ph.D. diss.
Agric. Univ., Wageningen, the Netherlands.
Christensen, T.H. 1984a. Cadmium soil sorption at low concentrations:I. Effect of time, cadmium load, pH and calcium. Water Air Soil
Pollut. 21:105–114.
Christensen, T.H. 1984b. Cadmium soil sorption at low concentrations:
II. Reversibility, effect of changes in solute composition and effectof soil aging. Water Air Soil Pollut. 21:115–125.
Gaber, H.M., W.P. Inskeep, S.D. Comfort, and J.M. Wraith. 1995.
Nonequilibrium transport of atrazine through large intact soil cores.
Soil Sci. Soc. Am. J. 59:60–67.
Gaudet, J.P., H. Jégat, G. Vachaud, and P.J. Wierenga. 1977. Solutetransfer with exchange between mobile and stagnant water throughunsaturated sand. Soil Sci. Soc. Am. J. 41:665–671.
Gilman, G.P. 1979. A proposed method for the measurement of exchange prpoerties of highly weathered soils. Aust. J. Soil Res.
17:129–139.
Glass, R.J., T.S. Steenhuis, and J.Y. Parlange. 1989. Mechanism forfinger persistence in homogeneous unsaturated porous media: Theory and verification. Soil Sci. 148:60–70.
Hinz., C., and H.M. Selim. 1994. Transport of zinc and cadmium insoils: Experimental evidence and modeling approaches. Soil Sci.
Soc. Am. J. 58:1316–1327.
Koopal, L.K., W.H. van Riemsdijk, J.C.M. de Wit, and M.F. Benedetti.
1994. Analytical isotherm equations for multicomponent adsrption
to heterogeneous surfaces. J. Colloid Interface Sci. 166:51–60.
Kung, K.J.S. 1993. Laboratory observation of funnel flow mechanism
and its influence on solute transport. J. Environ. Qual. 22:91–102.
Lapidus, L., and N.R. Amundson. 1952. Mathematics of adsorption
in beds. VI. The effect of longitudinal diffusion in ion exchange
and chromatographic columns. J. Phys. Chem. 56:984–988.
Nkedi-Kizza, P., J.W. Biggar, H.M. Selim, M.Th. van Genuchten,
P.J. Wierenga, J.M. Davidson, and D.R. Nielsen. 1984. On theequivalence of two conceptual models for describing ion exchange
during transport through an aggregated oxisol. Water Resour.
Res. 20:1123–1130.
Nkedi-Kizza, P., J.W. Biggar, M.Th. van Genuchten, P.J. Wierenga,
H.M. Selim, J.M. Davidson, and D.R. Nielsen. 1983. Modeling
tritium and chloride-36 transport an aggregated Oxisol. Water Resour. Res. 19:691–700.
Pang, L., and M. Close. 1999. Field-scale physical nonequilibriumtransport in an alluvial gravel aquifer. J. Contam. Hydrol. 38:
447–464.
Reedy, O.C., P.M. Jardine, G.V. Wilson, and H.M. Selim. 1996. Quantifying the diffusive mass transfer of non-reactive solutes in columnsof fractured sprolite using flow interruption. Soil Sci. Soc. Am.
J. 60:1376–1384.
Ritsema, C.J., L.W. Dekker, J.M.H. Hendrickx, and W. Hamminga.
1993. Preferential flow mechanism in a water repellent sandy soil.
Water Resour. Res. 29:2183–2193.
Selim, H.M., B. Büchter, C. Hinz, and L. Ma. 1992. Modeling thetransport and retention of cadmium in soils: Multireaction andmulticomponent approaches. Soil Sci. Soc. Am. J. 56:1004–1015.
Selim, H.M., J.M. Davidson, and R.S. Mansell. 1976. Evaluation ofa two-site adsorption–desorption model for describing solute transport in soil. p. 444–448. In Proc. of the Computer Simulation Conf.,
Washington, DC. Am. Inst. of Chem. Eng., Washington, DC.
Seuntjens, P., C. Cornelis, N. De Brucker, and P. Geuzens. 1999a.
Derivation of functional layers in a podzol toposequence for simulating cadmium transport. Phys. Chem. Earth 24:869–873.

Seuntjens, P., D. Mallants, N. Toride, C. Cornelis, and P. Geuzens.
1999b. Grid lysimeter study of nonreactive solute transport throughtwo distinct multi-layered Spodosols using TDR and wick samplers.
p. 284–293. In J. Feyen and K. Wiyo (ed.) Modelling of transport
processes in soils at various scales in time and space. Int. Workshopof EurAgEng’s Field of Interest on Soil and Water, Leuven, Belgium. 24–26 Nov. 1999. Wageningen Pers, Wageningen, the Netherlands.
Sillen, L.G., and A.E. Martell. 1971. Stability constants of metal ion
complexes. Suppl. no. 1. The Chemical Society, Burlington House,
London.
Šimůnek, J., K. Huang, and M.Th. van Genuchten. 1998. The HYDRUS code for simulating the one-dimensional movement of water, heat and multiple solutes in variably saturated media. Res.
Rep. no. 144. U.S. Salinity Lab., USDA, Riverside, CA.
Šimůnek, J., M.Th. van Genuchten, M. Šejna, N. Toride, and F.J. Leij.
1999. The STANMOD computer software for evaluating solutetransport in porous media using analytical solutions of convection–
dispersion equation. Versions 1.0 and 2.0. IGWMC-TPS-71. Int.
Ground Water Modeling Center, Colorado School of Mines,
Golden, CO.
Smith, L.G., and A.E. Martell. 1981. Critical stability constants. Volume 4: Inorganic complexes. Plenum Press, New York.
Smolders, E., K. Brans, A. Földi, and R. Merckx. 1999. Cadmium
fixation in soils measured by isotopic dilution. Soil Sci. Soc. Am.
J. 63:78–85.
Soil Survey Staff. 1998. Keys to soil taxonomy. 7th ed. USDA Natural
Resour. Conserv. Serv.
Statsoft. 1996. Statistica for Windows. Computer program manual.
Statsoft, Tulsa, OK.
Temminghoff, E., S.E.A.T.M. van der Zee, and F.A.M. de Haan. 1995.
Speciation and calcium competition effects on cadmium sorption by
sandy soils at various pH levels. Eur. J. Soil Sci. 46:649–655.
Tirez, K., P. Seuntjens, and N. De Brucker. 1999. Full uncertainty
calculation on quantitative determination of tracer (111Cd) cadmium and natural cadmium in soil column effluents with ICP–MS.
J. Anal. Atomic Spectrom. 14:1475–1484.
Toride, N., F.J. Leij, and M.Th. van Genuchten. 1995. The CXTFIT
code for estimating transport parameters from laboratory or field
tracer experiments. Version 2.0. Res. Rep. no. 137. U.S. Salinity
Lab., USDA, Riverside, CA.
Valocchi, A.J. 1985. Validity of the local equilibrium assumption formodeling sorbing solute transport through homogeneous soils. Water Resour. Res. 21:808–820.
Vanclooster, M., D. Mallants, J. Diels, and J. Feyen. 1993. Determining
local-scale solute transport parameters using time domain reflectometry. J. Hydrol. 148:93–107.
van der Zee, S.E.A.T.M., and W.H. Van Riemsdijk. 1987. Transportof reactive solutes in spatially variable soil systems. Water Resour.
Res. 23:2059–2069.
van Genuchten, M.Th., and R.W. Cleary. 1982. Movement of solutes
in soil: Computer-simulated and laboratory results. p. 349–386.
In G.H. Bolt (ed.) Soil chemistry. B. Physico–chemical models.
Elsevier Scientific Publ. Company, Amsterdam, the Netherlands.
van Genuchten, M.Th., and P.J. Wierenga. 1976. Mass transfer studies
in sorbing porous media: I. Analytical solutions. Soil Sci. Soc. Am.
J. 40:473–481.
Wilkens, B.J. 1997. Evidence for groundwater contaimnation by heavy
metals through soil passage under acidifying conditions. Ph.D. thesis. Universiteit Utrecht, the Netherlands.