OPEN REPORTSCK•CEN-BLG-998
05/DJa/P-36

User Manual of the Multicomponent
Variably- Saturated Flow andTransport Model HP1

Description, Verification, and Examples
Version 1.0

Diederik Jacques* and Jirka Šimůnek#*Waste and Disposal
SCK•CEN, Mol, Belgium#
Department of Environmental Sciences
University of California, Riverside, USA

June, 2005

SCK•CEN
Boeretang 200
2400 Mol
Belgium

Waste & Disposal Department

OPEN REPORT OF THE BELGIAN NUCLEAR RESEARCH CENTRESCK•CEN-BLG-998

User Manual of the Multicomponent
Variably- Saturated Flow andTransport Model HP1

Description, Verification, and Examples
Version 1.0

Diederik Jacques* and Jirka Šimůnek#*Waste and Disposal
SCK•CEN, Mol, Belgium#
Department of Environmental Sciences
University of California, Riverside, USA

June, 2005
Status: Unclassified
ISSN 1379-2407

SCK•CEN
Boeretang 200
2400 Mol
Belgium

Waste & Disposal Department

© SCK•CEN
Belgian Nuclear Research Centre
Boeretang 200
2400 Mol
Belgium
Phone +32 14 33 21 11
Fax +32 14 31 50 21
http://www.sckcen.beContact:
Knowledge Centre
library@sckcen.be

RESTRICTED
All property rights and copyright are reserved. Any communication or reproduction of this document, andany communication or use of its content without explicit authorization is prohibited. Any infringement tothis rule is illegal and entitles to claim damages from the infringer, without prejudice to any other right incase of granting a patent or registration in the field of intellectual property.
SCK•CEN, Studiecentrum voor Kernenergie/Centre d'Etude de l'Energie Nucléaire
Stichting van Openbaar Nut – Fondation d'Utilité Publique - Foundation of Public Utility
Registered Office: Avenue Herrmann Debroux 40 – B-1160 Brussel
Operational Office Boeretang 200, 2400 Mol, Belgium

i

TABLE OF CONTENTS
LIST OF TABLES.................................................................................................................... v
LIST OF FIGURES................................................................................................................vii
1

Introduction ....................................................................................................................... 1
1.1

Background.................................................................................................................. 1

1.2 Features and limitations of the coupled HP1 model.................................................. 3
1.2.1
Features ............................................................................................................. 3
1.2.2
Limitations ........................................................................................................ 3
2

Description of the model ................................................................................................... 5
2.1 Water flow in the vadose zone..................................................................................... 5
2.1.1
The water flow equation .................................................................................. 5
2.1.2
Initial and boundary conditions for water flow............................................. 6
2.2 Solute transport in the vadose zone ............................................................................ 7
2.2.1
The solute transport equation ......................................................................... 7
2.2.2
Initial and boundary conditions for the solute transport equation ............. 9
2.3 Heat transport in the vadose zone............................................................................... 9
2.3.1
The heat transport equation ............................................................................ 9
2.3.2
Initial and boundary conditions for heat transport .................................... 10
2.4 Geochemical reactions .............................................................................................. 10
2.4.1
Homogeneous aqueous reactions................................................................... 11
2.4.2
Heterogeneous ion exchange reactions ......................................................... 11
2.4.3
Heterogeneous mineral dissolution/precipitation ........................................ 12
2.4.4
Kinetic reactions ............................................................................................. 13
2.5 Multicomponent reactive transport........................................................................... 13
2.6

Coupling procedures ................................................................................................. 15

2.7 Model structure.......................................................................................................... 18
2.7.1
PHREEQC ...................................................................................................... 18
2.7.2
HYDRUS-1D ................................................................................................... 19
3

4

Description of data input ................................................................................................ 21
3.1

Input data................................................................................................................... 21

3.2

Running the model and output ................................................................................. 23

Verification problems...................................................................................................... 25
4.1 Modelling the transport of a single component or decay chain .............................. 25
4.1.1
Physical equilibrium and nonequilibrium transport of chloride ............... 25
4.1.2
Transport of a nonlinearly adsorbing solute subject to first-order decay 28
4.1.3
Transport of first-order decay chain of nonlinearly adsorbing
contaminants during transient flow ............................................................................ 35

ii

4.2 Multicomponent transport during steady-state flow ................................................ 39
4.2.1
Cation exchange reactions ............................................................................. 39
4.2.2
Mineral dissolution......................................................................................... 42
4.3 More complicated verification problems of HP1 model .......................................... 44
4.3.1
Heavy metal transport in a medium with a pH-dependent cation exchange
complex (Verification Problem 8 – MCATEXCH) .................................................... 44
4.3.2
Infiltration of a hyperalkaline solution in a clay sample (Verification
Problem 9 – ALKALINE) ............................................................................................ 50
5

Example ............................................................................................................................ 59
5.1 Long term transient flow and transport of major cations and heavy metals in a soil
profile (HEAVMET) ........................................................................................................... 59
5.1.1
Problem definition and governing chemical reactions................................ 59
5.1.2
Analysis of the simulation results.................................................................. 61

6

Concluding Remarks....................................................................................................... 67

References................................................................................................................................ 71

Acknowledgements ................................................................................................................. 69
References................................................................................................................................ 71
APPENDIX A – LISTING OF NOTATION........................................................................ 75
APPENDIX B – RELEASED VERSIONS AND BENCHMARKING.............................. 79

iii

ABSTRACT
Jacques, D., and J. Šimůnek, 2004. User manual of the Multicomponent variably-saturated
transport model HP1 (Version 1.0): Description, Verification and Examples. SCK•CEN, Mol,
Belgium, BLG-998, 79 p.
This report describes a new comprehensive simulation tool HP1 (HYDRUS1D-PHREEQC)
that was obtained by coupling the HYDRUS-1D one-dimensional variably-saturated waterflow and solute transport model with the PHREEQC geochemical code. The HP1 code
incorporates modules simulating (1) transient water flow in variably-saturated media, (2)transport of multiple components, and (3) mixed equilibrium/kinetic geochemical reactions.
The program numerically solves the Richards equation for variably-saturated water flow andadvection-dispersion type equations for heat and solute transport. The flow equation incorporates
a sink term to account for water uptake by plant roots. The heat transport equation considers
transport due to conduction and convection with flowing water. The solute transport equations
consider advective-dispersive transport in the liquid phase. The program can simulate a broad
range of low-temperature biogeochemical reactions in water, soil and ground water systems
including interactions with minerals, gases, exchangers, and sorption surfaces, based onthermodynamic equilibrium, kinetics, or mixed equilibrium-kinetic reactions. The program
may be used to analyze water and solute movement in unsaturated, partially saturated, or fully
saturated porous media. The flow region may be composed of nonuniform soils or sediments.
Flow and transport can occur in the vertical, horizontal, or a generally inclined direction. The
water flow part of the model can deal with prescribed head and flux boundaries, boundaries
controlled by atmospheric conditions, as well as free drainage boundary conditions. The
governing flow and transport equations were solved numerically using Galerkin-type linear finite
element schemes.
To test the accuracy of the coupling procedures implemented in HP1, simulation results werecompared with (i) HYDRUS-1D for transport problems of multiple components subject tosequential first-order decay, (ii) PHREEQC for steady-state flow conditions, and (iii)
calculations obtained from an independent geochemical transport model (CRUNCH) for
several relatively complex problems. Nine verification examples of increasing complexity aredescribed in this report.
This report serves as both a user manual and reference document. Detailed instructions for input
data preparation and interpretation of output data are given in the manuals of the original
HYDRUS-1D and PHREEQC codes. The graphical user interfaces of both HYDRUS-1D and
PHREEQC can be used for easy input data preparation and output display in the MS Windows
environment.

KEYWORDS
Biogeochemical model, variably-saturated water flow, multicomponent solute transport,
vadose zone, dissolution/precipitation, cation exchange, aqueous complexation, benchmark
calculations

iv

The software has been verified against selected test cases. However, no warranty is given that the
program is completely error-free. If you do encounter problems with the code, find errors, or
have suggestions for improvement, please contact one of the authors at
Diederik Jacques
Tel: +32-14-333209
Fax: +32-14-323553
E-mail: djacques@sckcen.be
Jirka Šimůnek
Tel: 1-951-827-7854
Fax: 1-951-827-3993
E-mail: jiri.simunek@ucr.edu

v

LIST OF TABLES
Table 2.1 Description of variables that are used in both HYDRUS-1D and PHREEQC,
or newly defined in HYDRUS-1D...........................................................................20
Table 4.1 Transport and simulation parameters for Verification Problem 3.........................30
Table 4.2 Soil hydraulic, transport, and reaction parameters for Verification Problem
5. .............................................................................................................................36
Table 4.3 Input simulation parameters and numerical results for Verification Problem
5. .............................................................................................................................37
Table 4.4 Initial and inflow concentrations for Verification Problem 6. ...............................39
Table 4.5 Soil hydraulic properties and cation exchange capacities of five soil layers.........45
Table 4.6 pH and solution concentrations used in the simulation (µmol l-1). ........................45
Table 4.7 Overview of aqueous equilibrium reactions and corresponding equilibrium
constants (data from phreeqc.dat database, Parkhurst and Appelo, 1999). ..........45
Table 4.8 Log K parameters for multi-site exchange complex. ..............................................46
Table 4.9 Mineralogical reactions, Log(K), molar volume, volume percent and molespresent in Opalinus Clay, and secondary minerals considered. ............................51
Table 4.10 Parameters for the kinetic dissolution reactions (Eq. (4.7)). ...............................51
Table 4.11 Overview of aqueous equilibrium reactions, and cation exchange half
reactions with corresponding equilibrium constant.............................................52
Table 4.12 Chemical composition of the initial pore water and the inflowing alkaline
solution. ................................................................................................................53
Table 5.1 Soil hydraulic and other properties of six soil horizons.........................................59
Table 5.2 Chemical composition of the rain water (from Stolk, 2001) and in the initialsoil solution. ...........................................................................................................60

vi

vii

LIST OF FIGURES
Figure 2.1 Schematic of the modelling approach of the coupled HP1 model. .......................17
Figure 4.1 Outflow concentrations for Verification Problem 1. Full and dashed lines
are results for physical equilibrium and nonequilibrium transport,
respectively, obtained with HYDRUS-1D. Dots are results obtained with
HP1. ......................................................................................................................26
Figure 4.2 Outflow concentrations for Verification Problem 2. The full line was
generated with HYDRUS-1D, while dots were obtained with HP1. ....................27
Figure 4.3 Depth profiles of Cont after 250, 500, and 750 d for different simulations
defined in Table 4.1 for Verification Problem 3. Tests consider effects ofgrid size for PHREEQC (a), and effects of maximum time steps for HP1 (b)......31
Figure 4.4 Depth profiles of Cont after 250, 500, and 750 d for different simulations
defined in Table 4.1 (Verification Problem 4 with a first-order decaycoefficient of 0.2 d-1). ............................................................................................34
Figure 4.5 Distribution of Conta, Contb, and Contc versus depth after 250, 500, and1000 d for different simulations as defined in Table 4.3 (Verification
Problem 5). ...........................................................................................................38
Figure 4.6 Distribution of pH, dissolved Ca, Cl, Zn, and Cd concentrations,exchangeable Ca, Na, Zn, and Cd concentrations versus depth after threedays of infiltration for Verification Problem 6. ....................................................41
Figure 4.7 pH, Ca, Cl, Zn, and Cd concentrations in the effluent for Verification
Problem 6 (Full line: HP1, dashed line: PHREEQC)...........................................42
Figure 4.8 Distribution of pH, Si, Al, amorphous SiO2, and Gibbsite versus depth after
150 days of infiltration for Verification Problem 7 (Full line: HP1, dashedline: PHREEQC). .................................................................................................43
Figure 4.9 Selected results for cation and heavy metal transport with a pH-dependent
cation exchange complex. (a) and (b) pH and Cd concentrations in outflow
at 50 cm depth, respectively; (c) and (d) distributions of pH and Cd
concentrations versus depth after 0.3, 0.5, and 0.7 y, respectively; (e) and(f) distribution of the fraction of deprotonated sites (1 - H-sites) and sites
with Cd versus depth after 0.3, 0.5, and 0.7 y, respectively. ................................48
Figure 4.10 Effect of pH and Cl concentration on Cd concentrations of leaching waterat the 50-cm depth. .............................................................................................49
Figure 4.11 Distribution of pH and selected elements versus depth (mol/l) after 0.3,
0.7, and 1.1 y of infiltration of a high pH solution in an Opalinus Claycore. ....................................................................................................................56
Figure 4.12 Outflow curves of selected elements (mol/l) during infiltration of a highpH-solution in an Opalinus Clay core................................................................57
Figure 5.1 Cumulative precipitation, evaporation and surface flux into the soil between
1972 and 1982. .....................................................................................................60
Figure 5.2 Time series of inflow minus ouftlow (dashed line) and change in the total
amount present in the soil profile (dots) for Cl, Na, Ca, and Cd. ........................62

viii

Figure 5.3 Relative mass balance errors, εr, as a function of time for Na, Ca, and Cd. ........63
Figure 5.4 Space-time plots of (a) water content, (b) amount of Cl (mmol / 1000 cm³soil), (c) pH, and (d) amount of aqueous Cd (mmol / 1000 cm³ soil) in the Ahorizon (0 – 7 cm depth)........................................................................................64
Figure 5.5 Time series of water content, Cl concentration (mmol / 1000 cm³ soil), pH,
and Cd concentration (mmol / 1000 cm³ soil) in the third horizon (Bh1horizon) at 22 cm depth. .......................................................................................65.

1

1 Introduction

1.1 Background
The migration of many naturally occurring elements and contaminants in the subsurface isaffected by a multitude of complex, interactive physical, chemical, mineralogical, geological,
and biological processes. Cycles of precipitation and evaporation largely determine if
contaminants remain near the soil surface. Changes in the chemical composition or pH of thesoil solution may impact the retention of heavy metals on organic matter or iron oxides.
Dissolution and precipitation of minerals generally buffer the transport of a solution with adifferent pH through the soil profile. Simulation of these and related processes requires acoupled reactive transport code that integrates the physical processes of water flow andadvective-dispersive solute transport with a range of biogeochemical processes.
In this report, we present a new code that resulted from the coupling of two existing codes:
HYDRUS-1D (Šimůnek et al., 1998) and PHREEQC-2 (Parkhurst and Appelo, 1999). The
new code should significantly expand applicability of the individual codes by preserving mostof their original features and capabilities.
HYDRUS-1D simulates water, heat, and solute movement in one-dimensional variably
saturated porous media for various boundary conditions, including precipitation andevaporation. A sink term accounting for water uptake by roots is also included in the model.
Solutes can be exchanged between the water and gas phase and may interact linearly or
nonlineary with the solid phase assuming either equilibrium or nonequilibrium reactions
between the dissolved and adsorbed solutes. The only possible interaction between the
different solutes is a consecutive chain reaction in which the solutes are sequentiallytransformed along the chain by means of first-order reactions (which hence depend only on
the concentration of the first solute) (van Genuchten, 1985). These chain reactions can beused to describe the degradation of pesticides (e.g., Wagenet and Hutson, 1987) or chlorinatedvolatile organic compounds (e.g., Schaerlaekens et al., 1999; Casey and Šimůnek, 2001), andconsecutive decay chains involving radionuclides (e.g., Mallants et al., 2003), endocrine
disrupting chemicals (Casey et al., 2003, 2004), and many other chemicals. No otherinteractions between different species or components are currently considered in the
HYDRUS-1D model.
The only attempts to include geochemistry into HYDRUS-1D-related models were those byŠimůnek and Suarez (1994, 1996) and Suarez and Šimůnek (1997), resulting in
UNSATCHEM-2D (which was based on the SWMS-2D code (Šimůnek et al., 1992), a twodimensional precursor of HYDRUS-2D) and UNSATCHEM. These two models consideredinteractions only between major ions whose possible components (Ca, Mg, Na, K, SO4, Cl,2alkalinity, and CO2) and geochemical reactions (speciation, cation exchange, andprecipitation/dissolution of minerals (and their kinetics)) are predefined and hence can not bechanged by the user. Although these two codes can be applied to a wide range of important
problems such as salinization of agricultural soils under irrigation or reclamation of sodic
soils (Šimůnek and Suarez, 1997), they can not be used to simulate the transport and reactionsof other chemical species, such as trace elements, radionuclides, and other chemicals.
PHREEQC-2 can simulate a large number of low-temperature geochemical reactions in water,
soil and ground water systems, including interactions with minerals, gases, solid solutions,exchangers, and sorption surfaces, based on thermodynamic equilibrium, kinetics, or mixed
equilibrium-kinetic reactions (e.g., Appelo et al., 2002). PHREEQC-2 also allows one tosimulate one-dimensional reactive transport using a mixing cell solution approach (seeAppelo and Postma, 1999, for details; e.g., Appelo et al., 1998; Postma and Appelo, 2000).
This model hence can be used to simulate reactive transport during steady-state flow,
including a wide variety of geochemical reactions. However, the model can not deal withsolute transport during transient water flow conditions.
In recent years, the PHREEQC geochemical code has been coupled to various (groundwater)water flow and solute transport models. For example, the PHT-3D model (Prommer, 2002)
couples PHREEQC-2 with MT3DMS (Zheng and Wang, 1998), the latter being an extensionof the single-species transport simulator by Zheng (1990). However, this coupled program
deals only with solute transport and reactions, while the groundwater flow field needs to becomputed using a separated simulation model. Xu (1996) developed a coupled model TRANPHREEQE that links PHREEQE with a two-dimensional finite element code for flow andtransport in aquifers. PHREEQE is a precursor of PHREEQC without cation exchange,
surface complexation, or kinetic reactions. Another example of coupling between PHREEQC
and a two-dimensional transport model for water-saturated conditions and constanttemperature was presented recently by Källvenius and Ekberg (2003) in the TACK model(Transport And Chemical Kinetics). However, this model does not simulate water flow. A
recent detailed review of available numerical multicomponent transport models is given byŠimůnek and Valocchi (2002), including an overview of the mathematical equations
representing the major chemical reactions and governing transport processes, a brief discussion
how these equations can be implemented in reaction multispecies transport models, and a
description of several possible applications.
This report describes a new comprehensive simulation tool HP1 (HYDRUS1D-PHREEQC)
that results from coupling the HYDRUS-1D one-dimensional variably-saturated water flow
and solute transport model with the PHREEQC geochemical code. The model incorporates
modules simulating (1) transient water flow in variably-saturated media, (2) transport ofmultiple components, and (3) mixed equilibrium/kinetic geochemical reactions. The accuracyof the coupled HP1 model was evaluated by comparing simulation results with HYDRUS-1D,

3

PHREEQC for steady-state flow conditions, and with calculations obtained with anindependent geochemical transport model CRUNCH (Steefel, 2002) for several more
complicated problems.

1.2 Features and limitations of the coupled HP1 model
1.2.1

FeaturesAny combination of the following features can be described with the HP1 model:• One-dimensional transient water flow for different boundary conditions including
atmospheric conditions (precipitation, evaporation, transpiration)•

Root water uptake as a sink for water•

Root growth•One-dimensional transient convective and conductive heat transport for time-variabletemperatures at the soil surface••

One-dimensional advective, dispersive and diffusive transport of multiple solutes
Effect of temperature on transport parameters, thermodynamic constants, and rate
parameters•

Different functional forms for the soil hydraulic properties, including hysteresis•

Physical non-equilibrium solute transport•

Physical and chemical spatial heterogeneity•

Equilibrium aqueous speciation reactions and kinetically controlled aqueous reactions
such as radioactive decay••

Multi-site cation exchange related to type and amount of minerals present
Equilibrium and kinetic dissolution/precipitation of primary and secondary minerals•

User-defined kinetic reactions•

Simultaneous presence of different reactions (sequential and parallel kinetic reactions,equilibrium and kinetic reactions, homogeneous and heterogeneous reactions,biogeochemical reactions)

1.2.2

Limitations

Specific limitations of the PHREEQC model for various geochemical calculations are
discussed in detail by Parkhurst and Appelo (1999, p. 4-6), and are not further discussed here.
Of these limitations, those related to flow and transport modelling are no longer of concern
here, since HYDRUS-1D is used to simulate transient water flow and solute transport.

4

One possible limitation involves the invoked method of coupling HYDRUS-1D and
PHREEQC, i.e. a non-iterative sequential coupling method (SNIA). This method is still
being discussed in the literature (e.g., Steefel and MacQuarrie, 1996, Mayer, 1999), withsome authors claiming that mass balance errors may occur when this coupling procedure isused. However, by using appropriate time steps, accurate results can be obtained as we willshow with examples later in this manual (see also Mayer, 1999, for some guidelines). In
addition, we believe that uncertainty in the assumed processes and its parameters likely willcontribute much more to uncertainty in the model simulations than possible (limited)
numerical errors caused by the coupling procedure.
This manual documents Version 1.0 of the coupled HP1 model. The following features of
PHREEQC are not yet operational in the current version: surface complexation, solid
solutions, and redox reactions. Diffusion and advection of components in the gas phase arealso not considered. We further do not account for changes in the volume of minerals andcorresponding changes in porosity, hydraulic properties, and solute transport parameters.

5

2 Description of the model

The HP1 model is the result of coupling the variably-saturated water flow and solute transportmodel HYDRUS-1D (Šimůnek et al., 1998) with the geochemical model PHREEQC(Parkhurst and Appelo, 1999). Details about the governing equations, initial and boundary
conditions, parameterization, and adopted numerical methods are given in manuals of theoriginal HYDRUS-1D and PHREEQC models. In this chapter we give only a very concise
overview of these topics, mainly to provide a better understanding of the verification
problems described in chapter 4. For more detail, users are encouraged to examine theoriginal manuals.

2.1 Water flow in the vadose zone

2.1.1 The water flow equation
Combination of the mass balance equation with the Darcy-Buckingham law results in the
Richards equation that describes water flow in variably-saturated porous media. The onedimensional form of the Richards equation can be written as)(∂θ (h)= ∂ ⎡ K (h) ∂h + cos α ⎤ − S (h)
∂t
∂x ⎣⎢
∂x
⎦⎥

(2.1)

where h is the water pressure head [L], θ is the volumetric water content [L3L-3], t time [T], x
is the spatial coordinate [L] (positive upward), S is the sink term [L3L-3T-1], α is the angle
between the flow direction and the vertical axis, and K is the unsaturated hydraulicconductivity [LT-1].
Both the water content and the unsaturated hydraulic conductivity are nonlinear functions of
the pressure head. Three analytical models are available in HYDRUS-1D to describe thesesoil hydraulic properties (Brooks and Corey, 1994; van Genuchten, 1980; and Vogel andCislerova, 1988). Only van Genuchten’s functions will be used in the verification problems
documented in this manualθ ( h ) = θr +θs − θr(2.2)(1 + α h )n m(

K ( h ) = K s S el ⎡1 − 1 − S e1 / m
⎢⎣

) ⎤⎥⎦
m

2

(2.3)

6

where θr is the residual water content [L3L-3], θs is the saturated water content [L3L-3], α [L-1],
n [-] and m (= 1 – 1/n) [-] are shape parameters, l is the pore connectivity parameter [-], Ks is
the saturated hydraulic conductivity [LT-1], and Se=(θ-θr)/(θs-θr) is effective saturation. These
parameters can be obtained by direct measurement of the θ-h or K-h relationships (e.g.,chapters 3.3 and 3.4 in Dane and Topp, 2002), by inverse optimization (Hopmans et al., 2002)
or indirect estimation from basic soil textural properties using pedotransfer functions. The
latter approach in HYDRUS-1D uses neural network predictions from textural data as
implemented in the ROSETTA program of Schaap et al. (1998). The following additionalfeatures related to the soil hydraulic properties are also included in HYDRUS-1D: (i)hysteresis in the water retention characteristic and the hydraulic conductivity function, (ii)description of small-scale spatial variability in the soil hydraulic properties by means ofscaling factors, and (iii) temperature dependence of the soil hydraulic functions. These
additional features are described in sections 2.6, 2.4, and 2.5 of Šimůnek et al. (1998),
respectively.
The sink term S in (2.1) is defined as the volume of water extracted from the soil by the roots.
The potential root water uptake rate Sp(x) is obtained by multiplying a normalized water
uptake distribution b(x) [L-1] with the potential transpiration Tp [LT-1]. The b(x) function is
obtained from the root distribution with depth, whereas Tp depends on climate conditions andvegetation. Both b(x) and Tp are input to the HYDRUS-1D model. The actual root water
uptake rate S(x) is obtained by multiplying Sp(x) with a root water stress response function(e.g., Feddes et al., 1978; van Genuchten, 1987) to account for a possible reduction in root
water uptake due to water and salinity stress conditions in the soil profile. Soil water uptake
reduction due to the salinity stress can be included using an osmotic head reduction functionthat can be either additive or multiplicative to water stress. The actual water uptake
distribution is then of the formS ( h, hφ , x ) = α ( h, hφ )b ( x )T p

(2.4)

where S is the root water uptake as a function of the pressure head h (related to water stress),the osmotic head hφ [L] (related to salinity stress), and depth x (related to the root spatial
distribution), while α(h,hφ) defines the reduction in root water uptake due to the water andsalinity stress. Root growth can be included using the Verhulst-Pearl logistic growth function,with the assumption of an exponential root distribution with depth.
2.1.2 Initial and boundary conditions for water flow
To solve Eq. (2.1), initial and boundary conditions must be specified. Initial conditions can be
defined in terms of pressure heads or water contents. Possible system-independent boundary
conditions are time series of pressure heads or soil water fluxes at the soil surface and/or the

7bottom of the soil profile, and a zero-gradient (free drainage) pressure head boundary
condition at the bottom of the soil profile. In addition, system-dependent boundary conditions
that depend on the status of the system are also available. When atmospheric conditions(precipitation, evaporation, and transpiration) defining the potential water flux across the top
boundary are specified, the actual water flux across this boundary depends also on the soil
moisture conditions. When the potential surface flux (precipitation) is larger than theinfiltration capacity, any excess water on the soil surface is either assumed to be immediately
removed by surface runoff, or is permitted to build up on the soil surface. At the bottom of thesoil profile, the following boundary conditions can be implemented: (i) a seepage faceboundary condition that assumes a zero-flux when the bottom of the soil profile is unsaturated
and a zero pressure head when it is saturated, (ii) a tile drain boundary condition thatapproximates flow to horizontal subsurface drains using selected analytical solutions, and (iii)
a deep drainage boundary condition that uses a functional relationship that relates the watertable depth with the deep recharge from the soil profile. Mathematical descriptions of theseboundary conditions can be found in section 2.7 of Šimůnek et al. (1998).

2.2 Solute transport in the vadose zone

2.2.1

The solute transport equation

The HYDRUS-1D code allows simulations of the transport of multiple solutes involved in asequential first-order decay chain in three phases (liquid, solid, and gaseous) using thephysical nonequilibrium advection-dispersion equations. However, many solute transportfeatures of HYDRUS-1D are not used in the coupled HP1 code since they are considered in
the PHREEQC module. These include interactions between the liquid and solid phases,
degradation/production, and the presence of sequential decay chain reactions. These
interactions are defined in PHREEQC using equilibrium or kinetic reactions. In fact, solutetransport in the HYDRUS module is modelled as the transport of inert tracers (i.e., nointeraction with the solid phase, and no solute sink terms) since reactions are considered in the
PHREEQC module. Note, however, that it is still possible to simulate reactive transport with
HYDRUS-1D in the coupled model when particular parameters are entered with non zero
values. In the present coupled model, interactions with (and thus also diffusion in) the gas
phase are not considered.
The general solute transport equations, as given by (3.1) and (3.2) in Šimůnek et al. (1998),
for conditions described above reduce to∂θ ci∂c ⎞ ∂qci⎛= ∂ ⎜ θ Diw i ⎟ −− Scr , i + Ri∂t∂x ⎝∂x ⎠ ∂x(2.5)

8

where i (= 1,…,Nm) is the aqueous species number (Nm is the total number of aqueous
species), ci is the aqueous concentration of the ith species [ML-3], q is the volumetric flux
density [LT-1], S the sink term in the water flow equation (Eq. (2.1)), cr the concentration of
the sink term [ML-3], Dw is the dispersion coefficient in the liquid phase [L2T-1], and Ri is thegeneral source/sink term due to geochemical reactions [ML-3T-1]. This sink/source term
contains heterogeneous equilibrium reactions, and homogeneous and heterogeneous kinetic
reactions (see section 2.4).
Physical nonequilibrium solute transport is modeled using a two region model that assumes
that the liquid phase can be divided into a mobile (flowing) region, θm [L3L-3], and an
immobile (stagnant) region, θim [L3L-3]. Solutes are exchanged between the mobile andimmobile regions by means of a first-order exchange process. The mathematical formulationof this nonequilibrium model is given by (van Genuchten and Wierenga, 1976):∂θ m ci , m ∂θ im ci , im∂c ⎞ ∂qci , m⎛+= ∂ ⎜ θ m Diw i , m ⎟ −− Scr , i + Ri , m
∂t
∂t
∂x ⎝
∂x ⎠
∂x
∂θ im ci ,im= ωi ( ci , m − ci ,im ) + Ri , im
∂t

(2.6)

where ci,m and ci,im are concentrations of the ith aqueous species in the mobile and immobile
regions [ML-3], respectively, Ri,m and Ri,im are the source/sink terms due to geochemical
reactions in the mobile and immobile regions [ML-3T-1], respectively, and ωi is the masstransfer coefficient for the ith aqueous species [T-1]. Note that HYDRUS-1D can also consider
chemical nonequilibrium transport (kinetic adsorption/desorption reactions). However, we
strongly suggest not to use this option of HYDRUS-1D, and to simulate chemical
nonequilibrium reactions using options in the PHREEQC module.
The dispersion coefficient is given byθ Diw = DL q + θ Di ,w τ w(2.7)

where Di,w is the molecular diffusion of the ith aqueous species in free water [L2T-1], DL is thelongitudinal dispersivity [L], and τw is a tortuosity factor in the liquid phase [-] that is related
to the water content as follows (Millington and Quirk, 1961)

τw =

θ 7/3
θ s2

(2.8)

The dependence of the diffusion parameter Diw on temperature can be described using the theArrhenius equation (see section 3.4 of Šimůnek et al., 1998).

9

2.2.2

Initial and boundary conditions for the solute transport equation

To solve Eqs. (2.5) or (2.6), the initial and boundary conditions must be specified. Initial total
aqueous concentrations of all aqueous species as a function of depth at time zero in both themobile and immobile regions must be defined. Concentrations of adsorbed secondary or
precipitated species must also be specified at time zero when kinetic adsorption andprecipitation/dissolution reactions are considered.
Possible boundary conditions include first-type (or Dirichlet type) boundary conditions
defining a prescribed boundary concentration, and third-type (or Cauchy type) boundary
conditions defining a prescribed boundary solute flux. At an impermeable boundary (i.e.,
where q=0) or at a boundary where water flows out of the domain, the third-type boundary
condition reduces to a second-type (Neumann type) boundary condition.

2.3 Heat transport in the vadose zone

2.3.1

The heat transport equation

The one-dimensional heat transport equation (neglecting water vapour diffusion) is given by∂C p (θ )T∂qT= ∂ ⎡ λ (θ ) ∂T ⎤ − C w− C w ST
∂t
∂x ⎣
∂x ⎦
∂x

(2.9)

where λ(θ) is the apparent thermal conductivity of the soil [MLT-3K-1], and Cp(θ) and Cw arevolumetric heat capacities of the porous medium and the liquid phase, respectively, [ML-1T2 -1
K ]. The volumetric heat capacity of the porous medium is estimated based on its
constituents (de Vries, 1963) as follows

C p (θ ) = C nθ n + C oθ o + C wθ + C aθ v

(2.10)

where Cn, Co, and Ca are the volumetric heat capacities of the solid phase, the organic matter,and the gas phase, respectively, [ML-1T-2K-1], and θn, θo, and θv are the volumetric fractions of
the solid phase, the organic matter, and the gas phase, respectively [L3L-3]. The apparent
thermal conductivity is defined as (de Marsily, 1986)λ (θ ) = λ0 (θ ) + β t C w q(2.11)

10

where βt is the thermal dispersivity [L] and λ0(θ) is the thermal conductivity of the soil
defined as (Chung and Horton, 1987)

λ0 (θ ) = b1 + b2θ + b3θ 0.5(2.12)

where b1, b2, and b3 are empirical parameters [ML-1T-2K-1].

2.3.2

Initial and boundary conditions for heat transport

To solve Eq. (2.9), initial and boundary conditions must be specified. The initial temperature
as a function of space at time zero must be defined. Possible boundary conditions again
include first-type (or Dirichlet type) boundary conditions defining a prescribed boundary
temperature and third-type (or Cauchy type) boundary conditions defining a prescribed heat
flux through the boundary. At an impermeable boundary (q = 0) or at a boundary where waterflows out of the domain, a third-type boundary condition reduces to a second-type (Neumanntype) boundary condition.
At the soil surface, the temperature can be represented using a sine function (Kirkham andPowers, 1972)⎛⎞
T = Tavg + A sin ⎜ 2π t − 7π ⎟
p
12
⎝ t
⎠

(2.13)

where pt is the period of time necessary to complete one temperature cycle [T], Tavg is theaverage temperature at the soil surface [K], and A is the amplitude of the sine wave [K].

2.4 Geochemical reactions

In general, species present in a system are related to each other by a set of reaction equations.
It is possible to write the various reaction equations in terms of a limited set of independent
components. The latter group permits one to define the stoichiometry of the system, and are
independent of each other (Morel and Herring, 1993). The species are thus divided in two
groups: (i) the components, primary species or master species (in the terminology of
PHREEQC), and (ii) the secondary species. The number of master species, Nm, equals thetotal number of species minus the number of reactions if the reactions are written
stoichiometrically independent (none of the reactions is a linear combination of the others).
The number of secondary species, Ns, is then defined as the number of species minus thenumber of master species. Thus, each reaction can be written in the canonical form

11
Nm

∑ν A = A
ji

j =1

m
j

(2.14)iwhere i = 1, …, Ns, Ns is the number of secondary species, Amj and Ai are the chemical formulaof the master and secondary species, respectively, and νji are the stoichiometric coefficients in
the reaction.
In the remainder of this section we present the mass-actions relations for different reactions(aqueous speciation, ion exchange, mineral precipitation/dissolution), and the activities of theaqueous and exchange species.
2.4.1

Homogeneous aqueous reactions

For an aqueous complexation reaction, Eq. (2.14) is written as
Nm∑ν A = A
l
ji

j =1

m
j

l
i

(2.15)

where the superscript l indicates the liquid phase, i=1, …, Nsa, where Nsa is the number ofaqueous secondary species. For equilibrium conditions, the mass-action equation isNm(

K = γ c ∏ γ mj c mj
l
i

l l
i i

)

−ν lji(2.16)

j =1

where Kli is the equilibrium constant for reaction (2.15), γli is the activity coefficient of the ith
aqueous complex, and γmj is the activity coefficient of the jth master species in solution. The
activity coefficients are defined with the Davies equation or the extended Debye-Hückel
equation (Langmuir, 1997; Parkhurst and Appelo, 1999).

2.4.2

Heterogeneous ion exchange reactions

In PHREEQC, the ion exchange sites are defined by exchange primary (or master) species(Xm) and ion exchange reactions are defined as half reactions. For the Gaines and Thomas
convention (Gaines and Thomas, 1953), the half reaction is written as

Nm∑ν A + ν
j =1

e
jie

m
j

e
jeie

X mje = Aiee

(2.17)

12

where je=1, …, NX (NX is the number of master exchangers), ie=1,…, Nse (Nse is the number ofthe secondary exchange species), and the superscript e refers to exchange reactions. Let theactivity of an exchange species be defined as

aie = γ iee β ie , je

(2.18)

where a is the activity, γei e is the activity coefficient of the ieth exchange species, and β ie , je isthe equivalent fraction of the ieth exchange species on the jeth exchanger, defined asbie je nie je / T je where bie je is the number of equivalents of exchanger je occupied by the ieth

exchange species, nie je are the moles of the ieth exchange species on exchanger je, and T je isthe total number of exchange sites for the jeth exchanger (in equivalents). Then the massaction equation for equilibrium conditions can be written as(γ β )

K iee = γ iee β ie , je ∏ (γ mj c mj )Nm−ν eji e

e
je

j =1

−ν ej e ieje , je(2.19)

The activity coefficients for the exchange species are calculated with the WATEQ DebyeHückel equation (Treusdell and Jones, 1974), the Davies equation, or are equal to one. Notethat in PHREEQC it is also possible to express ion exchange reactions using the Gapon
convention (Gapon, 1933; see also Appelo and Postma, 1999).

2.4.3

Heterogeneous mineral dissolution/precipitation

For equilibrium precipitation/dissolution reactions of a mineral, Eq. (2.14) is written as
Nm∑ν A = A
j =1

p
ji

m
j

p

i

(2.20)

where i=1, …, Np (Np is the number of minerals), and Api is the formula of the mineral, while
the superscript p refers to pure phases (minerals). Note that in the database of PHREEQC,
dissolution/precipitation can be written in terms of any of the aqueous species. These
reactions can always be transformed to a canonical form of Eq. (2.20). For equilibrium
conditions, the mass-action equation isNm(

K ip = ∏ γ mj c mj)

−ν jipj =1

since the activity of a pure phase (mineral) is assumed to be 1.

(2.21)

13

2.4.4

Kinetic reactions

Both homogeneous and heterogeneous reactions can be treated as kinetic reactions.
Homogeneous reactions define the production or consumption of a master species from othermaster species in the aqueous phase
Nm∑ν A = 0
j =1

lk
ji

m
j

(2.22)

where i=1,…, Nlk (Nlk is the number of kinetically-controlled homogeneous reactions in the
aqueous phase), the superscript lk indicates the homogeneous kinetic reaction, and ν lkji are thestoichiometric coefficients involved in the ith homogeneous kinetic reaction. The rateequation itself can be of any form and be dependent upon the total concentrations of a given
master species, on concentrations or activities of the master species, on concentrations or
activities of secondary aqueous species, or on external factors such as temperature andbiomass. It is possible to define rather complex sequential and parallel kinetic pathways (e.g.,
Steefel, 2000). This flexibility in PHREEQC is possible due to an embedded BASIC
interpreter which permits one to define rate expressions in a general way in the input file (seethe section ‘Numerical method and rate expressions for chemical kinetics’ in Parkhurst andAppelo, 1999). Typical examples of homogeneous kinetic reactions are oxidation-reduction
reactions, radioactive decay, and degradation reactions. Another type of homogeneous kinetic
reaction is the production or consumption of a particular master species (theconsumption/production of the specific secondary species is also possible). Only mineraldissolution/precipitation can be considered as a heterogeneous kinetic reaction. The BASIC
interpreter in PHREEQC allows one to use a broad range of reaction rate expressions.

2.5 Multicomponent reactive transport

Multicomponent reactive transport system may be viewed as consisting of Nm transport
equations for the master species⎛∂θ c j∂c j ⎞ ∂qc j= ∂ ⎜ θ D wj−
− Scr , j + R j
∂t
∂x
∂x ⎟
∂x
⎝
⎠

(2.23)

14(j=1,…, Nm), and Nsa transport equations for the secondary aqueous species⎛∂θ ci∂c ⎞ ∂qci= ∂ ⎜ θ Diw i ⎟ −− Scr , i + Ri
∂t
∂x ⎝
∂x ⎠
∂x

(2.24)(i=1,…, Nsa). The terms Rj and Ri include both rapid and slow reactions involving the given
species. However, some reactions can be so fast that the rate of the reaction is controlled by
the rate of transport of the species to or from the site of the reaction rather than by the reaction
itself. For these reactions, equilibrium can be assumed. Consequently, the sink/source term of
the geochemical reactions can be divided as
,hom
, het
,hom
, het
R j = R eq
+ R eq+ R kin+ R kinjjjjRi = Rieq ,hom + Rieq , het + Rikin ,hom + Rikin , het(2.25)

where the superscripts eq, kin, hom, and het refer to local equilibrium reactions, kinetic
reactions, homogeneous reactions, and heterogeneous reactions, respectively. A system
involving both local equilibrium and kinetic reactions is in a state of local partial equilibrium(Lichter, 1996). A system of Np+Nsa transport equations can be reduced to the number ofprimary species. The Req,hom
can be expressed in terms of the reaction stoichiometry definedjby Eq. (2.15)

Req , homjN sa= − ∑ ν ji Rieq , hom(2.26)

i =1

Substituting Eq. (2.26) in Eq. (2.23) and replacing Ri by Eq. (2.24) allows the global transport
equation for the jth master species to be written as⎛∂θ c j N sa
∂c j ⎞ ∂qc j∂θ ci+ ∑ν ji= ∂ ⎜ θ D wj
−
+
∂t
∂t
∂x
∂x ⎟
∂x
i =1
⎝
⎠
N sa

N

N

sa
⎛
∂c ⎞ sa∂qci+ ∑ν ji ∂ ⎜ θ Diw i ⎟ − ∑ν ji− Scr ,i − ∑ν ji Scr , j + Ro , j∂x ⎝∂x ⎠ i = 1
∂x
i =1
i =1

(2.27)

where Ro,j is the term that includes all other heterogeneous equilibrium and kinetic reactions.
Defining Cj as the total concentration of a master species
N sa

C j = c j + ∑ν ji ci
i =1(2.28)

15

and Crj the total concentration in the sink term
N sa

C rj = crj + ∑ν ji cri(2.29)

i =1

and assuming that the diffusion coefficients Dwj and Dwi are all equal (=Dw; i.e., speciesindependent diffusion) and that concentrations of the sink term cr,j and cr,i are equal to eachother, then the transport equation for the master species becomes⎛∂θ C j∂C j ⎞ ∂qC j= ∂ ⎜θ Dw− SC r ,i + Ro , j
⎟−
∂t
∂x ⎜
∂x ⎟
∂x
⎝
⎠

(2.30)

Note that it is possible to introduce heterogeneous equilibrium reactions in Eq. (2.30) in asimilar way, and thus to define transport equations for the total concentration of the masterspecies using only kinetic reactions as source/sink terms (Lichtner, 1996; Mayer, 1999).
However, this is not done here since the solution method for solving the reactive transport
equations is based on a sequential non-iterative approach (Steefel and MacQuarrie, 1996).
This means that the transport equation (2.30) is solved without the reaction term Ro,i, whereas
the mass-action equations (2.16), (2.19), and (2.21) are solved sequentially (see section 2.6).
Note that when cr is not equal to zero, a component will be taken up by the roots, but not a
species. When cr is zero, no solute is removed from the soil solution due to the root uptake. Itis, however, still possible to define some active uptake mechanism in PHREEQC, i.e., onethat is independent of water uptake (contained in the Ro,i term of Eq. (2.30)).

2.6 Coupling procedures

Reactive transport systems as defined in the above sections involve many processes that areinterrelated and contain parameters that are dependent upon the state of the system. Without
attempting to be complete, the following interactions may occur in natural systems•

Effect of concentration and temperature on flow properties (by affecting the water
density and the viscosity, and the surface tension at the air-water interface),•

Effect of temperature on diffusion coefficients,•

Effect of temperature on the equilibrium constants and rate coefficients,•

Effect of water flow on solute and heat transport,•

Effect of mineralogical changes on water flow and solute transport parameters.

Not all of these interactions are included in the present model. The effect of concentration on
the flow properties and the effect of mineralogical changes on the water flow, solute and heattransport equations are neglected in the present version of HP1.

16

Yeh and Cheng (1999) discriminate between strong coupling (in which the governing water
flow, solute transport and heat transport equations, as well as the geochemical reactions are
solved simultaneously) and weak coupling (in which the governing equations are solved
sequentially). In the latter method, state variables obtained after solving a given equation areused to calculate properties and state variables in the next equation. Yeh and Cheng (1999)
used the following sequence for this purpose: first solve the water flow equation, then the heattransport equation, and finally the solute transport equation. The approach assumes that (i) thetemperature effect on chemical reactions is important so that updated temperature information
should be used for the geochemical equilibrium and kinetic calculations, and that (ii) the
effect of advection (water flow) in the solute and heat transport equations is larger than theeffect of concentration and temperature on the water flow equation. Therefore, the water flow
equation should be solved before the solute and heat transport equations. The weak coupling
method is also invoked in our modelling approach. The same solutions sequence was used in
the original HYDRUS-1D model (Šimůnek et al., 1998), and is followed in HP1 as well.
Different approaches also exist to solve the multicomponent reactive transport equation(2.30). This equation contains terms describing the physical transport of the component (thefirst three terms on the right hand-side) and a term describing the geochemical reactions (thefourth term on the right hand-side). The physical transport part and the geochemical reactions
can be solved either simultaneously (a global-implicit or one-step approach) or sequentially(an operator-splitting, two-step, or sequential approach). For a discussion of these two
approaches, the reader is referred to Steefel and MacQuarrie (1996) and Mayer (1999). In ourmodel we use the sequential approach. Following Walter et al. (1994), the solution space has
three degrees of freedom: spatial, temporal, and chemical. Physical transport is connected in
the spatial and temporal domains, and the geochemical reactions are only connected in thechemical domain. The physical part (coupled in space, uncoupled over the components) is
obtained by solving Eq. (2.30) without the reaction term:⎛∂θ C j∂C j ⎞ ∂qC j= ∂ ⎜θ D w−− SC r ,i
∂t
∂x
∂x ⎟
∂x
⎝
⎠

(2.31)and the chemical part (uncoupled in space, i.e., no transport, but coupled over thecomponents) by simultaneously solving the equilibrium and kinetic geochemical reactions.
An overview of the coupled multicomponent reactive transport calculations is shown in
Figure 2.1. The symbols used in this figure are:n
the nth time stepHvariables related to water flow (pressure heads, fluxes)
Tvariables concerning heat transport (temperature)Cvariables dealing with components and species in the system

17

G
pw
ph
ps
pgc

variables concerning the solid phase (mineralogical composition, exchange site,
surface site, reactive surfaces)
vector of parameters needed to solve the water flow equation (Eq. (2.1))
vector of parameters needed to solve the heat transport equation (Eq. (2.9))
vector of parameters needed to solve the solute transport equation (Eq. (2.31))
vector of parameters needed to solve the geochemical reactions

Initial conditions
with state variables
H0, T0, C0, G0

Start a new time step nwith state variables
Hn, Tn, Cn, Gn

Compute new water flow state variables
Hn+1

Calculate new heat transport properties
ph(Hn+1)

Compute new heat transport state variables
Tn+1 based on Hn+1

n=n+1
Transport step

Calculate new solute transport properties
ps(Tn+1,Hn+1)

Compute new chemical state variables(no chemical reactions)
Ctransport based on Hn+1

Compute new geochemical state variables(no transport)
Cn+1, Gn+1 based on Ctransport

End of simulation
No
Yes

Figure 2.1 Schematic of the modelling approach of the coupled HP1 model.

PHREEQC

Geochemical step

Calculate new geochemical properties
pgc(Tn+1)

HYDRUS1D

Calculate new water flow properties
pw(Tn)

18

2.7 Model structure

In this section we summarize the changes implemented in both codes. While the Fortranroutines of HYDRUS-1D are compiled into Hydrus.dll, the c functions of PHREEQC arecompiled into PHREEQC.dll.

2.7.1

PHREEQC

The original PHREEQC code contained the following files: Advect.c, Basic.c, Cl1.c, Integr.c,
Inverse.c, Kinetics.c, Main.c, Mainsubs.c, Model.c, P2clib.c, Parse.c, Phqalloc.c, Prep.c,
Print.c, Read.c, Readtr.c, Spread.c, Step.c, Structur.c, Tidy.c, Transp.c, and Utility.c. In the
coupled HP1 we modified mostly only the main.c routine, and added a new file Hydr_tr.c

main.cFrom the main routine (main_Phreeqc) we call first the “read_text_file” function that reads
from the Species.in file the names of components to be transported by HYDRUS-1D, andstores these in the vector of strings cRows (nRows). Then we call the“check_hydrus_species” function that checks whether these components are defined in thegeochemical database. Then instead of the original “Transport()” function, we call
“HYDRUS_MAIN” with three parameters, which refer to c functions called from the
HYDRUS Fortran code. At the end we call function “free_string_array”, which deallocatesmemory from cRows.
void __stdcall HYDRUS_MAIN(void * run_HYDRUS_reactions, void *Get_Concentrations, void * Initialization);
int read_text_file(const char *file_path_name, int *nRows, char ***cRows);
void free_string_array(int nRows, char **cRows);
void check_hydrus_species(int nRows, char **cRows);

Hydr_tr.c
This file contains four new functions: Initialization, Get_Concentrations,
Run_Hydrus_reactions, and Get_Names that are called from HYDRUS. They are described
below:
void Initialization(int node_number, int lMobil, float * Temperature, float * Theta, float
ThImob, int max_species_number, int species_number, float * Concetrations,float * Im_Concetrations, float TimeInit);
void Get_Concentrations(int node_number, int max_species_number, int species_number,float * Concetrations, float * Im_Concetrations, int lMobil, float ThImob);

19

void run_HYDRUS_reactions(int node_number, int max_species_number, int
species_number, float * Concetrations, float time, float *Temperature, float * Theta,
float * Im_Concetrations, float TimeStep, int iPrint, int Step_No, int lMobil, float
ThImob);char * Get_Names(int species_number);

2.7.2

HYDRUS-1D

The original Hydrus-1D code contained the following files (see section 12 of Šimůnek et al.,1998): Hydrus.for, Input.for, Material.for, Output.for, Sink.for, Solute.for, Temper.for,Time.for, and Watflow.for. We added one additional file (Exports.h) that interfaces Fortranwith c. Most of the changes were done in the main HYDRUS file, where the main program
was converted into a HYDRUS_MAIN subroutine. From this subroutine we call thefollowing c functions that are located in the Hydr_tr.c file:call Initialization(NumNP,lMobil,TempN,ThOld,ThImob,NSD,NS,Conc,!
Sorb,tInit)
call Get_Concentrations(NumNP,NSD,NS,Conc,Sorb,lMobil,ThImob)
call run_HYDRUS_reactions(NumNP,NSD,NS,Conc,t,TempN,ThNew,Sorb,!
dt,iPrint,TLevel,lMobil,ThImob)and one new Fortran subroutine
call PhreeqcMB(NSD,NS,NumNP,Conc,Conc1,Sorb,Sorb1,x,ThNew,!
ThImob,PhrExch)
The run_HYDRUS_reactions function calls PHREEQC to carry out all chemical andbiological reactions and transfers information from HYDRUS to PHREEQC, the
Get_Concentrations function transfers information back to HYDRUS from PHREEQC, theInitialization brings initial information from the PHREEQC part, and PhreeqcMB calculates
mass balances for the main components.

20
Table 2.1 Description of variables that are used in both HYDRUS-1D and PHREEQC, or newly defined in
HYDRUS-1D.

Fortranvariables
Initialization
NumNP
lMobilTempN
ThOld
ThImob
NSD
NS
Conc
Sorb
tInit

c variablesType

node_number
int
lMobilint
Temperature
float*
Theta
float*
ThImobfloat
max_species_number int
species_number
int
Concetrations
float*
Im_Concetrations
float*TimeInit
float

Get_Concentrations (see definitions above)Conc
Concetrations
float*

Sorb

Im_Concetrations

float*

Conc1

Concetrations

float*

Sorb1

Im_Concetrations

float*

Description
Number of nodes
Mobile-immobile water model
Temperatures
Initial water contents
Immobile water content
Max number of components
Number of components
Concentrations
Immobile concentration
Initial time
Concentrations before chemical
reactions
Immobile concentrations before
chemical reactions
Concentrations after chemical
reactions
Immobile concentrations after
chemical reactions

run_HYDRUS_reactions
t
time
ThNew
Theta
dt
TimeStep
iPrintiPrintTLevelStep_Nofloat
float*
float
int
int

Time
Water content
Time step
Print time flag
Step Number

PhreeqcMB
Conc

real*

Conc1

real*

Sorb

real*

Sorb1

real*

x
PhrExch

real*
real*

Component concentrations in themobile phase before chemical
reactionsComponent concentrations in themobile phase after chemical reactions
Component concentrations in theimmobile phase before chemical
reactionsComponent concentrations in theimmobile phase after chemical
reactionsNodal coordinates
Change in the mass of a component

21

3 Description of data input

3.1 Input data

The following separate input files are required to run the coupled HP1 model:
HYDRUS-1D
selector.in

profile.dat

atmosph.in

PHREEQC
xxxxxxx.xxx

contains the following blocks:
A. Basic Information
B. Water Flow Information
C. Time Information
D. Root Growth Information
E. Heat Transport Information
F. Solute Transport Information
G. Root Water Uptake Information
contains the following block:
H. Nodal Information (except the initial total aqueousconcentrations)contains the following block:I. Atmospheric Information

phreeqc.in

is the database file containing thermodynamic data for aqueous species,
pure phases and exchange reactions (e.g., phreeqc.dat, minteq.dat).
contains information about the geochemical reactions.

species.in

contains a list of master species to be transported.

HP1

The data and format of the different input files and blocks are described in section 10 of the
HYDRUS-1D manual (Šimůnek et al., 1998) for input files related to HYDRUS-1D(selector.in, profile.dat, and atmosph.in), and in Parkhurst and Appelo (1999) for input files
related to PHREEQC. We recommend that users refer to original manuals in order tocomplete the data input. Specific guidelines for the coupled HP1 code are:•

The species.in file contains on the first line the path to the geochemical database, andthen a list of master species to be transported. An example of the species.in file is22

c:\program files\phreeqc\phreeqc.dat
Na
Cl
Ca
Mg
where phreeqc.dat is the database file. This file forms a link between HYDRUS-1D
and PHREEQC. The master species in the file must be written in the same way as they
are defined in the SOLUTION_MASTER_SPECIES block of the phreeqc.dat input
file (see Parkhurst and Appelo, 1999). The order of the master species in the species.in
file refers to solute numbers in the selector.in files for the HYDRUS-1D model. Notethat a check is made whether or not names from the species.in file correspond tomaster species from the geochemical database.
•

When the graphical user interface of HYDRUS-1D (discussed later) is used, thenumber of master species is limited to 10. After the HYDRUS-1D input files are
created with the graphical interface, the number of master species can be increased
manually by making several changes in the Selector.in file described in Table 10.6 ofŠimůnek et al. (1998): increase the variable NS to the number of master species, add
records 9 and 12 for each additional solute, and expand record 16. The Profile.dat andAtmosph.in files (Table 10.8 and 10.9, respectively) must also be expanded for the
specified number of species.

•

We suggest to set solute transport parameters describing exchange with the solid or
gas phases, or describing degradation, equal to zero. It is important to use the same
molecular diffusion coefficient for all master species (see section 2.5). Note that this is
not checked in the current version of the model.
All boundary conditions are specified in the input files of HYDRUS-1D.

••

Initial chemical conditions are to be defined in the phreeqc.in file. This includes
concentrations of all master species. If the initial concentration of a given masterspecies is zero, a very small value (e.g., 1E-20 mol/l) should be given; otherwise it ispossible that the master species will be neglected during the simulation. There is noneed to define the initial conditions of the total aqueous concentrations in the input
files of HYDRUS-1D. Initial conditions for water flow and heat transport are to be
defined in profile.dat.

•

For each node in the HYDRUS-1D finite element mesh, a SOLUTION must be
defined in the phreeqc.in file. The node number at the soil surface is 1 and numbering
increases with depth. Solutions at different depths with the same initial composition
can be grouped in the phreeqc.in input file.

•

Since the phreeqc.in data file is read prior to the profile.dat data file, the amount ofinitial water for each SOLUTION (defined in the profile.dat file) is not known when
the initial solutions are initialized. Therefore, the water content at each node (and

23•

hence for each solution) should be specified for each SOLUTION key word in thephreeqc.in file. Alternatively, the SOLUTION_SPREAD key word can be used todefine solutions with varying amounts of water and temperatures as a function of
depth.
The input file phreeqc.in must contain the keyword TRANSPORT together with theoption -cells indicating the number of cells. No additional information is required
after the keyword (e.g., dispersivity and molecular diffusion are all defined in the
input files of HYDRUS-1D), although –punch_cells and –punch_frequency
can be used to control the output.

The different input files are most conveniently created using the interactive graphic-baseduser interfaces of the original models. These interfaces are public domain and can be
downloaded from the World Wide Web for both the HYDRUS-1D and PHREEQC models.
The HYDRUS-1D model is located at http://www.pc-progress.cz/Fr_Hydrus.htm. The mainpart of the input-files required for the water flow, heat and solute transport parts of the
problem can be constructed using the HYDRUS-1D interface (section B of Šimůnek et al.,1998).
A
graphical
interfacefor
PHREEQC
can
be
found
at
or
at
http://www.geo.vu.nl/users/posv/phreeqc/index.html
http://wwwbrr.cr.usgs.gov/projects/GWC_coupled/phreeqci/. This interface is helpful for
constructing the phreeqc.in input file for the geochemical reactions in the problem.

3.2 Running the model and output

An extra file, path.dat, must be placed in the same directory as the executable of the HP1model (i.e., HP1.exe) and two dynamically linked libraries with HYDRUS-1D (i.e.,
Hydrus.dll) and PHREEQC (i.e., Phreeqc.dll) subroutines and functions. The path.dat filespecifies the path to the input and output file folder. All input files should be placed in thesame folder (except the database file for which the path and name is defined in the species.infile). The simulation is started by activating the executable. Output files are created in thesame folder as the input files. The content of the output files created by HYDRUS-1D aredescribed in section 11 of Šimůnek et al. (1998). The output file of PHREEQC contains userdefined output data (the SELECTED_OUTPUT key word in the input file, see Parkhurst andAppelo, 1999). Part of the created output can be viewed through the interfaces.

24

25

4 Verification problems
In this section we describe different test examples of the coupled HP1 model. We present
examples having different levels of complexity. The following problems were solved with
HP1 for verifying the numerical correctness of the coupling procedure: (i) transport problems
with single or multiple components, subject to sequential first-order decay, will be comparedwith simulations using HYDRUS separately, (ii) multicomponent transport problems forsteady-state flow will be compared with simulations using PHREEQC separately, and (iii)several more complicated problems will be compared with simulations using the CRUNCH
model. The CRUNCH model was a result of further development of the GIMRT/OS3D codesof Steefel and Yabusaki (1996, see http://www.csteefle.com). Several of the verification
examples as well as additional applications for transient flow were previously reported by
Jacques et al. (2002, 2003).

4.1 Modelling the transport of a single component or decay chain

This first group of verification problems concerns the transport of one solute (either inert oradsorbing according to the Freundlich isotherm), or of a sequential first-order decay chain.
Although these problems can be solved directly with the HYDRUS-1D program, we use themhere to verify the coupling of HYDRUS-1D and PHREEQC, with HYDRUS-1D solving thetransport part and PHREEQC the reaction part of the problem. First, we will test if thetransfer of information (i.e., solute concentrations and water contents) between two modulesof the coupled model (i.e., HYDRUS-1D and PHREEQC) is done correctly (section 4.1.1).
Second, we will compare the numerical accuracy of both the coupled HP1 model and
PHREEQC against the stand-alone HYDRUS-1D using a single-component transport
problem involving either equilibrium (using the Freundlich adsorption isotherm) or kinetic(first-order decay) reactions. These geochemical reactions are coupled with transport within
HYDRUS-1D using a one-step or global implicit method (i.e., various reactions are directly
included into the governing transport equations). In contrast, a non-iterative sequential
approach is used in both the PHREEQC and HP1 codes. In section 4.1.2 we will compare
results of these three models (i.e., HYDRUS-1D, PHREEQC, and HP1) for different Peclet
and Courant numbers. Finally, the transport of three sequential first-order decaying
contaminants is simulated for transient flow with both HYDRUS-1D and HP1 (section 4.1.3).
4.1.1

Physical equilibrium and nonequilibrium transport of chloride

In this first section we simply test if the transport of an inert solute (Cl) is correctly simulatedusing HP1 for the following conditions: (i) steady-state flow with no immobile water, (ii)steady-state flow with immobile water, and (iii) transient water flow. Since this example does26

not consider any geochemical reactions, all transport processes in the coupled model (HP1)
are simulated only with the HYDRUS-1D module. However, since solute concentration andwater content values still pass between the two modules, we will evaluate in this example if
the transfer of information between the two components of the coupled model is done
correctly.
Verification Problem 1: Steady-state physical (non)equilibrium transport of chloride(EQCL – NEQCL)
This problem simulates the transport of chloride (i.e., a geochemically inert tracer) during
saturated steady-state flow in a 20-cm long soil core. The saturated hydraulic conductivity Ks
is 1 cm d-1 and the saturated water content is 0.5 cm3cm-3. The following solute parameters
were used: a dispersivity (DL) of 8 cm for both the equilibrium and nonequilibrium cases, andfor the latter case an immobile water content (θim) of 0.1 cm3cm-3 and a first-order exchange
coefficient (ω) of 0.01 d-1. Simulation results obtained with HP1 were found to be identical tothose obtained with HYDRUS-1D, as illustrated by the outflow concentrations in Figure 4.1.

1

HYDRUS-1D
equilibrium
HYDRUS-1D
nonequilibrium
HP1

Concentration (mol/l)

0.8

0.6

0.4

0.2

0
0

10

20

30

40

50

Time (d)Figure 4.1 Outflow concentrations for Verification Problem 1. Full and dashedlines are results for physical equilibrium and nonequilibrium transport,
respectively, obtained with HYDRUS-1D. Dots are results obtained with HP1.

27

Verification Problem 2: Transient physical nonequilibrium transport of chloride(TRANSCL)
In this second problem we simulate the transport of chloride through a 1-m deep soil profile
subject to a transient upper boundary condition given by daily values of precipitation andevaporation over a 300-d period. Physical nonequilibrium (i.e., the presence of immobile
water in the soil profile) was considered in this problem. Note again that all transport
calculations were done by HYDRUS-1D, which means that the test again applies only to thetransfer of concentrations and water contents between HYDRUS-1D and PHREEQC. We

used parameters of the soil hydraulic properties typical for a loamy soil (θr = 0.078 cm3cm-3,

θs = 0.43 cm3cm-3, α = 0.036 cm-1, n = 1.56, and Ks = 24.96 cm d-1 from Carsel and Parish,1988). Solute transport parameters were as follows: a dispersivity DL of 8 cm, an immobile
water content θim of 0.05 cm3cm-3, and a first-order exchange coefficient ω of 0.0125 d-1.
Precipitation and evaporation rates were typical for the Campine region in Belgium. The soil
profile was discretized into 100 elements of 1 cm each. Chloride was applied during the first
53 days of the simulation with a concentration of 0.1 mmol l-1. Results in Figure 4.2 show aperfect match between HYDRUS-1D and HP1. From these results we conclude that thetransfer of water contents and concentrations between the transport and reaction modules wasdone correctly.

Concentration (mol/l)

1E-005

1E-006HYDRUS-1D
HP1

1E-007

1E-008
0

100

200

300

Time (d)Figure 4.2 Outflow concentrations for Verification Problem 2. The full line was generated with HYDRUS-1D,
while dots were obtained with HP1.

28

4.1.2

Transport of a nonlinearly adsorbing solute subject to first-order decay

This section considers numerical simulations of a nonlinearly sorbing solute undergoing firstorder decay during steady-state water flow. We assume only heterogeneous reactions of acontaminant Cont with a sorbing surface Sor. Adsorption is assumed to be instantaneous anddescribed with the Freundlich equation:
S s = K d C nF

(4.1)

where Ss is the adsorbed concentration [MaMs-1], C is the aqueous concentration [MaL-3], Kd isan empirical adsorption coefficient [Ma(1-n)L3nMs-1] and nF is the empirical Freundlich
coefficient [-]. The contaminant Cont is additionally assumed to be subject to first-orderdegradation:

R = µw C(4.2)

where R is the decay rate [MaL-3T-1] and µw is the first-order rate constant for solutes in theliquid phase [T-1].
Modelling nonlinear Freundlich adsorption with PHREEQC
To model instantaneous adsorption using the Freundlich adsorption isotherm, we rewrote Eq.
(4.1) in terms of the amount adsorbed per unit volume of water:()

S w = K d ρ C nF

(4.3)

where Sw is the adsorbed concentration per unit volume of water [MaL-3] and ρ is the bulk
density [MsL-3]. Equation (4.3) corresponds with the following mass action equation:[ SorCont ] = K s ρ = K w( d ) dn
[ Sor ][Cont ] F(4.4a)[SorCont ] = K wd [Sor ][Cont ]n(4.4b)

F

where Kdw is the adsorption constant in mass per unit volume of water [Ma(1-n)L3nL-3]. In
PHREEQC, this sorption reaction is modelled as a surface complexation reaction. The amountof adsorption sites is chosen very large so that [Sor] in Eq. (4.4b) does not change
significantly when the amount of adsorbed species [SorCont] remains small. An outline of the
PHREEQC input file is given in Box 4.1. In this example, the new solution master species

29
Box 4.1 PHREEQC input for the Freundlich adsorption isotherm
1
2

SOLUTION_MASTER_SPECIESCont Cont 0.0 Cont 1.0

3
4

SOLUTION_SPECIESCont = Cont; log_k 0.0

56

SURFACE_MASTER_SPECIES
Sor Sor

7
8
9

SURFACE_SPECIES
Sor = Sor
log_k 010
1112

Sor + nF Cont = SorCont-no_check; -mole_balance SorCont
log_k K

13 SURFACE 1
14 equilibrate 1
15 -no_edl true
16 Sor S 1 S#nF Freundlich term
#K = LOG(KWD /S)#S Total amount of adsorption sites(Cont), its aqueous speciation, and the surface species (Sor) are defined on lines 1-6 with the
PHREEQC-keywords SOLUTION_MASTER_SPECIES, SOLUTION_SPECIES, and
SURFACE_MASTER_SPECIES. The definition of the mass action equation of the Freundlich
isotherm, Eq. (4.4b), is given on lines 10-12 with the PHREEQC-keyword
SURFACE_SPECIES. Two options are included: (1) -no_check since the mole balance ofthe reaction equation (line 10) is not fulfilled due to nonlinearity of the Freundlich isotherm(nF≠1), and (2) -mole_balance to ensure balancing of 1 Sor and 1 Cont. The association
coefficient K is defined by Kdw/[Sor], which is approximately equal to the total amount ofadsorption sites S, if S is chosen very large, e.g. 10100 mol l-1 (put S = 10100 in line 16). As an
example, for Kds = 5, ρ = 1.5, and S = 10100 we obtained a K-value (=log_k) of -99.1249.
Verification Problem 3: Steady-state transport of nonlinearly adsorbed contaminant(STADS)In this problem we consider saturated steady-state water flow and single-component (Cont)transport through a soil column of 1 m length for a period of 1000 d. Very low initialconcentrations of Cont are assumed to be present (10-15 mol l-1) in all cases. Transport andsimulation parameters as well as Peclet numbers (Pe = (q ∆x)/(θ Dw) with ∆x the characteristic
lenght of a finite element), and Courant numbers (Cr = (q ∆t) / (θ ∆x), Steefel and
MacQuarrie, 1996) are given in Table 4.1. HYDRUS-1D, PHREEQC and HP1 simulations
with several different combinations of the Peclet and Courant numbers are compared in30
Table 4.1 Transport and simulation parameters for Verification Problem 3.
Parameters

HYDRUSH-11
Ks (cm d )330.5
θs (cm cm )-11q (cm d )1
2
v (cm d )DL (cm)1-3
1.5
ρ (g cm )5
Kds ()
0.8
nF (-)
Kdw ()
Log(Sor)Log_k∆x (cm)

P-1
0.5
2
1

PHREEQC
P-2
0.5
2
1

1
1
0.8
0.8
7.5
7.5
100
100
-99.125 -99.125

HP1
P-3
0.5
2
1

HP1-1
1
0.5
1
2

HP1-2
1
0.5
1
2

HP1-3
1
0.5
1
2

HP1-4
1
0.5
1
2

1
0.8
7.5
100
-99.125

1
1.5
0.8
7.5
100
-99.125

1
1.5
0.8
7.5
100
-99.125

1
1.5
0.8
7.5
100
-99.125

1
1.5
0.8
7.5
100
-99.125∆tmax (d)

1
1

1
-

0.5
-

0.25
-

1
1

1
0.25

1
0.08

1
0.04

Pe
Cr

1
2

1
1.33

0.5
0.66

0.25
0.33

1
2

1
0.5

1
0.16

1
0.08

Figure 4.3. Depth profiles of aqueous Cont-concentrations are plotted at 3 different times(250, 500, and 750 d). Figure 4.3a compares three simulations carried out with PHREEQC(the solution method is equivalent to SNIA (sequential non-iterative approach) with an
explicit time weighting scheme) and HYDRUS-1D (the global implicit method with implicit
time weighing).
The concentration front for the P-1 space discretization scheme displayed significant
numerical dispersion as compared to the HYDRUS-1D solution. Increasing the number ofcells in PHREEQC, and thus decreasing the Peclet number, produced less numericaldispersion (e.g., P-2 and P-3). However, the simulation time became then very large becauseof the large number of nodal points in the discretization scheme. Results of simulations using
HP1 are compared with HYDRUS-1D in Figure 4.3b. While the same spatial discretization
was used for all simulations, the maximum time step was decreased. HP1 simulation with amaximum time step of 1 d (corresponding to Courant numbers larger than one) provided lessaccurate results as compared to HYDRUS-1D. However, the agreement gradually improved
with lower maximum time steps and Courant numbers, with the lowest two Cr-values
providing almost identical results. HP1 results were found to be better than those of
PHREEQC for the same Peclet and Courant numbers (compare HP-2 and P-1 in Figure 4.3a).
To obtain a similar degree of accuracy for both HP1 and PHREEQC, at least 200 cells wereneeded for the PHREEQC simulation (Figure 4.3a) compared to 100 for HP1. Also thecomputational time for HP1 was significantly smaller as compared to PHREEQC.

31

When the accuracy of the HP1-code is evaluated based on the simulation with HYDRUS-1D,
HP1 will need smaller time steps (up to 25 times smaller) with the same spatial discretization.
Note that a similar analysis can be done when ∆x is decreased in the HP1 simulations.
1

H
HP1-2
P-1
(∆x = 1 cm)
P-2(∆x = 0.5 cm)
P-3(∆x = 0.25 cm)

Concentration (mol/l)

0.8

0.6

0.4

0.2

0

(a)
0

20

40

60

80

100

Depth [cm]
1

H
HP1-1
(∆tmax = 1 d)Concentration (mol/l)

0.8

HP1-2
(∆tmax = 0.25 d)
HP1-3(∆tmax = 0.08 d)

0.6

HP1-4
(∆tmax = 0.04 d)
0.4

0.2

0

(b)
0

20

40

60

80

100

Depth [cm]Figure 4.3 Depth profiles of Cont after 250, 500, and 750 d for different simulations defined in Table 4.1 for
Verification Problem 3. Tests consider effects of grid size for PHREEQC (a), and effects of maximum time steps
for HP1 (b).

32Modelling first-order decay with PHREEQC
First-order decay (Eq. (4.2)) in PHREEQC is modelled using the keywords RATES andKINETICS. The implementation in PHREEQC is shown in Box 4.2. The kinetic reaction
itself is defined under the keyword RATES (lines 1-8). In this example, the reaction is called'degradation' (line 2). Between -start and -end a Basic-program is written for the kinetic
reaction of the 'phase' degradation consisting of standard Basic-statements (e.g., here onlyrem) and special Basic-statements for PHREEQC (e.g., MOL, SAVE, TOT, and TIME). The
first statement of the Basic-program (line 4) is only a comment indicating the meaning of thefirst parameter. The second statement evaluates the rate equation (Eq. (4.2)) with parm(1)being the value of µw, TOT("water") the amount of water in the cell (which is included here tohave an equivalent problem as in HYDRUS-1D), and MOL("Cont") the molality of the solute.
The third statement (line 6) integrates the rate over the subinterval with the special variable
TIME. Finally, the moles of reaction during the time interval are saved with the last special
statement SAVE. Note the negative sign on line 6 that results in a negative amount of moles
saved in the last statement. In general, a positive sign represents decreasing amounts of aphase (i.e., dissolution), whereas a negative sign results in precipitation of that phase.
Consequently, elements will enter the solution in the former case (dissolution) and will be
removed from the solution in the latter case (i.e., precipitation, degradation, or decay). In thisexample the imaginary phase 'degradation' is precipitating. This is done to prevent thecessation of the kinetic reaction (i.e., when the phase 'degradation' is completely removed
from the system).
The second keyword in Box 4.2 (KINETICS) defines the names of the rate expressions
related to a specific cell. In this example we have one rate expression called 'degradation'.
Since 'degradation' is used here as an imaginary phase (and not a phase defined in the
database), the option -formula is used to define the elements produced (i.e., when the
Box 4.2 PHREEQC input for first-order decay reactions

1
2
3
4
5
6
7
8

RATES
degradation-start10 rem parm(1) first-order degradation coefficient (sec-1)
20 rate = parm(1)*TOT("water")*MOL("Cont")
30 moles = - rate * TIME
40 SAVE moles-end

9 KINETICS 1
10 degradation
11 -formula Cont 1; -parms 2.3148E-7

33

product of the stoichiometric element coefficient (i.e., 1 in the formula option) multiplied by
the moles of reaction during a particular time step is positive) or consumed (i.e., when the
product is negative) during the kinetic reaction. Since in this example the coefficient for theelement Cont is 1 (line 11) and the reaction progress is negative, the concentration of Cont
will decrease with the formation of the imaginary phase 'degradation'. Note that we could
write the input also with a negative coefficient and a positive reaction progress (i.e.,
dissolution of the phase). However, in that case we could reach complete dissolution of thephase and, consequently, termination of the decay reactions. The last option under KINETICS
in Box 4.2 is -parms for the purpose of defining parameters in the rate expression.
Verification Problem 4: Steady-state transport of nonlinearly adsorbing contaminant with
first-order decay (STDECAY)
This verification problem is the same as the previous one (i.e., Verification Problem 3) but
with the inclusion of first-order decay. All transport and simulation parameters are equal tothose of Verification Problem 3 (see Table 4.1). However, we did not perform the P-3 andHP-1 simulations because the former (P-3) was computationally inefficient and the latter (HP1) produced relatively high amount of numerical dispersion. The first-order rate constant waschosen equal to 0.2 d-1 (or 2.31 10-6 sec-1). We again compare HYDRUS-1D results with
PHREEQC simulations using different spatial discretizations, and HP1 simulations using
different maximum time steps (Figure 4.4). The latter code is the reference.

Results of simulations with mixed equilibrium (Freundlich adsorption) and kinetic (first-orderdecay) reactions were very similar for all three models, thus supporting the same conclusions
as in the previous case. The transport of a nonlinearly adsorbing, first-order decaying
contaminant was accurately modelled with HP1 when Courant numbers were reasonably
small (Cr smaller than 0.2, i.e. 10 times smaller than Cr from HYDRUS-1D). Furthermore,
results obtained with HP1 were somewhat better compared to results obtained with
PHREEQC for the same Peclet and Courant numbers.

34

1

H
HP1-2
P-1
(∆x = 1 cm)
P-2(∆x = 0.5 cm)

Concentration (mol/l)

0.8

0.6

0.4

0.2

0

(a)
0

20

40

60

80

100

Depth [cm]
1

H
HP1-2
(∆tmax = 0.25 d)Concentration (mol/l)

0.8

HP1-3
(∆tmax = 0.08 d)
HP1-4(∆tmax = 0.04 d)

0.6

0.4

0.2

0

(b)
0

20

40

60

80

100

Depth [cm]Figure 4.4 Depth profiles of Cont after 250, 500, and 750 d for different simulations defined in Table 4.1(Verification Problem 4 with a first-order decay coefficient of 0.2 d-1).

35

4.1.3

Transport of first-order decay chain of nonlinearly adsorbing contaminants
during transient flow

Problem definition
In this example we consider the transport of three (non)linearly adsorbing contaminants,
Conta, Contb, and Contc that are involved in a sequential first-order decay chain defined as:µ'

w,k ,conta
⎯⎯⎯⎯
⎯
→

Conta
↓( Kd ,nF )conta

SoraContaµ'

w, k ,contb
⎯⎯⎯⎯
⎯
→

Contb
↓( Kd ,nF )contb

SoraContbµ

w, k ,contc
⎯⎯⎯⎯
⎯
→

Contc
↓( Kd ,nF )

(4.5)

contc

SoraContc

where µw',k are the first-order rate constants connecting two contaminants, µw,k is the firstorder rate constant for Contc, Kd and nF are the Freundlich isotherm parameters for the threecontaminants, and SoraConta, SorbContb, and SorcContc are the three surface species related
to Conta, Contb, and Contc on three surfaces Sora, Sorb, and Sorc, respectively. Reaction
parameters for the three contaminants are given in Table 4.2. Model simulations were carriedout for a 1-m deep homogeneous soil profile during 1000 d assuming transient flow. Upper
boundary conditions were taken as daily precipitation rates representative of the Campine
Region in Belgium. Evaporation was neglected during the simulations. The lower boundary
condition was defined as free drainage. A uniform initial pressure head of –60 cm wasassumed for the entire soil profile. For solute transport, the following initial and boundary
conditions were considered: (1) low initial concentrations (10-15 mole l-1) for all threecontaminants, (2) third-type solute fluxes as the top boundary conditions with 1, 0.1, and 0
mol l-1 for Conta, Contb, and Contc, respectively, and (3) zero-gradient bottom boundary
condition.
Modelling the decay chain with HP1
The decay chain was modelled by defining the rates of three kinetic reactions (degrad_conta,
degrad_contb, and degrad_contc) to model the degradation of solutes Conta, Contb, andContc, respectively. The definition of degradation rates is similar as in Verification Problem
4. To model the degradation of Conta into Contb, the option -formula for keyword
KINETICS (line 11) is changed to -formula Conta 1 Contb -1, which indicates that when
1 mole of the imaginary phase degrad_conta precipitates (a negative rate), 1 mole of Conta
disappears and 1 mole of Contb appears. A similar kinetic rate reaction is written for thetransformations of Contb into Contc.

36

Verification Problem 5: First-order decay chain of nonlinearly adsorbing contaminants
during unsteady flow (SEASONCHAIN)
We now compare HYDRUS-1D results with four HP1 runs using different time-stepping
schemes. The first two simulations were carried out with maximum time steps of 30,000 s(HP1-1) and 10,000 s (HP1-2) in order to obtain Courant numbers smaller than 1 (based onflow velocities at the top boundary). For the other two simulations we used the stabilisation
criterion defined by Perrochet and Berod (1993) and implemented in HYDRUS-1D:

Pe • Crv ≤ ω s(4.6)

where ωs is the performance index [-]. Simulations were performed with ωs = 0.4 (HP1-3) and
0.2 (HP1-4). The HYDRUS-1D module of HP1 automatically adjust the time steppeing tofulfil this stabilization criterion. Transport parameters are given in Table 4.2 and simulation
parameters (i.e., temporal and spatial steps, and performance index) in Table 4.3.
Distributions of Conta, Contb, and Contc versus depth are shown in Figure 4.5 at three times(250, 500, and 1000 d). Simulations obtained with HP1-1 were identical to the HYDRUS-1D
results for Conta and Contb, except at the leading edge of the concentration profile for the
first print time for Conta. Results for HP1-3 and HP1-4 shown a better agreement between thesimulations than those of HP1-1 and HP1-2. Concentration profiles for Contc were also invery good agreement between different runs, except for the peak concentrations, which were
smaller than those obtained with HYDRUS-1D. Note again that simulations for HP1-3 andHP1-4 runs shows a better agreement than those for HP1-1 and HP1-2 runs.

Table 4.2 Soil hydraulic, transport, and reaction parameters for Verification Problem 5.
Parameter
Hydraulic parametersθr (cm cm )θs (cm3cm-3)
α (cm-1)
3

-3

n (-)
Ks (cm d-1)
Transport parameters

DL (cm)ρ (g cm )-3

Reaction parameters

Kd
nF (-)
Log_k *

µw' ,k (d-1)µw,k (d-1)
*

Value
0.078
0.43
0.036
1.56
24.96

Conta

Contb

ContC

0.5
1
-100.12
0.005
-

2.5
0.9
-99.43
0.06
-

5
0.8
-99.120.02

11.5

Input for HP1 assuming the amount of adsorption sites S = 10100.

37
Table 4.3 Input simulation parameters and numerical results for Verification Problem 5.
Parameter∆tinit (sec)∆tmin (sec)∆tmax (sec)ωs∆x (cm)N**∆tmin (sec)

HYDRUS1D*

864
0.0864
86400
2
1

HP1
HP1-1
HP1-2
HP1-3
Input simulation parameters
864
864
864
0.08640.08640.0864
30000
10000
86400
2
2
0.4
1
1
1

9070
864
86400
1.09
0.47

9960
864
30000
1.20
0.99

Simulation results
13585
864
10000
1.17
0.89

9761
864
86400
1.08
0.47

HP1-4
864
0.0864
86400
0.2
1
11775
864
57600
1.04
0.24

∆tmax (sec)
Pemax
Crmax
*
Note that Courant numbers are calculated differently (R is included in the denominator).
**
Number of time steps for the complete simulation.

An overview of different numerical criteria is given in Table 4.3. The HP1 simulation
parameters (temporal and spatial steps) for this particular problem had no significant
influence on the accuracy of the simulation results. For purposes of computational efficiency,
it is however always better to use the stability criterion based on the performance index ωs.
This stability criterion ensures that small time steps are used when flow is rapid, but allows
for larger time steps when flow is slow. This is done by modifying time steps such that the
Courant number obeys the inequality given by Eq. (4.6). Note that HP1 requires a lower
performance index (five times smaller) than HYDRUS-1D to obtain similar results.

38
1

1

Conta

Conta
0.8

Concentration (mol/l)Concentration (mol/l)

0.8

500 d
0.6

1000 d
0.4

250 d
0.2

0

0.6

0.4

0.2

(a)
0

0
20

40

60

80

100

(d)
0

20

Depth (cm)Contb

1000 d

0.06

0.04

0.02

100

0.06

0.04

(b)
0

0
20

40

60

Contb

0.02

500 d

250 d

80

100

(e)
0

20

Depth (cm)

40

60

80

100

Depth (cm)

0.1

0.1

Contc
1000 d

0.06

0.04

500 d

0.02

(c)
0

0.06

0.04

0.02

250 d
0
20

Contc

0.08

Concentration (mol/l)

0.08

Concentration (mol/l)

80

0.08

Concentration (mol/l)Concentration (mol/l)

0.08

0

60

0.1

0.1

0

40

Depth (cm)

40

60

80

100

(f)
0

20

Depth (cm)

40

60

80

100

Depth (cm)HYDRUS-1D

HP1-1(∆tmax = 30000 s, ωs = 2, left)/ HP1-3 (∆tmax = 86400 s, ωs = 0.4, right)

HP1-2 (∆tmax = 10000 s, ωs = 2, left) / HP1-4 (∆tmax = 86400 s, ωs = 0.2, right)
Figure 4.5 Distribution of Conta, Contb, and Contc versus depth after 250, 500, and 1000 d for different
simulations as defined in Table 4.3 (Verification Problem 5).

39

4.2 Multicomponent transport during steady-state flow

In this section we test different keywords of the PHREEQC program while using the coupled
HP1 model. The following specific reactions are tested: (1) cation exchange (EXCHANGE)
and (2) equilibrium and kinetic dissolution / precipitation of minerals (EQUILIBRIUMPHASES. The key words TRANSPORT and KINETICS were tested in the previous section.

4.2.1

Cation exchange reactions

Verification Problem 6: Transport of heavy metals subject to multiple cation exchange(CATEXCH)
In this problem, the transport of ten ions (Al, Br, Ca, Cd, Cl, K, Mg, Na, Pb, Zn) through asoil column is modelled. Initial and inflow concentrations of the ions are given in Table 4.4.
The cation exchange capacity is equal to 0.011 mol / cell. The soil core has a length of 8 cm
and is discretized into 40 cells of 0.2 cm. The flow velocity is 2 cm d-1 and the dispersivity is
2 cm. Simulations were performed for 15 days. The maximum time step used in HP1 was
0.015 d.

Distribution versus depth after three days and outflow concentrations during 15 days forselected output variables (ions and pH) are shown in Figure 4.6 and Figure 4.7, respectively,
for simulations carried out with PHREEQC and HP1. Only small differences between the twosimulations were present in the concentration profiles (Figure 4.6). Deviations may be due toincreased numerical dispersion in PHREEQC, as noted in the Cl concentration profile. As
discussed in the previous section, simulations with PHREEQC showed a larger dispersion
compared to the simulations with HYDRUS-1D and HP1 if the same spatial discretization
was used (e.g., Figure 4.3a and Figure 4.4a). If the spatial discretization in PHREEQC was
Table 4.4 Initial and inflow concentrations for Verification Problem 6.
mmol l-1

Initial pore water
Initial concentrations of Inflow concentrations
composition
exchangeable cations*
Al
0.5
0.92
0.1
-4
-4
Ca
1 10
2.88 10
5
Cd
0.09
0.17
0
K
2
1.06
0
Mg
0.75
1.36
1
Na
6
0.62
0
Pb
0.1
0.34
0
Zn
0.25
0.76
0
+
Br
11
3.7
-4
Cl
1 10
10
pH
5.5
2.9
*
Calculated in equilibrium with the initial pore water composition.
+
Br is used to impose a charge balance at pH of 5.5.

40

reduced to 0.0005 m (i.e., 160 cells), PHREEQC converged towards the numericalsimulations obtained with HP1 (Figure 4.6). Small differences between HP1 and PHREEQC
exist when a spatial discretization of ∆x = 0.002 m is used (Figure 4.7). If ∆x is decreased to
0.0005 m, no significant differences between the two simulations were present at the end of
the column and in the outflow concentrations. We conclude that the keyword EXCHANGE is
correctly coupled in HP1.

41
6

HP1 (∆x = 0.002 m)

5pH

PHREEQC (∆x = 0.002 m)
PHREEQC (∆x = 0.0005 m)

4

3

2
0

2

4

6

8

0.01
0.008

Pore water concentration (mol/l)Pore water concentration (mol/l)

Depth (cm)

Cl

0.006
0.004

Ca
0.002
0
0

2

4

6

0.0008Zn

0.00060.0004Cd
0.0002

0

8

0

2

0.01
0.008

CaX2
0.006
0.004
0.002

NaX

0
0

2

4

Depth (cm)

4

6

8

Depth (cm)Concentration on exchanger (mol/l)Concentration on exchanger (mol/l)

Depth (cm)

6

8

0.002

ZnX2

0.00160.00120.0008CdX2
0.0004
0
0

2

4

6

8

Depth (cm)Figure 4.6 Distribution of pH, dissolved Ca, Cl, Zn, and Cd concentrations, exchangeable Ca, Na, Zn, and Cd
concentrations versus depth after three days of infiltration for Verification Problem 6.

426

HP1 (∆x = 0.002 m)
HP1 (∆x = 0.0005 m)

5pH

PHREEQC (∆x = 0.002 m)
PHREEQC (∆x = 0.0005 m)

4

3

2
0

3

6

9

12

15

Time (days)
Cl

0.01

0.001

Concentration (mol/l)Concentration (mol/l)

Zn
0.008
0.006
0.004

Ca

0.002
0

0.00080.0006Cd

0.00040.000200369Time (days)

12

15

0

3

6

9

12

15

Time (days)

Figure 4.7 pH, Ca, Cl, Zn, and Cd concentrations in the effluent for Verification Problem 6 (Full line: HP1,dashed line: PHREEQC).

4.2.2

Mineral dissolution

Verification Problem 7: Transport with mineral dissolution (MINDIS)
A 100-cm long soil column, consisting of amorphous SiO2 and gibbsite (Al(OH)3), is leachedwith a solution containing 5.10-7 mol l-1 Si, 1.10-6 mol l-1 Al, and 1.10-3 mol l-1 Na (to obtainan inflow pH of 11.15). Initial concentrations were 1.76 10-3 mol l-1 Si, 8.87 10-9 mol l-1 Al,
and 1 10-12 mol l-1 Na, corresponding to a pH of 6.33. In each 1-cm cell, 0.015 mol
amorphous SiO2 and 0.005 mol gibbsite is present. The flow velocity is 2 cm/day and thedispersivity is 1 cm. Simulations are again carried out with both PHREEQC and HP1.

Distribution of pH, Si, Al, amorphous SiO2, and gibbsite versus depth after 150 days ofsimulation are presented in Figure 4.8. No significant differences between PHREEQC andPH1 were apparent for pH, Si, Al, amorphous SiO2 and gibbsite. The keyword
EQUILIBRIUM-PHASES hence was coupled correctly in HP1.

43

12

pH

11

HP1
PHREEQC
10

9

8
0

20

40

60

80

100

0.003

Dissolved concentration (mol/l)Dissolved concentration (mol/l)

Depth (cm)

Si

0.002

0.001

0
0

20

40

60

80

0.0002Al

0.000160.000128E-0054E-005
0

100

0

20

Depth (cm)
0.016

60

80

100

0.006

Amount (mol/l)

0.012

Amount (mol)40

Depth (cm)

0.008

0.004

0.004

0.002

amorf SiO2
0

Gibbsite
0

0

20

40

60

Depth (cm)

80

100

0

20

40

60

80

100

Depth (cm)Figure 4.8 Distribution of pH, Si, Al, amorphous SiO2, and Gibbsite versus depth after 150 days of infiltration
for Verification Problem 7 (Full line: HP1, dashed line: PHREEQC).

44

4.3

More complicated verification problems of HP1 model

Examples presented in this chapter compare the combined HP1 model against a different
computer program, CRUNCH, that also simulates multicomponent reactive transport inporous media, with the limitation that only steady-state flow can be invoked. CRUNCH isbased on the GIMRT/OS3D package (Steefel and Yabuski, 1996; Steefel, 2000). The
geochemical reactions and transport in CRUNCH are coupled in one of two ways: (1) a global
implicit approach (GIMRT) that simultaneously solves the transport and reaction equations,or (2) SNIA. GIMRT generally leads to smaller numerical errors. A comparison between
HP1 and CRUNCH-GIMRT allows one to assess numerical discretization errors of the SNIA
coupling as a function of the maximum Courant number, Cr.

4.3.1

Heavy metal transport in a medium with a pH-dependent cation exchange
complex (Verification Problem 8 – MCATEXCH)

4.3.1.1 Problem definition and governing chemical reactions
Cation adsorption to negatively charged soil solid phases can greatly affect the migration of
cations in soils. In HYDRUS-1D, equilibrium isotherms, such as the linear, Freundlich or
Langmuir isotherms, describe the adsorption/desorption of cations. the use of such isotherms
assumes that the adsorption of a particular cation is independent of the presence of otherelements in the soil solution or on the soil solid phases. Their coefficients are constant andindependent of pH, other cations, complexing elements, and ligands in the soil solution.
Unlike HYDRUS-1D, the coupled HP1 code can include the effect of these factors on
adsorption, and consequently, on the migration of multiple cations, by using the ion-exchange
reactions of PHREEQC.

In this example we consider the transport of several major cations (Na, K, Ca, and Mg) andthree heavy metals (Cd, Zn, Pb) through a 50-cm-deep multi-layered soil profile during
unsaturated steady-state flow. Each soil layer has different soil hydraulic properties andcation exchange capacities (CEC) (Table 4.5). The top 28 cm of the soil is assumed to becontaminated with the three heavy metals (initial pH 8.5), while an acid metal-free solution(pH 3) infiltrates into the soil (Table 4.6). Assuming that the cation exchange complex isassociated solely with organic matter, CEC increases significantly with increasing pH due tothe acid-base properties of its functional groups. The higher the pH, the more functionalgroups of the organic matter are deprotonated and thus the higher the cation exchange
capacity. This behaviour is represented by a multi-site cation exchange complex consistingof six sites, each having a different selectivity coefficient for the exchange of protons (seeAppelo et al., 1998). Finally, chloride is present in the soil solutions resulting in the formationof aqueous complexes with the heavy metals.

45

Table 4.5 Soil hydraulic properties and cation exchange capacities of five soil layers.
Horizon
A
E
Bh1
Bh2
Bh/C

Layer
thickness
(cm)
13
10
5
5
17

θr

θs

0.065
0.035
0.042
0.044
0.039

0.476
0.416
0.472
0.455
0.464

α

n

Ks
(cm day-1)

l

0.016
0.015
0.016
0.028
0.023

1.94
3.21
1.52
2.01
2.99

93
311
39
860
1198

0.5
0.5
0.5
0.5
0.5

(cm-1)

Cation exchange
capacity(eq / 1000 cm³ soil)
0.01830.01140.06640.05420.0116

Table 4.6 pH and solution concentrations used in the simulation (µmol l-1).
Solution
pH
Na*
K
Ca
Mg
0-28 cm depth
8.5
401.9
120
98
5
28-50 cm depth
8.5
454.0
120
98
5
Applied water
3.5
127.5
120
98
5
*
Concentration of Na is adjusted to obtain the desired pH.

Br
780
780
780

Cl
0
0
0

Cd
0.8
0.0
0.0

Pb
2.5
0
0

Zn
50
0
0

The soil profile is assumed to contain five distinct layers with different soil hydraulicproperties and cation exchange capacities. Table 4.5 gives thicknesses of the different
horizons, parameters of van Genuchten’s equations for the water retention and hydraulicconductivity functions, and the total cation exchange capacities. The higher exchange
capacities of the Bh1 and Bh2 horizons reflect their enrichment with immobilized organic
matter. Flow is assumed to be steady at a constant flux of 0.05 m day-1 (18.25 m year-1),
which causes the soil profile to be unsaturated (water contents vary between 0.37 and 0.15 asa function of depth). The bottom boundary condition for water flow is free drainage.
HYDRUS-1D was used to calculate the steady-state water content profile corresponding tothese boundary conditions. The dispersivity and diffusion coefficient were taken to be 0.05 m
and 9.2 10-10 m2 s-1, respectively.
Table 4.7 Overview of aqueous equilibrium reactions and corresponding equilibrium constants (data fromphreeqc.dat database, Parkhurst and Appelo, 1999).
Nr
(1)
(2)
(3)
(4)
(5)

Aqueous speciation reaction
H2O = OH- + H++
Na + H2O = NaOH + H+
K+ + H2O = KOH + H+
Ca2+ + H2O = CaOH+ + H+
Mg2+ + H2O = MgOH+ + H+

(6)
(7)
(8)
(9)
(10)
(11)
(12)
(13)

X2+ + H2O = XOH+ + H+
2+
X + 2 H2O = X(OH)2 + 2 H+
X2+ + 3 H2O = X(OH)3- + 3 H+
X2+ + 4 H2O = X(OH)42- + 4 H+
X2+ + Cl- = XCl+
X2+ + 2 Cl- = XCl2
X2+ + 3 Cl- = XCl3X2+ + 4 Cl- = XCl42-

Cd-10.08-20.35-33.30-47.35
1.98
2.60
2.40
-

Log_K
-14
-14.18-14.46-12.78-11.44Pb-7.71-17.12-28.06-39.70
1.60
1.80
1.70
1.38

Zn
-8.96
-16.90-28.40-41.20
0.43
0.45
0.5
0.2

46
Table 4.8 Log K parameters for multi-site exchange complex.
Yexchanger(1)
HY(2)

NaY
KY
MgY2
CaY2
CdY2
PbY2
ZnY2
-1.0
-0.3
-0.4
-0.2
-0.2
0.05
-0.2
HYa
HYb
HYc
HYd
HYe
HYf
1.65
3.3
4.95
6.85
9.6
12.35
(1)
The value for NaY is taken from Appelo et al. (1998). Values for the other complexes are taken from thephreeqc.dat database (Parkhurst and Appelo, 1999) and adapted relative to the K for NaY.
(2)
Values taken from Appelo et al. (1998).

An overview of the considered aqueous equilibrium reactions is given in Table 4.7. The roleof chloride as a complexing agent is described by reactions (10) through (13). Other
geochemical reactions that are considered are the heterogeneous multi-site ion-exchange
reactions. The exchange coefficients for major cations and heavy metals were assumed to be
the same for all exchange sites. Table 4.8 gives parameters for this multi-site exchange
complex.

4.3.1.2 Database and model input
HP1
Chemical reactions and their equilibrium constants are stored in the database phreeqc.dat.
The following keywords are used for this problem: SOLUTION_MASTER_SPECIES,
SOLUTION_SPECIES, EXCHANGE_MASTER_SPECIES, and EXCHANGE_SPECIES. The
input file for water flow and solute transport is relatively straightforward and can be easily
created the standard way with the graphical interface of HYDRUS-1D. The problem involves
9 components that are transported, i.e., Na, K, Mg, Ca, Cl, Br, Cd, Zn, and Pb. Onlyparameters related to solute transport (i.e., the dispersivity and the aqueous diffusion
coefficient) are defined in HYDRUS-1D. All other solute transport parameters, except for theFreundlich exponent which is equal to one, are set to zero. This information is stored in the
HYDRUS-1D input files Selector.in and Profile.dat. Note that these files also contain
information about the steady-state flux (0.05 cm day-1), the boundary concentrations of thenine components in the incoming water, and spatial and temporal discretization parameters of
the numerical problem. The link between HYDRUS-1D and PHREEQC is defined in theSpecies.in input file that contains names of nine elements, the names of which must be thesame as in the Phreeqc.dat database. Finally, the initial solutions and exchange complexes aredefined in the Phreeqc.in input file. Since the transport problem involves variable watercontents with depth, a SOLUTION-keyword needs to be defined for nearly each cell. In
addition, the cation exchange complex for each soil layer must be defined using the keyword
EXCHANGE, with each layer containing six exchange sites. Also included in the input file are
the keywords TRANSPORT to indicate that HYDRUS-1D will be used for transportmodelling, and SELECTED_OUTPUT to specify the desired output. Details about the
keywords used in Phreeqc.in can be found in the PHREEQC-2 manual (Parkhurst andAppelo, 1999).

47

CRUNCHThe database contains chemical reactions and equilibrium constants for the aqueous
speciation reactions (between end of primary and end of secondary) and theexchange reactions (between begin of exchange and end of exchange). The input
file reactive.in defines the initial conditions for the five layers and the inflowing water(CONDITION), primary, secondary, and exchange species (PRIMARY_SPECIES,SECONDARY_SPECIES, and ION_EXCHANGE), the discretization, and the initial andboundary conditions (DISCRETIZATION, INITIAL_CONDITIONS, and BOUNDARYCONDITIONS), flow and transport properties (FLOW, and TRANSPORT), and output(OUTPUT). The keyword POROSITY defines the porosity. Since the porosity is different foreach soil layer, porosities are read from the porosity1.dat input file. Numerical issues forsolving the equations are defined in RUNTIME. The read_saturation identifies the filecontaining the saturation degree for each cell. Information for this file (saturation1.dat) was
obtained from an initial HYDRUS-1D run to find the steady-state water content profile.

4.3.1.3 Comparison between CRUNCH-GIMRT, CRUNCH-OS, and HP1
Simulations were performed with three different models: (1) the global implicit option in the
CRUNCH-model (CRUNCH-GIMRT), (2) the sequential non-iterative option in the
CRUNCH-model (CRUNCH-OS), and (3) the coupled HP1 model which also runs in asequential non-iterative mode. The maximum time step for the latter two simulations waschosen such that both were run at a similar Courant number (0.5). Figure 4.9 shows selected
simulation results.

Infiltration of the low-pH solution causes a decrease in the cation exchange capacity (i.e., anincrease in the number of protonated sites on the six cation exchange complexes) (see Figure4.9e). This relation between pH in the soil solution and the deprotonated sites on the multisite cation exchange complex is apparent when comparing Figure 4.9c and Figure 4.9e. This
leads to a desorption of Cd (and other heavy metals) and their subsequent leaching from thesoil profile. Cd leaching peaks after about 0.3 y (Figure 4.9b), with most Cd leached from theprofile after 1 y.
Results obtained with CRUNCH-OS and HP1 using Cr = 0.5 showed very good agreement,
especially for the outflow curves. This indicates that the coupling between HYDRUS-1D and
PHREEQC was done correctly since the coupled model was tested in this relatively
complicated example with a completely independent reactive transport model (in contrast to
the verification problems discussed earlier in sections 4.1 and 4.2). However, the pH-outflow(Figure 4.9a) and the concentration distributions versus depth (Figure 4.9c-f) showed a small

48

increase in numerical dispersion as compared to simulations using the global implicit
approach (CRUNCH-GIMRT).
9

1.2

(a)

(b)

Cd (µmol kgw -1 )

8

pH

7
6
5

0.8

0.4

4
3

0
0

0.2

0.4

0.6

0.8

1

0

0.2

Time (year)
9

1.2

(c)

0.6

0.8

1

(d)

8

Cd (µmol kgw-1)

0.3 y

7

pH

0.4

Time (year)

6
5

0.5 y
0.8

0.7 y
0.4

4
3

0
0

0.1

0.2

0.3

0.4

0.5

0

0.1

Depth (m)2.5x10-3(e)

0.5

Cd-sites / total CEC

1 - H-sites / total CEC

0.6

0.4
0.3
0.2

0.0x100
0.3

0.5

0.4

0.5

Depth (m)

0.4

0.5

CRUNCH-GIMRTCRUNCH-OS(f)

1.0x10-3

0
0.2

0.4

1.5x10-3

5.0x10-4

0.1

0.3

2.0x10-3

0.1

0

0.2

Depth (m)

0

0.1

0.2

0.3

Depth (m)HP1 (Cr 0.5)
HP1 (Cr 0.1)

Figure 4.9 Selected results for cation and heavy metal transport with a pH-dependent cation exchange
complex. (a) and (b) pH and Cd concentrations in outflow at 50 cm depth, respectively; (c) and (d)distributions of pH and Cd concentrations versus depth after 0.3, 0.5, and 0.7 y, respectively; (e) and (f)distribution of the fraction of deprotonated sites (1 - H-sites) and sites with Cd versus depth after 0.3, 0.5,
and 0.7 y, respectively.

49

Furthermore, at the boundary between the fourth and fifth layer, the pH shows a small
increase for both models using the sequential non-iterative solution approach. This increase
seems to be an artefact of the implemented time step. Reducing Cr to 0.1 in HP1 produced
very good agreement with CRUNCH-GIMRT, with the artefact being no longer present(Figure 4.9). Overall, differences in results obtained with Courant numbers of 0.1 (at least oneorder of magnitude smaller than typical values for Cr) and 0.5 were quite small (especially forbreakthrough curves), with both runs providing acceptable results.
To illustrate the effect of a decreasing cation exchange capacity with decreasing pH on Cd
transport, we changed concentrations of the infiltrating water to those found in the 28–50 cm
horizon, i.e., with a pH of 8.5 (see Table 4.6). Simulations were carried out with CRUNCHGIMRT and HP1 using a Courant number of 0.5; this since no considerable differences wereobserved for the Cd breakthrough curve in previous calculations. As shown in Figure 4.10,
the Cd breakthrough is significantly retarded compared to that produced with the low-pH
infiltration water. Finally, to evaluate the effect of heavy metals complexation with Cl- ions,
we increased Cl concentration in the infiltrating water to 780 µmol per kilogram of water (anddecreased Br concentration to 0 µmol per kilogram of water). Only a small increase in the Cd
breakthrough was obtained (Figure 4.10) compared with the previous case. All three cases
produced essentially the same results using both CRUNCH-GIMRT and HP1.

1.2

Cd (µmol kgw -1)

pH 3.5CRUNCH-GIMRT
HP1 (Cr 0.5)

0.8

pH 8.5 + Cl
0.4

pH 8.5

0
0

0.2

0.4

0.6

0.8

1

Time (y)Figure 4.10 Effect of pH and Cl concentration on Cd concentrations of leaching water at the 50cm depth.

50

4.3.2

Infiltration of a hyperalkaline solution in a clay sample (Verification Problem 9 –ALKALINE)

4.3.2.1 Problem definition and governing chemical reactions

Adler (2001) described a column experiment, in which a high-pH-Na-Ca-OH solution
infiltrated into a core containing Opalinus Clay. This clay is investigated in Switzerland as apotential host formation for long term disposal of high level radioactive waste (Thury andBossart, 1999). Adler (2001) also defined a conceptual geochemical model to describe theoutflow concentrations of this experiment. A slightly adapted version of this conceptual
geochemical model is used here as a benchmark comparison between CRUNCH and HP11.
Infiltration of a high pH plume into the clay core leads to different reactions. First, the cationexchange complex associated with the clay minerals will interact with the compositionalchange in the aqueous solution. Due to the high pH, primary minerals in the clay will become
unstable and may dissolve. Dissolution of primary minerals is described with kinetic
dissolution rate equations. The infiltration of Na-Ca-rich water and the increase of Al and Si
from mineral dissolution will cause the precipitation of secondary minerals.
The infiltration of a hyperalkaline solution in a 7.4-cm-long clay core was simulated for aperiod of 1.1 year. The flow domain was discretized in 100 cells of 7.4 10-2 cm each. The
effective porosity of the Opalinus Clay was kept constant at a value of 0.13 during thesimulations. A constant flux of 2.403 10-9 m sec-1 (7.58 cm y-1) was applied, resulting in asolute transport velocity of 1.85 10-8 m sec-1. The diffusion coefficient was 0.5 10-9 m sec-2
and the dispersivity zero.
The mineralogical composition of Opalinus Clay is given in Table 4.9. Due to the intrusion of
the high pH plume, primary minerals of the clay will become unstable and gradually dissolve.
The dissolution rate depends on the chemical composition of the aqueous phase and is
described with a ‘transition state’ based dissolution model (Aagaard and Helgeson, 1981;
Parkhurst and Appelo, 1999):⎛ M ⎞Rm = Am,0 ⎜ m ⎟⎝ M m ,0 ⎠

2/3{OH } k [1 − Q / K ]
− n

0,m

m

m

(4.7)

where Rm is the dissolution rate of mineral m (mol sec-1), Am,0 is the initial reactive surface
area of mineral m (m2 for HP1, m2 m-3 for CRUNCH), (Mm / Mm,0)2/3 is a factor accounting for
1

The conceptual geochemical model used here is somewhat simplified compared to the model described byAdler (2001, p. 55-60). In this report, the same model is used in both CRUNCH and HP1.

51

Table 4.9 Mineralogical reactions, Log(K), molar volume, volume percent and moles present in Opalinus
Clay, and secondary minerals considered.
mineralReaction

Log K

Primary minerals
Al2Si2O5(OH)4 + 6 H+ = 2 Al+3 + 2 SiO2 + 5 H2O
6.36
K0.6Mg0.25Al2.3Si3.5O10(OH)2 + 8 H+ = 0.6 K+ +0.25 Mg2+ + 2.3 Al+3 + 3.5 SiO2 + 5 H2O
8.51
Quartz
SiO2 = SiO2
-3.91
CalciteCaCO3 = Ca+2 + CO32-8.51
Dolomite
CaMg(CO3)2 = Ca+2 + Mg+2 + 2 CO32-16.71
GypsumCaSO4.2H2O = Ca+2 + SO4-2 + 2 H2O-4.48
Secondary minerals
Hydrotalci
Mg4Al2O4(OH)6 + 14 H+ = 4 Mg2+ + 2 Al+3 + 10
75.44
te
H2O
Mg4Si6O15(OH)2.6H2O + 8 H+ = 4 Mg2+ + 6 SiO2
29.96
Sepiolite+ 11 H2O(1)
Used as input in CRUNCH.
(2)
Used as input in HP1.
Kaolinite
IlliteMolar
volumeVolume%(1)

Moles(2)99.52023

2.31

22.68036.93464.365
74.69

19
11.5
13.5
2.0
0.02

0.38
5.06
3.655
0.31
0.00268207.570.00

0.0

285.600.00

0.00

a lower reactive surface area during dissolution, Mm,0 and Mm are the initial and current
numbers of moles (in HP1) and volume percentages (in CRUNCH) of mineral m, n is thereaction order of the OH- dependency of the reaction rate, k0,m is the intrinsic rate constant ofmineral m (mol m-2 sec-1), Km is the equilibrium constant of the reaction, and Qm is the ion
activity product. Table 4.10 lists the parameters of Eq. (4.7). Due to the stiffness of thesystem (a combination of high and low rate constants), not all dissolution processes of the
primary minerals were treated as kinetic reactions in HP1 (equilibrium was assumed for
calcite). Eq. (4.7) was also used to describe the kinetic precipitation of two secondary
minerals (see Table 4.9) in CRUNCH. These minerals were treated as being in equilibriumwith the solution in HP1. Quartz, kaolinite, and illite were not allowed to precipitate in either
model.

Table 4.10 Parameters for the kinetic dissolution reactions (Eq. (4.7)).
MineralLog k0nA0(mol m-2 sec-1)(m2 m-3)Kaolinite1-10.880.54
1.00 1071
Illite-13.910.22
0.33 107
1
Quartz-10.200.30
2.00 1052
Calcite-1.00
0.00
2.00 1031
Dolomite-7.70
0.00
2.00 1031
Gypsum-8.00
0.00
2.00 1032
Hydrotalcite
-1.00
0.00
1.00 102
2
Sepiolite
-1.00
0.00
1.00 102
1
Data from Table 5, Part III, Adler (2001)
2
Minerals are assumed to be in equilibrium in HP1. High dissolution
constants are used in CRUNCH.

52
Table 4.11 Overview of aqueous equilibrium reactions, and cation exchange half reactions with corresponding
equilibrium constant.
Nr
(1)
(2)
(3)
(4)
(5)
(6)
(7)
(8)
(9)
(10)
(11)
(12)
(13)
(14)
(15)
(16)
(17)
(18)
(19)
(20)
(21)
(22)
(23)
(24)
(25)
(26)
(27)
(28)
(29)
(30)
(31)
(32)
(33)
(34)
(35)
(36)

Reaction
Aqueous speciation reactions
H2O = OH-+ H+
HCO3- + H+ = CO2 + H2O
HCO3- = CO3-2 + H+
3+
Al + H2O = AlOH+2 + H+
Al3+ + 2 H2O = Al(OH)2 + 2 H+
Al3+ + 2 H2O =AlO2- + 4 H+
Al3+ + 2H2O = HAlO2 + 3 H+
Al+3 + SO42- = AlSO4++3
Al + Na+ + 2 H2O = NaAlO2 + 4 H+
Al+3 + 2 SO42- = Al(SO4)2Ca2+ + H2O = CaOH+ + H+
Ca2+ + HCO3- = CaCO3 + H+
Ca2+ + HCO3- = CaHCO3+
Ca2+ + Cl- = CaCl+
Ca2+ + 2 Cl- = CaCl2
Ca2+ = SO42- = CaSO4
Cl- + H+ = HCl+
K + H2O = KOH + H+
K+ + Cl- = KCl
K+ + SO42- = KSO4+
K + H+ + SO42- = KHSO4
Mg2+ + H2O = MgOH+ + H+
Na+ + H2O = NaOH + H+
Na+ + HCO3- = NaHCO3
Na+ + HCO3- = NaCO3- + H+
Na+ + Cl- NaCl
Na+ + SO42- = NaSO4SO42- + H+ = HSO4SO42- + 2 H+ = H2SO4
SiO2 + 2 H2O = H2SiO42- + 2 H+
SiO2 + H2O = HSiO3- + H+
SiO2 + Na+ + H2O + NaHSiO3 + H+
Cation exchange half reactions (X- is the exchange site)
Na+ + X- = NaX
K+ + X- = KX
Ca2+ + 2 X- = CaX2
Mg2+ + 2 X- = MgX2

Log K
-14.0
6.3447-10.3288
4.9571-10.5945-22.8833-16.4329
3.01-23.6266
4.9-12.78-7.00171.0467-0.6956-0.6436
2.111
0.67-14.46-1.4946
0.87960.8136-11.44-14.180.1541-9.8144-0.7770.82
1.9794-1.0209-22.96-9.9525-8.3040
0.00
0.757
0.604
0.505

An overview of the considered aqueous equilibrium reactions and the ion exchange reactions
is given in Table 4.11. The total cation exchange capacity of Opalinus Clay is 120 meq per kgof rock. The initial concentrations and concentrations in the infiltrating water are given in
Table 4.12.
4.3.2.2 Database and model input
HP1
The chemical reactions and their equilibrium constants are stored in the databasephreeqc.dat.
This
problemuses
thefollowing
keywords:
SOLUTION

53
Table 4.12 Chemical composition of the initial pore water andthe inflowing alkaline solution.
Element

Initial composition
Composition of theof the pore water
inflowing alkaline solution(mol kgw-1)1(mol kgw-1)2-3
K
5.42 10
1.61 10-1-1Na
3.28 10
6.52 10-2-2Ca
2.16 102.24 10-3-2
Mg
3.57 101.00 10-30-1Cl
2.64 101.00 10-30-2S
9.20 101.00 10-30-6Al
1.00 101.00 10-30-5Si
4.40 101.00 10-30-4C
3.96 10
4.21 10-5
pH
7.90
13.20
1
: From Table 2; Part III, Adler (2001)2: From Table 6, Part III, Adler (2001)_MASTER_SPECIES,
SOLUTION_SPECIES,
EXCHANGE_MASTER_SPECIES,
EXCHANGE_SPECIES, PHASES, and RATES. The last keyword is used to define thespecific rate reactions for the dissolution of primary minerals (except calcite) with Eq. (4.7)
translated into the BASIC-programming language for each mineral (for details, see the
PHREEQC-2 manual; Parkhurst and Appelo, 1999).
The input file for water flow and solute transport can again be easily created with thegraphical interface of HYDRUS-1D. Nine components (i.e., Ca, K, Na, C, S, Cl, Mg, Si, andAl) are transported in this problem. Only parameters related to transport are defined in
HYDRUS-1D, i.e., the dispersivity, and the aqueous diffusion coefficient. This information isstored in input files selector.in and profile.dat. Note that these files also contain information
for the steady-state inflow flux and concentrations of the nine components in the infiltratingwater.
The link between HYDRUS-1D and PHREEQC is defined in the input file species.in, whichcontains names of the 9 components. The phreeqc.in file contains chemical information for
the reactive transport problem. Initial mineralogical and chemical conditions for each cell aredefined with the following keywords: SOLUTION (initial chemical composition of theaqueous phase), EXCHANGE (size and initial composition of the cation exchange complex),
EQUILIBRIUM_PHASES (identifies the minerals in equilibrium with the aqeous phase andtheir initial amount), KINETICS (identifies the kinetic reactions with their parameters andtheir initial amount), TRANSPORT (to indicate that HYDRUS-1D is used for transportmodelling), and SELECTED_OUTPUT (to specify the desired output). Details on the
keywords used in phreeqc.in can be found in the PHREEQC-2 manual (Parkhurst andAppelo, 1999).

54

CRUNCHThe database contains chemical reactions and equilibrium constants for the aqueous
speciation reactions (between end of primary and end of secondary), the minerals(between end of secondary and end of minerals), the kinetic dissolution or
precipitation reactions of the minerals (between begin of mineral kinetics andend of mineral kinetics), and the exchange reactions (between begin of
exchange and end of exchange). The input file reactive.in defines the initialconditions for the clay and the inflowing water (CONDITION), the primary, secondary,
and/or exchange species, and the minerals (PRIMARY_SPECIES, SECONDARY_SPECIES,
ION_EXCHANGE, and MINERALS), the discretization, and the initial and boundary
INITIAL_CONDITIONS,andBOUNDARY
conditions(DISCRETIZATION,
CONDITIONS), the flow and transport properties (FLOW, and TRANSPORT), and the output(OUTPUT). The keyword POROSITY defines the porosity which is fixed in this simulation.
Numerical information for solving the equations are defined in RUNTIME.

4.3.2.3 Comparison between CRUNCH-GIMRT, CRUNCH-OS, and HP1
Simulations were performed with two different codes (HP1 and CRUNCH) where CRUNCH
was applied in two different modes: the global implicit option in the CRUNCH model(CRUNCH-GIMRT), and the sequential non-iterative option in the CRUNCH model(CRUNCH-OS). The coupled HP1 model runs in a sequential non-iterative mode. For thelatter model, different maximum time steps (corresponding to different Courant numbers:
0.125 and 0.065) were used. The maximum allowed Courant number of the CRUNCH-OS
simulation was set to 0.125.

Figure 4.11 shows distributions of pH and selected elements after 0.3, 0.7, and 1.1 y. Outflow
curves are shown in Figure 4.12. In general, correspondence between the different modelling
results is relatively good. Very good agreement was obtained for the propagation of both highpH and elements that are only influenced by cation exchange (Na and K, the latter not shown)(Figure 4.11a, b and Figure 4.12b). The pH near the inlet, as obtained with HP1 was slightly
larger than the pH of the CRUNCH simulations. The positions of the propagation waves of Si,
Al, and Mg (Figure 4.11c-e) and Ca (similar to Mg, not shown) were identical for thedifferent models, although a slight difference occurred in the concentration of these elements.
The reason for this is not clear; it might be due to difference in the precipitation kinetics of thesecondary minerals (hydrotalcite and sepiolite). Precipitation of these minerals in HP1 isgoverned by thermodynamic equilibrium, whereas kinetic precipitation is assumed in the
CRUNCH simulations. Due to this equilibrium assumption, Mg concentrations as calculatedwith HP1 should be slightly smaller than those obtained with CRUNCH. The Si or Al
concentrations are consequently also slightly larger. Other minor differences between the twomodels (e.g., the exact activity coefficients) may have contributed also to this difference.

55

Compared to cation exchange problems (see section 4.3.1), problems involving (kinetic)
dissolution/precipitation of minerals generally require somewhat smaller time steps if asequential non-iterative coupling method is used. Nevertheless, acceptable results wereobtained compared to those using the global implicit method.

56
0.25

13
12

1.1 y

0.2

11
10

Na

pH

0.15

0.7 y

0.1

0.3 t

9

0.05

8
7

(a)
0

0
0.02

0.04

0.06

0.08

(b)
0

0.04

0.016

0.03

0.012

0.02

0.01

0

0.02

0.04

0.06

0.08

Depth in clay core (m)AlSi

Depth in clay core (m)

0.008

0.004

(c)
0

0
0.02

0.04

0.06

0.08

(d)
0

Depth in clay core (m)

0.02

0.04

0.06

0.08

Depth in clay core (m)5x10-48x10-34x10-46x10-3S

Mg

3x10-42x10

4x10-3-4

2x10

1x10-4
0x100

(e)
0

-3

0x100
0.02

0.04

0.06

0.08

Depth in clay core (m)CRUNCH-GIMRTCRUNCH-OS

(f)
0

0.02

0.04

0.06

0.08

Depth in clay core (m)HP1 (Cr 0.125)
HP1 (Cr 0.065)

Figure 4.11 Distribution of pH and selected elements versus depth (mol/l) after 0.3, 0.7, and 1.1 y of infiltrationof a high pH solution in an Opalinus Clay core.

57
0.4

1x10-1

(a)

(b)

-2

1x10

0.3

Ca

Na

1x10-3
0.2

-4

1x10
0.1

1x10-5
1x10-6

0
0

0.2

0.4

0.6

0.8

1

0

0.2

Time (y)0.0008

0.4

0.6

0.8

1

0.8

1

Time (y)1x10-4(c)(d)

0.0006AlSi

8x10-50.00044x10-50.0002

0x100

0
0

0.2

0.4

0.6

0.8

1

0

0.2

Time (y)0.04

0.6Time (y)
11

(e)

(f)

10

pH

0.03

Mg

0.4

0.02

0.01

9

8

0

7
0

0.2

0.4

0.6

0.8

1

Time (y)

0

0.2

0.4

0.6

0.8

1

Time (y)CRUNCH-GIMRTCRUNC-OS

HP1 (Cr 0.125)
HP1 (Cr 0.065)

Figure 4.12 Outflow curves of selected elements (mol/l) during infiltration of a high pH-solution in an Opalinus
Clay core.

58

59

5 Example
This section considers an example that combines transient variably-saturated water flow,
solute transport, and geochemical reactions. The example is used to illustrate some of thecapabilities of HP1.

5.1 Long term transient flow and transport of major cations and heavy metals in a soil
profile (HEAVMET)

5.1.1

Problem definition and governing chemical reactions

This example considers similar processes as those discussed in Verification Problem 8. The
problem involves the transport of major cations (Na, K, Ca, and Mg) and heavy metals (Cd,Zn, and Pb) in a 1-m deep multi-layered soil profile subject to atmospheric boundary
conditions. While steady-state water flow was considered in Verification Problem 8, in thisexample the upper boundary condition varies with time and consists of daily values ofprecipitation and evaporation (no plants are assumed to be present, and transpiration and root
water uptake are not considered).
Soil hydraulic and physical parameters (Table 5.1) of the dry Spodosol located at the“Kattenbos” site near Lommel, Belgium were taken from Seuntjens (2000, Tables 3.1 and
7.1). The cation exchange complex was assumed to be associated solely with organic matter.
The cation exchange capacity hence is directly related to the amount of exchangeable protons
on the organic matter, taken to be 6 meq/g of the organic matter (proton dissociating groups
on fulvic acids are 6 – 10 meq/g and 4 – 6 meq/g on humic acids, Tipping, 2002). The cationexchange complex was assumed to consist of six different exchange sites as defined in Table
4.8.
The upper boundary condition (precipitation and evaporation) was defined using
meteorological data from the Brogel station (Belgium) for a period of 10 years (1972-1981).
Table 5.1 Soil hydraulic and other properties of six soil horizons.
Horizon

Depth
(cm)

ρ
(g cm³)

A
E
Bh1
Bh2
BC
C1
C2

0–7
7 – 1919 – 24
24 – 28
28 – 50
50 – 75
75 – 100

1.31
1.59
1.3
1.38
1.41
1.52
1.56

Organic
Carbon
%
2.75
0.75
4.92
3.77
0.89
0.12
0.08

θr

θs

α
(m-1)

n

Ks
(ms-1)

0.065
0.035
0.042
0.044
0.039
0.03
0.021

0.48
0.42
0.47
0.46
0.46
0.42
0.39

1.6
1.5
1.6
2.8
2.3
2.1
2.1

1.94
3.21
1.52
2.01
2.99
3.72
4.33

1.1 10-5
3.6 10-5
4.5 10-6
1.0 10-4
1.0 10-4
1.0 10-4
1.0 10-4

60

Cumulative P, E, P-E (cm)

800

600

400

200

0
1972

1974

1976

1978

1980

1982

Time (year)
Precipitation P
Evaporation E
Potential surface flux (P-E)
Actual surface flux
Figure 5.1 Cumulative precipitation, evaporation andsurface flux into the soil between 1972 and 1982.

Figure 5.1 shows precipitation P, potential evaporation E, and potential flow across the upperboundary (P – E). Since the amount of water available in the soil profile for evaporation was
smaller than that required for potential evaporation, the actual evaporation rate will be smaller
than the potential evaporation rate. Consequently, the actual surface flux was larger than thepotential surface flux (see Figure 5.1). For the considered 10-year period, the actual surface
flux (around 45 cm/year) was somewhat larger than expected for this region of Belgium. Thisis because no water uptake by plant roots was considered in this example. Element
concentrations in the precipitation (incoming water) were obtained from Stolk (2001) for
station 231 located in Gilze-Rijen (the Netherlands) close to the Belgian border. The
composition was based on the average of 13 measurements during 1999 (Table A5 in Stolk,
2001) (Table 5.2).

Table 5.2 Chemical composition of the rain water (from Stolk, 2001) and in the initial soil solution.
ConcentrationpHCl(µmol/l)Rain water
5.25
69
Initial soil solution
A
3.4
69
E
3.5
69
Bh1
3.6
69
Bh2
3.8
69
BC
4.4
69
C1
4.4
69
C2
4.5
69
*: Br is used as charge balance

Br*

Na

K

Ca

Mg

Cd

Zn

Pb

32

64

4

6

8

0

0

0

727
521
398
287
235
252
131

64
64
64
64
64
64
64

4
4
4
4
4
4
4

97
65
39
33
55
76
27

8
8
8
8
8
8
8

0.80
0.38
0.43
0.41
0.63
0.33
0.20

50
24
23
21
33
21
14

2.5
1.2
1.1
1.0
1.6
1.0
0.7

61

The initial water content was estimated by assuming a constant flux of 0.1206 cm/day(corresponding to the long term actual infiltration rate). Initial concentrations of Na, K, Mg
were assumed to be equal to the concentration in the precipitation, while initial concentrationsof Ca, Cd, and Zn were obtained from Seuntjens (2000; Table 5.2). Pb concentrations werearbitrarily set to be 20 times smaller than the Zn concentrations (no data available).
The simulation was carried out using the following numerical guidelines: a spatial
discretization of 1 cm, a performance index of 0.1, and a maximum time step of 0.125.

5.1.2

Analysis of the simulation results

5.1.2.1 Accuracy analysis
To assess the accuracy of the simulation, the mass balance was evaluated in terms of thefollowing two variables as a function of time:ttj
Qinj − Qout
= ⎛⎜ ∫ q sj0 dt − ∫ q sNj dt ⎞⎟
0
⎝ 0
⎠

[

(5.1)

]

NX
NX
⎡⎛
⎞
⎛
⎞ ⎤
m
m
M t − M = ∑ ⎢⎜⎜ (C j )i ,t + ∑ A j ( X j )ν j ⎟⎟θ i ,t − ⎜⎜ (C j )i ,0 + ∑ A mj ( X mj )ν j ⎟⎟θ i ,0 ⎥ ∆x
k =1
i =1 ⎣⎝
k =1
⎠
⎝
⎠ ⎦
j

j
0

N cells

[

]

(5.2)

where j = 1,…, Nm, Qijn and Qojut are the cumulative inflow and outflow rates of the jth masterspecies [moles L-2], respectively; qsj0 and qsjN are solute flux boundary conditions at thebottom and surface of the soil profile [moles T-1 L-2], respectively; Mtj, M0j are total amounts of
the jth master species in the flow region at time t and at the beginning of the simulation[moles L-2], respectively, and νj is the charge of the jth master species. All other symbols were
defined in section 2. Mtj and M0j in Eq. (5.2) consist of the total amount of the jth masterspecies in the aqueous phase, θCj (see Eq. 2.28), and the sum of the amounts adsorbed on six
cation exchange sites. Since an accurate simulation should conserve mass, the difference
between Eqs. (5.1) and (5.2), defined as the absolute error εa (Šimůnek et al., 1998), should besmall. Qijn and Qojut were calculated using information from the soluteN.out output files
generated by HYDRUS1D, while Mtj - M0j were calculated using values from the output filefrom PHREEQC defined by the keyword SELECTED_OUTPUT every 15 days. Figure 5.2
compares Eqs. (5.1) and (5.2) for Cl, Na, Ca and Cd for the 10-year simulation. Since
differences between values obtained with Eqs. (5.1) and (5.2) are very small, absolute errors
are small as well and no mass balance errors are expected.
The accuracy of the numerical solution can be evaluated also by using the relative massbalance error, εr (Eq. 6.52 in Šimůnek et al., 1998):

62

1.2x10-3

1.5x10-3
1.2x10-3∆Na (mol m-2)

∆Cl (mol m-2)

8.0x10-4

4.0x10-4

0.0x109.0x10-4

6.0x10-4
3.0x10-4

0

0.0x10-4.0x10-4
1972

1974

1976

1978

1980

1982

0

-3.0x10-4
1972

1974

Time (year)

1976

1978

1980

1982

1980

1982

Time (year)

2x10-40x100

Qin - Qout
Mt - M00∆Cd (mol m-2)

∆Ca (mol m-2)

0x10

-2x10-4
-4x10

-4

-1x10

-5

-2x10-5-6x10-4-8x10-4
1972

1974

1976

1978

1980

1982

-3x10-5
1972

1974

Time (year)

1976

1978

Time (year)

Figure 5.2 Time series of inflow minus ouftlow (dashed line) and change in the total amount present in the soil
profile (dots) for Cl, Na, Ca, and Cd.

ε rj =

100 ε a⎛ N cells⎜⎜⎜⎜ i =1⎝max ∑ M

t

j
t ,i

−M

j
t ,i

(

, M +∫ q
j
0

0

j
s0

+q

j
sN

)dt

⎞
⎟
⎟
⎟
⎟
⎠

(5.3)

where Mtj,i and M0j,i are the total amount of the jth master species in the ith cell at time t and at
the beginning of the simulation, respectively. Figure 5.3 shows relative mass balance errors
for Na, Ca, and Cd. The errors were consistently smaller than 1%, indicating that massbalance errors were small in this HP1 simulation using the non-iterative sequential approachfor solving the coupled transport-reaction equations.
5.1.2.2 Simulation results
As an illustration, two-dimensional time-depth plots are shown in Figure 5.4 for the top 7 cm
(the A horizon) for the water content, pH, mass of Cl (mmol / 1000 cm³ soil), and mass ofaqueous Cd (mmol / 1000 cm³ soil). The alternation between precipitation (wet conditions)and evaporation (dry conditions) as dictated by the atmospheric boundary conditions clearlyaffected the dynamics of the water content (Figure 5.4a), including upward flow of water
during the dry periods (not shown). The flow dynamics in turn significantly influenced the

63
0.6

0.2

0.4

ε Ca
r

ε Na
r

0.3

0.1

0.2

0
1972

0
1974

1976

1978

1980

1982

Time (year)

1972

1974

1976

1978

1980

1982

Time (year)

0.5
0.4

ε Cd
r

0.3
0.2
0.1
0
1972

1974

1976

1978

1980

1982

Time (year)
Figure 5.3 Relative mass balance errors, εr, as a function of time for Na, Ca, and Cd.

geochemistry near the soil surface. The most mobile elements, such as Cl and the monovalentcations (H+, Na+) moved upwards during the evaporation periods, leading to higher
concentrations near the soil surface (e.g., Cl in Figure 5.4b). The increased proton
concentration (induced by the lower water content and the supply from deeper horizons in thesoil profile) produced lower pH values near the soil surface during the summer periods with
their lower water contents (Figure 5.4c). The pH in actuality may have been affected also by
soil carbon dioxide concentrations, which usually change in response to the biological activity
and moisture status of the soil (Šimůnek and Suarez, 1993). This process, however, was notconsidered in our simulations here. Although upward flow during the summer had almost no
effect on the total amount of heavy metals near the soil surface (results not shown) due to thelow mobility of these elements, the concentrations of these elements did vary significantly
during the various seasons.
Because of lower water contents, concentrations of all aqueous species and elements
increased significantly during the summer periods. The changing aqueous concentrations inturn also caused changes in the cation exchange equilibrium by promoting monovalent cations
to sorb on the cation exchange complex. The lower cation exchange capacity (defined here as
the amount of deprotonated sites) during summer (lower pH) implied the presence of fewer

64(a) Water content(b) Log(Cl) (mmol / 1000 cm³ soil)

0

0

-3.8

-2

0.22
0.18

-4

0.14
0.1

-6
1972

-4

Depth (cm)Depth (cm)

0.26

1974

1976

1978

1980

1982

-2

-4.2
-4.4

-4

-4.6
-4.8

-6

-5
-5.2

0.06

1972

1974

1976

1978

1980

1982

Time (year)

Time (year)(c) pH(d) Log(Aqueous Cd) (mmol / 1000 cm³ soil)

0

0
-7

4.3

-2

4.1
3.9

-4

3.7

Depth (cm)Depth (cm)

4.5

-7.4

-2

-7.8
-8.2

-4

-8.6

3.5
-6
1972

3.3
1974

1976

1978

Time (year)

1980

1982

3.1

-6

-9
-9.4

1972

1974

1976

1978

1980

1982

Time (year)

Figure 5.4 Space-time plots of (a) water content, (b) amount of Cl (mmol / 1000 cm³ soil), (c) pH, and (d)amount of aqueous Cd (mmol / 1000 cm³ soil) in the A-horizon (0 – 7 cm depth).

sorption sites for the other cations. The increased supply of monovalent cations due to upward
water flow during summer further increased competition for the sorption sites, leading torelatively more sorption of the monovalent cations and higher concentrations of divalentcations and heavy metals in the soil solution. Another important process is the complexationof the heavy metals with Cl. Because of the presence of relative stable aqueous complexes(e.g., CdCl+), the amounts of Cd, Pb, and Zn will increase when Cl concentrations are high(i.e., during summer periods). Figure 5.4d shows the amount of Cd in the aqueous phase.
Differences in the amount of aqueous Cd between the summer and winter periods are morethan two orders of magnitude. Since water contents are lower during summer, the differencesin the concentrations will be even larger. Similar, but less distinct, patterns can be observed
also deeper in the soil profile (Figure 5.5).
This example shows that atmospheric boundary conditions can have a significant effect on theamounts of Cd and other heavy metals in the solution, and hence also on their bioavailability,
since uptake processes by plants and soil micro-organisms are often concentration dependent(i.e., passive uptake of solutes together with uptake of water by roots, as well as active uptake
which can be simulated with HP1 using Monod or Michealis-Menten kinetics). Moreover, thehigh concentrations of the heavy metals occur during periods with the highest

65

Water content

0.4
0.36
0.32
0.28

Cl (mmol / 1000 cm³ soil)

0.24
0.2

1x10-1

1x10-2

pH

4.4
4.2

Cd (mmol / 1000 cm3 soil)

4
3.8

1x10

-5

1x10-6
1972

1974

1976

1978

1980

1982

Time (year)
Figure 5.5 Time series of water content, Cl concentration (mmol / 1000 cm³ soil), pH, andCd concentration (mmol / 1000 cm³ soil) in the third horizon (Bh1-horizon) at 22 cm depth.

(micro)biological activity during the year (the summer months). HP1 hence seems to be avery attractive tool for evaluating the effects of atmospheric boundary conditions on the longterm mobility and leaching of the heavy metals in and from soils. While not further addressed
here, such evaluations possibly could also account for the transport of colloids such as microorganisms, humics or fulvics.

66

67

6 Concluding Remarks
The HP1-code, presented in this report, couples HYDRUS-1D (a one-dimensional variablysaturated water flow and solute transport model) with PHREEQC (a general biogeochemicalcode). The new code is able to handle (1) transient water flow in variably-saturated media, (2)transport of multiple components, (3) heat flow, and (4) mixed equilibrium (aqueous
speciation, mineral dissolution / precipitation, cation exchange) and kinetic geochemical
reactions.
The non-iterative sequential approach used to solve the coupled transport-reaction equations
generally worked very well. The accuracy of the coupling approach was tested using nine
verification examples with different degrees of complexity. Overall, compared to othermodels / codes, HP1 was found to be accurate when time steps were not too large. The easiest
and most efficient way to control the time step is by defining a relatively small maximum
time step and/or a small performance index (both in HYDRUS-1D). Our analyses indicatethat the performance index should be about five time smaller than the value previously
suggested for HYDRUS-1D simulations (i.e., 0.4 rather than 2).
The nine verification problems and the application example serve as guidelines to create theinput files using the graphical user interfaces of HYDRUS-1D and PHREEQC. The userinterfaces also display the output created with the HP1 simulations.

68

69

ACKNOWLEDGEMENTS
The work of Dr. J. Šimůnek was 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 andModel Integration). Additional support was obtained through the bilateral agreement project“Development and evaluation of a coupled geochemical transport model” (Agreement No. 585310-0-F105) between SCK•CEN – USDA, ARS and (Agreement No. C0-90001412.01)
between SCK•CEN – University of California, Riverside.
The work of Dr. D. Jacques was supported by SCK•CEN R&D program on “Environmental
Remediation” (project E022031 and CO91002).
The project coordinator at SCK•CEN was Dr. D. Mallants.

70

71

References
Aagaard, P., and H. C. Helgeson, 1982. Thermodynamic and kinetic constraints on reaction
rates among minerals and aqueous solutions, I. Theoretical considerations. Amer. J. Sci.,
282:237-285.
Adler, M., 2001. Interaction of claystone and hyperalkaline solutions at 30°C: A combined
experimental and modeling study. Ph. Dissertation, Universität Bern, 120 p.
Appelo, C.A.J., and D. Postma, 1993. Geochemistry, groundwater and pollution. A.A.
Balkema, Rotterdam, 536 p.
Appelo, C.A.J., E. Verweij, and H. Schäfer, 1998. A hydrogeochemical transport model for anoxidation experiment with pyrite/calcite/exchangers/organic matter containing sand.
Applied geochemistry, 13:257-268.
Appelo, C.A.J., M.J.J. Van der Weiden, C. Tournassat, and L. Charlet, 2002. Surface
complexation of ferrous iron and carbonate on ferrihydrite, and the mobilisation of arsenic.
Environ. Sci. Technol. 36:3096-3103.
Brooks, R.H., and A.T. Corey, 1966. Properties of porous media affecting fluid flow. J. Irrig.
Drainage Div., ASCE Proc. 72(IR2), 61-88.
Carsel, R.F., and R.S. Parish, 1988. Developing joint probability distributions of soil water
retention characteristics. Water Res. Res., 24:755-769.
Casey, F., and J. Šimůnek, 2001. Inverse analyses of transport of chlorinated hydrocarbons
subject to sequential transformation reactions. J. Envir. Qual., 30:1354-1360.
Casey, F. X. M., G. L. Larsen, H. Hakk, and J Šimůnek, 2003. Fate and transport of 17βEstradiol in soil-water systems. Environmental Science and Technology, 37(11):24002409.
Casey, F.X.M., G. L. Larsen, H. Hakk, and J. Šimůnek, 2004. Fate and transport of
testosterone in agriculturally significant soils. Environ. Sci. Technol., 38(3):790-798.
Chung, S.-O., and R. Horton, 1987. Soil heat and water flow with a partial surface mulch.
Water Resour. Res., 23:2175-2186.
Dane, J.H. and G.C. Topp, 2002. Methods of soil analysis, Part 4 Physical systems. SSSA
Book Series no. 5, SSSA, Inc. Madison, Wisconsin, USA, 1692 pp.
de Marsily, G., 1986. Quantitative Hydrogeology. Academic Press, London.
de Vries, D.A.. The thermal properties of soils. In ‘Physics of plant environment’ (Ed. R.W.
van Wijk), North Holland, Amsterdam, pp. 210-235.
Gaines, G.L., and H.C. Thomas, 1953. Adsorption studies on clay minerals. II. A formulationof the thermodynamics of exchange adsorption. J. Chem. Phys., 21:714-718.
Hopmans, J.W., J. Šimůnek, N. Romano, and W. Durner, 2002. Inverse methods. In ‘Methodsof soil analysis, Part 4’, (Co-eds. Dane and Topp), Book Series no. 5, SSSA, Inc. Madison,
Wisconsin, USA, 693-1008.
Jacques, D., J. Šimůnek, D. Mallants, M. Th. van Genuchten, 2002. Multicomponent transportmodel for variably-saturated porous media: Application to the transport of heavy metals in72soils, In: Computational Methods in Water Resources (Eds. Hassanizadeh, S.M., R.J.
Schotting, W.G. Gray, and G.F. Pinder), XIVth International Conference on
Computational Methods in Water Resources, June 23-28, 2002, Delft, The Netherlands.
Elsevier, Developments in Water Science, 47, 555-562.
Jacques, D., J. Šimůnek, D. Mallants, and M. Th. van Genuchten, 2003. The HYDRUSPHREEQC multicomponent transport model for variably-saturated porous media: Codeverification and application, MODFLOW and More 2003: Understanding throughModeling, Conference Proceedings, ed. E. Poeter, Ch. Zheng, M. Hill, and J. Doherty,International Ground Water Modeling Center, Colorado School of Mines, 23-27.
Jaynes, D.B., A.S. Rogowski, and H.B. Pionke, 1984. Acid mine drainage from reclaimed
coal strip mines. 1. Model description. Water Resources Research, 20: 233-242.
Källvenius, G., and C. Ekberg, 2003. TACK – a program coupling chemical kinetics with a
two-dimensional transport model in geochemical systems. Comp. and Geosci., 29:511-521.
Kirkham, D., and W.L. Powers, 1972. Advanced Soil Physics. John Wiley&Sons, New York,NY.
Lichtner, P.C., 1996. Continuum formulation of multicomponent-multiphase reactive
transport. In ‘Reactive transport in porous media’ (Eds. P.C. Lichtner, C.I. Steefel and
E.H. Oelkers), Reviews in Mineralogy, no. 34, Mineralogical Society of America,
Washington, DC, p. 1-81.
Mallants, D., D. Jacques, and T. Zeevaert, 2003. Modelling 226Ra, 222Rn, and 210Pb in a
proposed surface repository of very low-level long-lived radioactive waste. ICEM’03
Conference, Oxford, 21-25 Sept. 2003. ASME, Proc. On CD-ROM, session 25, paper 6.
Mayer, U.,1999. A numerical model for multicomponent reactive transport in variably
saturated porous media. Ph.D. thesis, Department of Earth Sciences, University ofWaterloo, 285p.
Millington, R.J., and J.M. Quirk, 1961. Permeability of porous solids. Trans. Faraday Soc.,
57:1200-1207.
Morel, F., and J. Hering, 1993. Principles and applications of aquatic chemistry. New York,John Wiley & Sons Inc.
Parkhurst, D.L., and C.A.J. Appelo, 1999. User's guide to PHREEQC (Version 2) – Acomputer program for speciation, batch-reaction, one-dimensional transport, and inverse
geochemical calculations. Water-Resources Investigations, Reprot 99-4259, Denver, Co,
USA, 312 pp.
Perrochet, P., and D. Berod, 1993. Stability and the standard Crank-Nicolson-Galerkinscheme applied to the diffusion-convection equation: some new insights. Water Resour.
Res., 29:3291-3297.
Postma, D. and C.A.J. Appelo, 2000. Reduction of Mn-oxides by ferrous iron in a flow
system: column experiments and reactive transport modeling. Geochim. Cosmochim. Acta,
64:1237-1247.

73

Prommer, H., 2002. A reactive multicomponent transport model for saturated porous media.
User's Manual Version 1.0. Contaminated Land Assessment and Remediation ResearchCentre, University of Edinburgh, UK, 132 p.
Schaerlaekens, J., D. Mallants, J. Šimůnek, M.Th. van Genuchten, and J. Feyen, 1999.
Numerical simulation of transport and sequential biodegradation of chlorinated aliphatichydrocarbons using CHAIN_2D. Hydrol. Process., 13:2847-2859.
Šimůnek, J., and D. L. Suarez, 1993. Modeling of carbon dioxide transport and production insoil: 1. Model development. Water Resour. Res., 29: 487-497.
Šimůnek, J., and D.L. Suarez, 1994. Two-dimensional transport model for variably saturated
porous media with major ion chemistry. Water Resour. Res., 30:1115-1133.
Šimůnek, J., and D. L. Suarez, 1997. Sodic soil reclamation using multicomponent transportmodelling. ASCE J. Irrig. and Drain. Engineering, 123(5):367-376.
Šimůnek, J. and A. J. Valocchi, 2002. Geochemical Transport, In: Methods of Soil Analysis, Part
1, Physical Methods, Chapter 6.9, Eds. J. H. Dane and G. C. Topp, Third edition, SSSA,
Madison, WI, 1511-1536.
Šimůnek, J., M. Šejna, and M.Th. van Genuchten, 1998. The HYDRUS-1D software package
for simulating the one-dimensional movement of water, heat, and multiple solutes invariably-saturated media. Version 2.0, IGWMC-TPS-70, International Ground Water
Modeling Center, Colorado School of Mines, Golden, CO, 202 p.
Šimůnek, J., D. L. Suarez, and M. Šejna, 1996. The UNSATCHEM software package forsimulating one-dimensional variably saturated water flow, heat transport, carbon dioxide
production and transport, and multicomponent solute transport with major ion equilibriumand kinetic chemistry, Version 2.0, Research Report No. 141, U.S. Salinity Laboratory,
USDA, ARS, Riverside, California, 186pp.
Šimůnek, J., T. Vogel, and M.Th. van Genuchten, 1992. The SWMS-2D code for simulatingwater flow and solute transport in two-dimensional variably saturated media, Version 1.1,
Res. Rep. 126, U.S. Salinity Lab., Agric. Res. Serv., U.S. Dept. of Agric., Riverside, Ca.
Singer, P.C., and W. Stumm, 1970. Acid mine drainage: the rate-determining step. Science,
167:1121-1123.
Steefel, C.I., 2000. New directions in hydrogeochemical transport modeling: Incorporating
multiple kinetic and equilibrium reaction pathways. In 'Computational Methods in Water
Resources XIII' (Eds. L.R. Bentley, J.F. Sykes, C.A. Brebbia, W.G. Gray, and G.F.
Pinder), A.A. Balkema, Rotterdam, 331-338.
Steefel, C.I., and K.T.B. MacQuarrie, 1996. Approaches to modeling of reactive transport inporous media. In ‘Reactive transport in porous media’ (Eds. P.C. Lichtner, C.I. Steefel and
E.H. Oelkers), Reviews in Mineralogy, no. 34, Mineralogical Society of America,
Washington, DC, p. 83-129.
Steefel, C.I., and S.B. Yabusaki, 1996. OS3D/GIMRT, Software for multicomponentmultidimensional reactive transport, user manual and programmer's guide. PNL-11166,
Pacific North west Laboratory, Richland, Washington.

74

Stolk, A.P., 2001. Landellijk Meetnet Regenwatersamenstelling, Meetresultaten 1999. RIVM,
The Netherlands, RIVM Rapport 723101 056, 61 p.
Suarez, D.L., and J. Šimůnek, 1997. UNSATCHEM: Unsaturated water and solute transportmodel with equilibrium and kinetic chemistry. Soil Sci. Soc. Am. J., 61:1633-1646.
Tipping, E., 2002. Cation binding by humic substances. Cambridge University Press,
Cambridge, UK, 434 p.
Truesdell, A.H., and B.F. Jones, 1974. WATEQ, A computer program for calculatingchemical equilibria of natural waters. J. Res. U.S. Geol. Surv., 2:233-274.
Thury, M., and Bossart, P., 1999. The Mont Terri rock laboratory, a new international
research project in a Mesozoic shale formation, in Switzerland. Eng. Geol., 52:347-359.
van Genuchten, M.Th., 1980. A closed-form equation for predicting the hydraulicconductivity of unsaturated soils. Soil Sci. Soc. Am. J., 44:892-898.
van Genuchten, M.Th., 1985. Convective-dispersive transport of solutes involved insequential first-order decay reactions. Comp. Geosci., 11:159-147.
van Genuchten, M.Th., and P.J. Wierenga, 1976. Mass transfer studies in sorbing porous
meida I. Analytical solutions. Soil Sci. Soc. Am. J., 40:473-480.
Vogel, T., and M. Cislerova, 1988. On the reliability of unsaturated hydraulic conductivity
calculated from the moisture retention curve. Transport in Porous Media, 3:1-15.
Wagenet R. J., and J. L. Hutson. 1987. LEACHM: Leaching Estimation And CHemistry Model,
A process-based model of water and solute movement, transformations, plant uptake andchemical reactions in the unsaturated zone, Continuum 2, Dept. of Agronomy, Cornell
University, Ithaca, New York, NY.
Walter, A.L., E.O. Frind, D.W. Blowes, C.J. Ptacek, and J.W. Molson, 1994. Modeling ofmulticomponent reactive transport in groundwater. 1. Model development and evaluation.
Water Res. Res., 30:3137-3148.
Xu, T., 1996. Modeling nonisothermal multi-component reactive solute transport throughvariably saturated porous media. Civil Engineering School, University of La Coruña,
Coruña, Spain, 310 p.
Yeh, G.-T., and H.-P. Cheng, 1999. 3DHYDROGEOCHEM: A 3-dimensional model ofdensity-dpeendent subsurface flow and thermal multispecies-multicomponent
HYDROGEOCHEMical transport. EPA/600/R-98/159, 150p.
Zheng, C., 1990. A modular three-dimensional transport model for simulation of advection,
dispersion and chemical reaction of contaminants in groundwater systems. Technicalreport, U.S. Enviromental Protection Agency.
Zheng, C., and P.P. Wang, 1998. MT3DMS, A modular three-dimensional transport model.
Technical report, U.S. Army Corps of Engineers Waterways Experiment Station,
Vicksburg, Miss.

75

APPENDIX A – LISTING OF NOTATION

aie

activity of ieth secondary exchange species[M L-3]

A
Ai
Ai
Am,0
Aiee

amplitude of temperature sine wave at soil surface
chemical formula of ith secondary species
chemical formula of ith secondary species
initial reactive surface area of a mineralchemical formula of ieth secondary exchange species

[K]
[-]
[-]
[L2]
[-]

Ali
Api
Amj
b
bi
bieje

chemical formula of ith aqueous secondary species
chemical formula of ith mineralchemical formula of jth master species
normalized water uptake distribution
parameters to calculate soil thermal conductivity (i=1,2,3)
number of equivalents of the ieth secondary exchange species on the jeth masterexchangeraqueous concentration of ith speciesaqueous concentration of ith species in the immobile region
aqueous concentration of ith species in the mobile region
aqueous concentration of ith aqueous secondary speciesaqueous concentration of jth master species
concentration of the sink term
concentration of ith secondary species in the sink term
concentration of jth master species in the sink term
aqueous concentration
volumetric heat capacity of the gas phasetotal concentration of jth master species
Courant number
total concentration of jth master species in the sink term
volumetric heat capacity of solid phasevolumetric heat capacity of organic mattervolumetric heat capacity of porous medium
volumetric heat capacity of liquid phasedispersion coefficient of ith species in liquid phasedispersion coefficient of ith species in free water
longitudinal dispersivity
evaporationsoil water pressure headosmotic headindex number of secondary species
index number of secondary exchange species
index number of master species
index number of master exchanger
index number of solute in a decay chain
intrinsic rate constant of mineral m
unsaturated hydraulic conductivity
empirical adsorption coefficient

[-]
[-]
[-]
[L-1]
[ML-1T-2K-1][equivalents]

ci
ci,im
ci,m
cli
cjm
cr
cr,il
crm,j
C
Ca
Cj
Cr
Crj
Cn
Co
Cp
Cw
Dwi
Di,w
DL
E
H
hφ
i
ie
j
je
k
k0,m
K
Kd

[M L-3][M L-3][M L-3][M L-3][M L-3][M L-3][M L-3][M L-3][L L-3][MLT-3K-1][M L-3][-][M L-3][MLT-3K-1][MLT-3K-1][MLT-3K-1][MLT-3K-1][L2 T-1][L2 T-1]
[L]
[L]
[L]
[L]
[-]
[-]
[-]
[-]
[-]
[M L-2 T-1]
[L T-1][Ma(1-n)L3nMs-1]

76
Kdw
Kieé

adsorption coefficient in mass per unit of volume water
equilibrium constant of reaction for ieth secondary exchange species[M(1-n)L3nL-3]
[-]

Kli
Kpi
Ks
m
m
Mjt
Mjt,i
M0j

[-]
[-]
[L T-1][-][-][M L-2][M L-2][M L-2]

l
n
n
nF
nieje

equilibrium constant of reaction for ith aqueous secondary species
equilibrium constant of reaction for ith mineralsaturated hydraulic conductivity
parameter in the water retention curve
index number of a mineralamount of the jth master species in the flow region at time t
amount of the jth master species in element i at time t
amount of the jth master species in the flow region at the beginning of thesimulation
amount of the jth master species in element i at the beginning of the simulationtpore connectivity parameterparameter in the water retention curve
reaction order of the OH- dependence of the reaction rateempirical Freundlich exponent
moles of ieth secondary exchange species on jeth master exchanger

[-]
[-]
[-]
[-]
[moles]

nm
nm,0
Nlk
Nm
Np
Ns
Nsa

number of moles of mineral minitial number of moles of mineral mnumber of kinetically-controlled homogeneous reactions
number of master species
number of minerals
number of secondary species
number of secondary aqueous species[moles][moles]
[-]
[-]
[-]
[-]
[-]

Nse

number of secondary exchange species

[-]

NX
P
Pe
Pt
q
qsj0
qsjN
Qjin
Qjin
Qmp
R
Ri
Ri,im

number of master exchange sites
precipitation
Peclet number
period of time necessary to complete one temperature cycle (one day)volumetric flux density
solute flux at the bottom of the soil profile of the jth master species
solute flux at soil surface of the jth master species
cumulative inflow of the jth master species
cumulative outflow of the jth master species
ion activity product of mineral mdecay rate
source/sink term due to geochemical reactions for ith species
source/sink term due to geochemical reactions for ith species in immobile
regionsource/sink term due to geochemical reactions for ith species in mobile region
dissolution rate of a mineralhomogeneous equilibrium reactions of source/sink termheterogeneous equilibrium reactions of source/sink term
homogeneous kinetic reactions of source/sink termheterogeneous kinetic reactions of source/sink termwater sink term
effective saturation
potential root water uptake

[-]
[L]
[-]
[T]
[L T-1][M L-2 T-1]
[M L-2 T-1]
[M L-2][M L-2][-][M L-3 T-1]
[M L-3 T-1]
[M L-3 T-1]

M0j,i

Ri,m
Rm
Req,hom
Req,hetRkin,homRkin,hetSSe
Sp[M L-2][M L-3 T-1]
[M T-1][M L-3 T-1]
[M L-3 T-1]
[M L-3 T-1]
[M L-3 T-1]
[L3 L-3 T-1]
[-]
[T-1]

77
Ss
Sw
t
T
Tavg
Tje

adsorbed concentration
adsorbed concentration per unit volume of watertime
temperature
average temperature at soil surface
total number of exchange sites for jeth master exchanger[MaMs-1][MaL-3][T][T][K][equivalents]

Tp
v
Vm
Vm,0
x
Xjme

potential transpiration
average pore water velocity
volume of mineral minitial volume of mineral m
spatial coordinatechemical formula of jeth master exchanger[L T-1][L T-1]
[L3]
[L3]
[L]
[-]

α
α
α
βie,jedimensionless water stress response function
parameter in the water retention curve
angle between flow direction and the vertical axis
equivalent fraction of ieth exchange species on jeth exchanger

[-]
[L-1]
[-]

βt

thermal dispersivity
absolute error in the solute mass balance[L][M L-2]

relative error in the solute mass balance of the jth master species[%]

activity coefficient of ieth secondary exchange species[-]

activity coefficient of ith aqueous secondary species
activity coefficient of jeth master exchanger

[-]
[-]

activity coefficient of jth master species
volumetric water content
immobile water content
mobile water content
volumetric fraction of the solid phasevolumetric fraction of the organic matter
residual water content
saturated water content
volumetric fraction of the air phaseapparent thermal conductivity of the soil
thermal conductivity of the soil
first-order rate constant for solutes in the liquid phasefirst-order rate constants for decay chain solutes in the liquid phasestoichiometric coefficient of jth master species in ith secondary species
stoichiometric coefficient of jth master species in ieth secondary exchange
speciesstoichiometric coefficient of jeth master exchanger in ieth secondary exchange
speciesstoichiometric coefficient of jth master species in ith aqueous secondary
speciesstoichiometric coefficient of jth master species in ith homogeneous kinetic
reactionstoichiometric coefficient of jth master species in ith mineralbulk density
tortuosity factor in liquid phasemass transfer coefficient for ith species[-][L3 L-3][L3 L-3][L3 L-3][L3 L-3][L3 L-3][L3 L-3][L3 L-3][L3 L-3][MLT-3K-1][MLT-3K-1]
[T-1]
[T-1]
[-]
[-]

εa
εrj
γiee
γli
γjee
γmj
θ
θim
θm
θn
θo
θr
θs
θv
λ
λ0
µw
µw’
νji
νej ieνjeeieνljiνlkj i
νpj i
ρ
τw
ωi

[-]
[-]
[-]
[-]
[MsL-3]
[-]
[T-1]

78
ωs

performance index

[-]

79

APPENDIX B – RELEASED VERSIONS AND BENCHMARKING
B.1 Version 1.0
Released date: November 2004
Status: beta version – controlled distribution
Versions used of the original codes:• HYDRUS-1D, Version (Version 2.0, July 1998)• PHREEQC (Version 2.4)
• Manual: Jacques, D., and J. Šimůnek, 2005. User manual of the Multicomponent
variably-saturated transport model HP1 (Version 1.0): Description, Verificationand Examples. SCK•CEN, Mol, Belgium, BLG-998, 79.
Overview of verification problems• Steady-state physical (non)equilibrium transport of chloride (Verification
Problem 1 - HYDRUS-1D – EQCL and NEQCL)• Transient physical nonequilibrium transport of chloride (Verification Problem 2- HYDRUS-1D – TRANSCL• Steady-state transport of a nonlinearly adsorbing contaminant (Verification
Problem 3 - HYDRUS-1D / PHREEQC – STADS)• Steady-state transport of a nonlinearly adsorbing contaminant subject to firstorder decay (Verification Problem 4 - HYDRUS-1D, PHREEQC – STDECAY)• Transport of first-order decay chain of nonlinearly adsorbing contaminants
during unsteady flow (Verification Problem 5 - HYDRUS-1D –SEASONCHAIN)• Transport of heavy metals subject to multiple cation exchange (Verification
Problem 6 – PHREEQC – CATEXCH)• Transport with mineral dissolution (Verification Problem 7 – PHREEQC MINDIS)• Heavy metal transport with a pH-dependent cation exchange complex(Verification Problem 8 – CRUNCH – MCATEXCH)• Infiltration of a hyperalkaline solution in a clay sample (Verification Problem 9– CRUNCH – ALKALINE)
Overview of examples• Long term transient flow and transport of major cations and heavy metals in asoil profile (HEAVMET)