HP2/3: Extensions of the HP1 Reactive Transport Code to Two
and Three DimensionsJiří Šimůnek1, Diederik Jacques2, and Miroslav Šejna3
1

Department of Environmental Sciences, University of California Riverside, Riverside, CA, USA,jiri.simunek@ucr.edu
2
Performance Assessment, Belgian Nuclear Research Institute, Mol, Belgium, djacques@sckcen.be3
PC-Progress, Ltd., Prague 2, Czech Republic, m.sejna@pc-progress.cz

Abstract
A large number of interacting physical, chemical, and biological processes determine the fate of major
cations and anions, contaminants (e.g., heavy metals, pesticides), and colloids in soil systems. Also,modeling of soil CO2 sequestration requires consideration of water flow, heat transport, gas diffusion, andmicrobiological soil respiration processes. The HP1 simulator (Jacques et al., 2008) is a state-of-the-art
model that was specifically developed to evaluate these and similar processes in the unsaturated zone.
HP1 couples the one-dimensional variably-saturated flow and transport model HYDRUS-1D (Šimůnek etal., 2008) with the generic geochemical model PHREEQC (Parkhurst and Appelo, 1999). Although HP1
is a versatile code for implementing different geochemical and transport conceptual models in variablysaturated porous media, one of its main limitations is that flow and transport is restricted to one
dimension. However, there are many applications that require consideration of various processes in either
two- or three-dimensions. Typical examples are flow and transport in soil systems with tiled drains, orwith drip and furrow irrigation, or in sloped, layered or heterogeneous systems. Therefore, PHREEQC has
recently been coupled also with HYDRUS (2D/3D) (Šimůnek et al., 2011) to handle flow and transport
problems, which require a higher dimensionality (HP2/3), and this manuscript describes the capabilitiesof this new tool. Typical two-dimensional flow and transport problems illustrating the capability of the
HP2 simulator are presented in this manuscript. The HP2/3 code uses the graphical user interface ofHYDRUS (2D/3D) (Šejna et al., 2011) for input and output processing, enabling definitions and finite
element discretization of very complex flow domains.

1. Introduction
The one-dimensional program HP1, which couples the PHREEQC geochemical code
(Parkhurst and Appelo 1999) with HYDRUS-1D, was first released in 2005 (Jacques andŠimůnek 2005). This comprehensive simulation tool (HP1 is an acronym for HYDRUSPHREEQC-1D) can simulate (1) transient water flow, (2) the transport of multiple components,(3) mixed equilibrium/kinetic biogeochemical reactions, and (4) heat transport in onedimensional variably-saturated porous media (soils). The code uses the Richards equation forsimulating one-dimensional variably-saturated water flow and advection-dispersion typeequations for heat and solute transport. However, the loosely coupled program can also simulatea broad range of low-temperature biogeochemical reactions in water, the vadose zone and inground water systems, including interactions with minerals, gases, exchangers and sorption
surfaces based on thermodynamic equilibrium, kinetic, or mixed equilibrium-kinetic reactions.
HP1 uses the operator-splitting approach with no iterations during one time step (a non-iterative
sequential modeling approach).
1Jacques and Šimůnek (2005), and Šimůnek et al. (2006), and Jacques et al. (2008ab, 2013),
demonstrated the versatility of HP1 in several examples, which included a) the transport ofheavy metals (Zn2+, Pb2+, and Cd2+) subject to multiple cation exchange reactions, b) transportwith mineral dissolution of amorphous SiO2 and gibbsite (Al(OH)3), c) heavy metal transport ina medium with a pH-dependent cation exchange complex, d) infiltration of a hyperalkaline
solution in a clay sample (this example considers kinetic precipitation-dissolution of kaolinite,
illite, quartz, calcite, dolomite, gypsum, hydrotalcite, and sepiolite), e) long-term transient flow
and transport of major cations (Na+, K+, Ca2+, and Mg2+) and heavy metals (Cd2+, Zn2+, and Pb2+)
in a soil profile, f) cadmium leaching in acid sandy soils, g) radionuclide transport, and h) longterm uranium migration in agricultural field soils following mineral P-fertilization.
Although HP1 is a versatile code for implementing different geochemical and transport
conceptual models in variably-saturated porous media, one of its main limitations is that flow
and transport is restricted to one dimension. However, there are many applications, such as flow
and transport in soil systems with tiled drains, or with drip and furrow irrigation, or in sloped,
layered or heterogeneous systems, that require consideration of various processes in either twoor three-dimensions. Therefore, PHREEQC has recently been coupled also with HYDRUS(2D/3D) (Šimůnek et al., 2011) to handle flow and transport problems, which require a higher
dimensionality (HP2/3), and this manuscript describes the capabilities of this new tool.
HP2/3 have, apart from the dimensionality (2D or 3D, respectively), the same capabilities as
HP1. The HP2/3 module may be used to analyze water flow, solute movement, and biogeochemical
reactions in unsaturated, partially-saturated, or fully saturated two- and three-dimensional porous
media. HP2/3 can handle flow domains delineated by irregular boundaries. The flow region itself
may be composed of nonuniform soils having an arbitrary degree of local anisotropy. Flow andtransport can occur in the vertical plane, the horizontal plane, or in a three-dimensional region
exhibiting radial symmetry about a vertical axis. The water flow part of the model considers
prescribed head and flux boundaries, as well as boundaries controlled by atmospheric conditions.
The HP2/3 code uses the graphical user interface of HYDRUS (2D/3D) (Šejna et al., 2011) for
input and output processing, enabling definitions and finite element discretization of very
complex flow domains. Typical two-dimensional flow and transport problems illustrating thecapability of the HP2 simulator are presented in this manuscript.

2. HYDRUS (2D/3D) Graphical User Interface for HP2/3
The HP2 code is fully incorporated into the HYDRUS (2D/3D) software package, and hence isinstalled automatically, together with selected examples, when one obtains HYDRUS (2D/3D)
and HP2 licenses and downloads HYDRUS from the Hydrus website. Here, we focus mostly on
the implementation of the HP2 module into HYDRUS (2D/3D). We will not discuss processes thatare related to geometry design, finite element discretization, variably-saturated water flow, heattransport, solute transport, or initial and boundary conditions, which are described in detail in the
HYDRUS (2D/3D) documentation (Šimůnek et al., 2011; Šejna et al., 2011). Nor will we discuss2

processes and reactions related to biogeochemical reactions that are described in detail in the
PHREEQC (Parkhurst and Apello, 1998] and HP1 (Jacques and Šimůnek, 2005, 2010) manuals.
2.1. Common Features with HP1
Similarly as HP1, the HP2 module is activated in the "Main Processes" window by selecting the"Solute Transport" check box and the "HP2 (Hydrus + Phreeqc)" radio button. The number of"Components" is again selected in the "Solute Transport - General Information" window. Largeparts of this window (e.g., iteration criteria) are disabled, since they are not relevant to HP2applications. Also the "HP2/3 Components and Database Pathway" window, in which users haveto select the thermodynamic database to be used with HP2 calculations and to list all main
components, is the same as in HP1.
2.2. HP2/3 Definitions Window
The main difference between HP1 and HP2 is in how PHREEQC-related inputs (i.e., those thatare printed into the PHREEQC.IN input file) are handled. The "HP2/3 Definitions" window has
been completely rewritten and reorganized. While in HP1, we had one "HP1 Definitions"window with four external text editors (Fig. 1), in HP2/3, the external text editors were converted
into four Internal pages (Fig. 2) (1 - Additions to Thermodynamic Database, 2 - Definitions of
Solution Compositions, 3 - Geochemical Model, and 4 - Addition Output), which can all beviewed directly in the "HP2/3 Definitions" window. One additional Page (5 - Database Viewer),in which one can view the selected thermodynamic database, has been added. The five pages canbe selected in the top left corner ("List of Pages") of the "HP2/3 Definition" window (Fig. 2).

Figure 1. HP1 Definitions window with four external edit windows.

This allowed us to include a section called, "Keywords (double-click to insert)" (bottom leftpart), which offers the most commonly used PHREEQC Keywords that are used in the first four
PHREEQC Pages. The Keywords are hierarchically organized in a tree-like structure in sevenmain groups: Solution Definition, Geochemical Model, Output, Chemical/Physical Reaction,
Database, Advanced, and Miscellaneous. A single click opens a particular tree subsection and adouble click copies the selected keyword to the Editor Window. This arrangement should
dramatically simplify the work of HP2/3 users, and prevent many common errors that they may
3have made with HP1. Note the most commonly used Solution Keywords in Figure 2 below density, -isotope, -pe, -pH, -redox, -temperature, -units, and -water.

Figure 2. HP2/3 Definitions window with an open page Definitions of Solution Compositions.

Additional new commands were included below the Editor window. For example, the "New
Solution" appears below the "Definitions of Solution Compositions" page in Fig. 2. Clicking thiscommand will result in including a template for a new solution (Fig. 3). This solution will have atemperature of 25oC, a pH of 7, be in equilibrium with the atmospheric content of oxygen andhave zero concentrations of Ca2+ and Cl-.
solution 1001
-temp 25
-units mmol/kgw-pH 7 charge
Ca 0
Cl 0
O(0) 1 O2(g) -0.68

Figure 3. A template for a New Solution.

Similar commands (e.g., Add Exchange, Surface, Equilibrium Phases, and Kinetics) appear
below the "Geochemical Model" page. These commands can be used to include templates fordefining cation exchange, surface complexation, and precipitation/dissolution reactions. For
more details about the HP2/3 related GUI, see the HP2 Manual (Šimůnek et al., 2012a).

4

2.3. HP2/3 Output
The output for the HP2 module is similar to the output for the standard HYDRUS module and forstandard variables, such as pressure head, water contents, and so on. Multiple variables can bedisplayed in the View window (Figure 4). Main components defined in the "HP2/3 Components
and Database Pathway" window, variables selected in the "HP2 Print and Punch Controls"
window, and variables specified in the "Additional Output" page of the "HP2/3 Definitions"window (Fig. 2) can be displayed this way.
Figure 4 below shows the "Results - Graphical Display" part of the Data Tab of the Navigator
Bar for the Leaching of the Uranium Tailings example given below (Section 3.2). From the HP2variables, the first seven are the main components (Total_H, Total_O, Ca, C, U, P, S, and Fe)
defined in the "HP2/3 Components and Database Pathway" window, the next two (pH and pe) arevariables selected in the "HP2 Print and Punch Controls" window, and finally the last four
(calcite, gypsum, ratherfordine and siderite) are mineral phases specified in the "Additional
Output" page of the "HP2/3 Definitions" window.

Figure 4. The "Results - Graphical Display" part of the Data Tab of the Navigator Bar for the HP2
module.

3. Example Problems with HP2
Most of the HP1 examples discussed above in the Introduction and provided with the HP1
installation have been rerun using HP2 to verify correct implementation of various componentsof the coupled program. Two additional two-dimensional examples are discussed below. All ofthese examples can be downloaded from the HYDRUS/HP2 website (the links are provided
below).

5

3.1. Furrow Irrigation with Cation Exchange
A furrow irrigation problem, similar to the one used in the UNSATCHEM manual (Šimůnek etal., 2012b), was used to simulate two-dimensional infiltration of gypsum saturated water into asodic soil. The example thus simulates a sodic soil reclamation problem and demonstrates thecation exchange feature of HP2. The schematic representation of the flow domain for theconsidered furrow irrigation is presented in Figure 5, together with the finite element mesh. It is
assumed that every other furrow is flooded with water and that the water level in the irrigatedfurrow is kept constant at a level of 6 cm. Due to symmetry, it is necessary to carry out thesimulation only for the domain between the axis of two neighboring furrows. Free drainage isused as the bottom boundary condition and zero flux is considered on the rest of the boundary.
The initial pressure head condition is -200 cm, and the soil hydraulic properties for silt are used.
15

15

40

15

15

15

85

100 cm100 cmFigure 5. Schematic representation and finite element mesh of the flow domain for the furrow
irrigation system.

The calculation was run at a constant temperature of 25 oC. The bulk density of the soil was
taken as 1.4 g cm-3 and molecular diffusion as 2 cm2day-1. Longitudinal and transverse
dispersivities were equal to 2 and 0.2 cm, respectively.
The solution composition of the water initially present in the soil profile is that of the following
highly sodic water: CaT=1.0 MgT=0.0, NaT=5.0, KT=0.0, SO4T=3.5, ClT=0.0 mmol L-1. The
cation exchange capacity is equal to 7.143 mmol kg-1 (10 mmol dm-3) and is divided between
exchangeable calcium and sodium (it is equilibrated with the solution). The solution compositionof the irrigation water was almost gypsum saturated: CaT=16.3, MgT=0.0, NaT=4.4, KT=0.0,
ClT=5.0, SO4T=16.0 mmolL-1. As a consequence of the reactions of the irrigation water with the6exchanger composition, cation exchange was the dominant chemical processes in the soil profile.
Cation exchange is treated as an instantaneous process in the model.
Figure 6 below presents the exchangeable concentrations of sodium to demonstrate typical output of
HP2. The exchange phase concentrations reflect the changes in aqueous Na and Ca concentrations.

Figure 6. Exchangeable concentrations of sodium (mol kgw-1) profiles at times: a) 0.5, b) 1, c) 3, andd) 5 days.

3.2. Leaching of the Uranium Tailings
This problem was inspired by one reported by Yeh and Tripathi (1991), and we have modified it
to make it more realistic. The problem considers the release and migration of uranium from a
simplified uranium mill tailings pile towards a river. The schematic of the transport domain isshown in Figure 7. The mill tailings pile is located adjacent to a surface that slopes down to ariver. The medium has the hydraulic properties of a loam with the saturated hydraulicconductivity of Ks = 3.78 m/day.
The horizontal bottom of the region is impermeable. The vertical left edge has the Dirichletboundary condition with a groundwater table 12 m above the bottom of the transport domain(exactly in the middle of the boundary). The top boundary (except for the mill tailing pile and theriver) has a flux boundary condition with a net rainfall rate of 0.139 cm/day. The horizontal
region on the top of the mill tailings pile is a Cauchy flow boundary with an infiltration rate of
1.39 cm/day. The nodes on the vertical line on the right side and the nodes on the river bottom
7

have the Dirichlet boundary condition reflecting the position of water in the river (4.5 m above
the bottom of the transport domain). A hypothetical pumping well with a withdrawal rate of
271.58 cm2/day is located at (x, z) = (400, 100). The region is discretized using a structured FE
mesh with 1564 elements and 852 nodes.
Cauchy Flow
= 0.139 dm/day
Mill Tailing
Pile

Infiltration Rate
= 0.0139 dm/day
Dirichlet Flow
Nodes

240 dmPumping Well

39 dm
1,050 dm

Figure 7. Problem description for the Uranium Tailing Problem (Yeh and Tripathi, 1991).

Chemical Reactions (35 Aqueous Complexation Reactions):
CaCO3(aq), CaHCO3+, CaSO4(aq), CaH2PO4+, CaPO4-, CaHPO4 (aq), CaOH+, FeSO4 (aq), FeOH+, Fe(OH)2 (aq),
Fe(OH)3-, Fe(OH)42-, (UO2)(OH)+, (UO2)2(OH)22+, (UO2)3(OH)42+, (UO2)3(OH)5+, (UO2)4(OH)7+, (UO2)3(OH)7-,(UO2)(CO3)(aq), (UO2)(CO3)22-, (UO2)(CO3)34-, (UO2)2(CO3)(OH)3-, (UO2)(SO4)(aq), (UO2)(SO4)22-,
H2(UO2)(PO4)+, H3(UO2)(PO4)2+, CaH4(UO2)(PO4)22+, CaH5(UO2)(PO4)23+, HCO3-, H2CO3(aq), HSO4-, HPO42-,
H2PO4-, and H3PO4

and Precipitation-Dissolution Reactions (14):
CaSO4(s), CaCO3(s), Ca5(OH)(PO4)3(s), FeCO3(s), Ca(UO2)2(PO4)2(s), Ca4H(PO4)3(s), CaH(PO4)(s), Ca(OH)2(s),Fe3(PO4)2(s), Fe(OH)2(s), (UO2)(OH)2(s), (UO2)(CO3)(s), Fe(UO2)2(PO4)2 (s), and H(UO2)(PO4)(s)and their thermodynamic equilibrium constants considered in this example are listed in Table 3of the HP2 manual. For reactive hydrogeochemical transport, the problem consists of eightcomponents: Total H, Total O, Ca, C, uranium, sulfate, phosphate, and Fe. A total of 35 species
and 14 minerals is defined for the problem; redox reactions were not considered. Table 1 lists
concentrations of the tailing and recharge waters.
Table 1. Initial and boundary compositions of recharge water and pore water in the tailings and regions
outside of the tailing for the uranium tailing problem.
Species

Inside the Tailing

Outside the Tailing

Recharge Water

Ca2+

1.0 x 10-2

1.0 x 10-2

1.0 x 10-3CO32-

1.0 x 10-2

1.5 x 10-3

1.5 x 10-3

UO22+

5.0 x 10-4

1.0 x 10-7

1.0 x 10-8PO43-

1.0 x 10-6

1.0 x 10-6

1.0 x 10-6SO42-

2.0 x 10-1

2.0 x 10-2

1.0 x 0-4H+

2.0 x 10-1

1.0 x 10-31.ox 10-3

Fe2+

3.5 x 10-2

1.0 x 10-7

1.0x 10-7

8The steady-state pressure heads and velocity fields are depicted in Figure 8. Uranium
concentration profiles are shown in Figure 9. pH, calcite, and gypsum profiles are shown in
Figure 10.

Figure 8. The steady-state pressure head (cm) (top) and flux (cm/d) (bottom) profiles for the Tailing
Pile Leaching example.

Figure 9. Uranium concentration profiles at time 0 (left), and 500 (right) d for the Tailing Pile
Leaching example.

Figure 10. pH (top), calcite (mol/L) (middle), and gypsum (mol/L) (bottom) profiles after 1000 d for
the Tailing Pile Leaching example.

9

4. Important Links
http://www.pc-progress.com/en/Default.aspx?h3d2-HP2: web page providing brief description of
the HP2 module and listing all main documents (reports and manuscripts) related to this module.
http://www.pc-progress.com/en/Default.aspx?H3D2-lib-HP2: Examples demonstrating the HP2
program using one-dimensional (for comparison with HP1) and two-dimensional (those
discussed above) problems.
http://www.pc-progress.com/Documents/HYDRUS3D_HP2_Manual.pdf: HP2 manual.

References

Jacques, D., and J. Šimůnek, User Manual of the Multicomponent Variably-Saturated Flow and Transport
Model HP1, Description, Verification and Examples, Version 1.0, SCK•CEN-BLG-998, Waste andDisposal, SCK•CEN, Mol, Belgium, 79 pp., 2005.
Jacques, D., J. Šimůnek, D. Mallants, and M. Th. van Genuchten, Modelling coupled hydrogeological andchemical processes in the vadose zone: a case study of long-term uranium transport following mineral
P-fertilization, Vadose Zone Journal, 7(2), 698-711, 2008a.
Jacques, D., J. Šimůnek, D. Mallants and M. Th. van Genuchten, Modelling coupled water flow, solutetransport and geochemical reactions affection heavy metal migration in a Podzol soil, Geoderma, 145,
449-461, 2008b.
Jacques, D., J. Šimůnek, D. Mallants, M. Th. van Genuchten, The HPx reactive transport models:
Summary of recent developments and applications, (this issue), 2013.
Parkhurst D. L., and C. A. J. Appelo, User’s guide to PHREEQ C (Version 2) – A computer program forspeciation, batch-reaction, one-dimensional transport and inverse geochemical calculations, WaterResources Investigations, Report 99–4259, Denver, Co, USA, 312 pp., 1999.
Šejna, M., J. Šimůnek, and M. Th. van Genuchten, The HYDRUS Software Package for Simulating Twoand Three-Dimensional Movement of Water, Heat, and Multiple Solutes in Variably-Saturated Media,
User Manual, Version 2.0, PC Progress, Prague, Czech Republic, pp. 280, 2011.
Šimůnek, J., D. Jacques, M. Th. van Genuchten, and D. Mallants, Multicomponent geochemical transportmodeling using the HYDRUS computer software packages, J. Am. Water Resour. Assoc., 42(6),1537-1547, 2006.
Šimůnek, J., M. Th. van Genuchten, and M. Šejna, Development and applications of the HYDRUS andSTANMOD software packages, and related codes, Vadose Zone Journal, 7(2), 587-600, 2008.
Šimůnek, J., M. Šejna, H. Saito, M. Sakai, and M. Th. van Genuchten, The HYDRUS-1D software
package for simulating the one-dimensional movement of water, heat and multiple solutes invariably-saturated media, Version 4.08, HYDRUS Software Series 3, Department of Environmental
Sciences, University of California Riverside, Riverside, California, USA, 2008.
Šimůnek, J., M. Th. van Genuchten, and M. Šejna, The HYDRUS Software Package for Simulating Twoand Three-Dimensional Movement of Water, Heat, and Multiple Solutes in Variably-Saturated Media,
Technical Manual, Version 2.0, PC Progress, Prague, Czech Republic, pp. 258, 2011.
Šimůnek, J., D. Jacques, M. Šejna, and M. Th. van Genuchten, The HP2 Program for HYDRUS (2D/3D): A
Coupled Code for Simulating Two-Dimensional Variably-Saturated Water Flow, Heat Transport, andBiogeochemistry in Porous Media, Version 1.0, PC Progress, Prague, Czech Republic, 76 pp., 2012a.
Šimůnek, J., M. Šejna, and M. Th. van Genuchten, The UNSATCHEM Module for HYDRUS (2D/3D)
Simulating Two-Dimensional Movement of and Reactions Between Major Ions in Soils, Version 1.0, PC
Progress, Prague, Czech Republic, 54 pp., 2012b.
Yeh, G. T., and V. S. Tripathi, A model for simulating transport of reactive multispecies components:
Model development and demonstration, Water Resour. Res., 27, 3075-3094, 1991.

10