Chapter 16

The Use of Multicomponent Solute Transport
Models in Environmental Analyses
Jiřı́ Šimůnek, Diederik Jacques, Tiago B. Ramos, and Bertrand Leterme

Abstract This chapter provides a brief overview of multicomponent solute transport models, which simulate the subsurface transport of multiple ions that may
mutually interact, can create various complex species, can compete with each otherfor sorption sites, and/or can precipitate or dissolve. These models are broadly
divided into two major groups, those with specific chemistry and general models,
and typical examples of these models are given. More detail is provided for theUnsatChem and HPx modules of HYDRUS. The applicability of the UnsatChem
module to simulate the multicomponent transport of major ions is demonstratedusing the data from a field experiment involving irrigation with waters of different
quality, carried out in Portugal. The second example illustrates the use of HP1 tosimulate the fate of mercury in a contaminated soil. Finally, the third example
demonstrates the versatility of HP2 to simulate the release and migration of
uranium from a simplified uranium mill tailings pile toward a river. Using thesethree examples, it is demonstrated that in spite of the considerable demand on input
data, the multicomponent solute transport models can be effective and versatile
tools for evaluating complex agricultural and environmental problems.

J. Šimůnek (*)
Department of Environmental Sciences, University of California, Riverside, CA 92521, USAe-mail: Jiri.Simunek@ucr.edu
D. Jacques • B. LetermeInstitute for Environment, Health and Safety, Belgian Nuclear Research Centre (SCK•CEN),
Boeretang 200, Mol 2400, Belgium
e-mail: djacques@sckcen.be
T.B. Ramos
CEER-Biosystems Engineering, Institute of Agronomy, Technical University of Lisbon,
Tapada da Ajuda, 1349-017 Lisbon, Portugal
e-mail: tiago_ramos@netcabo.pt
W.G. Teixeira et al. (eds.), Application of Soil Physics in Environmental
377
Analyses: Measuring, Modelling and Data Integration, Progress in Soil Science,
DOI 10.1007/978-3-319-06013-2_16, © Springer International Publishing Switzerland 2014

378

J. Šimůnek et al.

Keywords Variably-saturated • Water flow • Soil • Solute transport• Multicomponent models • Hydrus • HP1 • UnsatChem • Uranium transport• Salinity • SAR • Tailings pile

16.1

Introduction

In Chap. 15 of this volume, van Genuchten et al. provided a detailed description ofvarious approaches to modeling water flow and solute transport in soils and inunsaturated porous media in general. The focus of that chapter was on modeling thetransport of only one chemical species while assuming that the behavior of thissolute is independent of that of other species that may be present in the soil solution.
The interaction of this species with the environment (soil) is usually describedassuming linear reversible sorption (the Kd approach) and zero- or first-order
reactions (production and/or decay). The advantage of this approach is that alarge number of analytical solutions exist for one-, two-, and three-dimensional
solute transport problems with different types of boundary conditions and steadystate water flow conditions (e.g., van Genuchten 1981; Yeh 1981; Toride
et al. 1993; Leij and Bradford 1994; Šimůnek et al. 1999). These analytical models
are useful for a variety of applications, such as for “providing initial or approximate
analyses of alternative pollution scenarios, conducting sensitivity analyses toinvestigate the effects of various parameters or processes on contaminant transport,
extrapolating results over large times and spatial scales where numerical solutions
become impractical, serving as screening models, estimating transport parametersfrom laboratory or well-defined field experiments, providing benchmark solutions
for more complex transport processes that cannot be solved analytically, and forvalidating more comprehensive numerical solutions of the governing transport
equations” (van Genuchten et al. 2012).
A large number of numerical models (e.g., Šimůnek 2005; Šimůnek et al. 2008;
van Genuchten et al. 2014) exists for more complex solute transport situations,such as those involving one- or multi-dimensional transient water flow, heterogeneous systems, or nonequilibrium solute transport with nonlinear reactions, forwhich analytical solutions are generally not available and/or cannot be derived.
However, even such models represent a great simplification of real solute transport
problems, since these models usually consider the transport of a single solute
species, while the soil solution is always a mixture of many ions. These ions
may mutually interact, can create various complex species, can compete witheach other for sorption sites, and/or can precipitate or dissolve, depending on actual
conditions in the soil profile (e.g., Šimůnek and Valocchi 2002; Šimůnek and van
Genuchten 2006; Seaman et al. 2012). The transport codes accounting for thesemore complex reactions, usually called the hydrogeobiochemical transport models,offer the possibility to include processes such as aqueous complexation, cation andanion exchange, surface complexation, and/or equilibrium or kinetic precipitation/
dissolution of pure solid phases or solid-solutions. In addition, they may also

16

The Use of Multicomponent Solute Transport Models in Environmental Analyses

379

consider complex parallel and/or sequential geochemical and/or microbiologicalkinetic networks, including soil carbon cycles, which may involve more advanced
kinetic rate equations than simple zero-order or first-order reactions (e.g., Regnier
et al. 2005; Thullner et al. 2005; Centler et al. 2010; Langergraber and Šimůnek
2012; Jacques et al. 2013). There are many important agricultural and environmental problems that require models capable of analyzing transport and reactions ofmultiple chemical species involved in such mutual interactions. In this chapter, we
briefly review multicomponent solute transport models, give few typical examples(e.g., BioGeoChem, MIN3P, CrunchFlow, UnsatChem, HP1/2/3), and show a few
examples of their use in agricultural (i.e., modeling of soil salinity and the fate of
mercury in contaminated soil using HP1) and environmental (i.e., leaching of
uranium from U tailings pile) applications.

16.2

HydroGeoBioChemical Transport Models

Reviews of the development over the last two decades of hydrogeochemicaltransport models involving reactive multiple components are given by Mangold
and Chin-Fu Tsang (1991), Lichtner (1996), Steefel and MacQuarrie (1996),
Šimůnek and Valocchi (2002), Bell and Binning (2004), Steefel et al. (2005),
Seaman et al. (2012), and Šimůnek et al. (2013). Most coupled reactive transport
codes were developed for steady-state flow conditions for one, two, or three
dimensions, while the flow field (water contents and water fluxes) has to be
specified at the input or calculated outside of the reactive transport code using
other flow and transport numerical models and imported (Parkhurst and Appelo
1999; Steefel 2009; Parkhurst et al. 2004). Several multicomponent transportmodels have been published recently that also consider variably-saturated flow
(Zyvoloski et al. 1997; Mayer et al. 2002; Jacques et al. 2006, 2008a; Šimůnek
et al. 2013).

16.2.1 Models with Specific Chemistry
Geochemical models can be broadly divided into two major groups: those withspecific chemistry and general models (Šimůnek and Valocchi 2002). Models withspecific chemistry are usually constrained to very specific applications since they
are restricted to certain prescribed chemical systems. They are, however, mucheasier to use since the chemical system is already predefined by model developers
and can be more computationally efficient than general models because thenumerics can be optimized for this particular chemical system. Typical examples
of models with specified chemistry are those simulating the transport of major ions,such as LEACHM (Hutson and Wagenet 1992) and UnsatChem (Šimůnek andSuarez 1994, 1997), and various reclamation models (Šimůnek and Valocchi 2002).

380

J. Šimůnek et al.

These models typically consider the transport of major ions and their mutual
reactions such as complexation, cation exchange, and precipitation/dissolution.
The UnsatChem model, as an add-on module to HYDRUS-1D and HYDRUS(2D/3D) (Šimůnek et al. 2008), and its application to a field study will be discussed
below.
Another group of models with specified chemistry are those simulating carbon
and nitrogen cycles. Some examples are CENTURY (Parton et al. 1987),
LEACHM (Hutson and Wagenet 1992), RZWQM (Ahuja and Hebson 1992),
COUP (Jansson and Karlberg 2001), and the Wetland module of HYDRUS(2D/3D) (Langergraber and Šimůnek 2005, 2012). These models typically distribute organic matter, carbon, and organic and mineral nitrogen over multiple computational pools, while allowing organic matter to be decomposed by multiple
microbial biomass populations. They can account for most of the major reaction
pathways, such as inter-pool transfer of carbon and nitrogen, nitrification (ammonium to nitrate-N), denitrification (leading to the production of N2 and N2O),
volatilization losses of ammonia (NH3), and microbial biomass growth and death(e.g., Yuan et al. 2011).

16.2.2 General Models
Models with generalized chemistry, on the other hand, provide users with much
more flexibility in designing particular chemical systems, thus permitting a muchbroader range of applications. Users can then either select species and reactions
from large geochemical databases or are able to define their own species withparticular chemical properties and reactions. However, general models require
much better understanding of the chemical system involved from their users,
since they themselves are responsible for correctly defining it.
As discussed above, most codes with general geochemistry are limited tosolute transport and biogeochemical reactions, while water flow paths must be
calculated outside of the reactive transport code. Typical examples are PHREEQC(Parkhurst and Appelo 1999), CRUNCH (Steefel 2000), and PHAST (Parkhurst
et al. 2004). Only a few models allow the velocity field to be internally calculated(e.g., MIN3P, Mayer et al. 2002). Several codes for transient, variably-saturated
flow have also been coupled to general biogeochemistry models. These include
3DHYDROGEOCHEM (Yeh and Cheng 1999), TOUGHREACT (Xu et al. 2012),
MIN3P (Mayer et al. 2002), HYTEC (Van der Lee et al. 2003), OpenGeoSys-Gem(Kosakowski and Watanabe 2014), and HP1/2/3 (Jacques and Šimůnek 2005;
Šimůnek et al. 2013). HP1/2/3, an add-on module to HYDRUS-1D and HYDRUS(2D/3D), and its applications for site remediation and risk assessment of mercurycontaminated sites and for evaluating leaching of uranium from a uranium minetailings pile will be discussed below.

16

The Use of Multicomponent Solute Transport Models in Environmental Analyses

16.3

381

Applications of HydroGeoBioChemical
Transport Models

16.3.1 Literature Applications
There are many important agricultural and environmental applications for generalhydrogeobiochemical models, and many of these appeared in the literature during
the last two decades. For example, water leaching from various tailing piles (such as
from uranium mills) often contains many trace metals including molybdenum,
selenium, arsenic, and chromium (Brookins 1984). Using the HYDROGEOCHEM
model, Yeh and Tripathi (1991) simulated the release of trace metals and acidityfrom an acidic uranium mill tailings pile. Narasimhan et al. (1986) used the
DYNAMIX model to study groundwater contamination from an inactive uranium
mill tailings pile. Similar studies have been carried out by Walter et al. (1994),
Lichtner (1996), and Gerke et al. (1998).
Another significant application for hydrogeobiochemical models is to design andevaluate facilities for the safe disposal and long term isolation of hazardous andradioactive waste, especially high-level nuclear waste. The migration of neptunium
between the repository and the ground water table after a hypothetical repository
breach at the high-level nuclear waste repository at Yucca Mountain, Nevada, wasstudied using a multicomponent solute transport model by Viswanathan
et al. (1998). The fate of metal-organic mixed wastes was studied by Rittmann
and VanBriesen (1996) and VanBriesen (1998). Possible consequences of longterm interactions of engineered barriers (clay plugs, concrete components) with thehost rock of a disposal facility or various perturbations as a consequence ofconstruction (e.g. oxidation fronts) or waste content are also frequently studied
using geochemical transport models (Martens et al. 2010, 2011; Kosakowski and
Watanabe 2014). The analysis of redox zone formation in organic-contaminatedaquifers was carried out by Abrams et al. (1998) and Essaid et al. (1995), while
reactive permeable barriers for aquifer remediation were evaluated by Fryar and
Schwartz (1994). Deep geological CO2 storage (e.g., Gherardi et al. 2012; Walshet al. 2013) represents another environmental application of reactive transportmodels. In a study of Jacques et al. (2008a), the transport of major cations andanions and U associated with inorganic P-fertilization were computed for transient
flow conditions. This study clearly demonstrated that the distribution coefficient Kd,
which is often used to simulate the transport of reactive solutes, may vary significantly during the year as a consequence of transient flow conditions. Thus, the “Kdapproach” may represent an inadequate over-simplification of a complex environmental system and may produce erroneous predictions.
Generic multicomponent transport models are sufficiently versatile to allow
consideration of soil biological processes such as those involved in the carbon
cycle (Jacques et al. 2013), nitrogen cycle (Maggi et al. 2008; Gu et al. 2009), andmicrobial dynamics (e.g., Thullner et al. 2005; Lim et al. 2007; Centler et al. 2010).
For example, Thaysen et al. (2014a, b) implemented the CO2 module of UnsatChem

382

J. Šimůnek et al.

into a generic biogeochemical model HP1 and extended it to include the effects of
mixed cementitous materials in the soil profile on enhancing inorganic C-leachingto groundwater. Zhang et al. (2013) evaluated the bacterial transport using the
HYDRUS-PHREEQC framework and transport parameters depending on solution
chemistry.
Although the models with specific chemistry are limited to simulating only thetransport and reactions of a limited number of ions (e.g., major ions), there are many
important agricultural and environmental problems, such as salinization andsodification of irrigated soils, reclamation of saline or sodic soils, sustainability ofvarious irrigation practices, or the release of saline brine water during mining operations, that can be evaluated using these types of models (Šimůnek and Suarez 1997;
Gonçalves et al. 2006; Corwin et al. 2007; Ramos et al. 2011; Rasouli et al. 2013).
Several interesting applications of both the specific and general chemistrymodels have recently appeared in a special issue of Vadose Zone Journal entitled“Reactive Transport Modeling” (Seaman et al. 2012). The model with specificchemistry was used in a study of Reading et al. (2012), who used the
UNSATCHEM code to evaluate the use of gypsum to ameliorate a nonsaline
sodic clay soil in North Queensland, Australia. General reactive transport models
have been used in the recent studies of Bea et al. (2012), Chang et al. (2012), andLichtner and Hammond (2012). Bea et al. (2012) used the MIN3P-D
multicomponent reactive transport code to describe processes occurring within amine tailing impoundment and their impact on CO2 sequestration. Changet al. (2012) used PHREEQC, in combination with laboratory column experiments,
to evaluate the efficacy of Fe(II)-containing solutions as an in-situ means for
reducing toxic hexavalent chromium (Cr(VI)) to trivalent chromium (Cr(III)),
which subsequently coprecipitates with Fe(III). Lichtner and Hammond (2012)used PFLOTRAN to evaluate geochemical processes controlling the aqueous
phase distribution and solid phase speciation of a hexavalent uranium [U(VI)]plume below the Hanford 300 Area bordering the Columbia River in Hanford, WA.

16.3.2 Modeling of Soil Salinity Using the UnsatChem
Module of HYDRUS
The UnsatChem geochemical module (Šimůnek and Suarez 1994) simulates thetransport of major ions (i.e., Ca2+, Mg2+, Na+, K2+, SO42, CO32, and Cl) andtheir equilibrium and kinetic geochemical interactions, such as complexation,
cation exchange, and precipitation-dissolution (e.g., of calcite, gypsum and/or
dolomite). This module has been implemented into both the one- andtwo-dimensional computational versions of HYDRUS (Šimůnek et al. 2008). Possible applications of this module include studies of the salinization/reclamation ofagricultural soils (Šimůnek and Suarez 1997), sustainability of various irrigation
systems (Gonçalves et al. 2006; Ramos et al. 2011; Rasouli et al. 2013), the disposalof brine waters from mining operations, and many other applications (e.g., Schoups
et al. 2006; Skaggs et al. 2006; Corwin et al. 2007). Gonçalves et al. (2006) andThe Use of Multicomponent Solute Transport Models in Environmental Analyses
Rainfall

Water applied (mm)100

Irrigation

ET0

10

80

8

60

6

40

4

20

2

0
01/06/2004

18/12/2004

06/07/2005

22/01/2006

10/08/2006

0
26/02/2007

383

Reference evapotranspiration (mm)16

Fig. 16.1 Daily values of precipitation, irrigation and reference evapotranspiration from June
2004 until February 2007

Ramos et al. (2011) recently demonstrated the applicability of UnsatChem to thesimulation of the multicomponent transport of major ions both in soil lysimeters andin the field, respectively, irrigated with waters of different quality. They used theUnsatChem module of HYDRUS-1D to describe field measurements of the watercontent, overall salinity (expressed as Electric Conductivity, EC), the concentrationof individual soluble cations, the Sodium Adsorption Ratio (SAR), and the
Exchangeable Sodium Percentage (ESP). The experimental data and their analysis
using the UnsatChem module are described in detail by Gonçalves et al. (2006) andRamos et al. (2011) and only selected results from the latter study will be given here.

Problem Description
The field experiment was conducted from June 2004 to February 2007 in the
Alentejo region in Southern Portugal on a medium texture soil. The climate inthis region is mostly dry, sub-humid to semi-arid, with hot dry summers, and mild
winters with irregular rainfall. The experiment involved irrigating maize with either
locally available water or with synthetic (with added NaCl) saline water. The
experimental field was divided into three groups, A, B, and C (only results forplots A and C are reported below), each having a surface area of 6.75 m2 (2.25 m
wide  3 m long; 0.75 m between maize lines), with the Na+ gradient decreasingfrom A (saline waters with EC of 7.8 dS m1 during the first two seasons and15.6 dS m1 during the third season) to C (local water, EC  1.2 dS m1). A drip
irrigation system with an overall constant cumulative discharge of 18 L/h/m(24 mm/h) was used to deliver water, salts (NaCl), and fertilizer (N) to the crop.
The quality of irrigation water was thus the key variable between plots. The
reference evapotranspiration rates (ET0) were calculated using the PenmanMonteith method. Crop evapotranspiration rates (ETc) were then calculated using
the product of ET0 and Kc, where Kc is a crop coefficient accounting for both soil
evaporation and crop transpiration (Fig. 16.1).

384

J. Šimůnek et al.

Water content (cm3 cm-3)

0.5

Measured data (plot A)

HYDRUS-1D (plot A)

Measured data (plot C)

HYDRUS-1D (plot C)

0.4
0.3
0.2
0.1

I

0.0
01/06/04

I
18/12/04

I

06/07/05

22/01/06

10/08/06

26/02/07

Time

Fig. 16.2 Measured and simulated water contents at a depth of 40 cm. Irrigation seasons areindicated by the letter I

TDR probes and ceramic cups were installed at 20, 40, and 60 cm depths tomeasure soil water contents and collect soil solutions twice a week during irrigationseasons, generally 24 h after an irrigation event, and twice a month during the restof the year. The soil solution was monitored for the concentrations of soluble Na+,
Ca2+, Mg2+, ECsw (Electrical Conductivity of Soil Water), and SAR (Sodium
Adsorption Ratio). Soil hydraulic properties for three soil layers were measured
using suction tables, the pressure plate apparatus, and the evaporation method on
undisturbed soil samples, and fitted using the van Genuchten model. The Gapon
selectivity coefficients for four cations (Ca, Mg, K, and Na) were calculated fromthe initial measured soluble and exchangeable cations concentrations.

Results
Figure 16.2 shows measured and simulated water contents at a depth of 40 cm
during the three seasons. During irrigation periods (indicated by a letter I in
Fig. 16.2), the amount of water applied in the experimental field was considerably
higher than the usual amount of water used to irrigate maize in this region in order
to increase water contents above soil field capacity so that soil solution samples
could be collected using the installed ceramic cups. As a result, soil water content
rapidly increased in the beginning of each irrigation season, and then varied
between soil saturation and soil field capacity (Ramos et al. 2011). Between theend of the irrigation seasons (September) and the beginning of rainy seasons(October or November), soil water contents gradually decreased, allowing maize
to mature and be harvested. During rainy seasons, soil water contents reflected theoccurrence of individual rainfall events (Fig. 16.2).
Figure 16.3 shows that soil salinity increased considerably in the plot irrigatedwith saline waters. Measured ECsw in plot A reached values higher than 6.0 dS m1

16

The Use of Multicomponent Solute Transport Models in Environmental Analyses

385

25

ECsw (dS m-1)

20

I

15
10

I

I

Measured data (plot A)
Measured data (plot C)
HYDRUS-1D (plot A)
HYDRUS-1D (plot C)

50
01/06/04

18/12/04

06/07/05

22/01/06

10/08/06

26/02/07

Time

Fig. 16.3 Measured and simulated electrical conductivity of soil water (ECsw) at a depth of 40 cm.
Irrigation seasons are indicated by the letter I

during the first two irrigation seasons, and higher than 16.0 dS m1 during the thirdyear, as applied water became more saline. In the plot irrigated with the locallyavailable water (plot C), soil salinity remained below 5.0 dS m1 throughout all
irrigation seasons. Soil salinity decreased in all plots during rainfall periods, due tosoil leaching. Only in the rainy season of 2004–2005, rainfall was not sufficient tocompletely remove salts from the root zone in plot A, where ECsw values between
2.0 and 6.3 dS m1 were observed in the soil profile (Ramos et al. 2011).
Figure 16.4 shows measured and simulated concentrations of soluble Na+, Ca2+,and Mg2+ for plots A and C at a depth of 40 cm. As the only cation being added tosynthetic saline irrigation waters, the general behavior of sodium was similar to
ECsw. The highest concentrations were reached in plot A during the irrigationseasons. The highest measured Na+ concentrations were about 50 mmol(c)L1during the first 2 years and 130 mmol(c)L1 during the third year. During therainy seasons, Na+ concentrations decreased considerably due to soil leaching,similar to ECsw.
The general dynamics of calcium and magnesium concentrations were similar to
those of sodium and ECsw. The highest measured Ca2+ and Mg2+ concentrations
were about 20 and 10 mmol(c)L1, respectively. During rainfall events, these values
also decreased due to soil leaching. As Ca2+ and Mg2+ were not added to thesynthetic saline irrigation waters applied in plot A, concentrations of these cations
in saline waters were the same as in the fresh waters in plot C. Since Na+ wasapplied in large concentrations to plot A, both Na+ concentrations in the soil
solution and in the solid phase (not shown) increased, leading to soil sodification.
While the exchangeable Na+ concentration in the solid phase increased, the othercations, namely Ca2+ and Mg2+, were inevitably released to the soil solution.
For this reason, concentrations of Ca2+ and Mg2+ were higher in plot A than inplot C where low concentrations of Na+ were present in the irrigation waters. Sincecompared to simpler models based on adsorption isotherms, the UnsatChem module is able to consider cation exchange, it was able to simulate the release of Ca2+

386

J. Šimůnek et al.

Na+ concentration (mmol(c)L-1)

140
120

I

I

I

100
80
60
40
20
0
01/06/04

18/12/04

06/07/05

22/01/06

10/08/06

26/02/07

Mg2+ concentration (mmol(c) L-1)

Ca2+ concentration (mmol(c)L-1)

60
I

I

I

40

20

0
01/06/04

18/12/04

06/07/05

22/01/06

10/08/06

26/02/07

40
30
20

I

I

I

Measured data (plot A)

HYDRUS-1D (plot A)

Measured data (plot C)

HYDRUS-1D (plot C)

10
0
01/06/04

18/12/04

06/07/05

22/01/06

10/08/06

26/02/07

Time

Fig. 16.4 Measured and simulated Na+ (top), Ca2+ (middle), and Mg2+ (bottom) concentrations
at a depth of 40 cm

and Mg2+ from the solid phase as soil sodification increased. The UnsatChem
module thus proved to more adequately represent reality than linear models withadsorption isotherms.
Figure 16.5 presents the results for measured and simulated SAR values. SAR isan integral variable that characterizes salt-affected soils and provides informationon comparative concentrations of Na+, Ca2+, and Mg2+ in soil solutions. Thisvariable takes into consideration that the adverse effects of sodium are moderated

16

The Use of Multicomponent Solute Transport Models in Environmental Analyses

387

SAR (mmol(c) L-1)0.5

60
50

I

40
30
20

I

I

Measured data (plot A)
Measured data (plot C)
HYDRUS-1D (plot A)
HYDRUS-1D (plot C)

10
0
01/06/04

18/12/04

06/07/05

22/01/06

10/08/06

26/02/07

Time

Fig. 16.5 Measured and simulated SAR at a depth of 40 cm

by the presence of calcium and magnesium ions. This makes SAR an importantvariable to consider when managing saline irrigation waters. In plot A, SAR
increased rapidly in the surface layers after irrigation events, and then also gradually at deeper depths. The application of saline irrigation waters (SAR ranged from
38.8 to 80.2 (mmol(c)L1)0.5) led to very high SAR values in the soil solution,
reaching about 54 (mmol(c)L1)0.5 at the end of the experiments. In plot C, irrigatedwith fresh waters, SAR values did not vary significantly, having roughly the samevalues in the beginning and at the end of the experiments.
Various statistical measures evaluating the goodness of correspondence between
the model (the UnsatChem module of HYDRUS-1D) and measurements are given
by Ramos et al. (2011) and will not be given here. The correspondence between
measurements and model results is very good, considering that most model inputs
were independently measured in the laboratory and used in simulations without any
further adjustments and/or calibration. The correspondence between measurements
and model results obviously could have been better, had the input parameters beencalibrated (Ramos et al. 2011). An extensive discussion on the possible causes ofdeviations between measured data and simulation can also be found in Ramos
et al. (2011, 2012).

16.3.3 Modeling the Fate of Mercury in ContaminatedSoils Using HP1
The HP1 module, which is used in this example, couples the PHREEQC geochemical code (Parkhurst and Appelo 1999) with HYDRUS-1D (Šimůnek et al. 2008).
The one-dimensional version (HP1) was first released in 2005 (Jacques andŠimůnek 2005; Jacques et al. 2008a, b), and used successfully in many applications.
HP1, which is an acronym for HYDRUS-PHREEQC-1D, is a relatively comprehensive simulation module that can be used to simulate (1) transient water flow,

388

J. Šimůnek et al.

(2) the transport of multiple components, (3) mixed equilibrium/kinetic biogeochemical reactions, and (4) heat transport in one-dimensional variably-saturated
porous media. HP2 is a corresponding two-dimensional version of HP1, which
couples PHREEQC with HYDRUS (2D/3D). Both HP1 and HP2 modules cansimulate a broad range of low-temperature biogeochemical reactions in water, thevadose zone and/or ground water systems, including interactions with minerals,
gases, exchangers and sorption surfaces based on thermodynamic equilibrium, andkinetic or mixed equilibrium-kinetic reactions.
The example described below illustrates the possibility to simulate the fate of
mercury (Hg) in anthropogenic contaminated soil systems (Leterme and Jacques
2013; Leterme et al. 2014).

Problem Definition
Mercury (Hg) poses threats to human health and the environment, notably due to its
persistence in the environment and its ability to bioaccumulate in ecosystems (Liuet al. 2012). Mercury pollution can result from either direct contamination (e.g.,
spills, landfills, mine tailings) or indirect pathways, such as atmospheric deposition(Guédron et al. 2013). Mercury mining, gold and silver mining, manufacturing(chlor-alkali plants, manometer spill), wood preservation and cemeteries (throughthe release of Hg from dental amalgams) have been identified as the main sources ofsoil contamination (UNEP 2002). In addition, Hg geochemistry is quite complex(Fig. 16.6) and thus Hg contamination can occur in different forms: as aqueous
inorganic and organic HgII species, in elemental form (Hg0) in a (residual)
non-aqueous liquid phase (NAPL), or as a solid phase (e.g., cinnabar). Hg transport
is highly influenced by its strong binding to organic matter, especially to the thiolgroups, which decreases (sorption to soil organic matter) or increases (complexation with dissolved organic matter) the mobility of mercury. Reduction of HgII to
Hg0 enhances the Hg transport through the gaseous phase, which is an important
pathway in mercury-related risk assessment studies. Note that under reducing
conditions, organic forms of mercury (monomethylmercury and dimethylmercury)
may occur and may have a strong tendency for bioaccumulation in wetland
systems.
Leterme et al. (2014) implemented the conceptual model presented in Fig. 16.6
into HP1. This model addresses Hg-contaminated soils under oxidizing conditions(organic forms of Hg are excluded). They investigated seven scenarios with initialconditions representing different types of contamination (cinnabar, residual NAPL,
and aqueous HgII, and any combination of these three contamination sources)
polluting the first 10 cm of the soil (Fig. 16.7). Only the first three scenarios arepresented below. They also evaluated the sensitivity of four different indicators,
illustrated in Fig. 16.7, to different parameters and processes.
Daily climatic data (50 years) at Dessel, Belgium, were used to define surface
atmospheric boundary conditions and free drainage was considered as the bottom
boundary condition (representing deep groundwater and oxidising conditions).

16

The Use of Multicomponent Solute Transport Models in Environmental Analyses

389

Fig. 16.6 Conceptual model of Hg speciation and reactions in soil systems (DOM dissolved
organic matter, SOM soil organic matter, MMHg monomethylmercury, DMHg dimethylmercury)

#

1
3

4

2

Indicator

1

Percentage of Hg volatilized to atmosphere

2

Percentage of Hg leached out of the soil profile

3

Percentage of Hg present in the originally-polluted
soil horizon (top 10 cm)

4

Percentage of Hg present in the soil profile (10100 cm) below the originally-polluted soil horizon

Fig. 16.7 Problem definition and evaluated indicators. Initial Hg contamination is in the first
10 cm of the soil (red shaded). The top 30 cm of the soil contains soil organic matter (grey shaded),
the remaining part of the soil profile is inert

Incoming water has a typical rain water composition and contains dissolved organic
matter (DOM). Soil texture is sandy (this a priori corresponds to coarse-textured
soil types, likely to be encountered at polluted sites such as mine tailings) with agrass cover (30 cm rooting depth). Soil organic matter (SOM) and the initial Hg
contamination are uniform in the top 30 and 10 cm of the soil profile, respectively.
There is a background atmospheric Hg(g) concentration of 1 ng m3. Soil anddissolved organic matter are both represented by four different reactive sites (three
representing oxygen sites of fulvic and humic acids and one representing thiols).

390

J. Šimůnek et al.

Results
Figure 16.8 shows depth profiles of the initial contamination source and Hgconcentrations sorbed to SOM for three contamination scenarios and selected
times (initial, and 5, 10, 25, and 50 years). Since the dissolution rate of cinnabar
is relatively slow, about 50 % of the initial amount still remains after 50 years(Fig. 16.8, top). Note the non-uniform dissolution profile of cinnabar due to thedependence of the dissolution kinetics on DOM concentrations. Because of theslow dissolution rate of cinnabar and the presence of the strongly binding DOM
thiols, almost no Hg is retained on the SOM in the first few centimetres of the soil.
As DOM thiols move downward through the contaminated zone (top 10 cm), they
are increasingly more saturated with mercury due to continued cinnabar dissolution, and more Hg is retained by SOM, resulting in a peak retention of Hg at thebottom of the initially polluted soil layer (Fig. 16.8, top).
In the scenario with residual NAPL (Fig. 16.8, middle), all elemental Hg isdissolved after about 10 years. Because of this gradual release, mercury is mostly
retained within the top 10 cm during the first 25 years. Since the release of Hg from
NAPL is fast enough to saturate the thiol groups of DOM, a significant portion of
mercury is retained on SOM. After 25 years, the gradual leaching of mercury to thelower soil horizons is obvious.
In the scenario with initial aqueous contamination (Fig. 16.8, bottom), a large
portion of Hg leaches very quickly (mainly in the form of inorganic species such as
Hg(OH)2, HgCl2, HgOCl) from the top 10 cm and sorbs significantly to the SOM in
the underlying 20 cm. Note that the “Initial” concentration profile is drawn after the
first time step, which explains the “Initial” profile of Hg sorbed to SOM. Also fromthe start, all Hg is available for leaching via complexation on the thiols of DOM(although binding is much stronger on the thiols of DOM than of SOM, theircapacity, i.e., the number of thiol sites, is very limited), resulting in a gradualleaching of Hg from the top 10 cm.
Figure 16.9 shows HgII fluxes from the top 1 m to underlying layers for the first
25 years of the simulations. A quick initial release of mobile Hg in the scenario witha HgCl2(aq) contamination (a green line in Fig. 16.9) leads to earlier Hg leaching
than in scenarios with cinnabar and NAPL contaminations (red and blue lines,
respectively). In these two scenarios, Hg leaching occurs only after a few years, but
remains at only about half of the Hg flux in the scenario with the
HgCl2(aq) contamination. This is because in the latter case Hg was rapidly
redistributed in the top 30 cm (Fig. 16.8, bottom right) and less re-sorption occurs
when it complexes with DOM thiols. Hg fluxes for scenarios with cinnabar and Hg
NAPL are quite similar because the limiting factor for Hg mobility is not the Hgrelease through dissolution nor the sorption capacity of the soil solid phase, butrather the availability of DOM thiol groups. Figure 16.9 also shows that thetemporal dynamics (not its absolute values) of leaching is driven by the climate.
The sensitivity analysis indicated that the most sensitive parameters are related
to SOM and DOM. DOM concentration strongly influence the mobility and16

The Use of Multicomponent Solute Transport Models in Environmental Analyses

391

Fig. 16.8 Depth profiles at selected times for different forms of Hg for different types of initial
contamination: cinnabar (top), residual NAPL (middle), and aqueous Hg species (bottom). The left
column shows Hg concentrations of the original contamination, the right columns shows Hg
sorbed to SOM. Note different scales of the x-axis

leaching of mercury due to the high affinity of Hg to DOM thiol groups. On theother hand, the exchange capacities and the selectivity coefficients of the fulvic andhumic acids of SOM greatly influence the retention of mercury in the top soil. Notealso that the selectivity coefficient of the thiol groups in SOM has only a minor
effect because in anthropogenic mercury-contaminated sites, the thiol groups in
SOM are initially saturated with mercury. The dissolution rate of the initial

392

J. Šimůnek et al.

Fig. 16.9 Precipitation (left axis) and HgII fluxes (monthly moving averages) at the bottom of thesoil profile for three initial contamination scenarios (right axis)

contamination is also a sensitive parameter as illustrated above in scenarios withcinnabar and residual NAPL.

16.3.4 Modeling of Uranium Leaching from a U Tailings
Pile Using the HP2 Module of HYDRUS (2D/3D)
The HP2 is a two-dimensional equivalent of the one-dimensional HP1 program(Šimůnek et al. 2008). We will demonstrate below the versatility of HP2 using an
example involving the release and migration of uranium from a simplified uranium
mill tailings pile toward a river. This example includes the processes of water flow,
solute transport, precipitation/dissolution of the solid phase, cation exchange,
complexation, and many other reactions. This problem was inspired by, and is a
modification of (to make it more realistic), a problem reported by Yeh and Tripathi(1991). Detailed information about this project, as well as instructions on how to
implement the project using HP2, can be found in Šimůnek et al. (2012) and on the
HYDRUS website.

Problem Description
The schematic of the transport domain is shown in Fig. 16.10. The mill tailings pile
is located adjacent to a surface that slopes down to a river. The horizontal bottom of
the region is impermeable. The vertical left-edge has the Dirichlet boundary

16

The Use of Multicomponent Solute Transport Models in Environmental Analyses

393

Cauchy Flow
= 0.139 dm/day
Mill Tailing
Pile

Infiltration Rate
= 0.0139 dm/day
Dirichlet Flow
Nodes

240 dmPumping Well

39 dm
1,050 dm

Fig. 16.10 Problem description for the uranium tailings problem (Adopted from Yeh and Tripathi
1991)

condition with a groundwater table 12 m above the bottom of the transport domain.
The top boundary (except for the mill tailings pile and the river) has a flux boundary
condition with a net rainfall rate of 0.139 cm/day. The horizontal region on the topof the mill tailings pile is a flux boundary with an infiltration rate of 1.39 cm/day.
The nodes on the vertical line on the right side of the transport domain and thenodes on the river bottom have the Dirichlet boundary condition, reflecting theposition of water in the river (4.5 m above the bottom of the transport domain). A
hypothetical pumping well with a withdrawal rate of 271.58 cm2/day is located at(x, z) ¼ (400, 100). The region is discretized using a structured finite element meshwith 1,564 elements and 852 nodes.
The medium has the hydraulic properties of a loam (the residual water content,
θr ¼ 0.078, the saturated water content, θs ¼ 0.430, and the van Genuchten (1980)
parameters α ¼ 0.036 cm1, and n ¼ 1.56; Carsel and Parrish 1988) with thesaturated hydraulic conductivity of Ks ¼ 3.78 m/day. A longitudinal dispersivityof 2.5 m and a transverse dispersivity of 0.25 m are assumed.
With respect to solute transport, uranium was modeled as a complex
multicomponent system involving multiple interacting components. Table 16.1
lists chemical reactions and their thermodynamic equilibrium constants consideredin this example. For reactive hydrogeochemical transport, the problem consists ofeight components: Total H, Total O, Ca, C, uranium, sulfate, phosphate, and Fe. A
total of 35 aqueous species, 18 exchange species, 9 surface species and 2 minerals
are defined for the problem; redox reactions were not considered. Precipitation/
dissolution and sorption are the main chemical processes considered in this example. Gypsum (0.0037 mol/kg) and calcite (0.00047 mol/kg) are assumed to bepresent in the waste zone and the rest of the profile, respectively.
Sorption of U was described using a multi-site cation exchange model (see Jacqueset al. 2008a). This sorption model also buffers the acid pH due to proton exchange. In
the current model, protons will replace Ca in the exchange complex. Parameters of themulti-site cation exchange model are also given in Table 16.1 (adapted from Jacques

394

J. Šimůnek et al.

Table 16.1 Reaction network for the uranium tailings pile leaching problem
Reaction
Aqueous complexation reactions
H2O , H+ + OH
Ca2+ + CO32 , CaCO3 (aq)Ca2+ + H+ + CO32 , CaHCO3+
Ca2+ + SO42 , CaSO4 (aq)Ca2+ + 2H+ + PO43 , CaH2PO4+
Ca2+ + PO43 , CaPO4
Ca2+ + H+ + PO43 , CaHPO4 (aq)Ca2+ + H2O , H+ + CaOH+
Fe2+ + SO42 , FeSO4 (aq)Fe2+ + H2O , H+ + FeOH+
Fe2+ + 2H2O , 2H+ + Fe(OH)2 (aq)Fe2+ + 3H2O , 3H+ + Fe(OH)3
Fe2+ + 4H2O , 4H+ + Fe(OH)42
UO22+ + H2O , H+ + (UO2)(OH)+
2UO22+ + 2H2O , 2H+ + (UO2)2(OH)22+
3UO22+ + 4H2O , 4H+ + (UO2)3(OH)42+
3UO22+ + 5H2O , 5H+ + (UO2)3(OH)5+
4UO22+ + 7H2O , 7H+ + (UO2)4(OH)7+
3UO22+ + 7H2O , 7H+ + (UO2)3(OH)7
UO22+ + CO32 , (UO2)(CO3)(aq)
UO22+ + 2CO32 , (UO2)(CO3)22
UO22+ + 3CO32 , (UO2)(CO3)34
2UO22+ + CO32 + 3H2O , 3H+ + (UO2)2(CO3)(OH)3
UO22+ + SO42 , (UO2)(SO4) (aq)UO22+ + 2SO42 , (UO2)(SO4)22
2H+ + UO22+ + PO43 , H2(UO2)(PO4)+
3H+ + UO22+ + PO43 , H3(UO2)(PO4)2+
Ca2+ + 4H+ + UO22+ + 2PO43 , CaH4(UO2)(PO4)22+
Ca2+ + 5H+ + UO22+ + 2PO43 , CaH5(UO2)(PO4)23+
H+ + CO32 , HCO3
2H+ + CO32 , H2CO3(aq)H+ + SO42 , HSO4
H+ + PO43 , HPO42
2H+ + PO43 , H2PO4
3H+ + PO43 , H3PO4
Precipitation-dissolution reactions
Ca2+ + SO42 , CaSO4(s) (gypsum)Ca2+ + CO32 , CaCO3(s) (calcite)
Multi-site cation exchange modelb
Ca2+ + 2 Yi ¼ CaYi2
UO22+ + 2 Yi ¼ UO2Yi2
H+ + Ya ¼ HYa
H+ + Yb ¼ HYb
H+ + Yc ¼ HYc
H+ + Yd ¼ HYd
H+ + Ye ¼ HYe
H+ + Yf ¼ HYf

No.

Log Ka

(R1)
(R2)
(R3)
(R4)
(R5)
(R6)
(R7)
(R8)
(R9)
(R10)
(R11)
(R12)
(R13)
(R14)
(R15)
(R16)
(R17)
(R18)
(R19)
(R20)
(R21)
(R22)
(R23)
(R24)
(R25)
(R26)
(R27)
(R28)
(R29)
(R30)
(R31)
(R32)
(R33)
(R34)
(R35)

14.00
3.22
11.43
2.31
20.96
6.46
15.08
12.582.20
9.50
20.57
31.0043.00
5.30
5.68
11.88
15.82
21.9028.34
9.65
17.08
21.70
1.18
2.95
4.00
23.20
22.90
45.24
46.00
10.32
16.67
1.99
12.35
19.55
21.74

(R36)
(R37)

4.62
8.48

(R37)
(R38)
(R39)
(R40)
(R41)
(R42)
(R43)
(R44)

0.2
0.2
1.65
3.30
4.95
6.85
9.60
12.35
(continued)

16

The Use of Multicomponent Solute Transport Models in Environmental Analyses

395

Table 16.1 (continued)
ReactionNo.
Log KaSpecific sorption – surface complexation
(R45)
7.29
Hfo_OH + H+ ¼ Hfo_OH2+
(R46)
8.93
Hfo_OH ¼ Hfo_O + H+
(R47)
5.85
Hfo_OH + Ca2+ ¼ Hfo_OCa+ + H+
(R48)
2.8
Hfo_OH + UO22+ ¼ Hfo_OUO2+ + H+
(R49)
31.29
Hfo_OH + PO43+ + 3H+ ¼ Hfo_H2PO4 + H2O
Hfo_OH + PO43+ + 2H+ ¼ Hfo_HPO4 + H2O
(R50)
25.39
(R51)
17.72
Hfo_OH + PO43+ + H+ ¼ Hfo_PO42 + H2O
(R52)
7.78
Hfo_OH + SO42 + H+ ¼ Hfo_SO4 + H2O
(R53)
0.79
Hfo_OH + SO42 ¼ Hfo_OHSO42
a
K is a thermodynamic equilibrium constantb
Yi represents one of the six different exchange sites (Ya, Yb, Yc, Yd, Ye, Yf). Exchange
coefficients are the same for each exchange site for Ca2+ and UO22+, but are different for H+

Table 16.2 Initial andboundary compositions ofrecharge water and pore waterin the tailings and regions
outside of the tailings forthe uranium tailings problem(mol/l)Species
Inside the tailings
Outside the tailings
1.14  102a1.42  102a
Ca2+
2
2
CO3
1.0  10
1.5  103
2+
4
UO2
5.0  101.0  107PO43
1.0  106
7.5  107b
2c
2
SO4
5.33  10
1.35  102d2+
2
Fe
3.5  101.0  107
pH
2.3
7.14c
a
Calculated in equilibrium with calciteb
Calculated slightly under-saturated with hydroxyapatite (thesaturation index of 0.1)
c
Used as charge balancedCalculated slightly under-saturated with gypsum (the saturationindex of 0.1)

et al. 2008a). A relatively low (8.1  103 mol/kg) sorption capacity and a bulk
density of 1.5 g/cm3 were considered. Additionally, specific sorption of U and othercations and anions on Fe-oxides were described using a non-electrostatic surface
complexation model with a capacity of 8.1  104 mol/kg. Parameters of thisnon-electrostatic surface complexation model are also given in Table 16.1.
Two different water compositions were considered: (i) one solution composition
was considered for pore water and recharge water in the tailings, and (ii) the second
solution composition was used for pore water and recharge water outside of thetailings pile and for recharge water at the left hand side of the domain. Both watertypes are given in Table 16.2. Cauchy boundary conditions were considered on all
boundaries where boundary conditions are specified. A total of 1,000 days was
simulated.

396

J. Šimůnek et al.

Fig.16.11 The steady-state pressure head (cm) (top) and flux (cm/d) (bottom) profiles (at 1,000days) for the tailings pile leaching example

Fig. 16.12 Spatial distribution of pH initially (left) and at 1,000 days (right) for the tailings pileleaching example

Results
The steady-state pressure heads and velocity fields are depicted in Fig. 16.11.
Figure 16.11 shows that water fluxes in the groundwater (about 20-50 cm/day)were significantly larger than in the vadose zone (below 1.5 cm/day).
Figure 16.12 shows contour plots of the initial (left) and final (right) distributionsof pH. Acidity was initially confined to the tailings pile and was larger than 7 outsideof the tailings pile. The proton exchange on the cation exchanger buffered theacidity by replacing calcium with protons on the exchanger. The spatial extent of
the pH-altered region extended all the way to the river.
Figure 16.13 depicts contour plots of precipitated carbonate (Calcite, CaCO3(s))
and sulfate (Gypsum, CaSO4(s)) at 1,000 days. CaCO3(s), which was initially present
everywhere outside of the tailings, is dissolved completely in the zone where pH islower than about 5.5 (see Fig. 16.12). The influence of the acidity on calcite
precipitation is present almost in the entire domain, except for the upper right andbottom left parts of the domain. In these zones, the geochemistry was governed byincoming rain water (upper part) or ground water (bottom left). Dissolution of

16

The Use of Multicomponent Solute Transport Models in Environmental Analyses

397

Fig. 16.13 Spatial distribution of calcite (a) and gypsum (b) at 1,000 days for the tailings pileleaching example. Arrows indicate the color representing the initial conditions for calcite in thesoil region or for gypsum in the tailings pileFig. 16.14 Spatial distribution of aqueous U concentration (mol/l) at 1,000 days

Fig. 16.15 Spatial distribution of sorbed U concentration (mol/dm3) at 1,000 days

calcite and the abundance of sulphate in tailings water led to gypsum precipitation.
Gypsum, which was initially present only inside of the tailings, precipitated in thedomain where groundwater encountered the leachate from the tailings.
Figures 16.14 and 16.15 show dissolved and sorbed uranium concentration
profiles, respectively, as calculated using HP2 (i.e., using complex biogeochemical
calculations). U migration within the ground water system is considerably retarded.
U sorption seems to coincide quite well with the downstream edge of the aqueous U
plume.

398

16.4

J. Šimůnek et al.

Conclusions

In spite of the considerable demand on input data, the multicomponent solutetransport models such as HYDRUS, UnsatChem, or HP1/2, can be effective andversatile tools that are very useful for evaluating complex agricultural and environmental problems, such as irrigation management in regions with scarce water
resources, where suitable waters are not always available for irrigation or forevaluating the leaching of acidic waters from mine tailings. We have demonstrated
the use of the UnsatChem and HP2 modules of the HYDRUS software for these two
types of applications, respectively. These models, after proper calibration andvalidation, should be considered useful tools for establishing sound irrigation andenvironmental policies to mitigate various problems such as soil salinization/
sodification or point and non-point source pollution from mining operations. Additional examples of the use of the complex hydrobiogeochemical transport models
can be found at (and downloaded from) the HYDRUS website.

References
Abrams RH, Loague K, Kent DB (1998) Development and testing of a compartmentalized reaction
network model for redox zones in contaminated aquifers. Water Resour Res 34:1531–1541
Ahuja LR, Hebson C (1992) Root zone water quality model. GPSR Tech. Rep. No. 2. USDA,
ARS, Fort Collins
Bea SA, Wilson SA, Mayer KU, Dipple GM, Power IM, Gamazo P (2012) Reactive transportmodeling of natural carbon sequestration in ultramafic mine tailings. Vadose Zone J 11(2).
doi:10.2136/vzj2011.0053
Bell LSJ, Binning PJ (2004) A split operator approach to reactive transport with the forward
particle tracking Eulerian-Lagrangian localized adjoint method. Adv Water Resour
27:323–334
Brookins DG (1984) Geochemical aspects of radioactive waste disposal. Springer, New York
Carsel RF, Parrish RS (1988) Developing joint probability distributions of soil water retention
characteristics. Water Resour Res 24:755–769
Centler F, Shao H, De Biase C, Park C-H, Regnier P, Kolditz O, Thullner M (2010) GeoSysBRNS –a flexible multidimensional reactive transport model for simulating biogeochemical subsurface
processes. Comput Geosci 36:397–405
Chang H, Singer JH, Seaman JC (2012) In situ chromium(VI) reduction using iron(II) solutions:modeling dynamic geochemical gradients. Vadose Zone J 11(2). doi:10.2136/vzj2011.0172
Corwin DL, Rhoades JD, Šimůnek J (2007) Leaching requirement for soil salinity control: steadystate vs. transient-state models. Agric Water Manage 90(3):165–180
Essaid HI, Bekins BA, Godsy EM, Warren E, Baedecker MJ, Cozzarelli IM (1995) Simulation of
aerobic and anaerobic biodegradation processes at a crude oil spill site. Water Resour Res
31(12):3309–3327
Fryar AE, Schwartz FW (1994) Modeling the removal of metals from groundwater by a reactive
barrier: experimental results. Water Resour Res 30(12):3455–3469
Gerke HH, Molson JW, Frind EO (1998) Modelling the effect of chemical heterogeneityon acidification and solute leaching in overburden mine spoils. J Hydrol 209:166–185

16

The Use of Multicomponent Solute Transport Models in Environmental Analyses

399

Gherardi F, Audigane P, Gaucher EC (2012) Predicting long-term geochemical alteration of
wellbore cement in a generic geological CO2 confinement site: tackling a difficult reactive
transport modeling challenge. J Hydrol 420–421:340–359
Gonçalves MC, Šimůnek J, Ramos TB, Martins JC, Neves MJ, Pires FP (2006) Multicomponent
solute transport in soil lysimeters irrigated with waters of different quality. Water Resour Res
42:W08401, 17 pp. doi:10.1029/2006WR004802
Gu CH, Maggi F, Riley WJ, Hornberger GM, Xu T, Oldenburg CM, Spycher N, Miller NL,
Venterea RT, Steefel C (2009) Aqueous and gaseous nitrogen losses induced by fertilizer
application. J Geophys Res 114, G01006
Guédron S, Grangeon S, Jouravel G, Charlet L, Sarret G (2013) Atmospheric mercury incorporation in soils of an area impacted by a chlor-alkali plant (Grenoble, France): contribution of
canopy uptake. Sci Total Environ 445–446:356–364
Hutson JL, Wagenet RJ (1992) LEACHM: leaching estimation and chemistry model. Researchseries no. 92-3, Cornell University, Ithaca
Jacques D, Šimůnek J (2005) User manual of the multicomponent variably-saturated flow andtransport model HP1, description, verification and examples, version 1.0, SCK•CEN-BLG998, Waste and disposal. SCK•CEN, Mol, 79 pp
Jacques D, Šimůnek J, Mallants D, van Genuchten MT (2006) Operator-splitting errors in coupled
reactive transport codes for transient variably saturated flow and contaminant transport inlayered soil profiles. J Contam Hydrol 88:197–218
Jacques D, Šimůnek J, Mallants D, van Genuchten MT (2008a) Modeling coupled hydrologicaland chemical processes: long-term uranium transport following mineral phosphorus fertilization. Vadose Zone J 7(2):698–711. doi:10.2136/VZJ2007.0084
Jacques D, Šimůnek J, Mallants D, van Genuchten MT (2008b) Modelling coupled water flow,
solute transport and geochemical reactions affecting heavy metal migration in a podzol soil.
Geoderma 145:449–461
Jacques D, Šimůnek J, Mallants D, van Genuchten MTh (2013) The HPx reactive transportmodels: summary of recent developments and applications. In: Šimůnek J, van Genuchten
MTh, Kodešová R (eds) Proceedings of the 4th international conference HYDRUS software
applications to subsurface flow and contaminant transport problems. Invited paper, 21–22 Mar
2013. Department of Soil Science and Geology, Czech University of Life Sciences, Prague.
ISBN: 978-80-213-2380-3, pp 7–16
Jansson P-E, Karlberg L (2001) Coupled heat and mass transfer model for soil-plant-atmospheresystems. Royal Institute of Technology, Department of Civil and Environmental Engineering,
Stockholm, 325 pp
Kosakowski G, Watanabe N (2014) OpenGeoSys-Gem: a numerical tool for calculating geochemical and porosity changes in saturated and partially variable saturated media. Phys Chem Earth(in press). doi:http://dx.doi.org/10.1016/j.pce.2013.11.008
Langergraber G, Šimůnek J (2005) Modeling variably-saturated water flow and multi-component
reactive transport in constructed wetlands. Vadose Zone J 4(4):924–938
Langergraber G, Šimůnek J (2012) Reactive transport modeling of subsurface flow constructedwetlands. Vadose Zone J 11(2). doi:10.2136/vzj2011.0104, 14 pp
Leij FJ, Bradford SA (1994) 3DADE: a computer program for evaluating three-dimensional
equilibrium solute transport in porous media. Research report no. 134. U. S. Salinity Laboratory, USDA, ARS, Riverside
Leterme B, Jacques D (2013) Modeling Hg reactive transport in soil system using HP1. In:
Šimůnek J, van Genuchten MTh, Kodešová R (eds) Proceedings of the 4th international
conference HYDRUS software applications to subsurface flow and contaminant transport
problems. Invited paper, 21–22 Mar 2013. Department of Soil Science and Geology, Czech
University of Life Sciences, Prague. ISBN: 978-80-213-2380-3, pp 225–234
Leterme B, Blanc P, Jacques D (2014) A reactive transport model for mercury fate in soils –application to different anthropogenic pollution sources. Soil Pollut (in review)400

J. Šimůnek et al.

Lichtner PC (1996) Continuum formulation of multicomponent-multiphase reactive transport. In:
Lichtner PC, Steefel CI, Oelkers EH (eds) Reactive transport in porous media, reviews inmineralogy, Chapter 1, vol 34. Mineralogical Society of America, Washington, DC, pp 1–81
Lichtner PC, Hammond GE (2012) Using high performance computing to understand roles oflabile and nonlabile uranium(VI) on Hanford 300 area plume longevity. Vadose Zone J 11(2).
doi:10.2136/vzj2011.0097
Lim M-S, Yeo IW, Clement TP, Roh Y, Lee K-K (2007) Mathematical model for predicting
microbial reduction and transport of arsenic in groundwater systems. Water Res 41:2079–2088
Liu G, Cai Y, O’Driscoll N, Feng X, Jiang G (2012) Overview of mercury in the environment.
In: Liu G, Cai Y, O’Driscoll N (eds) Environmental chemistry and toxicology of mercury.
Chantilly, VA, pp 1–12
Maggi F, Gu C, Riley WJ, Venterea R, Hornberger GM, Venterea RT, Xu T, Spycher N, Steefel C,
Miller NL, Oldenburg CM (2008) A mechanistic treatment of the dominant soil nitrogen
cycling processes: model development, testing, and application. J Geophys Res Biogeosci
113, G02016. doi:10.1029/2007JG000578
Mangold DC, Chin-Fu Tsang (1991) A summary of subsurface hydrological and hydrochemical
models. Rev Geophys 29(1):51–79
Martens E, Jacques D, Van Gerven T, Wang L, Mallants D (2010) Geochemical modelling ofleaching of Ca, Mg, Al, and Pb from cementitious waste forms. Cement Concr Res
40(8):1298–1305
Martens E, Jacques D, Wang L, De Cannière P, Moors H, Van Gompel M, Mariën A, Valcke E(2011) Modelling of cation concentrations in the outflow of NaNO3 percolation experiments
through boom clay cores. Phys Chem Earth 36:1693–1699
Mayer KU, Frind EO, Blowes DW (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
Narasimhan TN, White AF, Tokunaga T (1986) Groundwater contamination from an inactive
uranium mill tailings pile, 2. Application of a dynamic mixing model. Water Resour Res
22(13):1820–1834
Parkhurst DL, Appelo CAJ (1999) User’s guide to PHREEQ C (version 2) – a computer programfor speciation, batch-reaction, one-dimensional transport and inverse geochemical calculations, water-resources investigations. Report 99–4259, Denver, 312 pp
Parkhurst DL, Kipp KL, Engesgaard P, Charlton SC (2004) PHAST – a program for simulating
ground-water flow, solute transport and multicomponent geochemical reactions. U.S. Geol
Surv Tech Methods 6-A8, 154 pp
Parton WJ, Schimel DS, Cole CV, Ojima DS (1987) Analysis of factors controlling soil organic
matter levels in Great Plains grasslands. Soil Sci Soc Am J 51:1173–1179
Ramos TB, Šimůnek J, Gonçalves MC, Martins JC, Prazeres A, Castanheira NL, Pereira LS (2011)
Field evaluation of a multicomponent solute transport model in soils irrigated with saline
waters. J Hydrol 407(1–4):129–144
Ramos TB, Šimůnek J, Gonçalves MC, Martins JC, Prazeres A, Pereira LS (2012)
Two-dimensional modeling of water and nitrogen fate from sweet sorghum irrigated withfresh and blended saline waters. Agric Water Manage 111:87–104
Rasouli F, Pouya AK, Šimůnek J (2013) Modeling the effects of saline water use in wheatcultivated lands using the UNSATCHEM model. Irrig Sci 31(5):1009–1024. doi:10.1007/
s00271-012-0383-8
Reading LP, Baumgartl T, Bristow KL, Lockington DA (2012) Applying HYDRUS to flow in a
sodic clay soil with solution composition–dependent hydraulic conductivity. Vadose Zone J 11(2). doi:10.2136/vzj2011.0137
Regnier P, Dale AW, Pallud C, van Lith Y, Bonneville S, Hyacinthe C, Thullner M, Laverman
AM, Van Cappellen P (2005) Incorporating geomicrobial processes in reactive transportmodels of subsurface environments. In: Nützmann G, Viotti P, Aagaard P (eds) Reactive
transport in soil and groundwater. Springer, Berlin Heidelberg, pp 109–125

16

The Use of Multicomponent Solute Transport Models in Environmental Analyses

401

Rittmann BE, VanBriesen JM (1996) Microbiological processes in reactive modeling. In: Lichtner
PC, Steefel CI, Oelkers EH (eds) Reactive transport in porous media, reviews in mineralogy,
Chapter 7, vol 34. Mineralogical Society of America, Chantilly, VA, pp 311–334
Schoups G, Hopmans JW, Tanji KK (2006) Evaluation of model complexity and space–timeresolution on the prediction of long-term soil salinity dynamics, western San Joaquin Valley,
California. Hydrol Process 20:2647–2668
Seaman JC, Chang H, Goldberg S, Šimůnek J (2012) Reactive transport modeling. Vadose Zone
J 11(2). doi:10.2136/vzj2012.0066, 7 pp
Šimůnek J (2005) Chapter 78: models of water flow and solute transport in the unsaturated zone.
In: Anderson MG, McDonnell JJ (eds) Encyclopedia of hydrological sciences. Wiley, Chichester,
pp 1171–1180
Šimůnek J, Suarez DL (1994) Two-dimensional transport model for variably saturated porous
media with major ion chemistry. Water Resour Res 30(4):1115–1133
Šimůnek J, Suarez DL (1997) Sodic soil reclamation using multicomponent transport modeling.
ASCE J Irrig Drain Eng 123(5):367–376
Šimůnek J, Valocchi AJ (2002) Geochemical transport. In: Dane JH, Topp GC (eds) Methods ofsoil analysis, part 4, physical methods, Chapter 6.9, 3rd edn. SSSA, Madison, pp 1511–1536
Šimůnek J, van Genuchten MT (2006) Contaminant transport in the unsaturated zone: theory andmodeling, Chapter 22. In: Delleur J (ed) The handbook of groundwater engineering, 2nd edn.
CRC Press, Boca Raton, pp 22.1–22.46
Šimůnek J, van Genuchten MTh, Šejna M, Toride N, Leij FJ (1999) The STANMOD computer
software for evaluating solute transport in porous media using analytical solutions ofconvection-dispersion equation. Versions 1.0 and 2.0, IGWMC – TPS – 71. International
Ground Water Modeling Center, Colorado School of Mines, Golden, 32 pp
Šimůnek J, van Genuchten MT, Sejna M (2008) Development and applications of the HYDRUS
and STANMOD software packages, and related codes. Vadose Zone J 7(2):587–600. doi:10.
2136/VZJ2007.0077
Šimůnek J, Jacques D, Šejna M, van Genuchten MTh (2012) The HP2 program for HYDRUS(2D/3D): a coupled code for simulating two-dimensional variably-saturated water flow, heattransport, and biogeochemistry in porous media, version 1.0, PC Progress, Prague, 76 pp
Šimůnek J, Jacques D, Langergraber G, Bradford SA, Šejna M, van Genuchten MTh (2013)
Numerical modeling of contaminant transport with HYDRUS and its specialized modules,
invited paper for the special issue “water management in changing environment”, Mohan
Kumar MS (ed) J Indian Inst Sci 93(2) 265–284. ISSN: 0970-4140 Coden-JIISAD
Skaggs TH, Shouse PJ, Poss JA (2006) Irrigation of forage crops with saline drainage waters:
2. Modeling drainage and root water uptake. Vadose Zone J 5:824–837
Steefel CI (2000) New directions in hydrogeochemical transport modeling: incorporating multiple
kinetic and equilibrium pathways. In: Bentley LR, Sykes JF, Brebbia CA, Gray WG, Pinder GF(eds) Computational methods in water resources XIII. A.A. Balkema, Rotterdam, pp 331–338
Steefel CI (2009) CrunchFlow: software for modeling multicomponent reactive flow and transport,
user’s manual, 12 Oct 2009
Steefel CI, MacQuarrie KTB (1996) Approaches to modeling of reactive transport in porous
media. In: Lichtner PC, Steefel CI, Oelkers EH (eds) Reactive transport in porous media.
Reviews in mineralogy, Chapter 2, vol 34. Mineralogical Society of America, Chantilly, VA,
pp 83–129
Steefel CI, DePaolo DJ, Lichtner PC (2005) Reactive transport modeling: an essential tool and a
new research approach for the earth sciences. Earth Planet Sci Lett 240:539–558
Thaysen EM, Jacques D, Jessen S, Andersen CE, Laloy E, Ambus P, Postma D, Jakobsen I (2014a)
Controls on carbon dioxide fluxes across the unsaturated zone of cropped and unplanted soil
mesocosms. Biogeosciences (in review)Thaysen EM, Jessen S, Postma D, Jakobsen R, Jacques D, Ambus P, Laloy E, Jakobsen I (2014b)
Plant-mediated transfer of CO2 to aquifers as influenced by crushed concrete waste and lime
amendments. Environ Sci Technol (in review)402

J. Šimůnek et al.

Thullner M, Van Cappellen P, Regnier P (2005) Modeling the impact of microbial activity on
redox dynamics in porous media. Geochim Cosmochim Acta 69:5005–5019
Toride N, Leij FJ, van Genuchten MT (1993) A comprehensive set of analytical solutions fornonequilibrium solute transport with first-order decay and zero-order production. Water
Resour Res 29:2167–2182
United Nations Environment Programme (UNEP) (2002) Global mercury assessment. United
Nations Environment Programme (UNEP), Geneva
VanBriesen J (1998) Modeling coupled biogeochemical processes in mixed waste systems, Ph.D.
thesis, Civil Engineering. Northwestern University
Van der Lee J, De Windt L, Lagneau V, Goblet P (2003) Module-oriented modeling of reactive
transport with HYTEC. Comput Geosci 29:265–275
van Genuchten MTh (1980) A closed-form equation for predicting the hydraulic conductivity ofunsaturated soils. Soil Sci Soc Am J 44:892-898
van Genuchten MTh (1981) Analytical solutions for chemical transport with simultaneous adsorption, zero-order production and first-order decay. J Hydrol 49:213-233
van Genuchten MT, Šimůnek J, Leij FL, Toride N, Sejna M (2012) STANMOD: model use,
calibration and validation, special issue standard/engineering procedures for model calibrationand validation. Trans ASABE 5(4):1353–1366
van Genuchten MTh et al (2014) The use of numerical flow and transport models in environmental
analyses. In: Teixeira WG, Ceddia MB, Ottoni MV, Donnagema GK (eds) Application of soil
physics in environmental analyses: measuring, modelling and data integration. Springer, Cham
Viswanathan HS, Robinsons BA, Valocchi AJ, Triay IR (1998) A reactive transport model of
neptunium migration from potential repository at Yucca Mountain. J Hydrol 209:251–280
Walsh SDC, Frane WLD, Mason HE, Carroll SA (2013) Permeability of wellbore-cement
fractures following degradation by carbonated brine. Rock Mech Rock Eng 46:455–464
Walter AL, Frind EO, Blowes DW, Ptacek CJ, Molson JW (1994) Modeling of multicomponent
reactive transport in groundwater 1. Model development and evaluation. Water Resour Res 30(11):3137–3148
Xu T, Spycher N, Sonnenthal E, Zheng L, Pruess K (2012) TOUGHREACT user’s guide: asimulation program for non-isothermal multiphase reactive geochemical transport in variably
saturated geologic media, version 2.0. Report LBNL-DRAFT. Lawrence Berkeley National
Laboratory, Berkeley
Yeh GT (1981) Analytical transient one-, two-, and three-dimensional simulation of wastetransport in the aquifer system, ORNL-5602. Oak Ridge National Laboratory, Oak Ridge
Yeh G-T, Cheng H-P (1999) 3DHYDROGEOCHEM: a 3-dimensional model of densitydependent subsurface flow and thermal multispecies-multicomponent HYDROGEOCHEMical transport. EPA/600/R-98/159, 150 pp
Yeh GT, Tripathi VS (1991) A model for simulating transport of reactive multispecies components: model development and demonstration. Water Resour Res 27(12):3075–3094
Yuan F, Meixner T, Fenn ME, Šimůnek J (2011) Impact of transient soil water simulation on
estimated nitrogen leaching and emission at high- and low-deposition forest sites in Southern
California. J Geophys Res 116:G03040. doi:10.1029/2011JG001644, 15 pp
Zhang H, Nordin NA, Olsen MS (2013) Evaluating the effects of variable water chemistry on
bacterial transport during infiltration. J Contam Hydrol 150:54–64
Zyvoloski GA, Robinson BA, Dash ZV, Trease LL (1997) Summary of the models and methods
for the FEHM application – a finite element heat- and mass-transfer code. Los Alamos National
Laboratory Rept. LA-13307-MS, Los Alamos