Available online at www.sciencedirect.com

Geoderma 145 (2008) 449 – 461
www.elsevier.com/locate/geoderma

Modelling coupled water flow, solute transport and geochemical reactions
affecting heavy metal migration in a podzol soil
D. Jacques a,⁎, J. Šimůnek b , D. Mallants a , M.Th. van Genuchten cabPerformance Assessment Unit, SCK•CEN, Boeretang 200, B-2400 Mol, Belgium
Department of Environmental Sciences, A135 Bourns Hall, University of California Riverside, 900 University Avenue, Riverside, CA 92521, USAcU.S. Salinity Laboratory, USDA, ARS, 450 W. Big Spring Rd, Riverside, CA, 92507, USA
Received 5 April 2007; received in revised form 5 September 2007; accepted 29 January 2008
Available online 10 March 2008

Abstract
Many or most subsurface pollution problems at the field scale involve such simultaneous processes as water flow, multicomponent solutetransport, heat transport and biogeochemical processes and reactions. Process-based models that integrate these various processes can be valuable
tools for investigating the mobility of a wide range of inorganic and organic contaminants subject to different hydrologic and geochemical
conditions. The HP1 reactive transport simulator, obtained by weak coupling of HYDRUS-1D and PHREEQC-2, was developed and designed toaddress multicomponent geochemical transport processes in the vadose zone. In this paper we discuss a hypothetical HP1 application involving
the transport of major cations and heavy metals in a soil during transient flow over a period of 30 years. Results show that variations in watercontent and water fluxes can significantly influence the speciation, and thus the mobility and availability, of elements. Decreasing water contents
near the soil surface lowered pH of the soil solution and produced new cation exchange equilibrium conditions. The upward transport of Cl during
summer due to increased evapotranspiration, and subsequent accumulation of Cl near the soil surface, caused an increase in the total aqueous Cd
concentration because of the formation of Cd–Cl complexes. Coupled reactive transport codes for the unsaturated zone, such as HP1, arepromising tools to unravel the complex interaction between soil physical and biogeochemical processes for all kinds of problems, including theimpact of natural processes and antropogenic activities on soil evolution.
© 2008 Elsevier B.V. All rights reserved.
Keywords: Reactive transport modelling; Unsaturated zone; Heavy metals

1. Introduction
Process-based reactive transport models can be valuable tools
for studying the subsurface fate and transport of organic or
inorganic contaminants, including radionuclides, nutrients andpesticides, that may be subject to a large number of often simultaneous interactive physical, chemical, mineralogical, geological,
and biological processes (e.g., Mayer et al., 2002; Metz et al.,
2002; Davis et al., 2004; Lichtner et al., 2004; Appelo andPostma, 2005; NRC, 2006). Another application involving longterm reactive transport is the genesis of a soil profile from parent
material which, besides weathering, includes aqueous complexa⁎ Corresponding author.
E-mail address: djacques@sckcen.be (D. Jacques).
0016-7061/$ - see front matter © 2008 Elsevier B.V. All rights reserved.
doi:10.1016/j.geoderma.2008.01.009

tion and the transport and sorption of chemical species andvarious clay or organic matter particles. These processes influencemineral precipitation and dissolution, which are the mainprocesses of pedogenesis in the long term (Jenny, 1941; Arkley,
1963; Hoosbeek and Bryant, 1992; Presley et al., 2004;
Rasmussen et al., 2005). The concentration of different elements
will determine whether parent materials dissolve and secondary
minerals form, depending on their saturation status. Someelements also play a role in the solubility or precipitation oforganic ligands, or the flocculation of clay particles. Changes in
the chemical composition or pH of the soil solution may impact
the sorption of elements on organic matter or iron oxides and also
affect the saturation index of minerals. Simulation of these andrelated processes requires a coupled reactive transport code that
integrates the physical processes of water flow and advective–

450

D. Jacques et al. / Geoderma 145 (2008) 449–461

dispersive solute transport with a range of biogeochemical
processes.
An excellent example of the interplay between geochemical
conditions and the mobility of elements was given by Jansen
et al. (2002) involving the binding of Al and Fe to (mobile)organic matter. Jansen et al. (2002) investigated factors
controlling soluble complexation between Fe/Al and dissolved
organic matter. Small pH changes (between 3.5 and 4.5) were
found to greatly affect soluble complexation of aluminium,
ferrous and ferric iron with dissolved organic matter. At a low
metal–organic matter ratio, complexation increased withincreasing pH due to proton dissociation of the functionalgroups on the organic matter. However, soluble complexation
decreased with increasing pH at a high metal–organic matterratio due to the formation of insoluble metal–organic mattercomplexes. Redox conditions are also important for Fe(II) andFe(III) reactions. As Fe(III)-dissolved organic matter complexes
are more stable but less soluble than Fe(II) complexes, Jansen
et al. (2002) were able to show that the Fe(III)-dissolved organic
matter complexes prevail because of their higher stability at
relatively low metal–organic matter ratios, whereas Fe(II)dissolved organic matter complexes dominate at higher ratios
due to the formation of insoluble Fe(III) complexes.
Since element speciation affected by the above and similarinteractions significantly influences the mobility and bioavailability of the elements, geochemical speciation calculations
should be integrated into the classical vadose zone flow andtransport models. At the same time, soils are inherently open
systems with changing air and water contents, fluctuating groundwater tables and space-time variations in carbon dioxide
concentrations. As will be illustrated in this paper, such transient
flow conditions also induce changes in pH and speciation (see also
Jacques et al., 2006b). To account for these effects, it is necessary
to couple all relevant physical and biogeochemical processes that
may occur in soil systems. Physical processes involve (i) transientwater flow in variably-saturated porous media, (ii) transport ofmultiple components, and (iii) heat transport. The conceptualmodels for each of these processes are typically written in terms ofpartial differential equations and then solved numerically subject
to certain initial and boundary conditions. The equations generally also involve a set of constitutive relationships for selectedvariables or soil properties (e.g., the unsaturated soil hydraulicproperties). The definition of a network of (bio)geochemical
reactions similarly allows to incorporate the current stage ofknowledge on relevant soil (bio)geochemical processes for a
specific problem. Multiple kinetic and equilibrium reaction
pathways are often a critical feature of such a reaction network(Steefel, 2000). One type of reactions in this approach proceeds
then relatively fast as compared to transport, may be pathindependent, and usually relates mostly to aqueous complexation
reactions and adsorption reactions. A second type of reactions
proceeding more slowly as compared to transport processes and toequilibrium reactions may involve oxidation/reduction, biologically mediated reactions, and mineral dissolution/precipitation
processes. Since these reactions are path-dependent, bothsequential and parallel reactions may be important (e.g., for
TNT-degradation pathways; Šimůnek et al., 2006b).

A quantitative approach based on integrating and coupling
multiple processes is also useful for pedogenesis modelling,
especially for mechanistic modelling at the mesoscopic scale(i.e., the profile, pedon or field-scale) within the hierarchy ofsoil process modelling as defined by Hoosbeek and Bryant(1992). Soil hydrology is clearly one of the major driving forcesof pedogenesis through its control on physical, biological andgeochemical processes (Lin, 2003). The term hydropedology(Lin, 2003; Lin et al., 2006) is now frequently used to explicitly
stress the interactive nature of pedologic and hydrologic
processes in the vadose zone. Geochemical processes in turn
may affect also the physical and hydraulic properties of soils,
including porosity and permeability (e.g., Lapen and Wang,
1999; Mecke et al., 2000). For saturated conditions, somemacroscopic (e.g., the Carmen–Kozeny, Brinkman or Fair–
Hartch models as discussed by Le Gallo et al., 1998) or grainscale (e.g., Panda and Lake, 1995) models are available for thispurpose. Similar models have also become available for theunsaturated hydraulic properties (Freedman et al., 2004).
Alternatively, pedotransfer functions relating basic soil physicalproperties to the unsaturated soil hydraulic functions could be
used to update the hydraulic properties (e.g., Finke, 2006).
The objective of this paper is to discuss the importance ofvariations in soil water contents and water fluxes on element
mobility as illustrated by a hypothetical numerical example
implemented in a coupled reactive transport code integrating thedifferent physical and biogeochemical processes.
2. Coupling geochemical reactions with water flow andsolute transport
2.1. Mathematical description of transport and geochemical
processes
This section defines some governing equations for the physical(transport) and the biogeochemical processes used in a mechanistic-based reactive coupled transport code. The overview islimited to the processes used in the numerical example discussed
further. A more complete overview of involved processes andtheir mathematical treatment is given e.g. in Lichtner (1996) andin the documentation of different coupled codes (e.g., 3DHYDROGEOCHEM, Yeh and Cheng, 1999; MIN3P, Mayer et al.,
2002; HP1, Jacques and Šimůnek, in preparation).
The geochemical processes used in the example are aqueous
complexation and ion exchange. The aqueous complexation
reactions are written in terms of components or master species
and their stoichiometric coefficients (Morel and Hering, 1993).
The general reaction equation for an aqueous species is:NmXmji Amj ¼ Aið1Þ

j¼1

where Nm is the number of aqueous master species, i = 1,…, Nsa,
Nsa is the number of aqueous secondary species, Ajm and Ai are
the chemical formula for the master and secondary species,
respectively, and νji are the stoichiometric coefficients in thereaction. At equilibrium (i.e., when the free energy of theD. Jacques et al. / Geoderma 145 (2008) 449–461

system is at a minimum), the mass law equations together withthe component conservation equations allow calculation of thespeciation. The mass law equations are:mljiNm 
Kil ¼ ali j am
ð2Þ
j
j¼1

where Kil is the thermodynamic equilibrium constant [−], ais the activity [−] defined as ai = γci / ci0 where c is the molal
concentration [mol kg− 1 water], ci0 is the standard state (i.e.
1 mol kg− 1 water, Appelo and Postma, 2005), γ is the activity
coefficient [−], and the superscript m and l refers to master andsecondary aqueous species. The activity coefficient is calculatedwith the extended Debye–Hückel equation (Langmuir, 1997;
Parkhurst and Appelo, 1999).
From the different models to describe sorption on the solid
soil phase (e.g., linear, Freundlich or Langmuir equilibrium
sorption isotherms, cation (or anion) exchange and surface
complexation models), the ion exchange model is used, written
in the Gaines and Thomas' (1953) convention:
Nm
X

e
m
e
mejie Amj þ mje ie Xje ¼ Aie

ð3Þ

j¼1

where je = 1,…, NX (NX is the number of master exchangers),
ie = 1,…, Nse (Nse is the number of the secondary exchange
species), and the superscript e refers to exchange reactions. Let
the activity of an exchange species be defined asaie ¼ geie bie ; je

ð4Þwhere a is the activity [−], γiee is the activity coefficient of theieth exchange species [−], and βie, je is the equivalent fraction of
the ieth exchange species on the jeth exchanger [−]. The massaction equations for equilibrium conditions can now be written asmeji mej iNm eeemKiee ¼ geie bie ; je j gm
geje bje ; je
:
j cj
j¼1

ð5Þ

The physical processes used in the example are water flow
and solute transport. The water flow equation is based onconservation of mass within a given soil volume and on thesummation of fluxes in and out of this volume. The fluxes aredescribed with the Darcy–Buckingham law relating the flux inthe unsaturated soil to the pressure head gradient in the soil anda moisture dependent proportionality factor. The resulting
Richards equation for water flow in soils is:
h

i
AhðhÞ
¼ A K ðhÞ Ah þ cos a  S ðhÞ
ð6Þ
At

Ax

Ax

where h is the soil water pressure head [L], θ is the volumetricwater content [L3L− 3], t is time [T], x is the spatial coordinate[L] (positive upward), S is the sink term [L3L− 3T− 1], α is theangle between the flow direction and the vertical axis, and K is
the unsaturated hydraulic conductivity [LT− 1]. To describe therelations between θ–h and K–h, the van Genuchten equations(van Genuchten, 1980) are used: hrhðhÞ ¼ hr þ ð1 hþs jahj
n m
Þ

ð7Þ

451

and
h

m i2
K ðhÞ ¼ Ks Sel 1  1  Se1=mð8Þ

where θr is the residual water content [L3L− 3], θs is the saturated water content [L3L− 3], α [L− 1], n [−] and m (= 1 − 1 / n)[−] are shape parameters, l is a pore connectivity parameter [−],
Ks is the saturated hydraulic conductivity [LT− 1], and Se =(θ − θr) / (θs − θr) is the effective saturation.
For equilibrium aqueous complexation reactions, it isconvenient to transport only total concentrations which aredefined as:
Ns a
X
Cj ¼ cj þ
mji ci :ð9Þ

i¼1

Assuming that all diffusion coefficients are equal (i.e.,
species-independent diffusion), the advection–dispersion equation that combines the mass conservation equation with Fick's
diffusion law can be written for each component as follows (seealso Lichtner, 1996; Mayer, 1999):
AhC¼ A hDw ACAxj  AqCj  SCr; j þ Ro; j
ð10Þ
At

Ax

Ax

where S the sink term in the water flow equation [L3L− 3T− 1],
Cr,j is the total concentration of the sink term [ML− 3], Dw is thedispersion coefficient in the liquid phase [L2T− 1], and Ro,j is thesource/sink term that represents various heterogeneous equilibrium and kinetic reactions (e.g., cation exchange, surface
complexation, mineral dissolution) and homogeneous kinetic
reactions (e.g., degradation reactions in the aqueous phase)[ML− 3T− 1].
The dispersion coefficient Dw in Eq. (10) is given by
hDw
i ¼ DL jqj þ hDi;w sw

ð11Þ

where Di,w is the molecular diffusion of the ith aqueous species
in free water [L2T− 1], DL is the longitudinal dispersivity [L], andτw is a tortuosity factor in the liquid phase [−] that is related to thewater content by the model of Millington and Quirk (1961).

2.2. HP1 simulator
All calculations discussed in this paper were carried out withversion 2.0 of HP1 (Jacques et al., 2006b). HP1 integrates a broad
range of physical and biogeochemical processes as described in
the original manuals of HYDRUS-1D (Šimůnek et al., 2005),
PHREEQC-2 (Parkhurst and Appelo, 1999), and HP1 (Jacquesand Šimůnek, in preparation). HYDRUS-1D was originallydesigned mostly for simulating water flow, solute transport andheat transport in soils from the soil column scale, to the lysimeterscale, and the field scale. Important to note is that the modelaccounts for soil heterogeneity at the macroscopic scale bydefining distinct soil horizons, at the small-scale by using scaling
factors (Vogel et al., 1991), and at the microscopic scale by using
the dual porosity or mobile/immobile flow models (Šimůnek

452

D. Jacques et al. / Geoderma 145 (2008) 449–461

Table 1
Soil hydraulic properties of the six soil horizons of a sandy dry Spodosol(Seuntjens, 2000)
Horizon Layer thickness θr
(cm)

θs

n
α
(cm− 1)

A
E
Bh1
Bh2
Bh/C
C

0.476
0.416
0.472
0.455
0.464
0.408

0.016
0.015
0.016
0.028
0.023
0.039

13
10
8
8
16
45

0.065
0.035
0.042
0.044
0.039
0.029

Ks
l
(cm day− 1)

2.411
93
2.446 311
2.197
392.010 860
2.534 1198
2.071 1212

0.5
0.5
0.5
0.5
0.5
0.5

et al., 2003). PHREEQC-2 stems mostly from the field of lowtemperature aquatic chemistry and geochemistry. The genericnature of PHREEQC-2 and its flexibility in defining the kinetic
reactions with respect to the reaction network, influencing factors
and mathematical formulation provide many opportunities todefine and include pedogenetic processes. The combination of theadvection–dispersion equation with the generic reaction term Ro, j
in Eq. (10) allows also the description of the movement of colloids(e.g. DeNovio et al., 2004; Šimůnek et al., 2006a). Version 2.0
couples version 3.0 of HYDRUS-1D to solve the water flow, heattransport and solute transport equations, and version 2.13 of
PHREEQC to calculate Ro,j in Eq. (10), aqueous complexationand the composition of the ion exchange complex.
HP1 implements a weak coupling method in which thegoverning equations for water flow, heat transport and solutetransport are solved sequentially. The multicomponent reactive
transport equations are solved using a non-iterative sequential
approach, meaning that the physical part is solved first without
any chemical interactions, while the chemical reactions that are
uncoupled in space and coupled over the components are solved
subsequently. A flow chart of the operator-splitting approach usedin HP1 was given in Jacques et al. (2006a). Although operatorsplitting errors may occur with our approach (e.g., Carrayrou
et al., 2004), Jacques et al. (2006a) showed that very accurate
results can be obtained by carefully selecting the spatial andtemporal discretizations for the numerical solution. We refer to
Jacques et al. (2006a) for specific discretization guidelines. HP1 isfreely available at www.pc-progress.cz.
3. Transient flow and transport of major cations and heavy
metals in a soil profile
As an illustrative example, HP1 is used in this section tosimulate the leaching of Cd and Zn in a dry Spodosol in a sandy
region of Northern Belgium using in-situ measured Cd and Zn(and additional elements) concentration profiles. Soils in theregion were contaminated by atmospheric deposition of Cd andZn from non-ferrous industry (Seuntjens, 2000). Water flow andthe transport of major cations (Na, K, Mg, and Ca), heavy
metals (Cd and Zn), anions (Cl, Br) and Al were simulated for a
1-m deep multi-layered soil profile subject to atmosphericboundary conditions for a period of 30 years. Our main focus
here is on how processes affecting water contents and water
fluxes also influence the geochemical conditions in the soil.
Specifically, the effect of cycles of evaporation and infiltration
on pH and Cd speciation will be discussed in detail.
3.1. Problem definition and model input
The experimental site “Kattenbos” is located in thecommunity of Lommel (Belgium) at a distance of about 2 km
east from a metallurgical plant. Various soil parameters and theinitial conditions were taken from Seuntjens (2000). The dry
Spodosol consists of six horizons, with a typical leached Ehorizon between 13–23 cm below the soil surface, and enriched(organic matter and Fe-oxides) Bh-horizons between 23–55 cm
below the surface. Table 1 gives depths, parameters of the water
retention characteristics, θ(h), and the saturated hydraulicconductivity, Ks, (both measured in the laboratory on coresamples of 100 cm3) for each soil layer.
Interactions between major cations and heavy metals withthe soil solid phase were simulated by means of cation exchange
processes assuming local equilibrium on a single type of exchange sites. Voegelin (2001) has shown that this approachadequately describes various features of Cd transport experiments. Although sorption on specific sites with a high affinity
for Cd may also occur in soils (Selim et al., 1992), this typeof binding is unlikely in acid sandy soils (Smolders et al.,
1999; Voegelin et al., 2001). Exchange parameters were
calibrated using concentrations measured in drainage waterfrom a steady-state flow experiment on large, undisturbed soil
lysimeters (1-m-long and 0.8-m-diameter, Seuntjens et al.,2001). The initial composition of the cation exchange site wasmeasured for each soil horizon (Table 2). Resulting exchange
parameters are given in Table 3. While the same log(K)parameters were used for all soil horizons, the size of the cationexchange complex was assumed to vary between horizons. We
note that Al in the soil profile was assumed to be in equilibriumwith Gibbsite (Al(OH)3). This resulted in a good description of
the measured pH and Al in the drainage water (results notshown).
The initial aqueous concentrations were obtained by equilibrating the soil solution with the cation exchange composition.

Table 2
Initial composition of the cation exchange sites (in mol/1000 cm3 soil) (Seuntjens, 2000)
Horizon
A
E
Bh1
Bh2
Bh/C
C

HX

NaX
−2

1.32 108.76 10− 3
1.78 10− 2
1.29 10− 2
3.61 10− 3
4.80 10− 3KX−42.73 102.90 10− 4
5.36 10− 4
2.41 10− 44.42 10− 4
9.05 10− 5

MgX2

CaX2
−4

2.99 102.03 10− 44.89 10− 42.81 10− 4
1.20 10− 42.01 10− 4−42.77 10
1.57 10− 4
5.00 10− 4
1.86 10− 4
1.24 10− 43.85 10− 5

CdX2
−4

1.44 10
7.26 10− 42.38 10− 4
1.24 10− 45.80 10− 56.34 10− 5

ZnX2
−5

1.52 101.70 10− 54.51 10− 5
1.84 10− 5
1.38 10− 5
1.25 10− 5

2.00 10− 4
8.20 10− 5
3.64 10− 42.45 10− 4
4.64 10− 5
1.84 10− 5

D. Jacques et al. / Geoderma 145 (2008) 449–461

453

Table 3Log K parameters for the cation exchange reactions
FormulaLog K

X− + H+ = HX
X− + Na+ = NaX
X− + K+ = KX
X− + Ca2+ = CaX2
X− + Mg2+ = MgX2
X− + Cd2+ = CdX2
X− + Zn2+ = ZnX2

4.2
1.8
2.85
2.8
3.1
4.0
2.85

Element concentrations in the rain water were obtained from
Stolk (2001) for station 231 located in Gilze-Rijen (theNetherlands) close to the site investigated here. The composition was based on the average of 13 measurements during 1999
(Table 4). The cation exchange sites were initially not in
equilibrium with the composition of the rain water. Consequently, long-term changes in the composition of the cationexchange sites were simulated for a period of 30 year, not only
for Cd and Zn (due to infiltration of heavy metal free water), butalso for the other major cations.
Atmospheric boundary conditions were taken for the soil
surface (see further), whereas a free drainage (a zero pressure
head gradient) was taken as the bottom boundary condition.
3.1.1. Precipitation and potential evaporation
Climatological data were taken from two nearby weather
stations (from Brogel for years 1969 through 1985, and from
Mol/Geel for years 1986 through 1998). The potential
evaporation, Ep [L T− 1], was calculated using the Penman
equation (Penmann, 1948).
The total precipitation during the 30-year simulation period was
24,933 mm with a mean of 831 mm/year and a standard deviationof 152 mm/year. The maximum and minimum yearly precipitation
was, respectively, 494 mm (1976) and 1169 mm (1998). The meanand standard deviation of the potential evaporation was 477 mm/
year and 91 mm/year, respectively, with 329 mm (1979) and
616 mm (1990) as minimum and maximum values, respectively.
This results in a mean yearly potential precipitation surplus Pp⁎,
defined as the difference between the precipitation P and thepotential evaporation Ep, of 353 mm/year, varying between 32(1976) and 621 (1988) mm. Fig. 1 shows the cumulative potential
precipitation surplus flux between 1972 and 1981. During thesummer months of 1975 and 1976, Pp⁎ was close to zero. An
overview of the monthly positive and negative Pp⁎ is given in
Fig. 2. On a monthly basis, at least 25% of months between April
and August had a negative potential precipitation surplus. The
actual evaporation rate, Ea [LT− 1], was lower than the potential
rate (10,473 mm versus 14,335 mm), resulting in a higher actual

Fig. 1. Cumulative potential and actual precipitation surplus for the 30-year
period (a) and a detail for the period 1972–1981 (b).

precipitation surplus (Fig. 1). The large difference is due to theabsence of plant roots in this example. The sandy soil with its low
hydraulic conductivity when dry was unable to provide enoughwater to satisfy the atmospheric demand, leading to lower Ea
fluxes.
3.1.2. Spatial and temporal discretization
Spatial and temporal discretizations leading to minimal
operator-splitting errors were determined using the guidelines

Table 4
Chemical composition of the rain water (from Stolk, 2001) (*: charge balance)
Concentration (µmol/l)

pH

Cl

Br⁎

Na

K

Ca

Mg

Cd

Zn

Rain water

5.25

69

32

64

4

6

8

0

0

Fig. 2. Statistics of the monthly potential precipitation surplus for the period
1969–1998.

454

D. Jacques et al. / Geoderma 145 (2008) 449–461

proposed by Jacques et al. (2006a). Simulations were firstcarried out for a shorter 3-year time period (1975–1977) involving alternating wet and dry periods, with at least one very
dry period having a relatively high potential evapotranspiration
rate, ETp, to optimise the maximum time step and the performance index ωs obtained by multiplying the Peclet and Courant
numbers. We found that smaller maximum time steps were mosteffective for reducing operator-splitting errors during periods
with larger evapotranspiration rates, whereas decreasing theperformance index was most effective during water infiltration
periods. A simulation with Δx = 0.5 cm, ωs = 0.125, and amaximum allowed time step of 0.125 day was used as thereference run. Nine other simulations were carried out with aspatial discretization of 1 cm and different maximum allowed
time steps (0.5, 0.25 and 0.125 day) and ωs values (0.4, 0.2
and 0.1). Since it is likely that the reference run may also
produce larger errors during relatively dry periods, the accuracyof different simulation runs was assessed by comparing
simulation results for periods with high precipitation rates andlow evaporation rates (i.e., during winter) to optimize ωs.
Evaluated over all chemical components, the lowest tested ωsvalue was selected (with a maximum time step of 0.125 day). To
optimize the maximum allowed time step, the maximum
difference between the reference run and simulations withdifferent maximum time steps was determined. Our current
example required a maximum time step of 0.125 day. For asimilar problem over a 10-year simulation period (1972–1982),
the relative error in the solute mass balance (using Δx = 1 cm,
ωs = 0.1, and a maximum time step of 0.125 day) was calculatedusing Eq. (8.55) of Šimůnek et al. (2005) for consecutive timeintervals of 15 days. The errors were consistently less than 1%,
which indicates that the HP1 mass balance errors were relatively
small for our example (see Jacques and Šimůnek, in
preparation).

Fig. 4. Water fluxes in the soil during year 1975 (12/27/1974–01/09/1976).

into three periods: period I (12/27/1974–04/20/1975) withmostly positive values of Pp⁎, period II (04/20/1975–09/20/
1975) with mostly negative values of Pp⁎, and period III (09/20/
1975–01/09/1976)) with mostly positive values of Pp⁎. Thisdivision was made to relate accumulation and leaching of thedifferent elements to time periods with either positive or
negative potential precipitation surpluses.

3.2. Simulation results
Detailed simulation results for the year 1975 are presented in
Figs. 3–6. On an annual basis, 1975 had the sixth lowest Pp⁎.
The year 1975 was chosen because it has relatively distinct
periods with mainly positive and negative potential precipitation surpluses, as shown in Fig. 3. This year could be divided

Fig. 3. Potential precipitation surplus of year 1975 (12/27/1974–01/09/1976).

3.2.1. Water fluxes during 1975
Water fluxes in the soil profile during the period 12/7/1974 to
01/09/1976 are shown in Fig. 4. As expected, the highest
variability in the water flux occurred in the top soil, ranging
from a 0.8 cm day− 1 downward flux to a 0.047 cm day− 1
upward flux. During evaporation periods, the zone with upward
flow propagated down to a depth of 50 cm in the soil profile.
3.2.2. Cl, Na, Ca, and Cd concentrations at different depths
during 1975
To circumvent the effect of changing water contents on theconcentrations, we expressed concentrations in terms of mass
per soil volume (mol/1000 cm3 of soil) instead of liquid
concentrations (mol/l). Fig. 5 shows the change (compared to thefirst date of 03/21/1975) in the total mass of Cl, Na, Ca and Cd(TX in mol/1000 cm3) in three soil compartments (0–9 cm, 20–
29 cm, and 40–49 cm) for a limited time period (03/21–11/28)
during 1975. For Na, Ca and Cd, TX represents the sum of theaqueous and adsorbed amounts. The cumulative inflow of Cl,
Na, and Ca from the rain water in the top compartment is alsoplotted.

D. Jacques et al. / Geoderma 145 (2008) 449–461

The results in Fig. 5 indicate a strong correlation between
the accumulation and leaching of Cl in the first compartment
with the precipitation surplus. Periods with a negative Pp⁎
caused an upward Cl flux and an increase in TCl in the firstcompartment. The contribution of Cl from rain water when Pp⁎
is negative was small. During rainstorms (positive Pp⁎), the Cl
amount first increased, but then decreased since Cl issubsequently leached from the first compartment. Althoughdeeper in the soil (e.g., 20–29 cm), where Cl fluxes are still
upwards during periods with negative Pp⁎ (not shown), TCl isdecreasing. Upward Cl fluxes are found down to a depth of
50 cm for the year 1975.
The total amount of Na, TNa, exhibits a similar behaviour.
However, Na is less mobile than Cl due to the exchange
processes. The upward Na fluxes and thus the Na accumulationin the first compartment are slightly smaller than those for Cl.
More importantly, leaching during positive Pp⁎ periods was less,455

resulting in a larger net increase in Na than in Cl during period I
in 1975. A slightly different pattern was simulated for Ca, a lessmobile element than Na. During period I, TCa increased in theupper compartment, whereas TNa (and TCl) were relatively
constant. During period II, changes in TCa when Pp⁎ was
positive were comparable or larger than when Pp⁎ was negative.
The opposite was simulated for Na. While the difference in Na
and Ca mobility during period I can be explained by theincreased retardation of Ca (and thus its lower leaching rate), thedifferences cannot be similarly explained during the second
period. This phenomenon is related to changing concentrationsof different cations due to variations in water contents andfluxes and the competition of cations with different affinities forthe exchange complex (further discussed below). Cd behaved
similarly as Ca, but without the accumulation in the firstcompartment for positive Pp⁎ values since the rain water did notcontain any Cd.

Fig. 5. Changes in amount of (a) Cl, (b) Na, (c) Ca, and (d) Cd (mol/1000 cm3 soil) in the soil between 0–9 cm (full line, circles), 20–29 cm (short dashes, squares) and
40–49 cm (long dashes, diamonds). For Cl, Na, and Ca the cumulative inflow in the rain water is shown as a full line (on the right axis, note the different scale). The
precipitation flux is shown on the top.

456

D. Jacques et al. / Geoderma 145 (2008) 449–461

Fig. 5 (continued ).

3.2.3. Distributions of Na, Ca, and Cd versus depth during
1975
Fig. 6 shows distributions versus depth of total Cl, Na, Ca
and Cd in the aqueous phase (TX,aq, X = Cl, Na, Ca, Cd). The
distributions of Cl illustrate the findings from the mass balance
analysis. During the first half of the year (from winter tosummer), chloride accumulated in the top soil due to limited
downward water flow (e.g., during period II for 1975) and only
slowly leached into the bottom half of the soil profile. TNa,aq
increased up to 2.5 times in the beginning of period II in
1975, but then remained more or less constant. On the other
hand, TCa,aq in the top soil increased more continuously up to 5
times at the end of period II compared to the beginning of
period II. TNa,cec increased during the summer months, whereas
TCa,cec remained relatively constant (results not shown; see
Fig. 7 for long-term variation). TCd,aq showed a strong correlation with the chloride dynamics (in contrast with Ca): on
11/28/1975, TCd,aq in the subsoil was as high as the highest
values in the top soil during the summer, whereas this was not

the case for Ca. During the summer months the correlation with
Cl was less clear in the top 15 cm: while TCl,aq increased
towards the soil surface, TCd,aq remained more or less constant.
TCd,cec did not change over one year (results not shown; see
Fig. 7). Since TCd,cec was much larger than TCd,aq (1.0 × 10− 5–
4.0 × 10− 5 and 1.0 × 10− 9 –1.0 × 10− 8 mol Cd, respectively), thetotal Cd content remained also constant during one year.
3.2.4. Long term changes in total mass (1969–1998)
The variability in the mass of a component in the aqueous
phase within a single year was found to be much larger than thelong-term variability over 30 years. A slight increase or decrease
occurred over the 30 year period, corresponding to an increase or
decrease in the adsorbed amounts. The total mass of Na, Ca, and
Cd in the different horizons over the period of 30 year is plottedin Fig. 7. For Na, variations in the top horizon (horizon A) within
a year were similar to the long-term variations, while for theother horizons the 30-year variability was larger than thevariations within one year. However, much of the variation over

D. Jacques et al. / Geoderma 145 (2008) 449–461

457

Fig. 6. Distributions versus depth of TCl,aq, TNa,aq, TCa,aq and TCd,aq (mol/1000 cm3 soil) at 10 times during year 1975 (12/27/1974–01/09/1976).

the 30 year period was due to disequilibrum between thecomposition of the rainwater and the initial composition of thecation exchange sites. It took almost 10 years to equilibrate eachhorizon with the inflowing water. After that, no long-term trend
was present anymore, and the within-year variability also
increased for the deeper horizons.
Almost no changes were found within one year in the total
amounts of Ca and Cd. Long-term changes occurred for Ca
during the 30-year simulation period. We found substantial
accumulation of Ca (and Mg, results not shown) in the top
horizon, accompanied by a decrease in proton-occupied cationexchange sites (results not shown). The change in TCd in the
Ap-horizon was small. Net leaching of Cd occurred in the E andBh1 horizon, resulting in Cd accumulation in the Bh2 horizon.
Cd leached also from the B/C horizon, while accumulating in
the upper part of the C-horizon.

3.3. Discussion
Fig. 8 shows time series between 1972 and 1982 for thewater content, Cl, pH and Cd in the aqueous phase at two
depths. The results illustrate the relation between water flow andgeochemical conditions in the soil. The alternation between
precipitation (wet conditions) and evaporation (dry conditions)as dictated by the atmospheric conditions clearly affected the
dynamics of the water content, with upward water flow during
dry periods. The flow dynamics in turn significantly influenced
the geochemistry near the soil surface. As illustrated, the mostmobile elements (anions such as Cl− and monovalent cations
such as Na+) move upwards during the evaporation periods,
thus causing these elements to accumulate near the soil surface.
The decrease in the water content near the soil surface due toevaporation produced higher concentrations and a lower pH.

458

D. Jacques et al. / Geoderma 145 (2008) 449–461

water contents, the concentration of all aqueous species increased during the summer period. Changes in aqueousconcentrations in turn caused changes in the cation exchange
equilibrium, thereby promoting monovalent cations to sorb on
the cation exchange complex and bivalent cations to desorb into
solution. This explains also the difference between TNa,aq andTCa,aq during the summer near the soil surface. The aqueousconcentration of Na was more controlled by the cation exchange
complex than that of Ca due to preferred adsorption of Na
during dry soil conditions. This (geochemical) process is further
amplified by the increased supply of monovalent cations due toupward water flow during summer, leading to relatively
increased sorption of the monovalent cations and higher
concentrations of divalent cations and heavy metals in the soil
solution. The complexation of Cd with Cl contributed also to
the increase in TCd,aq. At a depth of 3 cm, about 0.65 and up to
3.5% of the aqueous Cd was in the form of CdCl+ during winterand summer, respectively. Other Cd complexes of the form
CdCln(2 − n) were present in much smaller concentrations. Similar
to the pH changes, changing water contents and upward fluxes
both affected the amount of Cd in the aqueous phase.

Fig. 7. TNa,cec, TCa,cec and TCd,cec (mol/1000 cm3 soil) in each horizon during the
30-year simulation (the C horizon is divided in two parts: C1 — 75 cm, andC2 — 75–100 cm).

PHREEQC calculations showed that removing only water from
aqueous solutions results in a pH decrease (for a similargeochemical system, i.e., only cation exchange and equilibriumwith gibbsite). Another factor intensifying the pH decrease was
the greater mobility of anions (Cl− and Br−) compared tocations. Thus, the physical factors of decreasing water content
and upward water and solute fluxes caused the pH to decrease
when Pp⁎ is negative. It should be noted that the pH in reality isalso affected by other geochemical or biological processes notincluded in our conceptual model. For example, soil carbon
dioxide concentrations that usually change in response tobiological activity and moisture status of the soil can also affectsoil pH (Šimůnek and Suarez, 1993).
Although upward flow during the summer had almost no
effect on the total amount of heavy metals in the surface horizon
due to the low mobility of these elements, the aqueous concentrations of the metals did vary significantly during theseason. Several factors contribute to this. First, because of lower

Fig. 8. Time series of water content, Tcl,aq, pH and TCd,aq (mol/1000 cm3 soil)
between 1972 and 1982 at 5 cm (A-horizon) and 35 cm (Bh1-horizon).

D. Jacques et al. / Geoderma 145 (2008) 449–461

Our results above indicate that atmospheric boundary
conditions can have a significant effect on the amount andtransport of Cd in a soil profile, and on its bioavailability since
uptake processes by plants and soil micro-organisms are often
concentration-dependent. Passive root uptake of solutes together
with water uptake increases with increasing solute concentrations. Similarly, active uptake as described with Monod orMichealis–Menten kinetics will increase with increasing soluteconcentrations. Moreover, the high heavy metal concentrations
occurred during the summer months with the highest (micro)
biological activity. In addition, Cd speciation may also play arole in uptake. For example, Smolders and McLaughlin (1996)
observed more Cd uptake by chard when Cl concentrations
increased while the Cd2+ activity was kept constant, likely due tophytoavailability of CdCl+ and other CdCln2 − n-species.
4. Concluding remarks
The example in this paper on heavy metal mobility shows that
temporal climatic variations in precipitation and evaporation canlead to significant changes in the soil geochemical conditions,especially near the soil surface. The temporal variability in soil
physical variables (water contents and fluxes as a result ofchanging atmospheric conditions) causes temporal variability inthe geochemical variables and processes (e.g., pH, concentrations, exchange equilibria). In our example the pH decreased atlower water contents. Due to increasing concentrations at lower
water contents, monovalent cations are preferred on the cationexchange complex relative to bivalent cations. Upward water
fluxes during the summer, moreover, causes transport of chloride
to the soil surface, which in turn leads to higher Cd concentrationsin the aqueous phase due to Cd–Cl complex formation.
The example illustrates that simulators as HP1 are potentiallyattractive tools for studying reactive transport processes in thevadose zone during transient variably-saturated flow. Geochemical conditions are an important factor since they
determine the speciation (both in the soil solution and on thesolid phase) of the elements involved, and thus their mobility
and bioavailability. Small variations in prevailing geochemical
conditions may alter significantly the speciation and mobility ofheavy metals or other constituents. Soil systems, moreover, aresubject to large transient variations since they are open to theatmosphere. Hence, changes in the composition of rain water or
atmospheric deposition, among other external factors, canmaterially alter the geochemical conditions in a soil profile.
Reactive transport codes for soil systems can therefore behelpful for unravelling processes related to pedogenesis and tomodel and describe pedogenesis. Although the numericalexample discussed in this paper is not directly related tomodelling pedogenesis, the simulations do illustrate the interactions that may occur between transient water flow and a variety ofgeochemical reactions, including mass transfer processes between
the liquid and solid phases. The geochemical reactions together
with variably-saturated water flow and drying and wetting cycles
may explain many specific processes occurring during pedogenesis when extrapolated or evaluated over longer time periods. It
was beside the scope of the present paper to analyze the effect of459

short-term fluctuations on the long-term behaviour of contaminants in the soil system. The studies of Jacques et al. (2006b,
in press) revealed that there is a difference in U-leaching fromsoils between simulations using daily atmospheric values or longterm averaged ones. The breakthrough of U at 1 m depth occurredearlier in case of daily values.
The use of integrated reactive transport codes come with stillother challenges related to (i) computational demand, (ii) updating
physical properties, (iii) detail of process knowledge and model
input, and (iv) model validation.
Because of relatively high computational demand, coupledcodes appear mainly useful for short-time (10–1000 years)
simulations. Still, the model seems suitable for evaluating soil
genetic processes and genesis resulting from human impacts,
such as the effect of acid rain on forest soils or the developmentof technosoils.
A further limitation of most current codes is that soil horizons
and their properties (e.g., porosity, permeability and cationexchange properties) are considered to be static during thesimulations, which is the contrary to what is needed in modelling
pedogenesis. Further research should therefore address timedependency in flow and transport properties in the coupled codes.
Also, these codes combine relatively specialized informationfrom the field of soil physics and geochemistry. Their users thus
need a relatively thorough understanding of the basic physicaland biogeochemical processes involved. The integrated models
are also very parameter demanding, with many parameters
being subject to uncertainty, including especially those invarious chemical rate equations. For example, it is notstraightforward under saturated conditions to relate the reactive
surface area of a mineral to, for example, BET-measured surface
areas (Brantley, 2003). In addition, only a small fraction of thereactive surface area will be chemically active (i.e., in contact
with water) during unsaturated conditions (Kuechler andNoack, 2007). One way to address these uncertainties is byperforming uncertainty and sensitivity analysis to first identify
those parameters which most impact system response, thereby
prioritizing parameter estimation efforts.
For the time being, only a limited number of test cases andexperimental data exist for the unsaturated zone to test the fullcapabilities of a reactive transport code (see also Davis et al.,2004). One elaborate recent application is a study by Gonçalves
et al. (2006) who successfully applied the chemistry-specificmajor ion geochemistry module in HYDRUS-1D (version 3.0,Šimůnek et al., 2005) to a 4-year experimental data set
involving water flow and solute transport in lysimeters irrigatedwith waters of different quality and subjected to atmosphericconditions. Additional studies of this type should provide morecredibility to the use of coupled models for addressing flow andreactive transport problems in the vadose zone.
Acknowledgements
This work was partly supported by the National Science
Foundation, agreement no. EAR-9876800 and the Terrestrial
Sciences Program of the Army Research Office (Terrestrial
Processes and Landscape Dynamics and Terrestrial System

460

D. Jacques et al. / Geoderma 145 (2008) 449–461

Modeling and Model Integration). Additional support was
obtained through the bilateral agreement project “Development
and evaluation of a coupled geochemical transport model
between SCK•CEN and USDA-ARS (agreement no. 58-53100-F105), between SCK•CEN and the University of California,
Riverside (agreement no. C0-90001412.01), and by the
SCK•CEN R&D program on “Environmental Remediation”(project E022031 and CO91002).
References
Appelo, C.A.J., Postma, D., 2005. Geochemistry, Groundwater and Pollution,
Second edition. A.A. Balkema, Rotterdam, the Netherlands.
Arkley, R.J., 1963. Calculation of carbonate and water movement in soil fromclimatic data. Soil Sci. 96, 239–248.
Brantley, S., 2003. Reaction kinetics of primary rock-forming minerals underambient conditions. In: J.I. Drever (Editor), Fresh Water Geochemistry,
Weathering and Soils, K.K. Turekian and H.D. Holland (Editors), Vol. 5 ofTreatise on Geochemistry. Pergamon Press, Oxford, UK, pp. 73–118.
Carrayrou, J., Mosé, R., Behra, P., 2004. Operator-splitting procedures forreactive transport and comparison of mass balance errors. J. Contam.
Hydrol. 68, 239–268.
Davis, J.A., Yabusaki, S.B., Steefel, C.I., Zachara, J.M., Curtis, G.P., Redden, G.B.,
Criscenti, L.J., and Honeyman, B.D., 2004. Assessing conceptual models for
subsurface reactive transport of inorganic contaminants. EOS, 85: 449 and 455.
DeNovio, N.M., Saiers, J.E., Ryan, J.N., 2004. Colloid movement in unsaturated
porous media: recent advances and future directions. Vadose Zone J. 3, 338–351.
Finke, P., 2006. Modelling the genesis of Luvisols in late Weichsel loess. In:
Samouëlian, A., Cornu, S. (Eds.), Proceedings of the Workshop on
Modelling of Pedogenesis, pp. 29–30.
Freedman, V.L., Bacon, D.H., Saripalli, K.P., Meyer, P.D., 2004. A film
depositional model of permeability for mineral reactions in unsaturated
media. Vadose Zone J. 3, 1414–1424.
Gaines, G.L., Thomas, H.C., 1953. Adsorption studies on clay minerals. II. A
formulation of the thermodynamics of exchange adsorption. J. Chem. Phys.
21, 714–718.
Gonçalves, M.C., Šimůnek, J., Ramos, T.B., Martins, J.C., Neves, M.J., Pires, F.P.,
2006. Multicomponent solute transport in soil lysimeters irrigated withwaters of different quality. Water Resour. Res. 42, W08401. doi:10.1029/
2005WR004802.
Hoosbeek, M.R., Bryant, R.B., 1992. Towards the quantitative modeling of
pedogenesis — a review. Geoderma 55, 183–210.
Jacques, D., Šimůnek, J., (in preparation). User manual of the multicomponent
variably-saturated transport model HP1 (Version 2.0): description, verification and examples. SCK•CEN, Mol, Belgium.
Jacques, D., Šimůnek, J., Mallants, D., van Genuchten, M.Th., 2006a. Operatorsplitting errors in coupled reactive transport codes for transient variably
saturated flow and contaminant transport in layered soil profiles. J. Contam.
Hydrol. 88, 197–218.
Jacques, D., Šimůnek, J., Mallants, D., van Genuchten, M.Th., 2006b.
Modelling uranium leaching from agricultural soils to groundwater as acriterion for comparison with complementary safety indicators. In: Van
Iseghem, P. (Ed.), Scientific Basis for Nuclear Waste Management XXIX.
Mat. Res. Soc. Symp. Proc., vol. 932, pp. 1057–1064.
Jacques, D., Šimůnek, J., Mallants, D., van Genuchten, M.Th., (in press).
Modeling coupled hydrological and chemical processes in the vadose zone:
a case study of long-term uranium transport following mineral Pfertilization. Vadose Zone J.
Jansen, B., Nierop, K.G.J., Verstraten, J.M., 2002. Influence of pH and metal/
carbon ratios on soluble organic complexation of Fe(II), Fe(III) and Al in
soil solutions determined by diffusive gradients in thin films. Anal. Chim.
Acta 454, 259–270.
Jenny, H., 1941. Factors of Soil Formation. McGrawhill, New York, NY, USA.
Kuechler, R., Noack, K., 2007. Comparison of the solution behaviour of apyrite-calcite mixture in batch and unsaturated sand column. J. Contam.
Hydrol. 90, 203–220.

Langmuir, D., 1997. Aqueous Environmental Geochemistry. Prentice Hall, New
Jersey, USA.
Lapen, D.R., Wang, C., 1999. Placic and Ortstein horizon genesis and peatland
development, Southeastern Newfoundland. Soil Sci. Soc. Am. J. 63, 1472–1482.
Le Gallo, Y., Bildstein, O., Brosse, E., 1998. Coupled reaction-flow modellingof diagenetic changes in reservoir permeability, porosity and mineralcompositions. J. Hydrol. 209, 366–388.
Lichtner, P.C., 1996. Continuum formulation of multicomponent-multiphasereactive transport. In: Lichtner, P.C., Steefel, C.I., Oelkers, E.H. (Eds.),
Reactive Transport in Porous Media, Reviews in Mineralogy, no. 34.
Mineralogical Society of America, Washington, DC, USA, pp. 1–81.
Lichtner, P.C., Yabusaki, S., Pruess, K., Steefel, C.I., 2004. Role of competitive
cation exchange on chromatographic displacement of cesium in the vadosezone beneath the Hanford S/SX tank farm. Vadose Zone J. 3, 203–219.
Lin, H., 2003. Hydropedology: bridging disciplines, scales and data. Vadose
Zone J. 2, 1–11.
Lin, H., Bouma, J., Pachepsky, Y., Western, A., Thompson, J., van Genuchten,
M.Th., Vogel, H.-J., Lilly, A., 2006. Hydropedology: synergistic integrationof pedology and hydrology. Water Resour. Res. 42, W05301. doi:10.1029/
2005WR004085.
Mayer, U., 1999. A numerical model for multicomponent reactive transport invariably saturated porous media. Ph.D. thesis, Department of Earth
Sciences, University of Waterloo, Waterloo, Canada.
Mayer, K.U., Frind, E.O., Blowes, D.W., 2002. Multicomponent reactive
transport modelling in variably saturated porous media using a generalized
formulation for kinetically controlled reactions. Water Resour. Res. 38,
1174. doi:10.1029/2001WR000682.
Mecke, M., Westman, C.J., Ilvesniemi, H., 2000. Prediction of near-saturated
hydraulic conductivity in three Podzolic Boreal forest soils. Soil Sci. Soc.
Am. J. 64, 485–492.
Metz, V., Pfingsten, W., Lützenkirchen, J., Schüssler, W., 2002. TrePro 2002:
modelling of coupled transport reaction processes. Workshop of the
Forschungszentrum Karlsruhe, 20–21 March 2002, Karlsruhe, Germany
126 p. http://www.fzk.de/fzk/groups/ine/documents/internetdokument/
id_040291.pdf2002.
Millington, R.J., Quirk, J.M., 1961. Permeability of porous solids. Trans.
Faraday Soc. 57, 1200–1207002E.
Morel, F., Hering, J., 1993. Principles and Applications of Aquatic Chemistry.
John Wiley & Sons Inc., New York, USA.
NRC, 2006. Proceedings of the international conference on conceptual model
development for subsurface reactive transport modelling of inorganic
contaminants, radionuclides and nutrients, NUREG/CP-0193, Office of
Nuclear Regulatory Research, U.S. Nuclear Regulatory Commission,
Washington, DC, USA. http://www.iscmem.org/Documents/Proceedings_
03.pdf2006.
Panda, M.N., Lake, L.W., 1995. A physical model of cementation and its effects
on single-phase permeability. AAPG Bull. 79, 431–443.
Parkhurst, D.L., Appelo, C.A.J., 1999. User's guide to PHREEQC (Version 2) —a computer program for speciation, batch-reaction, one-dimensionaltransport, and inverse geochemical calculations. Water-Resources Investigations, Report 99-4259, Denver, Co, USA.
Penmann, H.L., 1948. Natural evaporation from open water, bare soils andgrass. Proc. Royal Soc. A. 190, 120–145.
Presley, D.R., Ransom, M.D., Kluitenberg, G.J., Finnell, P.R., 2004. Effects ofthirty years of irrigation on the genesis and morphology of two semi-aridsoils in Kansas. Soil Sci. Soc. Am. J. 68, 1916–1926.
Rasmussen, C., Southard, R.J., Horvath, W.R., 2005. Modeling energy inputs topredict pedogenic environments using regional environmental databases.
Soil Sci. Soc. Am. J. 69, 1266–1274.
Selim, H.M., Büchter, B., Hinz, C., Ma, L., 1992. Modeling the transport andretention of cadmium in soils: multireaction and multicomponent
approaches. Soil Sci. Soc. Am. J. 56, 1004–1015.
Seuntjens, P., 2000. Reactive Solute Transport in Heterogeneous Porous
Medium: Cadmium Leaching in Acid Sandy Soils. Ph.D., University of
Antwerp, Belgium.
Seuntjens, P., Mallants, D., Toride, N., Cornelis, C., Geuzens, P., 2001. Gridlysimeter study of steady state chloride transport in two Spodosol types
using TDR and wick samplers. J. Contam. Hydrol. 51, 13–39.

D. Jacques et al. / Geoderma 145 (2008) 449–461
Šimůnek, J., Suarez, D.L., 1993. Modeling of carbon dioxide transport andproduction in soil: 1. Model development. Water Resour. Res. 29, 487–497.
Šimůnek, J., Jarvis, N.J., van Genuchten, M.Th., Gärdenäs, A., 2003. Review
and comparison of models for describing non-equilibrium and preferential
flow and transport in the vadose zone. J. Hydrol. 272, 14–35.
Šimůnek, J., van Genuchten, M.Th., Šejna, M., 2005. The HYDRUS-1D
Software Package for Simulating the One-Dimensional Movement of Water,
Heat, and Multiple Solutes in Variably-Saturated Media. Version 3.0,
HYDRUS Software Series 1. Department of Environmental Sciences,
University of California Riverside, Riverside, California, USA.
Šimůnek, J., He, C., Pang, L., Bradford, S.A., 2006a. Colloid-facilitated solutetransport in variably saturated porous media: numerical model andexperimental verification. Vadose Zone J. 5, 1035–1047.
Šimůnek, J., Jacques, D., van Genuchten, M.Th., Mallants, D., 2006b.
Multicomponent geochemical transport modelling using HYDRUS-1Dand HP1. J. Am. Water Resour. Assoc. 46, 1537–1547.
Smolders, E., McLaughlin, M.J., 1996. Chloride increases cadmium uptake in
Swiss chard in a resin-buffered nutrient solution. Soil Sci. Soc. Am. J. 60,
1443–1447.
Smolders, E., Brans, K., Földi, A., Merckx, R., 1999. Cadmium fixation in soils
measured by isotopic dilution. Soil Sci. Soc. Am. J. 63, 78–85.

461

Steefel, C.I., 2000. New directions in hydrogeochemical transport modelling:
incorporating multiple kinetic and equilibrium reaction pathways. In:
Bentley, L.R., Sykes, J.F., Brebbia, C.A., Gray, W.G., Pinder, G.F. (Eds.),
Computational Methods in Water Resources XIII. A.A. Balkema,
Rotterdam, the Netherlands, pp. 331–338.
Stolk, A.P., 2001. Landelijk Meetnet Regenwatersamenstelling, Meetresultaten
1999. RIVM, The Netherlands. RIVM Rapport 723101 056.
van Genuchten, M.Th., 1980. A closed-form equation for predicting thehydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44,
892–898.
Voegelin, A., 2001. Competitive sorption and transport of heavy metal cations insoils. Diss. ETH. No. 14101, Zürich, Switzerland.
Voegelin, A., Vulava, V., Kretzschmar, R., 2001. Reaction-based model
describing competitive sorption and transport of Cd, Zn, and Ni in anacidic soil. Environ. Sci. Technol. 35, 1651–1657.
Vogel, T., Cislerova, M., Hopmans, J.W., 1991. Porous media with linearlyvariable hydraulic properties. Water Resour. Res. 27, 2735–2741.
Yeh, G.-T., Cheng, H.-P., 1999. 3DHYDROGEOCHEM: A 3-dimensionalmodel of density-dependent subsurface flow and thermal multispeciesmulticomponent HYDROGEOCHEMical transport. EPA/600/R-98/159.