MZ

Modeling Coupled Hydrologic and Chemical
Processes: Long-Term Uranium Transport
following Phosphorus FerƟlizaƟon

S

S

:V

D. Jacques,* J. Šimůnek, D. Mallants, and M.Th. van Genuchten
Coupling physical and geochemical processes within one integrated numerical simulator provides a process-based toolfor inves ga ng the mobility of contaminants as aﬀected by changing hydrologic regimes and geochemical condi ons.
We review interac ons between physical and biogeochemical processes in the vadose zone, and then present a casestudy demonstra ng these complex interac ons. A hypothe cal applica on is presented of the HP1 mul component
transport simulator to predict the transport of two major elements (Ca and P) and one trace element (U) applied annually for 200 yr to a field soil in the form of an inorganic P fer lizer. Interac ons of Ca, P, and U with the solid phase aredescribed using ca on exchange and surface complexa on reac ons. Simula ons assuming steady-state or transient
flow condi ons were analyzed in terms of temporal varia ons of the linear distribu on coeﬃcient, Kd, which depends
strongly on pH and the composi on of the aqueous phase. If the composi on of the aqueous phase is constant, adsorpon of Ca and U increases with increasing pH. Due to the annual addi on of Ca, P, and U, and compe on between P
and U for sorp on sites, the Kd of these elements decreased with me near the soil surface. Deeper in the soil, the Kdof U followed the pH increase because of a lack of compe on from P. Because of the combined eﬀects of changing
hydrologic and geochemical condi ons, the Ca and U distribu on coeﬃcients and solute fluxes during the transient
simula on exhibited large short- me varia ons of up to three orders of magnitude.

Fare complex three-dimensional heterogeneous
systems involving the transport, redistribution, and transformation of multiple components in diﬀerent phases under
nonisothermal conditions. The various transformations are to alarge extent initiated, directed, and catalyzed by (micro)biological agents. The transport of components is not limited to only
chemical elements such as major cations and anions, heavy metals,
radionuclides, or organic pollutants, but often also involves organicand inorganic colloids and microorganisms. Vegetation furthermore interacts with the soil, mainly through its roots, which extract
water and nutrients from the soil. Roots at the same time releasevarious organic ligands to mobilize nutrients and alter geochemical
conditions near the root–soil interface (e.g., Seuntjens et al., 2005;
Calba et al., 2004; Geelhoed et al., 1999).
An integrated approach is necessary for understanding andunraveling the complex interplay of different soil processes.
D. Jacques and D. Mallants, Performance Assessments Unit, Belgian
Nuclear Research Centre (SCK⋅CEN), Boeretang 200, B-2400 Mol, Belgium;J. Šimůnek, Dep. of Environmental Sciences, Univ. of California, Riverside,
CA 92521; M.Th. van Genuchten, USDA-ARS, U.S. Salinity Lab., 450 W. Big
Spring Rd., Riverside, CA 92507. Received 27 Apr. 2007. *Corresponding
author (djacques@sckcen.be).
Vadose Zone J. 7:698–711
doi:10.2136/vzj2007.0084© Soil Science Society of America
677 S. Segoe Rd. Madison, WI 53711 USA.
All rights reserved. No part of this periodical may be reproduced ortransmi ed in any form or by any means, electronic or mechanical,
including photocopying, recording, or any informa on storage andretrieval system, without permission in wri ng from the publisher.
www.vadosezonejournal.org · Vol. 7, No. 2, May 2008

698

Integration of all pertinent processes into a numerical transportmodeling tool is indispensable for evaluating the various coupledeﬀects involved (e.g., Mayer et al., 2002; Metz et al., 2002; Daviset al., 2004b; Lichtner et al., 2004; Appelo and Postma, 2005;
Nützmann et al., 2005; Nuclear Regulatory Commission, 2006;
Carrillo-González et al., 2006). Steefel et al. (2005) brieﬂy discussed some historical steps in the development of such tools.
Many reactive transport codes of diﬀerent levels of complexity
have been developed during the last few decades for saturated
porous media (e.g., van der Lee and De Windt, 2001). Morerecently, similar codes started to consider also ﬂow and transport
in the unsaturated zone. Examples are 3DHYDROGEOCHEM(Yeh and Cheng, 1999), CORE2D (Samper et al., 2000), MIN3P(Mayer et al., 2002), Pollutrans (Kuechler and Noack, 2002),
UNSATCHEM-2D (Šimůnek and Suarez, 1994), RETRASO(Saaltink et al., 2004), HYDRUS-1D (Šimůnek et al., 2005),
and HP1 (Jacques et al., 2006b). Thus far, most applicationsof these multicomponent transport models to both saturated
and unsaturated conditions have been limited to steady-state
ﬂow scenarios.
Transient ﬂow in soils can substantially aﬀect geochemical conditions for both equilibrium and kinetic reactions.
Experimental evidence shows that alternations in inﬁltration(due to precipitation and irrigation) and evaporation periods,
or changes in temperature, can strongly inﬂuence geochemical
conditions near the soil surface. Such eﬀects have been demonstrated during experimental cycles in laboratory setups (e.g.,Öztürk and Özkan, 2004), seasonal variations in the ﬁeld (e.g.,Berner et al., 1998), and for artiﬁcial recharge of groundwa-

ter (Greskowiak et al., 2005). The interaction of discharging
groundwater and shallow groundwater dynamics with temporalvariations in atmospheric conditions also aﬀects the soil watercomposition (e.g., de Mars et al., 1997; de Mars and Garritsen,
1997; Joris, 2005). Slattery and Ronnfeldt (1992) showed thatseasonal variations in pH were partly related to ionic strength.
Gonçalves et al. (2006) reported an experimental study to quantify salinization and alkalization risks in lysimeters exposed tonatural atmospheric conditions (rainfall and evapotranspiration)and manually irrigated with waters of diﬀerent quality. Overallsalinity and Na concentrations changed signiﬁcantly during theyear. Sodium concentrations near the soil surface peaked during
dry periods after irrigation events, but subsequently decreased
again to initial values during rainfall periods. The experimental data were successfully simulated using HYDRUS-1D, which
couples an ion chemistry module with a water ﬂow and solutetransport model.
A similar behavior for Na was found by Jacques et al. (2008)
in a hypothetical simulation of the 30-yr migration of heavy
metals under atmospheric conditions. During periods with a precipitation deﬁcit (i.e., when potential evaporation is larger than
precipitation), the pH near the soil decreased at the same timeas the water content. The pH decrease was caused by increased
ionic strength (i.e., higher concentrations) and the accumulationof anions such as Cl− near the soil surface due to upward ﬂow.
In addition, the total amount of Cd in soil water also increased
during the summer due to (i) increased competition of monovalent cations for the exchange sites (especially with Na, which is amore mobile element than K and thus accumulates near the soil
surface during upward ﬂow), and (ii) increased aqueous complexation of Cd and Cl.
Water content variations also inﬂuence the rate of kinetic
reactions due to their eﬀects on the wetted reactive area. The
reactive surface area is a critical parameter aﬀecting the rate ofmineral dissolution. It is not straightforward to link reactive surface area to experimentally (e.g., BET mineral surface areas) or
geometrically determined values. Brantley (2003) reported several
reasons for this discrepancy: (i) the diﬀerence in reactivity ofdiﬀerent surface sites (e.g., edge vs. basal plane sites for phyllosilicates); (ii) aging eﬀects; (iii) the occurrence of deep etch cavities
inducing transport-limited microenvironments for dissolution;and (iv) coatings of natural mineral surfaces with organic matter,
Fe, Al, and Si. The critical parameter during transient unsaturated
ﬂow is the speciﬁc solid–water interfacial area, which varies as afunction of the water content. In a theoretical study, Saripalli etal. (2006) showed that the interfacial area changed only slightly
near saturation since the larger pores entail only a small percentage of the interfacial area, whereas a much stronger decrease
in interfacial area occurs under dry conditions. Simulations forsteady-state ﬂow conditions involving diﬀerent degrees of saturation illustrated that changing interfacial areas can have an
important eﬀect on solute ﬂuxes. This eﬀect of the degree ofsaturation was recently illustrated by Kuechler and Noack (2007)
for pyrite–calcite dissolution during unsaturated ﬂow. Only 5%of the batch-measured reactive surface area for pyrite was required
to describe unsaturated column experiments at a degree of water
saturation of 11%.
These examples illustrate the important eﬀects transient waterﬂow can have on multicomponent solute transport. Our mainwww.vadosezonejournal.org · Vol. 7, No. 2, May 2008

objective is to further demonstrate these types of interactions
between transient ﬂow and biochemical processes in soils. We
ﬁrst discuss various feedbacks between transient ﬂow and otherprocesses, and then present a hypothetical example involving the
migration of several chemical species during a 200-yr time period.
The governing vadose zone ﬂow and chemical equations of thecoupled HP1 multicomponent reactive transport code (Jacquesand Šimůnek, 2005) are given next. Our focus is especially on
nonlinear feedbacks between diﬀerent processes relevant for soil
systems. Then we show the important eﬀects of the upper boundary condition (steady-state vs. transient ﬂow) on the long-term
migration of the main constituents (Ca and P) and a trace element (U) of an inorganic P fertilizer applied annually for 200 yr
to a soil proﬁle. The example demonstrates how transient watercontents and ﬂuxes aﬀect soil pH, and hence the retention, bioavailability, and ﬂuxes of various solute components.

Coupled ReacƟve Transport in the Vadose Zone
We will focus on interactions between various processes inthe vadose zone by deﬁning the governing equations (see also
Steefel et al., 2005). Without attempting to be complete, selected
feedbacks between processes important for unsaturated ﬂow andtransport problems will be discussed here. As a guide, we will
outline the equations used in the HP1 multicomponent reactive transport simulator. This simulator couples two previously
independent codes into one comprehensive tool. These are the
HYDRUS-1D ﬂow and transport model (Šimůnek et al., 2005)
and the PHREEQC-2 biogeochemical model (Parkhurst and Appelo,1999). More information about HP1, including the code itself, canbe found at www.sckcen.be\hp1 (veriﬁed 18 Apr. 2008).
Processes and Reac ons

Physical transport is described using the conservation
equations for momentum, ﬂuid mass, solute mass, and energy.
Typically for soils and other porous materials, momentum conservation is described with the Darcy–Buckingham equation in
combination with the conservation equation for water to yield
the Richards equation as given by Eq. [1.1] in Table 1 (symbols
are deﬁned in the Appendix). We assume that the gas phase iscontinuous and at atmospheric pressure, which simpliﬁes thedescription of ﬂuid ﬂow from a multiphase to a single-phase ﬂow
problem. To solve the Richards equation, the soil water content,
θ [L3 L−3], as a function of the soil water pressure head, h [L],and the unsaturated hydraulic conductivity, K [L T−1], as a function of θ or h must be deﬁned (e.g., Brooks and Corey, 1966;
van Genuchten, 1980). The Richards equation describes uniform
water ﬂow in soils, but not any preferential or nonequilibrium
ﬂow. Alternative formulations for such conditions include various dual-porosity or dual-permeability models (see Šimůnek etal., 2003; Köhne et al., unpublished data, 2008; Šimůnek andvan Genuchten, 2008).
When neglecting water vapor diﬀusion, the one-dimensional
heat transport equation takes the form of Eq. [1.2] in Table 1.
The volumetric heat capacity and thermal conductivity depend
on various soil constituents (Campbell, 1985; Šimůnek et al.,2005), including the soil water content. Heat and water are also
simultaneously transported as water vapor, which provides strong
coupling between soil water and heat dynamics (e.g., Saito et al.,2006). The latter process is not included in HP1.
699

T

1. Governing flow and transport equa ons, equilibrium geochemical mass-ac on equa ons, and examples of kine c reac on equaons used in HP1.
Process

Equa on†Equa on no.

Flow and Transport Equa ons
Water flow (Richards equa on)∂θ ( h ) ∂ ⎡⎤=K ( h ) ∂ h + cos α ⎥ − S ( h )∂t∂ x ⎢⎣∂x⎦[1.1]

Heat transport∂ C p ( θ )T∂ qT
= ∂ ⎡⎢ λ ( θ ) ∂ T ⎤⎥ − C w− C w ST
∂t
∂x ⎣
∂x ⎦
∂x

[1.2]

Solute transport (advec ondispersion equa on)⎛∂θ C j∂ C j ⎞⎟ ∂ qC j⎟−
= ∂ ⎜⎜ θ D w− SC r, j + Ro, j∂t∂ x ⎜⎝∂ x ⎠⎟⎟
∂x

[1.3]

(

)

Chemical Equilibrium Reac ons

Aqueous specia on

Reac on Equa on

Mass Ac on Law

Nm

K il = γ il c i ∏ (γ mj c mj )Nm

∑ ν lji A mj = Ai
j =1

Nm

j =1

e

e e

e

Nm

∑ ν sji A mj + ν sj i S mj = Si
j =1

s

s s

s

j =1

∑ ν pji A mj = M i
j =1p−ν eje ie(γ ej βej , j )
e

e

e

[1.5]

K is,int
=
e
Nm
⎡
e
⎛FΨj
⎞
−ν s ⎤
m m −ν jis⎢ s(β sjs , js ) jsis ⎥⎥ exp ⎜⎜⎜⎝⎜ RT s Δ z is ⎠⎟⎟⎟⎟⎢ β is , js ∏ (γ j c j )
j =1
⎢⎣
⎥⎦

s

Nm

Mineral dissolu on−ν ejie

K iee = γ eie β eie , je ∏ (γ mj c mj )eNm

Surface complexa on

[1.4]

j =1

∑ νeji A mj + νej i X mj = X i

Ion exchange−ν ljiNm−ν pjiK ipp = ∏ (γ mj c mj )

p

[1.6]

[1.7]

j =1

Chemical Kine c Reac ons
General rate equa on for
aqueous speciesNkdci= ∑ ν kjk ,i R jk
dt
j =1

[1.8]

k

Mineral dissolu on equa on(Lasaga, 1995)

n +R j = k0 Amin exp (−E a RT )a H+ g ( I ) ∏ ai i f (Δ Gr )
k

n

H

i

aﬃnity term based on transi on state theory (Aagaard and Helgeson, 1982; Lasaga, 1981, 1995, 1998):

f (Δ G r ) = 1− Q
Monod rate equa on (Schäfer
et al., 1998)

[1.9]

K ipp
N M , jk

K Ij ,iN I , jkc j ,m

R jk = v max, jk X mo,r ∏

∏

k

k

M
I
m=1 K j ,m + c j ,m i =1 K j ,i + c j ,mkkk

[1.10]k† Symbols are defined in the Appendix.

where cj is the concentration of the jth component, ci is the concentration of the ith secondary species, νji is the stoichiometriccoeﬃcient of the jth component in the reaction equation for the
ith species, and Nsa is the number of secondary species. Equation[1.3] is based on the assumption that diﬀusion coeﬃcients of all
species are equal (see also Lichtner, 1996; Mayer, 1999). The
Ro,j term in Eq. [1.3] represents sources–sinks for components
in the aqueous phase, accounting for various heterogeneous equilibrium and kinetic reactions and homogeneous kinetic reactions.
Depending on the type of component, diﬀerent reactions andprocesses will be involved.
Heterogeneous reactions involve chemical elements in otherthan the liquid phase and include adsorption reactions described

The conservation of solute mass combined with Fick’s lawof diﬀusion results in the advection–dispersion equation as given
by Eq. [1.3] in Table 1. This equation is not only applicable
to dissolved chemical elements, but also to colloids (e.g., vanGenuchten and Šimůnek, 2004; DeNovio et al., 2004; Šimůnek
et al., 2006) and microorganisms (e.g., Taylor and Jaﬀé, 1990; Jinet al., 2002; Rockhold et al., 2004, 2005; Gargiulo et al., 2007).
Note that Eq. [1.3] in Table 1 is expressed in terms of the total
concentration C [mol L−1] of a given component j:
N sa

C j = c j + ∑ ν ji c i

[1]

i =1

www.vadosezonejournal.org · Vol. 7, No. 2, May 2008

700

by adsorption isotherms (linear, Freundlich, or Langmuir isotherms), cation exchange (e.g., Voegelin et al., 2001; Bruggenwert
and Kamphorst, 1982), or surface complexation reactions (e.g.,
see Goldberg et al. [2007] for a review of soil systems), mineraldissolution–precipitation, and gas dissolution–exsolution. Except
for adsorption reactions described by isotherms, the other three
types of reactions, together with aqueous speciation, account for
multicomponent interactions. Their equilibrium state is determined by mass-action laws such as given by Eq. [1.4–1.7] inTable 1. See Morel and Hering (1993), Langmuir (1997), andAppelo and Postma (2005) for background on activities, exchange
conventions, and surface complexation models. Heterogeneous
equilibrium is obtained by inserting the number of moles ofthe species (derived from the mass-action equations) into themole- and charge-balance equations. Heterogeneous reactions
can also be treated kinetically (e.g., Dang et al., 1994; TebesStevens et al., 1998).
Kinetic homogeneous reactions are reactions that occur solely
in the liquid phase and include radioactive decay, biologicallymediated degradation, and solute uptake by roots (Barber, 1995;
Tinker and Nye, 2000; Hopmans and Bristow, 2002; Nowack etal., 2006). Examples of complex kinetic rate equations are mineral dissolution and precipitation rates based on transition-state
theory (e.g., Lasaga, 1995; Mayer et al., 2002), Monod-type rateequations (e.g., Schäfer et al., 1998; Salvage and Yeh, 1998), andbioﬁlm models (e.g., Cunningham and Mendoza-Sanchez, 2006)
for biodegradation and microbial growth such as those given by
Eq. [1.8–1.10] in Table 1.
Additional heterogeneous reactions may occur for colloids
and microorganisms, such as attachment–detachment processes atthe solid phase, straining, ﬁltration, partitioning of colloids to theair–water interface, and size exclusion (e.g., Wan and Tokunaga,
1997, 2002; Bradford et al., 2003, 2006; Šimůnek et al., 2006;
Gargiulo et al., 2007). Biomass production and decay are otherimportant processes that need to be taken into account when
predicting the fate of microorganisms in soil systems.

feedback (i.e., the eﬀect of temperature on solute transport andreactions) is, however, much more important. Solute transportand reactions additionally are strongly aﬀected by the temperaturedependence of the diﬀusion coeﬃcient and by thermodynamic
constants and reaction rates.
How and what processes are incorporated into the reaction
term may have a considerable eﬀect on the predicted mobilityof various components. Many unsaturated ﬂow and transportmodels coupling water ﬂow and solute and heat transport incorporate reaction terms for single components to account for suchreactions as equilibrium or nonequilibrium sorption by using
linear, Freundlich, or Langmuir isotherms and ﬁrst- or zero-orderkinetics. In these types of models, the concentration of one element has no eﬀect on the fate of other elements except when
considering a sequential decay reaction. By including interactions between diﬀerent components, one can calculate aqueous
speciation, competitive adsorption described by cation exchangeor surface complexation reactions, or complex kinetic reactions(see Table 1). Aqueous speciation can enhance the apparent solubility, adsorption potential, and mobility of some of the elements
due to cation–ligand complexation. Examples of such complexation reactions are U complexes (Waite et al., 1994; McKinley
et al., 1995; Morrison et al., 1995; Lenhart and Honeyman,
1999; Lenhart et al., 2000; Barnett et al., 2002) or the formation of soluble Al– or Fe–organic matter complexes (Jansen et al.,
2002; Carrillo-González et al., 2006). As discussed above, variably saturated ﬂow aﬀects many geochemical variables (e.g., pH,concentrations) that consequently inﬂuence geochemical equilibrium. A linear equilibrium single-component adsorption modelin which the solute distribution between the liquid and solid
phases is described using the distribution coeﬃcient, Kd, leads tocompletely diﬀerent concentration patterns than those predictedusing a competitive cation exchange model (e.g., Gonçalves et al.,2006) or a surface-complexation model (Zhu, 2003). This is dueto the inability of the Kd approach to account for spatial and temporal variations in the soil water chemistry and water contents.
These variations lead to changes in the ratio between exchanger
and ﬂuid content and consequently aﬀect competition between
cations for the exchange complex (Steefel et al., 2005).
An important feedback of geochemical reactions on transport
processes is the eﬀect of mineral dissolution and precipitation
on not only the overall porosity of the system, but also on the
pore-size distribution. This feedback aﬀects all transport-related
soil properties. These properties or parameters can be updated invarious ways. For example, diﬀusion coeﬃcients can be evaluatedusing Archie’s law by introducing a cementation factor (Soler andLasaga, 1998; van der Lee et al., 2003), while saturated hydraulicconductivities can be modiﬁed using macroscopic or grain-scalemodels (e.g., Le Gallo et al., 1998), and the unsaturated hydraulicproperties using various physical models (Freedman et al., 2004)
or pedotransfer functions (Finke, 2006).
In addition to geochemical reactions, microbial colonizationand microbially produced compounds can also signiﬁcantly impact
the hydraulic regime. Microbial processes, for example, may lower
the capillary fringe or induce locally dry zones (Yarwood et al.,2006). Microbial processes can aﬀect not only the physical properties of a porous medium by clogging pores, producing capsules,
or generating gases, but also the air–liquid interfacial tension or
the solid–liquid contact angle due to the production of surface-

Feedback

The heat transport equation (Eq. [1.2]) is coupled to the
Richards equation (Eq. [1.1]) through the eﬀects of the ﬂuid
ﬂux on heat advection and the apparent thermal conductivity.
Conversely, temperature aﬀects surface tension, viscosity, andﬂuid density and, consequently, the soil hydraulic properties (e.g.,
Zhang et al., 2003).
Integrating the water ﬂow and advection–dispersion equations into one numerical tool allows one to account for the eﬀects
of water ﬂow on solute fate and transport. Solute advection andhydrodynamic dispersion are directly inﬂuenced by water ﬂuxes,
while many chemical and biological reactions depend strongly
on water contents. Water contents also indirectly inﬂuence soil
O2 and CO2 concentrations, which have a direct eﬀect on pH
and consequently on many chemical and biological reactions andprocesses. Although CO2 transport and production processes are
considered in HYDRUS-1D, they have not yet been incorporated
into HP1. Solute concentrations aﬀect the density of the liquid
phase and can also inﬂuence root water uptake due to salinity
stress (van Genuchten, 1987; Hopmans and Bristow, 2002). Soil
temperatures may be aﬀected by heat produced or consumed
during various geochemical and biological reactions. The reverse
www.vadosezonejournal.org · Vol. 7, No. 2, May 2008701

active compounds (Yarwood et al., 2006). Rockhold et al. (2002)
modeled such changes by using pressure head scaling factors.
Several of the processes and feedbacks discussed here (notably
the eﬀects of chemical precipitation–dissolution and microbial
growth on porosity and pore size distribution and CO2 transportand production) are not yet included into the current version of
the HP1 simulator (Jacques and Šimůnek, 2005).

ally deposited U by fertilizers is transported to drainage channels
in the Kanovci area of the Republic of Croatia. Their data indicated that high U concentrations in surface waters and shallow
groundwater were caused by phosphate fertilizer applications. On
the other hand, Zielinski et al. (1997) found that leaching of U
from inorganic P fertilizers contributed much less to high U concentrations in drainage water than natural or irrigation-enhanced
leaching from soils. The radiological impact expressed as dose isalso quite limited (Saueia and Mazzilli, 2006).
Zielinski et al. (1997) pointed out that the mobility and thenature of P and U ﬁxation and precipitation depend on soil conditions, such as pH, moisture, mineralogy, and texture. Guzman
et al. (2002) provided a number of processes controlling U mobility. The list included sorption on solid surfaces, formation ofcomplexes with Ca, PO42−, and CO32−, and (co-)precipitationof U–P, Ca–U, Al–P, or Fe–P (acid soils) or Ca–P (alkaline soils)
minerals (see also Zielinski et al., 1997; Brady, 1990).
To model the subsurface transport of U, either naturally
occurring in soils or applied with mineral fertilizers, the interaction of U with other elements and the solid phase needs tobe considered. The amount of U adsorption to the solid phasedepends on the speciﬁc composition of both the solid phase (e.g.,clay minerals, Fe and Al oxides, and organic matter) and theaqueous phase (e.g., pH, complexing agents such as Cl−, SO42−,
PO42−, or the presence of competing cations for cation exchangeor surface complexation). The amount of Fe oxides in the soil
strongly aﬀects the pH dependence of U sorption (Barnett et al.,2000). Changing chemical conditions signiﬁcantly inﬂuences themobility of U in soils and sediments (Kohler et al., 1996). Here
we use HP1 to obtain insight into the complex system of interacting geochemical processes that govern Ca, P, and U mobility
in the soil when an inorganic P fertilizer is applied annually. Our
study is an update of preliminary calculations by Jacques et al.
(2006a). We focus especially on the eﬀects of the imposed upperboundary condition (constant vs. transient atmospheric boundary
conditions) on the migration of U in an undisturbed acid soil
proﬁle. The analysis of migration of Ca, P, and U is performedin terms of temporal variations in the distribution coeﬃcient Kd
to illustrate the limitation of such an empirical approach, which
was illustrated in previous studies by Bethke and Brady (2000),among others.

Case Study: Uranium Leaching from Soils
Following Long-Term Mineral
Phosphorus FerƟlizaƟon
Background

Common inorganic phosphate fertilizers such as triple
superphosphate, single superphosphate, monoammonium phosphate, and diammonium phosphate are produced from apatite,
Ca5(PO4)3(OH,F,Cl). Both PO42− and Ca during production
are replaced by other elements such as radionuclides from the
238U and 232Th decay series (mainly 226Ra and 210Pb) andrare earth elements originating from the original raw material(Saueia and Mazzilli, 2006). Phosphoric acid and the byproductphosphogypsum are produced using H2SO4 according to thefollowing reaction:
Ca10F2(PO4)6 + 10 H2SO4 + 10 H2O →10 CaSO4⋅2H2O + 6 H3PO4 + 2 HF

[2]Phosphogypsum (a calcium sulfate dihydrate) is often enrichedwith 226Ra because of its chemical similarity to Ca (Papastefanou
et al., 2006, Santos et al., 2006). Phosphoric acid produced
during the phosphate production process reacts with phosphaterocks to yield monocalcium phosphate or triple superphosphate,with overall reactions as follows (Saueia and Mazzilli, 2006):
Ca10F2(PO4)6 + 7 H2SO4 + 6.5 H2O →3 Ca(H2PO4)2⋅H2O + 7 CaSO4.0.5H2O+ 2 HF

[3]Ca10F2(PO4)6 + 10 H3PO4 + 10 H2O →10 Ca(H2PO4)2⋅H2O + 2 HF

[4]Problem Defini on

Compared with phosphate rock, the mineral P fertilizers areenriched in U since U remains in the liquid phase together with
H3PO4 during the acidiﬁcation of the phosphate rock. The
U concentration of the P fertilizer is related to the amount of
the phosphate in the fertilizer (Barišić et al., 1992). Uranium
concentrations in P-bearing fertilizers have been reported to bein the range of 300 to 3000 Bq kg−1 of fertilizer (European
Commission, 1999), or 1700 to 9200 Bq kg−1 of fertilizer forboth 238U and 234U (Cogné, 1993).
The environmental and radiological impacts of the application of inorganic P fertilizers are, in general, relatively modest.
Spalding and Sackett (1972) attributed the increase of U in North
American rivers to applications of large amounts of P fertilizers
to agricultural lands. Rothbaum et al. (1979) found that most U
originating from applied superphosphate was retained within theplow layer of arable soils or in the organic surface layers undergrassland. Barišić et al. (1992) concluded that >20% of the annuwww.vadosezonejournal.org · Vol. 7, No. 2, May 2008

Phosphorus fertilizers are often applied annually to agricultural ﬁelds, partly in inorganic form. In Flanders (northern
Belgium), the amount of applied P that originates from mineral
fertilizers decreased from 13.9 × 106 kg P in 1990 to 2.38 × 106
kg P in 2003 (MIRA Team, 2004). To simulate the long-term
migration of P and U through the soil, the following assumptions
were made: (i) an average P fertilization of 1 g P m−2 (or 6.35
× 106 kg P for all of Flanders) is applied each year on 1 May inthe form of Ca(H2PO4)2; and (ii) 1 kg of Ca(H2PO4)2 contains
10−3 mol of U. Based on a speciﬁc activity 2.96 × 106 Bq mol−1
U, 2960 Bq of alpha activity exclusively due to 238U is present
in 1 kg of Ca(H2PO4)2. Thus, a fertilization of 1 g P m−2 corresponds to 1.61 × 10−2 mol of Ca(H2PO4)2 containing 3.77
× 10−6 mol of U.
The soil proﬁle studied here is a dry Spodosol located atthe “Kattenbos” site (Lommel, Belgium) consisting of seven soil702

2. Soil hydraulic parameters† (van Genuchten, 1980) of the dry Spodosol (data fromhorizons in the top 1 m. The soil has a T
Seuntjens[2000] and Seuntjens et al. [2001a]).
typical leached E horizon between 7 and19 cm below the surface and enriched (inHorizon‡
Depth
n
Ks
θr
θs
α
organic matter and Fe oxides) Bh hori−1
cm
m
m s−1
zons between 19 and 28 cm. We note Top mineral layer, A
0–7
0.065
0.48
1.6
1.94 1.1 × 10−5that the top horizons are not typical for Elluviated layer, E
7–19
0.035
0.42
1.5
3.21 3.6 × 10−5
an agricultural soil, but rather reﬂect an Illuviated layer with organic ma er Bh1
19–24
0.042
0.47
1.6
1.52 4.5 × 10−6
undisturbed proﬁle covered by natural
Bh2
24–28
0.044
0.46
2.8
2.01 1.0 × 10−4
vegetation (i.e., heather). The proﬁle Transi on layer, BC
28–50
0.039
0.46
2.3
2.99 1.4 × 10−4
C1
50–75
0.030
0.42
2.1
3.72 1.4 × 10−4
was selected for our modeling pur- Unconsolidated material
C2
75–100 0.021
0.39
2.1
4.33 1.4 × 10−4
poses because it was well characterized
in terms of physical and geochemical † θ and θ are the residual and saturated water contents, respec vely; α and n are shape parameters; andrsproperties. Table 2 gives the locations
Ks is saturated hydraulic conduc vity.
of the diﬀerent soil horizons and their ‡ Descrip on from Brady (1990).
soil hydraulic parameters as determinedon undisturbed soil samples of 100 cm3(Seuntjens, 2000; Seuntjens et al., 2001a).
4 and 6 molc kg−1 organic matter, respectively (Tipping, 2002).
Our conceptual chemical model considers the transport ofWe used an average value of 6 molc kg−1 organic matter.
the following 14 components: the elements C, Ca, Cl, F, Mg,
Both cations and anions adsorb on Fe oxides. Adsorption
N(V), Na, P, S(VI), and U(VI), H2O, O (not included in water),
was described by speciﬁc binding using a nonelectrostatic surface
H (not included in water), and the charge of the soil solution.
complexation model. Surface complexation reactions included
Normally a solution is always charge balanced; however, in asurface hydrolysis reactions (protonation and deprotonation) andmodel that includes surface reactions without compensating thereactions involving the cations Ca2+, Mg2+, and UO2+ and thesurface charge in the diﬀuse double layer, both the aqueous phaseanions PO42−, SO42−, and Fe− (denoting ≡FeOH as a neutral
and the solid surface will have a charge of opposite sign. The total
Fe oxide surface) as follows: ≡FeOH2+, ≡FeO−, ≡FeH2PO4,
system must still be charge balanced.
≡FeHPO4−, ≡FePO42−, ≡FeSO4−, ≡FeOHSO42−, ≡FeF,
Three types of geochemical equilibrium reactions were
≡FeHF−, ≡FeOMg+, ≡FeOCa+, and ≡FeOUO2+. Equilibrium
taken into account: (i) aqueous speciation reactions; (ii) cation
constants for weak sites were taken from Dzombak and Morelexchange reactions on organic matter; and (iii) surface complex(1990) (see also Parkhurst and Appelo, 1999). The capacity of theation reactions on Fe oxides. The thermodynamic data for thesurface was derived from the amount of Fe2O3 as measured in the
aqueous species were taken from the WATEQ4F database (Ballsoil proﬁle assuming 0.875 reactive sites per mole of Fe (Waite etand Nordstrom, 1991), complemented with data for U fromal., 1994). The cation exchange capacity and the reactive sites ofLangmuir (1997).
each horizon are given in Table 4.
The cation exchange complex was assumed to be associatedAlthough CO2 in general is important for U speciation, thissolely with organic matter. Organic matter typically consists of an
type of complexation reaction is not important in our simulations
assemblage of diﬀerent functional groups and types of adsorption
since the pH varied only between 3.5 and 4.5. The partial pressites. A multisite cation exchange complex allows simulations ofsures of O2 and CO2 in the soil air phase were assumed to be
the acid–base properties of the organic matter and the correspondconstant and equal to 20 and 0.32 kPa, respectively. These values
ing increase in the cation exchange capacity with increasing pH
are smaller and larger, respectively, than the partial pressures in(e.g., Appelo et al., 1998). Appelo et al. (1998) included proton
the atmosphere. The assumption of having a constant CO2 partial
exchange on a multisite cation exchange complex to mimic thepressure represents a signiﬁcant simpliﬁcation since soil CO2 conpresence of diﬀerent functional groups. Six cation exchangers
centrations in general display strong seasonal variations related towere assumed to be present, each with a diﬀerent Ke value forthe microbiological activity in the soil (e.g., Šimůnek and Suarez,
the half reaction for H+:
1993; Suarez and Šimůnek, 1993).
Meteorological data from two weather stations (from Brogelfor years 1969–1985 and from Mol/Geel for 1986–1998) wereY i + H+ = H Y i[5]used to deﬁne precipitation, P [L T−1], and potential evaporawhere Yi is a cation exchange site (i = a, b, c, d, e, f ). Table 3 gives
tion, Ep [L T−1], rates, the latter based on the Penmann equationethe K values for the multisite exchange complex. Dependingon the pH, some sites are deprotonated such that cations can T
3. Logarithmic equilibrium constants (Ke) parameters for the mul adsorb. Other sites, however, will remain protonated and are site exchange complex.
thus not available for cation exchange. Besides H+, ﬁve cation
YNaY
KY
MgY2
CaY2
UO2Y2
species (Na+, K+, Mg2+, Ca2+, and UO2+) participated in the Exchanger†
−1.0
−0.3
−0.4
−0.2
−0.2
exchange reactions. The exchange coeﬃcients for the cations HY‡
HYa
HYb
HYc
HYd
HYe
HYf
were assumed to be the same on all six sites. The size of the
1.65
3.3
4.95
6.85
9.6
12.35
total cation exchange complex was estimated from the organic† The value for NaY was taken from Appelo et al. (1998). Values for the other commatter content and the amount of exchangeable protons on
plexes were taken from the phreeqc.dat database (Parkhurst and Appelo,
the organic matter. The proton dissociating groups on fulvic
1999) and adapted rela ve to the K for NaY.
and humic acids are between 6 and 10 molc kg−1 and between ‡ Values taken from Appelo et al. (1998).
www.vadosezonejournal.org · Vol. 7, No. 2, May 2008703

T
4. Bulk densi es, weight contents of organic ma er andFe2O3 of the dry Spodosol (Seuntjens et al., 2001a; Seuntjens,
2000), and calculated sizes of the mul site ca on exchange complex (CEC) and the Fe-oxide surfaces.
Horizon Bulk densityg cm−3
1.31
1.59
1.30
1.38
1.41
1.52
1.56

A
E
Bh1
Bh2
BC
C1
C2

Organic
ma er

Fe2O3

————%————
3.65
0.20
0.78
0.09
3.03
1.19
1.23
1.19
0.55
0.68
0.30
0.57
0.13
0.57

CEC

composition, the selected soil air partial pressures of O2 and CO2,
and an assumed aqueous U concentration of 8.0 × 10−9 mol
L−1. This value corresponds with the average U concentration at
5-cm depth obtained after 50 yr in the steady-state ﬂow simulations. When the background P concentration in rainwater is the
only P source, UO22+ is the dominant U species in the soil solution, followed by UO2F+, UO2SO4, UO2OH+, and UO2HPO4.
When the P concentration is much higher due to the additionof the P fertilizer, UO22+ is still the dominant U species, followed by UO2HPO4− (dominant at pH values >4), UO2F+,
and UO2OH+. For both cases, H2PO4− was the most abundantaqueous P species, followed by H3PO4.
Surface speciation was also calculated for the same pH range,
the rainwater composition, the soil air partial pressures of O2 andCO2, and the total aqueous U concentration of 10−9 mol L−1.
The Na concentration was adapted to obtain the required pH
and charge balance. Figure 1a shows the most abundant surface
species on the Fe oxides. The surface of the Fe oxides was mostly
protonated. Below pH 5.5, FeH2PO4 was the second most abundant surface species. The amounts of FeOH, FeHPO4−, andFeOUO2+ increased substantially across the considered pH range.
For example, FeOUO2+ increased up to seven orders of magnitude. The amount of adsorbed U consequently also increased

Sites onsurface——molc m−3 soil ——
287
28.7
74
15.7
236
169.5
102
180.0
46
105.1
27
94.9
12
97.4

(Penmann, 1948). Statistical information about the imposed climatological time series is given in Jacques et al. (2008). To obtain
a 200-yr time series for modeling purposes, the 30-yr period was
simply repeated several times. A free-drainage boundary condition was assumed at the bottom of the soil proﬁle.
Element concentrations in the rainwater were obtained
from Stolk (2001) for Station 231 located in Gilze-Rijen (theNetherlands), which is sufficiently close to the investigatedsite. The solution composition was based on the average of 13
measurements during 1999 (Table 5). The U concentration of
the rainwater was set equal to an arbitrary low value of 10−24
mol L−1 to have the required non-zero initial concentration
for PHREEQC simulations. The P fertilizer was assumed to be
applied every year on 1 May in the form of Ca(H2PO4)2 in 1 cmof irrigation water with the same composition as rainwater. The
fertilizer was assumed to be completely dissolved in the irrigationwater, thus neglecting any possible kinetics in the dissolution
process of the fertilizer grains. Before the 200-yr simulation, a5-yr steady-state ﬂow simulation was performed without fertilizer
applications to establish chemical initial conditions in topsoil horizons in equilibrium with the chemistry of the inﬂowing rainwater.
Simula on Results
Specia on of Phosphate, Uranium, and the Solid-Phase Surfaces

Aqueous speciation of U and P in the soil solution was calculated with PHREEQC between pH 2 and 6 for the rainwaterT
5. Composi on of the rainwater (Stolk, 2001) and rainwaterwith P fer lizer.
Element

Rainwater

Rainwater + P fer lizer————————mol L−1 ————————
C(4)
Ca
Cl
F
K
Mg
N(V)
Na
P
S(VI)
U(VI)
pO2(g), Pa
pCO2(g), PapH

1.45 × 10−5
6.00 × 10−6
6.95 × 10−5
1.00 × 10−6
4.00 × 10−6
8.00 × 10−6
3.70 × 10−5
6.40 × 10−5
7.00 × 10−7
3.10 × 10−5
1.00 × 10−24
2.07 × 10 4
3.20 × 101
4.1

www.vadosezonejournal.org · Vol. 7, No. 2, May 2008

1.45 × 10−5
1.62 × 10−3
6.95 × 10−5
1.00 × 10−6
4.00 × 10−6
8.00 × 10−6
3.70 × 10−5
6.40 × 10−5
3.22 × 10−3
3.10 × 10−5
3.77 × 10−7
2.07 × 10 4
3.20 × 101
4.2

F . 1. (a) Surface specia on (Fe oxides) as a func on of pH, (b) pH
dependency of the adsorp on of U(VI) in the A horizon on both the
ca on exchange complex (CEC) and the surface complex (SC), and(c) pH dependency of the size of the ca on exchange complex onorganic ma er for the A, E, and Bh1 horizons of the dry Spodosol.
704

drastically. Figure 1b shows the ratio of total adsorbed U to total
U in the system together with the distribution of U between thetwo sorption processes involved (i.e., cation exchange on organic
matter and surface complexation on Fe oxides). Almost all U wasin solution at low pH since both the cation exchange complex
and the surface sites were mostly protonated. While about 50%of the total U was adsorbed around pH 3, almost all U wasadsorbed at pH 4. The relative contribution of U adsorbed bycation exchange compared with surface complexation decreased
with increasing pH. This is consistent with the experimental ﬁndings of McKinley et al. (1995). We found that U adsorption on
the cation exchange complex increased up to pH 4.7 (results notfurther shown). Since the Na concentration in solution increases
with increasing pH owing to the Na addition for the purpose ofobtaining charge balance, Na will increasingly replace the othercations on the exchange complex at higher pH. Note that above
pH 6, U adsorption will decrease again when dissolved CO2 ispresent due to the formation of U(VI)–CO3 complexes (Daviset al., 2004a; Waite et al., 1994). Almost all P was adsorbed in
the studied pH range, with 99.3 and 99.8% P being adsorbed at
pH 2 and 6, respectively (results not shown).
The soil horizons diﬀered in their adsorption capacity, as
illustrated in Fig. 1c and Fig. 2. Figure 1c shows the increasing
capacity of the exchange complex to retain cations other thanprotons with increasing pH. This capacity was lowest for the E
horizon due to the low organic matter content. Figure 2 shows Kd
values deﬁned as the total adsorbed Ca, P, and U concentrations(mol kg−1) across the aqueous Ca, P, and U concentrations (molL−1), respectively, for the A, E, and Bh2 horizons. The horizons
had diﬀerent sorption types and capacities: the A horizon had ahigh capacity of the exchange complex and a low capacity of thesurface complex, the E horizon had low capacities of both theexchange and surface complexes, while the Bh2 horizon had thehighest capacity of the surface complex. The Kd values are shownfor two aqueous-phase compositions: rainwater without the P
fertilizer (having an aqueous U concentration of 8 × 10−9 mol
L−1 for the A and E horizons and 4 x10−9 mol L−1 for the Bh2
horizon representative of the long-term averaged U concentrations during the steady-state ﬂow simulation) and rainwater with
the P fertilizer (Table 5). The pH was controlled by changing
the Na and Cl concentrations. The Kd values of Ca and U wereboth strongly pH dependent, whereas the Kd of P was relatively
independent of pH.
The composition of soil water had a signiﬁcant eﬀect on theKd value. Since the concentrations of Ca, P, and U were muchlarger when rainwater inﬁltrated following P-fertilizer application,competition for the same number of sites on the cation exchange
and surface complexes increased and more elements remained in
the aqueous phase. At pH 4, this resulted in a decrease in Kd ofalmost two orders of magnitude for Ca and P, and a decrease ofabout three orders of magnitude for U. Accordingly, the percentages of U and P adsorbed in the A horizon ranged from 1.5 to
41% for U and from 89.6 to 89.8% for P (results not shown).
Note that almost 100% of U was adsorbed when rainwater in theabsence of P fertilizer inﬁltrated into the soil (Fig. 1b). A slightdecrease in Kd for U occurred between pH 3.5 and 4.5, especiallyin the A and E horizons for rainwater with the P fertilizer. For
the A horizon, the ﬁrst maximum in the percentage of adsorbed
U was around pH 3.4 (33%), with mainly sorption on the cationwww.vadosezonejournal.org · Vol. 7, No. 2, May 2008

exchange complex (results not shown). This is diﬀerent from thecase for rainwater without P fertilizer when already >10% of thetotal amount of U was sorbed onto the surface complex (Fig. 1b).
The Kd reached a minimum at pH 4.7 in the A horizon (Fig. 2a),with only 13% of the total U amount adsorbed on the exchange
and surface complexes (6.5% on both the cation and surface
complexes). From pH 4.7 to 6, the percentage of U adsorbedon the surface complex increased to 42%. The small capacity ofthe surface complex and the high concentration of competing
species (especially P) resulted in a shift in U adsorption on thesurface complex toward higher pH values. This eﬀect was lesspronounced in the Bh horizon because of its higher capacity ofthe surface complex.

F . 2. The linear distribu on coeﬃcient, Kd, of Ca, P, and U as afunc on of pH for rainwater with a U concentra on of 8.0 × 10−9
mol L−1 for the A and E horizons and 4.0 × 10−9 mol L−1 for the Bh2
horizons (lines without dots) and for rainwater with the P fer lizer(lines with dots) for the (a) A, (b) E, and (c) Bh2 horizons.
705

Diﬀerent solution compositions aﬀect the relative
mobility between Ca, P, and U. For rainwater without
P fertilizer, U was less mobile than P at pH values
above 4.3 to 4.5, depending on the soil horizon; however, U was always more mobile when the rainwatercontained P fertilizer. The relative mobility of P and
Ca did not depend only on the water composition,
but also on the soil horizon. In the A horizon witha high exchange capacity and low surface complex
capacity, Ca was less mobile than P at pH values >4.3
for rainwater without P fertilizer and 3.5 for rainwaterwith P fertilizer. In the Bh2 horizon with a relatively
high surface complex capacity, the Kd of P was always
larger that the Kd of Ca.
Calcium, Phosphorus, and Uranium Depth Profiles

Figure 3 shows distributions vs. depth of total
Ca, P, and U per cubic meter of soil at selected times
for both the steady-state and transient water ﬂow
simulations. Total concentrations were obtained bysumming the amounts in the aqueous phase, on the
cation exchange complex (for Ca and U), and on thesurface complex (for Ca, U, and P). Note the diﬀerent F . 3. Concentra on profiles vs. depth of total (a and b) Ca, (c and d) P, and (econcentration scales for P in plots for the ﬁrst 50 and and f) U at (a, c, and e) 10, 25, and 50 yr and (b, d, and f) 100, 150, and 200 yr
the last 150 yr. Calcium accumulated mainly in the for steady-state (thin blue lines) and transient (red thick lines) flow simula ons.
upper horizon during the ﬁrst 50 yr, although more No ce the diﬀerent concentra on scales for P in the top and bo om plots.
accumulated in the Bh1 horizon under transient-ﬂow
conditions. During the next 100 yr, the amount of
Ca in the two Bh horizons gradually increased. This occurred
long-time changes in the amount of Ca, P, and U in the soil
somewhat faster for the transient than the steady-state simulation.
solution due to application of the P fertilizer. Steady-state solute
After 200 yr, however, the amount of Ca in the top 50 cm wastransport conditions are reached when Kd becomes constant forquite similar for the two cases.
the steady-state ﬂow simulation. This condition was reached faster
Phosphorus migrated faster through the A and E horizons
for P and U (at about 40–50 yr) than for Ca (at about 100 yr) at
than did Ca (Fig. 3a vs. 3c) and was more mobile in these two
5-cm depth. This is similar to the conclusion drawn from Fig. 3.
horizons since the amount of Fe2O3 was very small. On the other
Since the total concentrations increased from the initial condition,
hand, Ca transport was retarded by cation exchange reactions in
the Kd at 5-cm depth decreased by almost one order of magnitude
the organic-matter-rich A horizon. The migration of
P through the Bh horizons was slower than that of
Ca due to strong sorption of P on the Fe oxide surface. Phosphorus had not yet reached the BC horizon
after 200 yr. Note that the concentration front for thetransient simulation moved slightly faster than for thesteady-state ﬂow simulation.
Uranium showed the largest diﬀerences between
the two simulations, especially during the ﬁrst 50 yr.
Uranium migrated relatively fast, similar to P, throughthe A and E horizons. As discussed with Fig. 2, U was
more mobile than P when the pH was <4.3, especiallyat relatively high Ca and P concentrations. Much U
accumulated in the Bh horizon during the last 150
yr. Uranium transport was again slightly faster for thetransient simulation than during steady-state ﬂow.
Temporal Varia ons in the Linear Distribu on Coeﬃcient

Figures 4 and 5 show time series of pH and theKd values for Ca, P, and U at depths of 5 (A horizon)and 25 (Bh2 horizon) cm, respectively. A long-time
trend in Kd was observed for both steady-state andtransient ﬂow simulations. This trend is related towww.vadosezonejournal.org · Vol. 7, No. 2, May 2008

F . 4. Time series of (a) pH and of the linear distribu on coeﬃcient, Kd, for (b)Ca, (c) P, and (d) U at a depth of 5 cm for the steady-state (dashed blue line) andtransient (solid black line) flow simula ons during the first 100 yr.
706

for Ca, by one-third for P, and slightly more than one
order of magnitude for U. The increase in competition between elements for the surface sites had a larger
eﬀect on the long-time sorption potential than thelong-time pH increase on the Kd values. A diﬀerent
picture emerged for the Bh2 horizon at a depth of 25
cm (Fig. 5). The increase in pH during the 200-yr simulation caused here an increase in adsorption potential
for U (increasing Kd). The competition with P for thesurface complex was limited since P migrated more
slowly through the Bh horizons than U (see Fig. 3).
After 200 yr, the peak in U occurred at a depth ofabout 25 cm, whereas the peak in P was at 20 cm.
The Kd of P at the 25-cm depth did not show a pronounced long-time trend until 175 yr, which further
reﬂects the limited transport of P to this depth.
The transient flow simulation showed significant short-time variations in pH and Kd for all threeelements. As was discussed in Jacques et al. (2008),
F . 5. Time series of (a) pH and of the linear distribu on coeﬃcient, Kd, for (b)transient weather conditions in precipitation and Ca, (c) P, and (d) U at a depth of 25 cm for the steady-state (dashed blueline) andevaporation resulted in variable water contents and transient (solid black line) flow simula ons during the first 200 yr.
water ﬂuxes. These variations subsequently induced
variations in geochemical conditions. Decreasing
by the Kd time series, steady-state solute transport was obtainedwater contents due to evaporation coincided with decreasing pH
after 40 yr for U and P and after about 100 yr for Ca at the
because of an increase in ionic strength (Jacques et al., 2008). In
bottom of the A horizon. For the transient ﬂow simulation, solute
this study, the decrease in pH caused, in general, the adsorption
ﬂuxes continued to vary by almost three orders of magnitude
potential for Ca and U to decrease. Short-time variations in P
after steady-state solute transport conditions were reached for theconcentrations were limited due to the pH independence of its
steady-state ﬂow simulation. Downward U ﬂuxes for the transient
Kd, as discussed above. Note that the repetition of a 30-yr time
ﬂow simulation were always larger than those for the steady-stateseries of meteorological data (up to 200 yr) may result in some
ﬂow simulation during the ﬁrst 20 yr.
repetitive pattern in the Kd time series, e.g., for P at a depth of
At the bottom of the Bh2 horizon, solute ﬂuxes for the
5 cm (Fig. 4c).
steady-state ﬂow simulation started to increase after 30, 80, andThe seasonal variations in pH and Kd are pronounced close
150 yr for U, Ca, and P, respectively. Uranium ﬂuxes reached ato the soil surface and for early simulation times when P and U
steady-state value after approximately 100 yr. Phosphorus never
are not yet at their long-time steady state (at 5-cm depth andreached the bottom of the soil proﬁle, even after 200 yr. For
for the ﬁrst 40 yr). The Kd values of Ca and U changed up to
the steady-state ﬂow conditions, U started to leach after 120
two orders of magnitude during the ﬁrst 20 yr. Although theyr, whereas Ca ﬂuxes increased 140 yr after an initial decrease.
diﬀerence in total amounts of the three elements between the
Uranium and Ca solute ﬂuxes for the transient ﬂow simulation
steady-state and transient ﬂow simulations was relatively small
were much higher than for the steady-state simulation. Uranium(see Fig. 3), the aqueous concentrations diﬀered dramatically. This
reached the bottom of the soil proﬁle after 80 yr, while Ca ﬂuxes
may have a profound eﬀect on processes that are dependent on
started to increase after about 100 yr. Similarly as for the simuaqueous concentrations, such as nutrient uptake by roots.
lated concentration proﬁles (Fig. 3), the U (and also Ca) ﬂuxes
increased earlier for the transient ﬂow simulation than for the
Fluxes in the Soil and at the Bo om of the Soil Profile
steady-state simulation. This was most obvious at the 100-cm
The solute mass ﬂux density, Js [M L−2 T−1], is the sum of thedepth. Faster leaching of all three elements was due to interacmass ﬂux due to advection of the dissolved solute and the mass
tions caused by short-time variations in water contents, waterﬂux due to hydrodynamic dispersion and molecular diﬀusion. In
ﬂuxes, and geochemical conditions.
terms of the resident ﬂuid concentration, Clr [M L−3], Js may be
written as (Jury and Roth, 1990)
Conclusions
Accurate
investigationof the mobility of contaminants in∂ C lrthe vadose zone requires consideration of all pertinent physicalJ s = −θ D
+ J wC lr[6]∂z
and biogeochemical processes. Coupling of ﬂuid ﬂow, advective–
where Jw is the water ﬂux density [L T−1]. Figure 6 comparesdispersive transport, and multicomponent geochemistry into one
downward ﬂuxes of Ca, P, and U for the steady-state and transient
integrated simulator provides a much-needed tool for studying
ﬂow simulations. Fluxes are shown for the bottom of the A, Bh2,
the interplay between diﬀerent processes, including the various
and C2 horizons at depths of 7, 28, and 100 cm, respectively.
nonlinear feedbacks that occur between the diﬀerent processes.
Although some upward ﬂuxes were calculated at the bottom of
Variable water contents and ﬂuid ﬂuxes during transient ﬂow in
the A horizon during periods of low precipitation or high potenthe vadose zone signiﬁcantly aﬀect prevailing geochemical contial evaporation, these were not included in Fig. 6. As indicatedditions such as soil pH. Consequently, the adsorption potential
www.vadosezonejournal.org · Vol. 7, No. 2, May 2008707

F . 6. Downward
solute fluxes of (a,d, and g) Ca, (b, e,
and f) P, and (c, f,
and i) U (mol d−1m−2) at the bottom of the (a–c) Ahorizon, the (d–f)
Bh2 horizon, and(g–i) the soil profile
for the steady-state(dashed blue line)
and transient (dots)flow simula ons.
Note that the me
scale for the A horizon (a–c) is limited
to 100 yr.

Appendix

and thus the element mobility may change with time. It is similarly important for reactive transport simulators to account forinteractions between diﬀerent species resulting from chemical
reactions such as aqueous speciation and competitive sorption
on the cation exchange or surface complex. The long-time surface application of some elements (e.g., by fertilization) leads
to temporal variations in the distribution coeﬃcient Kd and anonunique relation between Kd and pH.
We illustrated the feedbacks between the diﬀerent processes,
as well as the multiple interactions that may occur between thechemical species present, by using the HP1 reactive transport
simulator to predict the migration of Ca, P, and U following
inorganic P fertilization through a layered soil. The numericalsimulations showed that the distribution coeﬃcients for Ca, P,
and U depend strongly on pH, the composition of the aqueous phase, the mineralogical properties of the soil horizons, andthe amount and type of competing elements. Simulations also
showed that P is a stronger competitive element for U than Ca.
In horizons with large amounts of P (e.g., the A horizon), thesorption potential and Kd of U did not follow the long-term
increase in pH due to the stronger counter eﬀect of P sorption
on the surface complex. On the other hand, when U and Ca
reached a certain horizon (e.g., the Bh2 horizon in our example),while P did not, the Kd for U closely followed the increase in pH.
Furthermore, yearly additions of Ca, P, and U in the form of P
fertilizer resulted in a long-term decrease in Kd of the upper horizons. Short-term temporal variations in the water contents andﬂuxes, when transient ﬂow was considered, caused signiﬁcant pH
variations, which produced very signiﬁcant variations in the Kdvalues for Ca and U by up to two orders of magnitude. Leachingof Ca and U from the soil proﬁle occurred faster in the transient
ﬂow simulation than in the steady-state simulation as a resultof the interplay between changing hydrologic and geochemicalsoil conditions.
www.vadosezonejournal.org · Vol. 7, No. 2, May 2008

ai
Am
j

M ip

activity of secondary species i (dimensionless)chemical formula of aqueous master species j(dimensionless)chemical formula of aqueous secondary species i(dimensionless)reactive surface area [L2]concentration of ith aqueous species (mol kg−1
H2O)
total concentration of jth component (mol kg−1
H2O)
volumetric heat capacity of porous medium [M L−1
T−2 K−1]
total concentration in the sink term (mol kg −1
H2O)
volumetric heat capacity of the liquid phase [M L−1
T−2 K−1]
dispersion coeﬃcient in the liquid phase [L2 T−1]
apparent activation energy (kJ mol−1)
Faraday constant, 96,485 C mol−1
Gibbs free energy (kJ mol−1)
water pressure head [L]ionic strength (mol kg−1 H2O)
intrinsic rate constant [mol L−2 T−1]
unsaturated hydraulic conductivity [L T−1]
equilibrium constant: superscripts x = l, e, s, or p arefor aqueous, exchange, surface species, and minerals, respectively; subscript y = ix is an index for
speciﬁc species (dimensionless)Monod (x = M) or inhibition (x = I) constant in
the ith term of degradation reaction jk (mol kg−1
H2O)
chemical formula of mineral ip

ni
Nk
Nm

order of the reaction of secondary species i
number of kinetic reactions
number of aqueous components or master species

Ai
Amin
ci
Cj
Cp
Cr,j
Cw
Dw
Ea
F
G
h
I
k0
K
Kyx

K jxk , i708

N M, jk

number of inhibition terms in degradation reaction
jk
number of Monod terms in degradation reaction jk

q
Q
R
R jk

volumetric ﬂuid ﬂux density [L T−1]
ion activity product (dimensionless)universal gas constant (kJ mol−1 K−1)
rate of degradation reaction jk [mol kg−1 H2O T−1]

Ro,j
S
Sm
js

source–sink term [L3 kg−1 H2O T−1]
sink term for root water uptake [mol T−1]
chemical formula of surface master species js

S ischemical formula of secondary surface species istTv max, jkXm
jetime [T]temperature (K)maximum speciﬁc rate of substrate utilization indegradation reaction jk [mol kg−1 H2O T−1 (molL−1)−1]
spatial coordinate (positive upward) [L]chemical formula of exchange master species je

X iechemical formula of secondary exchange species ie

Xmo,rconcentration of the microbial population r [M kg−1
H2O]
charge on surface species is

N I, jkxz is
αβ eie , je
β iss , js
γjm
γyx
λ
θ
x
νy,z

Ψ js

Appelo, C.A.J., and D. Postma. 2005. Geochemistry, groundwater, and pollution. 2nd ed. A.A. Balkema, Rotterdam, the Netherlands.
Appelo, C.A.J., E. Verweij, and H. Schäfer. 1998. A hydrogeochemical transportmodel for an oxidation experiment with pyrite/calcite/exchangers/organic
matter containing sand. Appl. Geochem. 13:257–268.
Ball, J.W., and D.K. Nordstrom. 1991. WATEQ4F: User’s manual with revised
thermodynamic data base and test cases for calculating speciation of major,
trace and redox elements in natural waters. Open-File Rep. 90-129. USGS,
Denver, CO.
Barber, S.A. 1995. Soil nutrient bioavailability: A mechanistic approach. 2nd ed.
John Wiley & Sons, New York.
Barišić, D., S. Lulić, and P. Miletić. 1992. Radium and uranium in phosphate
fertilizers and their impact on the radioactivity of waters. Water Res.
26:607–611.
Barnett, M.O., P.M. Jardine, and S.C. Brooks. 2002. U(VI) adsorption toheterogeneous subsurface media: Application of a surface complexation
model. Environ. Sci. Technol. 36:937–942.
Barnett, M.O., P.M. Jardine, S.C. Brooks, and H.M. Selim. 2000. Adsorption
and transport of uranium(VI) in subsurface media. Soil Sci. Soc. Am. J.
64:908–917.
Berner, R.A., J.-L. Rao, S. Chang, R. O’Brien, and K. Keller. 1998. Seasonalvariability of adsorption and exchange equilibria in soil waters. Aquat.
Geochem. 4:273–290.
Bethke, C.M., and P.V. Brady. 2000. How the Kd approach undermines groundwater cleanup. Ground Water 38:435–443.
Bradford, S.A., J. Šimůnek, M. Bettahar, M.Th. van Genuchten, and S.R. Yates.
2003. Modeling colloid attachment, straining, and exclusion in saturated
porous media. Environ. Sci. Technol. 37:2242–2250.
Bradford, S.A., J. Šimůnek, M. Bettahar, M.Th. van Genuchten, and S.R. Yates.
2006. Signiﬁcance of straining in colloid deposition: Evidence and implications. Water Resour. Res. 42:W12S26, doi:10.1029/2005WR004791.
Brady, N.C. 1990. The nature and properties of soils. 10th ed. Macmillan,
New York.
Brantley, S. 2003. Reaction kinetics of primary rock-forming minerals underambient conditions. p. 73–118. In J.I. Drever (ed.) Fresh water geochemistry, weathering, and soils. Treatise on geochemistry. Vol 5. Pergamon Press,
Oxford, UK.
Brooks, R.H., and A.T. Corey. 1966. Properties of porous media aﬀecting ﬂuid
ﬂow. J. Irrig. Drain. Div. Am. Soc. Civ. Eng. 92:61–88.
Bruggenwert, M.G.M., and A. Kamphorst. 1982. Survey of experimental information on cation exchange in soil systems. p. 141–203. In G.H. Bolt (ed.)
Soil chemistry: B. Physico-chemical models. Elsevier, Amsterdam.
Calba, H., Firdaus, P. Cazevieille, C. Thée, R. Poss, and B. Jaillard. 2004. The
dynamics of protons, aluminium, and calcium in the rhizosphere of maize
cultivated in tropical acid soils: Experimental study and modelling. Plant
Soil 260:33–46.
Campbell, G.S. 1985. Soil physics with BASIC. Elsevier, Amsterdam.
Carrillo-González, R., J. Šimůnek, S. Sauvé, and D. Adriano. 2006. Mechanisms
and pathways of trace elements in soils. Adv. Agron. 91:111–178.
Cogné, F. 1993. L’énergie nucléaire dans notre environnement radioactive: Principales données chiﬀreés. Rev. Gen. Nucl. 2:97–103.
Cunningham, J.A., and I. Mendoza-Sanchez. 2006. Equivalence of two models
for biodegradation during contaminant transport in groundwater. Water
Resour. Res. 42:W02416, doi:10.1029/2005WR004205.
Dang, Y.P., R.C. Dalal, D.G. Edwards, and K.G. Tiller. 1994. Kinetics of zincdesorption from Vertisols. Soil Sci. Soc. Am. J. 58:1392–1399.
Davis, J.A., D.E. Meece, M. Kohler, and G.P. Curtis. 2004a. Approaches to surface complexation modeling of uranium(VI) adsorption on aquifer sediments. Geochim. Cosmochim. Acta 68:3621–3641.
Davis, J.A., S.B. Yabusaki, C.I. Steefel, J.M. Zachara, G.P. Curtis, G.B. Redden, L.J.
Criscenti, and B.D. Honeyman. 2004b. Assessing conceptual models for subsurface reactive transport of inorganic contaminants. EOS 85:449, 455.
de Mars, H., and G. Garritsen. 1997. Interrelationship between water quality
and groundwater ﬂow dynamics in a small wetland system along a sandy
hill ridge. Hydrol. Processes 11:335–351.
de Mars, H., M.J. Wassen, and H. Olde Venterink. 1997. Flooding and groundwater dynamics in fens in eastern Poland. J. Veg. Sci. 8:319–328.
DeNovio, N.M., J.E. Saiers, and J.N. Ryan. 2004. Colloid movement in unsaturated porous media: Recent advances and future directions. Vadose Zone
J. 3:338–351.
Dzombak, D.A., and F.M.M. Morel. 1990. Surface complexation modeling:

angle between ﬂow direction and vertical axis
equivalent fraction of exchange species i e onexchanger je (dimensionless)mole fraction of surface species is on surface js
activity coeﬃcient for the jth aqueous master species(kg H2O mol−1)
activity coeﬃcient for species y: x = l or e for aqueous or exchange species; y = ix is the index for thesecondary species (kg H2O mol−1)
apparent thermal conductivity of the soil [M L T−3
K−1]
volumetric water content [L3 L−3]
stoichiometric coeﬃcient of master species y in thereaction for secondary species z: superscripts x =l, e, s, p, and k are for aqueous species, exchange
species, surface species, minerals, and kinetic reaction, respectively; subscripts y = jx is an index forthe master species or kinetic reaction; z = ix is anindex for the secondary species (dimensionless)surface potential for surface js [V]A
This work was partly supported by the National Science Foundation, Agreement no. EAR-9876800 and the Terrestrial Sciences Program of the Army Research Oﬃce (Terrestrial Processes and Landscape
Dynamics and Terrestrial System Modeling and Model Integration).
Additional support was obtained through the bilateral agreement project “Development and evaluation of a coupled geochemical transportmodel” 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
Project CO91002.

References
Aagaard, P., and J.C. Helgeson. 1982. Thermodynamic and kinetic constraints
on reaction rates among minerals and aqueous solutions: I. Theoreticalconsiderations. Am. J. Sci. 282:237–285.
www.vadosezonejournal.org · Vol. 7, No. 2, May 2008709

Lasaga, A.C. 1981. Transition state theory. p. 169–195. In A.C. Lasaga et al.
(ed.) Kinetics of geochemical process. Rev. Mineral. 8. Min. Soc. Am.,
Washington, DC.
Lasaga, A.C. 1995. Fundamental approaches in describing mineral dissolution and precipitation rates. p. 23–86. In A.F. White et al. (ed.) Chemical weathering rates of silicate minerals. Rev. Mineral. 31. Min. Soc. Am.,
Washington, DC.
Lasaga, A.C. 1998. Kinetic theory in the earth sciences. Princeton Univ. Press,
Princeton, NJ.
Le Gallo, Y., O. Bildstein, and E. Brosse. 1998. Coupled reaction–ﬂow modelling of diagenetic changes in reservoir permeability, porosity and mineralcompositions. J. Hydrol. 209:366–388.
Lenhart, J.J., S.E. Cabaniss, P. MacCarthy, and B.D. Honeyman. 2000.
Uranium(VI) complexation with citric, humic and fulvic acids. Radiochim. Acta 88:345–353.
Lenhart, J.J., and B.D. Honeyman. 1999. Uranium(VI) sorption to hematite in
the presence of humic acid. Geochim. Cosmochim. Acta 63:2891–2901.
Lichtner, P.C. 1996. Continuum formulation of multicomponent–multiphasereactive transport. p. 1–81. In P.C. Lichtner et al. (ed.) Reactive transportin porous media. Rev. Mineral. 34. Min. Soc. Am., Washington, DC.
Lichtner, P.C., S. Yabusaki, K. Pruess, and C.I. Steefel. 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.
Mayer, U. 1999. A numerical model for multicomponent reactive transport invariably saturated porous media. Ph.D. diss. Dep. of Earth Sciences, Univ.
of Waterloo, Waterloo, ON, Canada.
Mayer, K.U., E.O. Frind, and D.W. Blowes. 2002. Multicomponent reactive
transport modeling in variably saturated porous media using a generalized formulation for kinetically controlled reactions. Water Resour. Res.
38:1174, doi:10.1029/2001WR000862.
McKinley, J.P., J.M. Zachara, S.C. Smith, and G.D. Turner. 1995. The inﬂuenceof uranyl hydrolysis and multiple site-binding reactions on adsorption of
U(VI) to montmorillonite. Clays Clay Miner. 43:586–598.
Metz, V., W. Pﬁngsten, J. Lützenkirchen, and W. Schüssler. 2002. TrePro 2002:
Modelling of coupled transport reaction processes. Worksh. of the Forschungszentrum Karlsruhe, 3rd, Karlsruhe, Germany. 20–21 Mar. 2002.
Available at www.fzk.de/fzk/groups/ine/documents/internetdokument/id_040291.pdf (verified 20 Mar. 2008). Forschungszentrum
Karlsruhe, Karlsruhe.
MIRA Team. 2004. Milieu: en natuurrapport Vlaanderen. Achtergronddocument 2004 Vermesting. Available at www.milieurapport.be (veriﬁed 17
Apr. 2008). Vlaamse Milieumaatschappij (Flemish Environment Agency),
Mechelen, Belgium.
Morel, F., and J. Hering. 1993. Principles and applications of aquatic chemistry.
John Wiley & Sons, New York.
Morrison, S.J., R.R. Spangler, and V.S. Tripathi. 1995. Adsorption of
uranium(VI) on amorphous ferric oxyhydroxide at high concentrations of
dissolved carbon(IV) and sulfur(VI). J. Contam. Hydrol. 17:333–346.
Nowack, B., K.U. Mayer, S.E. Oswald, W. Van Beinum, C.A.J. Appelo, D.
Jacques, P. Seuntjens, F. Gérard, B. Jaillard, A. Schnepf, and T. Roose.
2006. Veriﬁcation and intercomparison of reactive transport codes to describe root-uptake. Plant Soil 285:305–321.
Nuclear Regulatory Commission. 2006. Proc. Int. Conf. on Conceptual Model
Dev. for Subsurface Reactive Transport Modeling of Inorganic Contaminants, Radionuclides, and Nutrients. NUREG/CP-0193. Available at www.
nrc.gov/reading-rm/doc-collections/nuregs/conference/cp0193/cp0193.
pdf (veriﬁed 17 Apr. 2008). Oﬃce of Nuclear Regulatory Research, U.S.
Nuclear Regulatory Commission, Washington, DC.
Nützmann, G., P. Viotti, and P. Aagaard. 2005. Reactive transport in soil andgroundwater: Processes and models. Springer-Verlag, Berlin, Germany.
Öztürk, H.S., and L. Özkan. 2004. Eﬀects of evaporation and diﬀerent ﬂow
regimes on solute distribution in soil. Transp. Porous Media 56:245–255.
Papastefanou, C., S. Stoulos, A. Iaonnidou, and M. Manolopoulou. 2006. The
application of phosphogypsum in agriculture and the radiological impact.
J. Environ. Radioact. 89:188–198.
Parkhurst, D.L., and C.A.J. Appelo. 1999. User’s guide to PHREEQC (Version
2): A computer program for speciation, batch-reaction, one-dimensionaltransport, and inverse geochemical calculations. Water-Resour. Invest. Rep.
99-4259. USGS, Denver, CO.
Penmann, H.L. 1948. Natural evaporation from open water, bare soils and grass.
Proc. R. Soc. London Ser. A 190:120–145.

Hydrous ferric oxide. John Wiley & Sons, New York.
European Commission. 1999. Radiation protection 95: Reference levels forworkplaces processing materials with enhanced levels of naturally occurring radionuclides. A guide to assist implementation of Title VII of theEuropean Basic Safety Standards Directive (BSS) concerning natural radiation sources. Eur. Commission, Luxembourg.
Finke, P. 2006. Modelling the genesis of Luvisols in late Weichsel loess. p. 29–30.
In A. Samouëlian and S. Cornu (ed.) Proc. Worksh. on Modelling of Pedogenesis, Nancy, France. 2–4 Oct. 2006. INRA, Orléans.
Freedman, V.L., D.H. Bacon, K.P. Saripalli, and P.D. Meyer. 2004. A ﬁlm depositional model of permeability for mineral reactions in unsaturated media.
Vadose Zone J. 3:1414–1424.
Gargiulo, G., S. Bradford, J. Šimůnek, P. Ustohal, H. Vereecken, and E. Klumpp.
2007. Bacterial transport and deposition under unsaturated conditions:
The role of the matrix grain size and the bacteria surface protein. J. Contam. Hydrol. 92:255–273.
Geelhoed, J.S., W.H. Van Riemsdijk, and G.R. Findenegg. 1999. Simulation of
the eﬀect of citrate exudation from roots on the plant availability of phosphate adsorbed on goethite. Eur. J. Soil Sci. 50:379–390.
Goldberg, S., L.J. Criscenti, D.R. Turner, J.A. Davis, and K.J. Cantrell. 2007.
Adsorption–desorption processes in subsurface reactive transport modeling. Vadose Zone J. 6:407–435.
Gonçalves, M.C., J. Šimůnek, T.B. Ramos, J.C. Martins, M.J. Neves, andF.P. Pires. 2006. Multicomponent solute transport in soil lysimeters irrigated with waters of diﬀerent quality. Water Resour. Res. 42:W08401,
doi:10.1029/2005WR004802.
Greskowiak, J., H. Prommer, G. Massmann, C.D. Johnston, G. Nützmann,and A. Pekdeger. 2005. The impact of variably saturated conditions onhydrogeochemical changes during artiﬁcial recharge of groundwater. Appl.
Geochem. 20:1409–1426.
Guzman, E.T.T., M.V.E. Alberich, and E.O. Regil. 2002. Uranium and phosphate behavior in the vadose zone of a fertilized corn ﬁeld. J. Radioanal.
Nucl. Chem. 254:509–517.
Hopmans, J.W., and K.L. Bristow. 2002. Current capabilities and future needsof root water and nutrient uptake modelling. Adv. Agron. 77:103–183.
Jacques, D., and J. Šimůnek. 2005. User manual of the multicomponent variably
saturated transport model HP1: Description, veriﬁcation and examples.
Version 1.0. BLG-998. SCK⋅CEN, Mol, Belgium.
Jacques, D., J. Šimůnek, D. Mallants, and M.Th. van Genuchten. 2006a. Modelling uranium leaching from agricultural soils to groundwater as a criterion for comparison with complementary safety indications. Mater. Res.
Soc. Symp. Proc. 932:1057–1064.
Jacques, D., J. Šimůnek, D. Mallants, and M.Th. van Genuchten. 2006b. Operator-splitting errors in coupled reactive transport codes for transient variably saturated ﬂow and contaminant transport in layered soil proﬁles. J.
Contam. Hydrol. 88:197–218.
Jacques, D., J. Šimůnek, D. Mallants, and M.Th. van Genuchten. 2008. Modelling coupled water ﬂow, solute transport and geochemical reactionsaﬀecting heavy metal migration in a podzol soil. Geoderma (in press),doi:10.1016/j.geoderma.2008.01.009.
Jansen, B., K.G.J. Nierop, and J.M. Verstraten. 2002. Inﬂuence of pH andmetal/carbon ratios on soluble organic complexation of Fe(II), Fe(III) andAl in soil solutions determined by diﬀusive gradients in thin ﬁlms. Anal.
Chim. Acta 454:259–270.
Jin, Y., M.V. Yates, and S.R. Yates. 2002. Microbial transport. p. 1481–1509. In
J.H. Dane and G.C. Topp (ed.) Methods of soil analysis. Part 4. Physicalmethods. SSSA Book Ser. 5. SSSA, Madison, WI.
Joris, I. 2005. Soil water dynamics in alluvial wetlands: Field and modellingstudy. Ph.D. diss. Univ. of Leuven, Leuven, Belgium.
Jury, W.A., and K. Roth. 1990. Transfer functions and solute movement throughsoil: Theory and applications. Birkhäuser, Berlin.
Kohler, M., G.P. Curtis, D.B. Kent, and J.A. Davis. 1996. Experimental investigation and modeling of uranium(VI) transport under variable chemical
conditions. Water Resour. Res. 32:3539–3551.
Kuechler, R., and K. Noack. 2002. Transport of reacting solutes through theunsaturated zone. Transp. Porous Media 49:361–375.
Kuechler, R., and K. Noack. 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.
www.vadosezonejournal.org · Vol. 7, No. 2, May 2008

710

Rockhold, M.L., R.R. Yarwood, M.R. Niemet, P.J. Bottomley, and J.S. Selker. 2002.
Considerations for modeling bacterial-induced changes in hydraulic propertiesof variably saturated porous media. Adv. Water Resour. 25:477–495.
Rockhold, M.L., R.R. Yarwood, M.R. Niemet, P.J. Bottomley, and J.S. Selker.
2005. Experimental observations and numerical modeling of coupled microbial and transport processes in variably saturated sand. Vadose Zone J.
4:407–417.
Rockhold, M.L., R.R. Yarwood, and J.S. Selker. 2004. Coupled microbial andtransport processes in soils. Vadose Zone J. 3:368–383.
Rothbaum, H.P., D.A. McGaveston, T. Wall, A.E. Johnston, and G.E.G. Mattingly. 1979. Uranium accumulation in soils from long-continued applications of superphosphate. J. Soil Sci. 30:147–153.
Saaltink, M.W., F. Battle, C. Ayora, J. Carrera, and S. Olivella. 2004. RETRASO,
a code for reactive transport in saturated and unsaturated porous media.
Geol. Acta 2:235–251.
Saito, H., J. Šimůnek, and B.-P. Mohanty. 2006. Numerical analysis of coupled water, vapor, and heat transport in the vadose zone. Vadose Zone J.
5:784–800.
Salvage, K.M., and G.T. Yeh. 1998. Development and application of a numerical model of kinetic and equilibrium microbiological and geochemical reactions (BIOKEMOD). J. Hydrol. 209:27–52.
Samper, J., R. Juncosa, J. Delgado, and L. Montenegro. 2000. CORE2D: A codefor nonisothermal water ﬂow and reactive solute transport. User manual
Version 2. Tech. Publ. 06/2000. Empresa Nacional de Residuos Radioactivos, Madrid.
Santos, A.J.G., B.P. Mazzilli, D.I.T. Favoro, and P.S.C. Silva. 2006. Partitioning of radionuclides and trace elements in phosphogypsum and its source
materials based on sequential extraction methods. J. Environ. Radioact.
87:52–61.
Saripalli, K.P., V.L. Freedman, B.P. McGrail, and P.D. Meyer. 2006. Characterization of the speciﬁc solid–water interfacial area–water saturation relationship and its import to reactive transport. Vadose Zone J. 5:777–783.
Saueia, C.H.R., and B.P. Mazzilli. 2006. Distribution of natural radionuclides
in the production and use of phosphate fertilizers in Brazil. J. Environ.
Radioact. 89:229–239.
Schäfer, D., W. Schäfer, and W. Kinzelbach. 1998. Simulation of reactive processes related to biodegradation in aquifers. 1. Structure of the three-dimensional reactive transport model. J. Contam. Hydol. 31:167–186.
Seuntjens, P. 2000. Reactive solute transport in heterogeneous porous media:
Cadmium leaching in acid sandy soils. Ph.D. diss. Univ. of Antwerp, Antwerp, Belgium.
Seuntjens, P., D. Mallants, N. Toride, C. Cornelis, and P. Geuzens. 2001a. Gridlysimeter study of steady-state chloride transport in two Spodosol types
using TDR and wick samplers. J. Contam. Hydrol. 51:13–39.
Seuntjens, P., B. Nowack, and R. Schulin. 2005. Root-zone modeling of heavy
metal uptake and leaching in the presence of organic ligands. Plant Soil
265:61–73.
Šimůnek, J., C. He, L. Pang, and S.A. Bradford. 2006. Colloid-facilitated solutetransport in variably saturated porous media: Numerical model and experimental veriﬁcation. Vadose Zone J. 5:1035–1047.
Šimůnek, J., N.J. Jarvis, M.Th. van Genuchten, and A. Gärdenäs. 2003. Review
and comparison of models for describing non-equilibrium and preferential
ﬂow and transport in the vadose zone. J. Hydrol. 272:14–35.
Šimůnek, J., and D.L. Suarez. 1993. Modeling of carbon dioxide transport andproduction in soil: 1. Model development. Water Resour. Res. 29:487–497.
Šimůnek, J., and D.L. Suarez. 1994. Two-dimensional transport model for variably saturated porous media with major ion chemistry. Water Resour. Res.
30:1115–1133.
Šimůnek, J., and M.Th. van Genuchten. 2008. Modeling nonequilibrium ﬂow andtransport processes using HYDRUS. Vadose Zone J. 7:782–797 (this issue).
Šimůnek, J., M.Th. van Genuchten, and M. Šejna. 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 Softw. Ser. 1. Dep. of Environmental Sciences, Univ. of California,
Riverside.
Slattery, W.J., and G.R. Ronnfeldt. 1992. Seasonal variation of pH, aluminium
and manganese in acid soils from north-eastern Victoria. Aust. J. Exp. Agric. 32:1105–1112.
Soler, J.M., and A.C. Lasaga. 1998. An advection–dispersion–reaction model ofbauxite formation. J. Hydrol. 209:311–330.

www.vadosezonejournal.org · Vol. 7, No. 2, May 2008

Spalding, R.F., and W.M. Sackett. 1972. Uranium in runoﬀ from the Gulf
of Mexico distributive province: Anomalous concentrations. Science
175:629–631.
Steefel, C.I., D.J. DePaolo, and P.C. Lichtner. 2005. Reactive transport modeling: An essential tool and new research approach for the Earth sciences.
Earth Planet. Sci. Lett. 240:539–558.
Stolk, A.P. 2001. Landelijk Meetnet Regenwatersamenstelling, Meetresultaten
1999. RIVM Rep. 723101056. Natl. Inst. for Public Health and the Environ. (RIVM), Bilthoven, the Netherlands.
Suarez, D.L., and J. Šimůnek. 1993. Modeling of carbon dioxide transport andproduction in soil: 2. Parameter selection, sensitivity analysis, and comparison of model predictions to ﬁeld data. Water Resour. Res. 29:499–513.
Taylor, S.W., and P.R. Jaﬀé. 1990. Substrate and biomass transport in a porous
medium. Water Resour. Res. 26:2181–2194.
Tebes-Stevens, C., A.J. Valocchi, J.M. VanBriesen, and B.E. Rittman. 1998. Multicomponent transport with coupled geochemical and microbiological reactions: Model description and example simulations. J. Hydrol. 209:8–26.
Tinker, P.B., and P.H. Nye. 2000. Solute movement in the rhizosphere. Oxford
Univ. Press, New York.
Tipping, E. 2002. Cation binding by humic substances. Cambridge Univ. Press,
Cambridge, UK.
van der Lee, J., and L. De Windt. 2001. Present state and future directions ofmodeling of geochemistry in hydrogeological systems. J. Contam. Hydrol.
47:265–282.
van der Lee, J., L. De Windt, V. Lagneau, and P. Goblet. 2003. Module-oriented
modeling of reactive transport with HYTEC. Comput. Geosci. 29:265–275.
van Genuchten, M.Th. 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892–898.
van Genuchten, M.Th. 1987. A numerical model for water and solute movement in and below the root zone. Res. Rep. 121. U.S. Salinity Lab., Riverside, CA.
van Genuchten, M.Th., and J. Šimůnek. 2004. Integrated modeling of vadosezone ﬂow and transport processes. p. 37–69. In R.A. Feddes et al. (ed.)
Unsaturated-zone modelling: Progress, challenges and applications. Frontis, Wageningen, the Netherlands.
Voegelin, A., V. Vulava, and R. Kretzschmar. 2001. Reaction-based model describing competitive sorption and transport of Cd, Zn, and Ni in an acidic
soil. Environ. Sci. Technol. 35:1651–1657.
Waite, T.D., J.A. Davis, T.E. Payne, G.A. Waychunas, and N. Xu. 1994.
Uranium(VI) adsorption to ferrihydrite: Application of a surface complexation model. Geochim. Cosmochim. Acta 58:5465–5478.
Wan, J.M., and T.K. Tokunaga. 1997. Film straining of colloids in unsaturated
porous media: Conceptual model and experimental testing. Environ. Sci.
Technol. 31:2413–2420.
Wan, J.M., and T.K. Tokunaga. 2002. Partitioning of clay colloids at air–waterinterfaces. J. Colloid Interface Sci. 247:54–61.
Yarwood, R.R., M.L. Rockhold, M.R. Niemet, J.S. Selker, and P.J. Bottomley.
2006. Impact of microbial growth on water ﬂow and solute transport inunsaturated porous media. Water Resour. Res. 42:W10405, doi:10.1029/
2005WR004550.
Yeh, G.-T., and H.-P. Cheng. 1999. 3DHYDROGEOCHEM: A 3-dimensional model of density-dependent subsurface ﬂow and thermal multispecies–multicomponent hydrogeochemical transport. EPA/600/R-98/159.
USEPA, Ada, OK.
Zhang, F., R. Zhang, and S. Kang. 2003. Estimating temperature eﬀects onwater ﬂow in variably saturated soils using activation energy. Soil Sci. Soc.
Am. J. 67:1327–1333.
Zhu, C. 2003. A case against Kd–based transport models: Natural attenuation ata mill tailings site. Comput. Geosci. 29:351–359.
Zielinski, R.A., S. Asher-Bolinder, A.L. Meier, C.A. Johnson, and B.J. Szabo.
1997. Natural or fertilizer-derived uranium in irrigation drainage: A casestudy in southeastern Colorado, U.S.A. Appl. Geochem. 12:9–21.

711