JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATIONDECEMBER

AMERICAN WATER RESOURCES ASSOCIATION

2006

MULTICOMPONENT GEOCHEMICAL TRANSPORTMODELING USING HYDRUS-1D AND HP11

Jirka S̆imůnek, Diederik Jacques, Martinus Th. van Genuchten, and Dirk Mallants2

ABSTRACT: The transport of reactive contaminants in the subsurface is generally affected by a large number of nonlinear andoften interactive physical, chemical, and biological processes.
Simulating these processes requires a comprehensive reactive
transport code that couples the physical processes of waterflow and advective-dispersive transport with a range of biogeochemical processes. Two recently developed coupled geochemical models that are both based on the HYDRUS-1D software
package for variably saturated flow and transport are summarized in this paper. One model resulted from coupling HYDRUS1D with the UNSATCHEM module. While restricted to major ion
chemistry, this program enables quantitative predictions ofsuch problems as analyzing the effects of salinity on plant
growth and the amount of water and amendments required toreclaim salt-affected soil profiles. The second model, HP1,
resulted from coupling HYDRUS-1D with the PHREEQC biogeochemical code. The latter program accounts for a wide range of
instantaneous or kinetic chemical and biological reactions,
including complexation, cation exchange, surface complexation, precipitation dissolution and/or redox reactions. The versatility of HP1 is illustrated in this paper by means of two
examples: the leaching of toxic trace elements and the transport of the explosive TNT and its degradation products.
(KEY TERMS: vadose zone; transport and fate; transport models; heavy metals; multicomponent geochemical transport;
reclamation models.)S̆imůnek, Jirka, Diederik Jacques, Martinus Th. van Genuchten,
and Dirk Mallants, 2006. Multicomponent Geochemical Transport Modeling Using HYDRUS-1D and HP1. Journal of theAmerican Water Resources Association (JAWRA) 42(6):15371547.

INTRODUCTION
A multitude of interactive physical, chemical, andbiological processes generally affects the fate andtransport of reactive agricultural and industrial contaminants, as well as many naturally occurring elements, in the subsurface. Physical processes in thesubsurface may include capillary and gravitational
forces that cause transient water flow; advection, dispersion and diffusion in the liquid and gaseous phasesthat lead to solute transport; and convection and conduction that result in heat transport. The soil solutionis always a mixture of many ions that may be
involved in mutually dependent chemical processes
such as complexation reactions, cation exchange, precipitation-dissolution, sorption-desorption, volatilization, redox reactions, and degradation (S̆imůnek andValocchi, 2002; Appelo and Postma, 2005). The transport and transformation of many contaminants is further mediated by subsurface aerobic or anaerobic
bacteria. Bacteria catalyze redox reactions in whichorganic compounds (e.g., hydrocarbons) act as theelectron donor and inorganic substances (oxygen,nitrate, sulfate, or metal oxides) as the electron acceptor. By catalyzing such reactions, bacteria gain energy
and organic carbon to produce new biomass. The
transport of bacteria is generally more complex thanthat of other solutes in that they are affected by such
additional processes as filtration, straining, sedimentation, adsorption, and desorption, growth, and inactivation (e.g., Bradford et al., 2003; van Genuchten andS̆imůnek, 2004). These and related processes can be

1Paper No. 05144 of the Journal of the American Water Resources Association (JAWRA) (Copyright © 2006). Discussions are open until
June 1, 2007.
2Respectively, Professor and Hydrologist, Department of Environmental Sciences, University of California, Riverside, California 92521;
Research Scientist, Waste and Disposal Department, SCK•CEN, Boeretang 200, B-2400 Mol, Belgium; Senior Research Scientist, George E.
Brown, Jr., Salinity Laboratory, 450 West Big Spring Road, Riverside, California 92507; and Research Scientist, Waste and Disposal Department, SCK•CEN, Boeretang 200, B-2400 Mol, Belgium (E-Mail/Simunek: Jiri.Simunek@ucr.edu).

JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATION

1537

JAWRA

S̆IMŮNEK, JACQUES, VAN GENUCHTEN, AND MALLANTS
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 interpool transfer of carbon andnitrogen, 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.
On the other hand, models with generalized chemistry provide users with much more freedom indesigning particular chemical systems and thus also
permit a much broader range of applications. Users
can either select species and reactions from large geochemical databases or define their own species withparticular chemical properties and reactions. Thesemodels generally invoke a variety of sequential iterative or noniterative operator splitting approaches inwhich transport and chemistry are numerically solved
in separate steps. Examples of models with generalized chemistry have been recently reviewed byS̆ imůnek and Valocchi (2002) and Jacques et al.
(2006). Most codes with general geochemistry are limited to solute transport and biogeochemical reactions,while water flow paths have to be calculated outsideof the reactive transport code (e.g., PHREEQC,
Parkhurst and Appelo, 1999; CRUNCH, Steefel, 2000;
PHAST, Parkhurst et al., 2004). Only a few models
allow the velocity field to be calculated internally.
Several codes for transient unsaturated flow havealso been coupled to general biogeochemistry models.
These include, for example, 3DHYDROGEOCHEM(Yeh and Cheng, 1999) and HP1 (Jacques andS̆imůnek, 2005).
The HYDRUS-1D (Version 3.0) (S̆imůnek et al.,2005) and HP1 models (Jacques and S̆imůnek, 2005)
that represent typical examples of models with specific and generalized chemistry, respectively, are discussed below.

simulated using integrated reactive transport codesthat couple the physical processes of water flow andadvective-dispersive solute transport with a range ofbiogeochemical processes.
In this paper two recently developed coupled geochemical models are summarized that are both basedon the HYDRUS-1D software package for simulating
variably saturated flow and transport (S̆imůnek etal., 1998). The first program, which resulted from coupling HYDRUS-1D with the UNSATCHEM module(S̆ imůnek et al., 1996), is restricted to major ions
(Ca 2+, Mg 2+, Na +, K +, Cl - , SO 42-, and HCO 3- , andtheir complexes). This program enables quantitative
predictions of such problems as analyzing the effectsof salinity on plant growth, water and salinity management practices and irrigation techniques underarid and semiarid conditions, evaluation of water suitability for irrigation, and estimating the amount ofwater and amendments required to reclaim saltaffected soil profiles. The other program (HP1) resulted from coupling HYDRUS-1D with the PHREEQC
geochemical code (Parkhurst and Appelo, 1999). This
program accounts for a wide range of instantaneousand kinetic chemical and biological reactions. Several
examples will be given to demonstrate the utility ofboth programs.

COUPLED GEOCHEMICAL MODELSAs summarized by S̆imůnek and Valocchi (2002),
geochemical transport models can be divided into two
major groups: those with specific chemistry and moregeneral models. Models with specific chemistry are
generally restricted to certain prescribed chemicalsystems and thus are usually constrained to very specific applications. However, they are often much easier to use and computationally much more efficientthan general models. Models simulating the transportof major ions, such as LEACHM (Wagenet and Hutson, 1987) and UNSATCHEM (S̆imůnek and Suarez,
1994; S̆imůnek et al., 1996), and various reclamation
models (S̆ imůnek and Valocchi, 2002) are typicalexamples of models with specific chemistry. Thesemodels typically consider the transport of major ionsand their mutual reactions such as complexation,
cation exchange, and precipitation/dissolution. Models
simulating carbon and nitrogen cycles are also becoming a standard feature of many environmental models, for example in DAISY (Hansen et al., 1990),
LEACHM (Hutson and Wagenet, 1992), RZWQM(Ahuja and Hebson, 1992), COUP (Jansson andKarlberg, 2001), and HYDRUS-2D coupled withCW2D (Langergraber and S̆ imůnek, 2005). Thesemodels typically distribute organic matter, carbon,
JAWRA

Overview of HYDRUS-1DHYDRUS-1D (S̆ imůnek et al., 2005) is a publicdomain computer software package (www.hydrus2d.
com) that simulates the one-dimensional movement ofwater, heat, and multiple solutes in variably saturated media. The program uses finite elements to numerically solve the Richards equation for saturated
unsaturated water flow and Fickian-based advectiondispersion equations for both heat and solute transport. The flow equation includes a sink term toaccount for water uptake by plant roots as a functionof both water and salinity stress. The heat transport
equation considers conduction as well as advection
with flowing water. The solute transport equations
1538

JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATIONMULTICOMPONENT GEOCHEMICAL TRANSPORT MODELING USING HYDRUS-1D AND HP1assume advective-dispersive transport in the liquid
phase and diffusion in the gaseous phase. The transport equations also include provisions for nonlinear
and/or nonequilibrium reactions between the solidand liquid phases, linear equilibrium reactions
between the liquid and gaseous phases, zero-orderproduction, and two first-order degradation reactions– one that is independent of other solutes and onethat provides the coupling between solutes involved in
the sequential first-order decay reactions. The program allows analysis of flow and transport in both
unsaturated, partially saturated, or fully saturated
media. The flow region itself may consist of nonuniform (layered) soils. The unsaturated soil hydraulicproperties (the constitutive relationships) aredescribed using Brooks and Corey (1964), vanGenuchten (1980), Durner (1994), and Kosugi (1996)
type analytical functions; Durner (1994) uses modified van Genuchten type functions that permit a better description of the hydraulic properties nearsaturation. The software optionally also considers
hysteresis by assuming that drying scanning curves
are scaled from the main drying curve and wetting
scanning curves from the main wetting curve.
Root growth in HYDRUS-1D is simulated by meansof a logistic growth function. Water and salinity
stress-response functions also can be considered.
Physical nonequilibrium solute transport can be
accounted for by assuming a two-region, dual-porosity
type formulation that partitions the liquid phase into
separate mobile and immobile regions. HYDRUS-1D
was recently modified to also allow simulations of
the transport of viruses, colloids, and bacteria (Schijven and S̆imůnek, 2002; Bradford et al., 2003). In
addition, HYDRUS-1D implements a MarquardtLevenberg type parameter estimation technique(Marquardt, 1963; S̆imůnek and Hopmans, 2002) for
inverse estimation of soil hydraulic (Hopmans et al.,2002) and/or solute transport and reaction (S̆imůnek
et al., 2002) parameters from measured transient or
steady state flow and/or transport data. The main
new features of the current version 3.0, as comparedto version 2.0, are new analytical models for the soil
hydraulic properties (i.e., the Kosugi and Durner analytical models), provisions for compensated root water
uptake, dual porosity features for water flow, attachment/detachment models for solute transport to allow
simulations of virus, colloid, and bacteria transport,
and modules for carbon dioxide dynamics and for
major ion chemistry.
HYDRUS-1D uses a Microsoft Windows-basedgraphics user interface (GUI) to manage the input
data required to run the program, as well as for nodal
discretization and editing, parameter allocation, problem execution, and visualization of results. All spatially distributed parameters – such as those for
JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATION

various soil horizons, the root water uptake distribution, and the initial conditions for water, heat, andsolute movement – are specified in a graphical environment. The program offers graphs of the distributions of the pressure head, water content, water andsolute fluxes, root water uptake, temperature, andsolute concentrations in the soil profile at pre-selected
times. Also included are a small catalog of unsaturated soil hydraulic properties (Carsel and Parish, 1988)
and a set of pedotransfer functions based on neuralnetworks (Schaap et al., 2001).
The general formulation of the transport equationsin the HYDRUS codes allows one to simulate not only
non-adsorbing and linearly or nonlinearly sorbing
chemicals but also a variety of other contaminants
such as viruses (Schijven and S̆imůnek, 2002), colloids (Bradford et al., 2003), and cadmium (Seuntjens
et al., 2001) or chemicals involved in the sequential
biodegradation reactions such as chlorinated aliphatichydrocarbons (Schaerlaekens et al., 1999; Casey andS̆imůnek, 2001) or hormones (Casey et al., 2003).
Specific (Major Ion) Chemistry Module of
HYDRUS-1DThe solute transport module of HYDRUS-1D as
described above is limited to single ions or ions subject to relatively simple first-order consecutive decay
reactions (e.g., nitrification-denitrification chains, orradionuclide decay chains). As an alternative to thisapproach, we recently implemented into HYDRUS-1D
a major ion chemistry module based on the
UNSATCHEM model (S̆ imůnek and Suarez, 1994,
1997; S̆imůnek et al., 1996, 2005). This module considers the transport of major ions and carbon dioxide
in soils. CO2 transport is assumed to be governed bydiffusion in both the liquid and gas phases, by advection in the liquid phase, and by respiration by bothsoil microorganisms and plant roots.
The major variables of the UNSATCHEM chemicalsystem are Ca, Mg, Na, K, SO 4, Cl, NO 3, H 4SiO 4,alkalinity, and CO2 (Table 1). The model accounts forequilibrium chemical reactions between these components such as complexation, cation exchange, andprecipitation-dissolution. For the precipitationdissolution of calcite (CaCO3) and the dissolution of
dolomite (MgCO 3 ·3H 2 O), either equilibriumor multicomponent kinetic expressions can beused, including both forward and backward reactions. Other precipitation-dissolution reactions considered involve gypsum (CaSO4·2H2O), hydromagnesite (Mg 5 (CO 3 ) 4 (OH) 2 ·4H 2 O), nesquehonite(MgCO3·3H2O), and sepiolite (Mg2Si3O7.5(OH)·3H2O).
Since the ionic strength of soil solutions can vary considerably in time and space and often reach high
1539

JAWRA

S̆IMŮNEK, JACQUES, VAN GENUCHTEN, AND MALLANTSTABLE 1. Chemical Species Included in the Major Ion Module of HYDRUS-1D.
Aqueous Components

7Ca2+, Mg2+, Na+, K+, SO42-, Cl-, NO3-

Complexed Species

10

CaCO3o, CaHCO3+, CaSO4o, MgCO3o, MgHCO3+, MgSO4o, NaCO3-, NaHCO3o, NaSO4-, KSO4-

Precipitated Species

6

CaCO3, CaSO4·2H2O, MgCO3·3H2O, Mg5(CO3)4(OH)2·4H2O, Mg2Si3O7.5(OH)·3H2O,
CaMg(CO3)2

Sorbed (exchangeable) Species

4Ca, Mg, Na, K

CO2-H2O Species

7

PCO2, H2CO3*, CO32-, HCO3-, H+, OH-, H2O

Silica Species

3

H4SiO4, H3SiO4-, H2SiO42-

values, both the modified Debye-Hückel (Truesdell
and Jones, 1974) and Pitzer (Pitzer, 1979) expressions
were incorporated into the model, thus providing
options for calculating single-ion activities. This module also considers the effects of solution compositionon the unsaturated soil hydraulic properties. The new
UNSATCHEM module of HYDRUS-1D enables quantitative predictions of processes involving major ions,such as simulations of the effects of salinity on root
water uptake and plant growth, evaluation of alternative irrigation, salinity and crop management practices, evaluation of water suitability for irrigation,and estimation of the amount of water and amendment required to reclaim soil profiles to desired levelsof salinity and exchangeable sodium percentage(ESP). A recent application of HYDRUS-1D and its
major ion chemistry module is given by Gonçalves etal. (2006), who simulated solute transport in threelysimeters irrigated with different-quality waters overa time period of three years. HYDRUS-1D successfully described field measurements not only of the overall salinity but also of individual soluble cations aswell as sodium adsorption ratio (SAR) and ESP(Gonçalves et al., 2006).
A typical application of the above UNSATCHEM
major ion chemistry module within HYDRUS-1D issodic soil reclamation involving the removal of excessive salts from the soil profile and replacement of
exchangeable Na+ with Ca2+. This problem is relatively complicated in that one must simulate not only
variably saturated water flow and solute transport
but also the dissolution of solid phases (i.e., of reclamation amendments such as gypsum or calcite),cation exchange (replacement of Na+ with Ca2+), complexation reactions, the effects of solution composition
on hydraulic conductivities and infiltration rates, andthe effects of elevated CO2 concentrations on pH andhence on calcite dissolution (S̆ imůnek and Suarez,
1997). We briefly summarize here results of a hypothetical sodic soil reclamation study by S̆imůnek andSuarez (1997), who considered six reclamation strategies. In the first two cases, (A) high quality (i.e., low
ion concentration) and (B) gypsum saturated waters
JAWRA

were applied to the soil without additional amendments. In the third example, (C) high quality water
was applied to a soil profile in which gypsum wasincorporated in the top 20 cm. In three other cases,(D) high quality water or water having one of two different acid concentrations at (E) pH = 2.05 and (F) at
pH = 1.09 were applied to the surface of a calcareous
soil. The solution composition of water initially present in the soil profile was assumed to be that of atypical highly sodic soil system. The cation exchange
capacity was set at 200 mmolc/kg and the ESP at 60(Ca = 40.0, Mg = 40.0, Na = 120.0 mmolc/kg). The
solution composition of four irrigation waters used forthe reclamation process is given in Table 2. Modelpredicted CO2 concentrations were used for the last
three cases, where CO2 has an important effect on
calcite dissolution. All calculations consideredhydraulic conductivity reductions due to changes inthe solution composition. The example assumes thatthe soil surface was flooded with water, with thewater kept at a constant level of 1 cm above the soil
surface.
Figure 1 shows calculated sodium adsorption ratio,
SAR {=Na/[(Ca+Mg)/2]0.5}, profiles for all six reclamation strategies. The various scenarios were evaluatedin terms of the time required for reclamation to becompleted and the amount of water needed for this.
Times to full profile saturation (saturation equal toporosity) for the six strategies above were calculated
to be 358, 15.3, 2.0, 10, 6.4, and 3.2 days, respectively.
Times to profile reclamation were 86,870, 195.6,
116.5, 91, 64, and 15.5 days, respectively, while theamounts of infiltrated water needed for reclamation
were 7,230, 559, 470, 600, 485, and 197 cm, respectively (S̆imůnek and Suarez, 1997). The first strategy
was included only as a reference to compare withother scenarios, mainly to emphasize the effects ofgypsum and calcite dissolution on reclamation. This
strategy is clearly inappropriate for sodic soil reclamation, since high quality water does not have a sufficient supply of calcium to replace the exchangeablesodium and as such can lead to soil dispersion andsubsequent reductions in the hydraulic conductivity.
1540

JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATIONMULTICOMPONENT GEOCHEMICAL TRANSPORT MODELING USING HYDRUS-1D AND HP1
TABLE 2. Solution Composition of Waters Used in the Reclamation Example (in mmolcL-3).

Water

Used in
Scenario

CaT

MgT

NaT

ClT

SO4T

Alkalinity

High Quality

1.5

0.5

2.0

1.0

2.5

0.5

A,C,D

Gypsum Saturated

32.0

0.5

2.0

1.0

33.0

0.5

B

Acid Water I

1.5

0.5

2.0

3.0

11.0

-10.0

E

Acid Water II

1.5

0.5

2.0

93.0

11.0

-100.0F

Initially Present in the Soil

0.2

0.2

4.8

4.8

0.0

0.4

A,B,C,D,E,F

Figure 1. Sodium Adsorption Profiles for Six Reclamation Strategies. (A) Irrigation with high quality water and no amendments;(B) Irrigation with gypsum saturated water; (C) Irrigation with high quality water and gypsum incorporated in
the top 20 cm; (D) Irrigation with high quality water and and with calcite present in the entire soil profile;(E) Irrigation with acid water at ph 2.05 and with calcite in the entire soil profile; and (F) Irrigation withacid water at pH = 1.09 and with calcite in the entire soil profile (after S̆imůnek and Suarez, 1997).

Since sodic soil reclamation is a relatively complex
problem, only models that quantitatively integrate allof the major processes involved should be used, rather
JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATION

than resorting to highly simplified models based on
single-ion or decay chain transport (S̆imůnek and Valocchi, 2002).
1541

JAWRA

S̆IMŮNEK, JACQUES, VAN GENUCHTEN, AND MALLANTS
The HP1 General Chemistry Modelcadmium leaching in acid sandy soils; radionuclidetransport; and long-term uranium migration in agricultural field soils following mineral P-fertilization.
The first example is briefly discussed below, together
with a new application that involves the fate and subsurface transport of explosives.

HYDRUS-1D was recently coupled also with the
PHREEQC geochemical code (Parkhurst and Appelo,1999) to create a new comprehensive simulation tool,
HP1 (acronym for HYDRUS1D-PHREEQC) (Jacqueset al., 2003; Jacques and S̆imůnek, 2005; S̆imůnek etal., 2005). This new code contains modules simulatingtransient water flow in variably saturated media, thetransport of multiple components, mixed equilibrium/kinetic biogeochemical reactions, and heat transport. HP1 becomes a significant expansion of theindividual HYDRUS-1D and PHREEQC programs bypreserving most of their original features and capabilities. The code still uses the Richards equation
for simulating variably saturated water flow andadvection-dispersion type equations for heat andsolute transport. However, the loosely coupled program can now simulate also a broad range of low-temperature biogeochemical reactions in water, thevadose zone and in ground water systems, includinginteractions with minerals, gases, exchangers, andsorption surfaces based on thermodynamic equilibrium, kinetic, or mixed equilibrium-kinetic reactions.
HP1 uses the operator splitting approach with no iterations during one time step (noniterative sequential
approach). Jacques et al. (2006) evaluated the accuracy of the operator splitting approach for a kinetic
reaction network (i.e., sequential and parallel kinetic
degradation reactions) by comparing HP1 with ananalytical solution for trichloroethylene (TCE) degradation and for mixed equilibrium and kinetic reactions involving different flow conditions (steady state
and transient). HP1 is available on request fromwww.sckcen.be/hp1. The current version of HP1 islimited to one-dimensional flow and transport problems and does not consider oxidation reduction reactions, surface complexation, and diffusion ofcomponents in the gas phase. These reactions will beavailable in upcoming versions currently under development (Jacques et al., 2005b).
Jacques et al. (2003, 2005a) and Jacques andS̆imůnek (2005) demonstrated the versatility of HP1
on several examples such as the transport of heavy
metals (Zn 2+ , Pb 2+ , and Cd 2+ ) subject to multiple
cation exchange reactions; transport with mineral dissolution of amorphous SiO2 and gibbsite, Al(OH)3;
heavy metal transport in a medium with a pH dependent cation exchange complex; infiltration of a hyperalkaline solution in a clay sample (this example
considers kinetic precipitation-dissolution of kaolinite,
illite, quartz, calcite, dolomite, gypsum, hydrotalcite,and sepiolite); long term transient flow and transport of major cations (Na+, K+, Ca2+, and Mg2+) andheavy metals (Cd2+, Zn2+, and Pb2+) in a soil profile;
JAWRA

Leaching of Heavy Metals From a Soil Column. As a first example application of HP1, Figure 2
shows calculated effluent concentrations of major ions
(Ca2+, Na+, Al3+, and Cl-) and heavy metals (Zn2+,
Pb2+, and Cd2+) leached from a soil column 8 cm long
having an initial solution as defined in Table 3 andwith its ion exchange complex in equilibrium withthis solution. Heavy metals initially on the exchange
complex are being replaced during the leaching process by major ions (Ca2+, Mg2+, and Al3+). Ca richwater with a chemical composition as given in Table 3
was applied to the top of the column at a steady flowrate of 2 cm day-1. The longitudinal dispersivity was
assumed to be 0.2 cm.

Figure 2. Concentrations of (a) Major Ions and (b) Heavy Metals(Zn, Pb, and Cd) in the Effluent From an 8 cm Long Soil Column.

The example considered the transport of ten ions
(Al, Br, Ca, Cd, Cl, K, Mg, Na, Pb, Zn) during a 15-day
period. Removal of Cd, Pb and Zn from the solid
1542

JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATIONMULTICOMPONENT GEOCHEMICAL TRANSPORT MODELING USING HYDRUS-1D AND HP1
TABLE 3. Chemical Components (mmol/L) and Species Considered, and Elemental Compositions of Initial
and Boundary Solutions Used in the Column Simulation (X refers to ion exchanger).

Components

Solutions
Boundary
InitialSpeciesAl

Al3+, Al(OH)2+, Al(OH)2+, Al(OH)30, Al(OH)4-

0.1

0.5

Br

Br-

3.7

11.9

Cl

Cl- (and Cd, Pb, and Zn-species, see below)

10

0

Ca

Ca2+, Ca(OH)+

5

0

K

K+, KOH00

2Na

Na+, NaOH0

0

6

Mg

Mg2+, Mg(OH)+

1

0.75

Cd

Cd2+, Cd(OH)+, Cd(OH)20, Cd(OH)3-, Cd(OH)42-, CdCl+, CdCl2, CdCl3-

0

0.9

Pb

Pb2+, Pb(OH)+, Pb(OH)20, Pb(OH)3-, Pb(OH)42-, PbCl+, PbCl20, PbCl3-, PbCl42-

0

0.1

Zn

Zn2+, Zn(OH)+, Zn(OH)20, Zn(OH)3-, Zn(OH)42-, ZnCl+, ZnCl20, ZnCl3-, ZnCl42-

0

0.25

X

AlX3, AlOHX2, CaX2, CdX2, KX, NaX, MgX2, PbX2, ZnX2 (mmol)NA11phase was calculated to occur by ion exchange withmostly Ca and to a lesser degree with Mg. The 8 cmlong soil column for the numerical calculations wasdiscretized into 40 elements of 0.2 cm each. Moredetails about numerical stability and its dependence
on the selected temporal and spatial discretization,
the cation exchange complex, and the solution precision are given by Jacques and S̆imůnek (2005).

among other places (Pennington et al., 2001; Brannon
and Pennington, 2002). They often leach through thevadose zone into ground water and thus threaten
underlying or downgradient water resources. The costto complete the current U.S. Department of Defense(DoD) cleanup program using traditional approaches
is estimated to be more than $35 billion (Pennington
et al., 1999a).
Previous efforts to model the transport of major
explosives have been based mostly on either Modflow(Harbaugh et al., 2000) combined with MT3D (Zheng
and Wang, 1999) for ground water systems or
FEMWATER (Yeh et al., 1992) for vadose zone type
investigations (Pennington et al., 1999b). Neither ofthese codes at present accounts for all of the processes
known or hypothesized to affect the fate and subsurface transport of explosives and their metabolites,such as nonlinear and kinetic dissolution and sorption, simultaneous transport of multiple metabolites,
simultaneous transport in both the liquid and vapor
phase, and the temperature dependence of transportand reaction parameters. Most or all of these processes can be simulated with HP1 and other biogeochemical models such as 3DHYDROGEOCHEM or
CRUNCH – notably the transport of the major explosives themselves (the parent products) such as TNT,
RDX, or HMX but also the various metabolites (thedaughter products) such as 2ADNT, 2,4DANT,
1,3,5TNB as they are being created sequentionally bydegradation of the parent compounds. HP1 further
allows the parent and daughter compounds to havedifferent mobilities in the subsurface as dictatedby their specific dissolution, sorption and transport

Transport of TNT and Its Daughter Products.
The second example deals with the fate and transportof explosive species. Concern is increasing about thepresence of explosives and energetics in the subsurface environment. Such chemicals are the result of
the manufacture, distribution, testing and/or unsafe
disposal of ammunition. Contamination of soil andwater resources with residues of major explosives andenergetics now constitutes a widespread environmental problem. These include TNT (2,4,6-trinitrotoluene); RDX (hexahydro-1,3,5- trinitro-1,3,5triazine); HMX (octahydro-1,3,5,7-tetranitro-1,3,5,7tetrazocine); picrid acid; tetryl, hexyl, nitroglycerin;
and PETN; as well as their principal transformation
products (notably 2ADNT, 4ADNT, 2,4DANT,
2,6DANT, 1,3,5TNB, 1,3DNB, NB, and 3,5DNA) andsome of their manufacturing impurities (2,4DNT and
2,6DNT) (Gorontzy et al., 1994; Brannon and Pennington, 2002). Many of these compounds are known
or suspected to be potentially toxic, mutagenic, or carcinogenic. The chemicals have been found not only atammunition production and handling facilities in the
United States and worldwide but also at many disposal sites and army testing, training, and firing ranges,
JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATION

1543

JAWRA

S̆IMŮNEK, JACQUES, VAN GENUCHTEN, AND MALLANTS
properties. In its most general case, the model permits contaminants to reside in all three phases, i.e.,
the liquid, solid (precipitated and sorbed), andgaseous phases. HP1 allows for multiple degradation
pathways that involve both diverging and converging
branches. Since users of HP1 can define any kinetic
reaction through its Basic interpreter, there is virtually no limitation on the complexity of processes calculated. An example of an initially diverging and thenconverging branch is TNT degradation, which may
proceed along two branches before reaching the compound 2,4,6-triminotoulene (TAT) by reductive biotransformation (Figure 3).

Figure 3. Degradation Pathway of TNT (2,4,6-trinitrotoluene)
Through Two Degradation Products 2ADNT (2-amino-4,
6-dinitrotoluene) and 4ADNT (4-amino-2,6-dinitrotoluene)
to TAT (2,4,6-triminotoulene).

The example here considers the fate and transportof four explosive species in a 100 cm deep fully saturated soil profile leached for 10 days at a water flux
rate of 1 cm/h and using a longitudinal dispersivity
equal to 5 cm. TNT, the parent product, was assumed
to be initially present in the top 5 cm of the soil at aconcentration of 1 mg/kg (6.61e-6 mol). The kinetic
dissolution rate (mol/h) of TNT is given as[TNT]
r = Ak 1 −K Figure 4. Breakthrough Curves at Depths of 50 and
100 cm for TNT, 2ADNT, 4ADNT, and TAT.

(1)simulated concentration distributions versus depth.
Although TAT is the last daughter product, it arrived
at the 50 and 100 cm depths before the other species
as a result of having the smallest distribution coefficient and due to the relatively short half-life of
4ADNT (T 1/2 = 17.3 h) (Figure 4). Because of themuch longer half-life of 2ADNT (T1/2 = 115.5 h), mostof the observed TAT present at early times is due todegradation of 4ADNT. And while concentrations of
TNT were initially highest (e.g., after 2 d in Figure 5)
and then decreased due to degradation, concentrations of all daughter species (2ADNT, 4ADNT, andTAT) always increased during the 10-day simulation.
Because of less sorption, TAT moved fastest throughthe soil profile, followed by TNT, and then 2ADNT
and 4ADNT at very similar velocities. This example
demonstrates that ground water may be more vulnerable to leaching of TNT daughter products (notably

where A is the reactive surface of TNT (assumed to be
1 cm2 and kept constant during the simulation), k is
the dissolution rate (equal to 1.8e-5 mol/cm2/h), andK is the equilibrium constant of TNT dissolution,log(K) = -5.356. When dissolved, TNT was assumed todegrade at a rate of 0.01 1/h to form the daughter
products 2ADNT (2-amino-4,6-dinitrotoluene) (66.7
percent) and 4ADNT (4-amino-2,6-dinitrotoluene)(33.3 percent), which further degrade at rates of 0.006
and 0.04 1/h, respectively, to produce the same final
daughter product TAT (Pennington et al., 1999b). The
distribution coefficients, KD, for linear sorption wereassumed to be equal to 3, 5, 6, and 0 L/kg for TNT,
2ADNT, 4ADNT, and TAT, respectively (Pennington etal., 1999b). Figure 4 shows calculated breakthroughcurves at 50 and 100 cm depths, and Figure 5 shows
JAWRA

1544

JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATIONMULTICOMPONENT GEOCHEMICAL TRANSPORT MODELING USING HYDRUS-1D AND HP1

Figure 5. Concentration Profiles for TNT, 2ADNT, 4ADNT, and TAT (0 to 10 days)

relatively easy to use, but rather restrictive in termsof its range of applications – HP1 (HYDRUS1DPHREEQC) accounts for a very wide range of
instantaneous and kinetic chemical and biologicalreactions, including complexation, cation exchange,
surface complexation, precipitation dissolution, and/or
redox reactions, thus requiring users to be morefamiliar with the chemistry involved. Capabilities ofthe major ion module of HYDRUS-1D were demonstrated by simulating sodic soil reclamation involving
the removal of excessive salts from a soil profile andreplacement of exchangeable Na+ with Ca2+. Capabilities of HP1 were shown by means of example problems involving the leaching of heavy metals subject toTAT) than of the parent compound itself and that sensors detecting daughter products may provide an
early warning of possible TNT leaching.

SUMMARY AND CONCLUSIONS
Two new comprehensive simulation tools based on
the HYDRUS-1D software package were presented inthis paper. While the UNSATCHEM module (directly
incorporated into HYDRUS-1D) is restricted to a specific chemical system involving the transport andreactions of major ions – thus making the program
JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATION

1545

JAWRA

S̆IMŮNEK, JACQUES, VAN GENUCHTEN, AND MALLANTS
multiple cation exchange reactions and the transportof TNT and its daughter products 2ADNT, 4ADNT,and TAT. This example indicates that ground water
may be more vulnerable to leaching of TNT daughter
products, notably TAT, than of the parent compound
itself and that monitoring for the daughter products
may provide an early warning of possible TNT leaching.

Hansen, S., H.E. Jensen, N.E. Nielsen, and H. Svendsen, 1990.
DAISY: Soil Plant Atmosphere System Model. NPO Report No.
A 10, The National Agency for Environmental Protection,
Copenhagen, 272 pp.
Harbaugh, A.W., E.R. Banta, M.C. Hill, and M.G. McDonald, 2000.
MODFLOW-2000: The U.S. Geological Survey Modular GroundWater Model – User Guide to Modularization Concepts and theGround-Water Flow Process. U.S. Geological Survey Open-File
Report 00-92, 121 pp.
Hopmans, J.W., J. S̆ imůnek, N. Romano, and W. Durner, 2002.
Inverse Modeling of Transient Water Flow (Third Edition). In:
Methods of Soil Analysis, Part 1, Physical Methods, J.H. Daneand G.C. Topp (Editors). SSSA, Madison, Wisconsin, Chapter
3.6.2, pp. 963-1008.
Hutson, J.L. and R.J. Wagenet, 1992. LEACHM: Leaching Estimation And Chemistry Model. Research Series No. 92-3. Cornell
Univ., Ithaca, New York.
Jacques, D. and J. S̆imůnek, 2005. Multicomponent – Variable Saturated Transport Model, Description, Manual, Verification and
Examples, Waste and Disposal. SCK•CEN, BLG-998, Mol, Belgium.
Jacques, D., J. S̆imůnek, D. Mallants, and M. Th. van Genuchten,
2003. The HYDRUS-PHREEQC Multicomponent Transport
Model for Variably-Saturated Porous Media: Code Verification
and Application. In: MODFLOW and More 2003: Understanding
Through Modeling, E. Poeter, Ch. Zheng, M. Hill, and J. Doherty(Editors). Int. Ground Water Modeling Center, Colorado Schoolof Mines, Fort Collins, Colorado, pp. 23-27.
Jacques, D., J. S̆imůnek, D. Mallants, and M. Th. van Genuchten,
2005a. Long Term Uranium Migration in Agricultural Field
Soils Following Mineral P-Fertilization, Proc. of The 10th Int.
Conf. on Environ. Remediation and Radioactive Waste Management, Scottish Exhibition and Conference Center, Glasgow, Scotland.
Jacques, D., J. S̆imůnek, D. Mallants and M. Th. van Genuchten,
2005b. Modelling Uranium Leaching From Agricultural Soils to
Groundwater as Yardstick in Comparison With Complementary
Safety Indicators. Material Research Society 2005, Symposiumon the Scientific Basis for Nuclear Waste Management, Ghent,
Belgium.
Jacques, D., J. S̆imůnek, D. Mallants, and M. Th. van Genuchten,
2006. HP1: A Coupled Numerical Code for Variably Saturated
Water Flow, Solute Transport, and Biogeochemistry in Soil Systems. J. Contam. Hydrology (in press).
Jansson, P-E. and L. Karlberg, 2001. Coupled Heat and Mass
Transfer Model for Soil-Plant-Atmosphere Systems. Royal Institute of Technology, Dept of Civil and Environmental Engineering, Stockholm, Sweden, 325 pp.
Kosugi K., 1996. Lognormal Distribution Model for Unsaturated
Soil Hydraulic Properties. Water Resour. Res. 32(9):2697-2703.
Langergraber, G. and J. S̆ imůnek, 2005. Modeling VariablySaturated Water Flow and Multi-Component Reactive Transport in Constructed Wetlands, Vadose Zone Journal 4:924-938.
Marquardt, D.W., 1963. An Algorithm for Least-Squares Estimationof Nonlinear Parameters. SIAM J. Appl. Math. 11:431-441.
Parkhurst, D.L. and C.A.J. Appelo, 1999. User’s Guide to PHREEQC (Version 2) – A Computer Program for Speciation, BatchReaction, One-Dimensional Transport, and Inverse Geochemical
Calculations. Water Resources Investigation, Report 99-4259,
Denver, Colorado, 312 pp.
Parkhurst, D.L., K.L. Kipp, P. Engesgaard, and S.C. Charlton, 2004.
PHAST – A Program for Simulating Ground-Water Flow, Solute
Transport and Multicomponent Geochemical Reactions. USGSTechniques and Methods 6-A8, 154 pp.

ACKNOWLEDGMENT
This paper is based on work supported in part by SAHRA (Sustainability of Semi-Arid Hydrology and Riparian Areas) under the
STC Program of the National Science Foundation, Agreement No.
EAR-9876800, and the Terrestrial Sciences Program of the Army
Research Office (Terrestrial Processes and Landscape Dynamics
and Terrestrial System Modeling and Model Integration). Support
was further obtained through the bilateral agreement project“Development and Evaluation of a Coupled Geochemical Transport
Model” (Agreement No. CO-90001412.01) between SCK.CEN andthe University of California, Riverside.

LITERATURE CITED
Ahuja, L.R. and C. Hebson, 1992. Root Zone Water Quality Model.
GPSR Tech. Rep. No. 2, USDA, ARS, Fort Collins, Colorado
Appelo, C.A.J. and D. Postma, 2005. Geochemistry, Groundwaterand Pollution (Second Edition). A.A. Balkema Publ., Leiden, The
Netherlands, 649 pp.
Bradford, S.A., J. S̆imůnek, M. Bettehar, M. Th. van Genuchten,
and S.R. Yates, 2003. Modeling Colloid Attachment, Straining,and Exclusion in Saturated Porous Media, Environ. Sci. andTechnology 37(10):2242-2250.
Brannon, J.M. and J.C. Pennington, 2002. Environmental Fate andTransport Process Descriptors for Explosives. ERDC/EL TR-0210, U.S. Army Engineer Research and Development Center,
Vicksburg, Mississippi, p. 64.
Brooks, R.H., and A.T. Corey, 1964. Hydraulic Properties of PorousMedia. Hydrol. Paper No. 3, Colorado State Univ., Fort Collins,
Colorado.
Carsel, R.F. and R.S. Parrish, 1988. Developing Joint Probability
Distributions of Soil Water Retention Characteristics, Water
Resour. Res. 24:755-769.
Casey, F.X.M. and J. S̆imůnek, 2001. Inverse Analyses of the Transport of Chlorinated Hydrocarbons Subject to Sequential Transformation Reactions. Journal of Environmental Quality 30(4):1354-1360.
Casey, F.X. M., G.L. Larsen, H. Hakk, and J. S̆imůnek, 2003. Fate
and Transport of 17β-Estradiol in Soil-Water Systems. Environmental Science and Technology 37(11):2400-2409.
Durner, W., 1994. Hydraulic Conductivity Estimation for Soils With
Heterogeneous Pore Structure. Water Resour. Res. 32(9):211223.
Gonçalves, M.C., J. S̆imůnek, T.B. Ramos, J.C. Martins, M.J. Neves,and F.P. Pires, 2006. Multicomponent Solute Transport in Soil
Lysimeters Irrigated With Waters of Different Quality. Water
Resour. Res. 42, W08401, doi:10.1029/2006WR004802, p. 17.
Gorontzy, T., O. Drzyzga, M. Kahl, D. Bruns-Nagel, J. Breitung,E. von Loew, and K.-H. Biotevogel, 1994. Microbial Degradationof Explosives and Related Compounds. Critical Reviews in
Microbiology 20(4):265-284.

JAWRA

1546

JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATIONMULTICOMPONENT GEOCHEMICAL TRANSPORT MODELING USING HYDRUS-1D AND HP1
Pennington, J.C., R. Bowen, J.M. Brannon, M. Zakikhani, D.W. Harrelson, D. Gunnison, J. Mahannah, J. Clarke, T.F. Jenkins, andS. Gnewuch, 1999a. Draft Protocol for Evaluating, Selecting,and Implementing Monitored Natural Attenuation at Explosives-Contaminated Sites. Technical Report EL-99-10, U.S.
Army Engineer Research and Development Center, Vicksburg,
Mississippi, pp. 156.
Pennington, J.C., D. Gunnison, D.W. Harrelson, J.M. Brannon,
M. Zakikhani, T.F. Jenkins, J.U. Clarke, C.A. Hayes, T. Myers,
E. Perkins, D. Ringelberg, D. Townsend, H. Fredrickson, andJ.H. May, 1999b. Natural Attenuation of Explosives in Soil andWater Systems at Department of Defense Sites: Interim Report.
Technical Report EL-99-8, U.S. Army Engineer Waterways
Experiment Station, Vicksburg, Mississippi, pp. 232.
Pennington, J.C., T.F. Jenkins, J.M. Brannon, J. Lynch, T.A. Ranney,
T.E. Berry, Jr., C.A. Hayes, P.H. Miyares, M.E. Walsh, A.D.
Hewitt, N. Perron, and J.J. Delfino, 2001. Distribution and Fateof Energetics on DoD, Test and Training Ranges: Interim Report
1. ERDC TR-01-13, U.S. Army Engineer Research and Development Center, Vicksburg, Mississippi, p. 78.
Pitzer, K.S., 1979. Activity Coefficients in Electrolyte Solutions.
Chapter 7, CRC Press, Boca Raton, Florida.
Seuntjens, P., K. Tirez, J. S̆imůnek, M. Th. van Genuchten, C. Cornelis, and P. Geuzens, 2001. Aging Effects on Cadmium Transport in Undisturbed Contaminated Sandy Soil Columns. Journalof Environmental Quality 30:1040-1050.
Schaerlaekens, J., D. Mallants, J. S̆imůnek, M. Th. van Genuchten,
and J. Feyen, 1999. Numerical Simulation of Transport andSequential Biodegradation of Chlorinated Aliphatic Hydrocarbons Using CHAIN_2D. J. of Hydrological Processes 13(17):2847-2859.
Schaap, M.G., F.J. Leij, and M. Th. van Genuchten, 2001. ROSETTA: A Computer Program for Estimating Soil Hydraulic Parameters With Hierarchical Pedotransfer Functions, J. Hydrol.
251:163-176.
Schijven, J. and J. S̆imůnek, 2002. Kinetic Modeling of Virus Transport at Field Scale, J. of Contam. Hydrology 55(1-2):113-135.
S̆imůnek, J. and J.W. Hopmans, 2002. Parameter Optimization and
Nonlinear Fitting (Third Edition). In: Methods of Soil Analysis,
Part 1, Physical Methods, J.H. Dane and G. C. Topp (Editors).
SSSA, Madison, Wisconsin, Chapter 1.7, pp. 139-157.
S̆imůnek, J., D. Jacques, J.W. Hopmans, M. Inoue, M. Flury, and
M. Th. van Genuchten, 2002. Solute Transport During VariablySaturated Flow – Inverse Methods (Third Edition). In: Methodsof Soil Analysis, Part 1, Physical Methods, J.H. Dane and G.C.
Topp (Editors). SSSA, Madison, Wisconsin, Chapter 6.6, pp.
1435-1449.
S̆imůnek, J., D. Jacques, and M. Th. van Genuchten, 2004. Development of General Conceptual-Numerical Model to Simulate theFate and Subsurface Transport of Explosives, and the Moisture
and Temperature Signatures Around Land Mines, ARO Review
Meeting of Basic Research in Surface and Groundwater Hydrology, Geotechnical Engineering, and Geomorphology, USACEERDC Waterways Experiment Station Vicksburg, Mississippi.
S̆ imůnek, J., M. S̆ ejna, and M. Th. van Genuchten, 1998. The
HYDRUS-1D Software Package for Simulating the One-Dimensional Movement of Water, Heat, and Multiple Solutes in Variably-Saturated Media, Version 2.0. IGWMC-TPS-70, Int.
Ground Water Modeling Center, Colorado School of Mines, Golden, Colorado, 202 pp.
S̆ imůnek, J., M. S̆ ejna, and M. Th. van Genuchten, 1999. The
HYDRUS-2D Software Package for Simulating Two-Dimensional Movement of Water, Heat, and Multiple Solutes in Variably
Saturated Media, Version 2.0. IGWMC-TPS-53, International
Ground Water Modeling Center, Colorado School of Mines, Golden, Colorado, 251 pp.

JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATION

S̆imůnek, J., and D. L. Suarez, 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., and D.L. Suarez, 1997. Sodic Soil Reclamation Using
Multicomponent Transport Modeling. ASCE J. Irrig. and Drain.
Engineering 123(5):367-376.
S̆imůnek, J., D.L. Suarez, and M. S̆ejna, 1996. The UNSATCHEM
Software Package for Simulating One-Dimensional Variably
Saturated Water Flow, Heat Transport, Carbon Dioxide Production and Transport, and Multicomponent Solute Transport With
Major Ion Equilibrium and Kinetic Chemistry, Version 2.0.
Research Report No. 141, U.S. Salinity Laboratory, USDA, ARS,
Riverside, California, 186 pp.
S̆imůnek, J. and A.J. Valocchi, 2002. Geochemical Transport (Third
Edition). In: Methods of Soil Analysis, Part 1, Physical Methods,J. H. Dane and G.C. Topp (Editors). SSSA, Madison, Wisconsin,
Chapter 6.9, pp. 1511-1536.
S̆ imůnek, J., M. Th. van Genuchten, and M. S̆ 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 Software Series 1.
Department of Environmental Sciences, University of California
Riverside, Riverside, California, 270 pp.
Steefel, C.I., 2000. New Directions in Hydrogeochemical Transport
Modeling: Incorporating Multiple Kinetic and Equilibrium Reaction Pathways. In: Computational Methods in Water Resources
XIII, L.R. Bentley, J.F. Sykes, C.A. Brebbia, W.G. Gray, and G.F.
Pinder (Editors). A.A. Balkema, Rotterdam, pp. 331-338.
Truesdell, A.H. and B.F. Jones, 1974. WATEQ: A Computer Program for Calculating Chemical Equilibria of Natural Waters.
J. Res.. U. S. Geol. Surv. 2(2): 233-248.
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. and J. S̆imůnek, 2004, Integrated Modelingof Vadose Zone Flow and Transport Processes. In: Unsaturated
Zone Modelling: Progress, Challenges and Applications, R.A.
Feddes, G.H. de Rooij, and J. C. van Dam (Editors). Wageningen,
The Netherlands, pp. 37-69.
Wagenet, R.J., and J.L. Hutson, 1987. LEACHM: Leaching Estimation And CHemistry Model: A Process-Based Model of Waterand Solute Movement, Transformations, Plant Uptake andChemical Reactions in the Unsaturated Zone. Continuum 2.
Water Resour. Inst., Cornell University, Ithaca, New York.
Yeh, G.-T. and H.-P. Cheng, 1999. 3DHYDROGEOCHEM: A
3-Dimensional Model of Density-Dependent Subsurface Flowand Thermal Multispecies-Multicomponent HYDROGEOCHEMical Transport. EPA/600/R-98/159, 150 pp.
Yeh G.T., S. Sharp-Hansen, B. Lester, R. Strobl and J. Scarbrough,1992. 3DFEMWATER/ 3DLEWASTE: Numerical Codes forDelineating Wellhead Protection Areas in Agricultural Regions
Based on the Assimilative Capacity Criterion. Office of Researchand Development, Environmental Protection Agency, Athens,
Georgia, EPA/600/R-92/223, 278 pp.
Zheng, Chunmiao and P. Patrick Wang, 1999. MT3DMS, A Modular
Three-Dimensional Multi-Species Transport Model for Simulation of Advection, Dispersion and Chemical Reactions of Contaminants in Groundwater Systems; Documentation and User’s
Guide, U.S. Army Engineer Research and Development Center
Contract Report SERDP-99-1, Vicksburg, Mississippi, 202 pp.

1547

JAWRA