Biologia 64/3: 465—469, 2009
Section Botany
DOI: 10.2478/s11756-009-0106-7

Selected HYDRUS modules for modeling subsurface flow
and contaminant transport as influenced by biological processes
at various scales
Jirka Šimůnek1, Diederik Jacques2, Navin K. C. Twarakavi3& Martinus Th. van Genuchten41Department of Environmental Sciences, University of California, Riverside, CA 92521, USA
Institute for Environment, Health and Safety, Belgian Nuclear Research Centre (SCK-CEN), Boeretang 200, B-2400 Mol,
Belgium
3
Agronomy and Soils, Auburn University, Auburn, AL 36849, USA
4
Department of Mechanical Engineering, COPPE/LTTC, Federal University of Rio de Janeiro, UFRJ, Rio de Janeiro,
RJ 21945-970, Brazil

A
ut
h
co o
py r's2

Abstract: A large number of modiﬁcations or special modules of the HYDRUS software packages have been developed
during the past several years to evaluate the eﬀects of a range of biohydrological processes on subsurface water ﬂow andthe transport of various chemicals and contaminants. The objective of this manuscript is to brieﬂy review the diﬀerent
modules that were included, and to present various applications illustrating the eﬀects of biological processes on water ﬂow
and solute transport and reactions in variably-saturated media.
Key words: biological processes; contaminant transport; CO2 production and transport; reactive transport modeling;
unsaturated

Introduction

There is increasing evidence that biological processes
can have major eﬀects on water ﬂow and the fate
and transport of contaminants in the subsurface atvarious temporal and spatial scales (Seuntjens et al.
2005; Rockwood et al. 2005; amongst others). Like
many other models, the HYDRUS software packages(Šimůnek et al. 2008) were originally developed to simulate only water ﬂow and heat and solute transport insoils, and as such ignored many important biologicalprocesses. During recent years a large number of modiﬁcations or special modules have been included in the
HYDRUS models to allow evaluation of possible eﬀects
of biological processes on ﬂuid ﬂow and contaminanttransport in the subsurface. The objective of this manuscript is to brieﬂy review these modules and to present
selected applications of the HYDRUS codes illustrating the eﬀects of biological processes on water ﬂow andsolute transport and reactions in soils. However, it isbeyond the scope of this paper to provide many details
about the modules. Readers are invited to ﬁnd thesedetails in referenced literature.
HP1 – Coupled HYDRUS-1D and PHREEQC
codes
A comprehensive modeling tool in terms of available
c 2009 Institute of Botany, Slovak Academy of Scienceschemical and biological reactions was recently developed by coupling HYDRUS-1D with the PHREEQC
geochemical code of Parkhurst & Appelo (1999). This
coupling resulted in the very ﬂexible simulator, HP1,
which is an acronym for HYDRUS1D-PHREEQC(Jacques & Šimůnek, 2005; Jacques et al. 2006). The
combined code contains modules simulating (1) transient water ﬂow in variably-saturated media, (2) thetransport of multiple components, (3) mixed equilibrium/kinetic biogeochemical reactions, and (4) heattransport. HP1 is a signiﬁcant expansion of the individual HYDRUS-1D and PHREEQC programs bycombining and preserving most or all of the featuresand capabilities of the two codes into a single numerical model. The code uses the Richards equation forvariably-saturated ﬂow and advection-dispersion typeequations for heat and solute transport. However, theprogram can now simulate also a broad range of lowtemperature biogeochemical reactions in water, the vadose zone and in ground water systems, including interactions with minerals, gases, exchangers, and sorption
surfaces, based on thermodynamic equilibrium, kinetics, or mixed equilibrium-kinetic reactions (Šimůnek et
al. 2008). Various applications of HP1 were presented
by Jacques & Šimůnek (2005) and Jacques et al. (2006,
2008a,b).
There is much evidence in the literature that microorganisms, present throughout the subsurface, can3.00E-04

1.60E-09
Sorbed Co
Sorbed CoNTA
Biomass

2.00E-04
1.00E-09
1.50E-04

8.00E-10
6.00E-10

1.00E-04
4.00E-10

Concentration [mol/kg]

2.50E-04

1.20E-09

Biomass [g/L]

Concentrations [mol/g]

1.40E-09

3.50E-06

6.30E-05

3.00E-06

6.25E-05

2.50E-06

6.20E-05
NTA
Co
O(0)

2.00E-06

6.15E-05

1.50E-06

6.10E-05

1.00E-06

6.05E-05

5.00E-07

6.00E-05

Concentration [mol/kg]

J. Šimůnek et al.

466

5.00E-05
2.00E-10
0.00E+00

0.00E+00
0

10

20

30

40

50

60

70

80

0.00E+00

5.95E-050Time [hours]

10

20

30

40

50

60

70

80

Time [hours]

A
ut
h
co o
py r'sFig. 1. Selected output of the HP1 program: a) sorbed concentrations and biomass, and b) total concentration of main components
and oxygen.

catalyze many chemical reactions and thus may affect the fate of both organic and inorganic contaminants. Many bacteria can beneﬁcially degrade or immobilize environmental contaminants (e.g., Ginn et al.
2002; Gargiulo et al. 2007b). As an illustration of theuse of HP1, we present here its application to a problem ﬁrst presented by Tebes-Stevens & Valocchi (1998),
and also used by Parkhurst & Appelo (1999), involving
bacterially mediated degradation of an organic contaminant. This example involves the transport of a pulse ofcobalt and nitrilotriacetate (NTA), an organic chelating agent that has been detected along with metals andradionuclides in some contaminated groundwater, into
a 10–m long soil column. The example accounts for several interacting chemical processes that are common
to many environmental problems: bacterially-mediateddegradation of an organic substrate (NTA), bacterialcell growth and decay, kinetic metal sorption (Co2+ ),and aqueous speciation including metal-ligand complexation (CoNTA− ). For this application we assumed
that only a single protonated form is biodegradable,
which leads to the following overall stoichiometry underaerobic conditions (Tebes-Stevens & Valocchi, 1998):
HNTA2− + 1.62 O2 + 1.272 H2 O + 2.424 H+ →→ 0.576 C5 H7 O2 N + 3.12 H2 CO∗3 + 0.424 NH+
4(1)

where C5 H7 O2 N is the simple elemental formula forbiomass (Rittmann & VanBriesen, 1996) and H2 CO∗3
represents both aqueous CO2 and H2 CO3 . Equation(1) assumes that electrons are transferred from thedonor substrate, HNTA2− , to the acceptor, O2 , in order to produce the energy required for synthesis of new
biomass. The overall rate of the above biological reaction is described using a Monod-type kinetic model asfollows:AS
r = rr Xm(2)
KS + S KA + A
where r [mol/sec/L] and rr [mol/sec/g biomass] are theactual and maximum reaction rates, respectively, Xm
is the biomass [g cells/L], S and A are the concentrations of the electron donor substrate (HNTA2− ) and theelectron acceptor (O2 ), respectively, and KS and KA are

the half-saturation constants for the substrate and theelectron acceptor, respectively. The change in biomass
is related to HNTA2− degradation rate by a yield coeﬃcient and to the actual biomass by a biomass decay
constant.
This example demonstrates several unique capabilities of PHREEQC, and hence also of HP1 which
preserved most of the PHREEQC geochemical features. For example, the HP1 code allows users to deﬁne species that may not appear in existing geochemical databases (e.g., NTA), as well as include additional
kinetic equations, which may be speciﬁed using short
Basic programs that are read as input (such as in thiscase the kinetic degradation of NTA, the corresponding growth of biomass, and kinetic sorption of metals and metal-ligand complexes). These features oﬀer
users of the program unprecedented ﬂexibility to deﬁnetheir own chemical problems and reactions. Figure 1
presents selected output from HP1 for the ﬁnal node of
the 10m long soil column. Only dissolved oxygen enters
the soil column in this saturated ﬂow problem. During the breakthrough of the NTA pulse, oxygen is consumed by the NTA degradation reaction. Although thisexample involves some sophisticated coupling between
chemical and biological processes, it is still limited to
one oxygen consuming reaction (NTA degradation), one
electron donor (NTA) and one bacterial population.
The particle (bacteria or viruses) transportmoduleBioremediation and/or bio-augmentation strategies toclean up recalcitrant chemicals in subsurface environments could be greatly enhanced by the eﬃcientdelivery of specialized bacteria to targeted locationsin the subsurface (Gargiulo et al. 2007b). At thesame time, bacterial contamination of groundwater can
be a serious problem that may lead to large outbreaks of waterborne diseases (e.g., Ginn et al. 2002;
Reynolds et al. 2004). An ability to accurately predict the fate and transport of pathogenic microorganisms in unsaturated soil environments is therefore needed in eﬀorts to protect surface and ground-HYDRUS modules for modeling subsurface ﬂow and contaminant transport

467

Fig. 2. Measured and simulated volumetric water contents (VWC) and CO2 ﬂuxes at the experimental ﬁeld in Davis, California(Buchner et al. 2008).

matrix grain sizes, at diﬀerent water saturations, and at
diﬀerent growth stages (metabolically active or stationary phases). With further modiﬁcations to account forbacteria aggregation, the model was additionally usedby Bradford et al. (2006) to simulate the transport andstraining of E. coli bacteria at the column scale. Yet
further modiﬁcations of this module allow considerationof the bacteria- (or colloid-) facilitated solute transport(Šimůnek et al. 2006; Pang & Šimůnek 2007).

A
ut
h
co o
py r'swater resources from contamination (Gargiulo et al.
2007b).
The fate and transport of bacteria in porous media is more complex than that of solutes since bacteria transport is aﬀected by several additional processes,
such as ﬁltration, straining, sedimentation, adsorption,and desorption, growth, and decay/death of bacteria.
Various modiﬁcations of the HYDRUS models havebeen made to allow the code to simulate these processes involved in the transport of viruses and bacteria.
The governing convection-dispersion equation for virus
or bacteria transport is similar as the general equation
for colloid transport assuming applicability of the attachment/detachment concept:
∂θc
∂s1
∂s2
+ρ
+ρ
=
∂t
∂t
∂t


∂
∂c
∂qc
− µw θc − µs ρ(s1 + s2 )
=
θD
−
∂x
∂x
∂x

(3)

where c is the (virus or bacteria) concentration in the
aqueous phase, s is the solid phase (virus or bacteria) concentration, the subscripts 1 and 2 represent
two kinetic sorption sites, respectively, θ is the watercontent, ρ is the bulk density, D is the hydrodynamic
dispersion coeﬃcient, t is time, x is the distance, andµw and µs represent inactivation and degradation coeﬃcients in the liquid and solid phases, respectively.
The latter two coeﬃcients, when assigned negative values, can also be used to account for bacteria growth(Gargiulo et al. 2007a). Two kinetic sorption sites canbe used in Equation (3) to represent various chemical or physical processes. One sorption site, for example, may be used to represent the chemical processes ofattachment/detachment, while the other sorption site
may represent the physical process of straining. Alternatively, one of the sorption sites could be used to represent attachment/detachment to/from air-water interfaces.
The transport model described above has beensuccessfully used by, among many others, Schijven &Šimůnek (2002) to simulate the transport of viruses
at the ﬁeld scale, and Gargiulo et al. (2007a,b; 2008)
to simulate the column-scale transport of hydrophobic(Rhodococcus rhodochrous) and hydrophilic (Deinococcus radiodurans) bacteria strains in sands of diﬀerent

The SOILCO2 module of HYDRUS-1D
Carbon dioxide production and oxygen consumption
represent important processes resulting from microbiological activity in the soil environment and in the unsaturated zone. Microbiological respiration, combinedwith plant root respiration and high water contents,
can change CO2 concentrations several orders of magnitude from values that are at equilibrium with the CO2
content of the atmosphere (about 0.036% at present),
up to extreme values of about 20% (Šimůnek & Suarez,
1993). The solubility of many solid phases, such as carbonates and oxihydroxides, signiﬁcantly changes within
this range of CO2 , primarily because changes in soil
CO2 lead to changes in the soil pH for all but acid soils(Šimůnek & Suarez 1994).
The SOILCO2 and UNSATCHEM modules (Šimůnek & Suarez 1993, 1994) of HYDRUS-1D consider
the eﬀects of CO2 producing microbiological activity
and CO2 transport in soils on geochemical transport.
We refer to the original codes for detailed descriptionsof the many processes that were considered. Suarez &Šimůnek (1993) discussed several examples of seasonal
simulations of soil CO2 concentrations and CO2 ﬂuxes
to the atmosphere for diﬀerent crops, and compared
their results with published ﬁeld data. The coupled program was obtained by combining various features from
SOILCO2 (i.e., water ﬂow, heat transport, and CO2transport) and UNSATCHEM (major ion chemistry)
within a new module within HYDRUS-1D. The resulting module of HYDRUS-1D has seen several applications recently to lysimeter and ﬁeld data (e.g., Baueret al. 2008; Herbst et al. 2008; Buchner et al. 2008).
Results from one application by Buchner et al. (2008)
are shown in Fig. 2.

J. Šimůnek et al.

468
Average watertable depths
MODFLOW Sub-model

Bottom fluxes arerecharge at thewater table for the
next MODFLOW
time stepHYDRUS Sub-model

Solve for the bottom fluxes
in each profile using
atmospheric data and the
1D Richards equation

A
ut
h
co o
py r'sWater table

Zone 1

Zone 2

Fig. 3. Schematic description of the coupling procedure for water ﬂow in HYDRUS-MODFLOW.

The HYDRUS module for MODFLOW

Water ﬂuxes in the subsurface of forested and agricultural ecosystems at the ﬁeld or catchment scale canbe highly variable spatially because of several interconnected processes. In addition to natural subsurface
heterogeneity, inﬁltrating water is often distributed unevenly due to above-ground interception and redistribution of rainfall by plant canopies (Sansoulet et al. 2008).
Transpiration and evaporation often also vary widelyover the landscape as a function of geographical location, slope exposure, plant species, and many other factors. These various processes can have important ecohydrological consequences since they signiﬁcantly aﬀect
groundwater recharge and nutrient leaching.
To be able to address spatially variable saturatedunsaturated hydrological processes at the larger scale,
Seo et al. (2007) and Twarakawi et al. (2008) recently developed the HYDRUS package for Modﬂow
2000 (Harbaugh et al. 2000). By fully incorporating
HYDRUS-1D within the widely used MODFLOW program, the resulting HYDRUS package provides MODFLOW with recharge ﬂuxes at the water table, while
MODFLOW provides HYDRUS with the position of
the groundwater table that is used as its bottom boundary. While the HYDRUS package considers the mainprocesses and factors aﬀecting ﬂuxes in the vadose zone,such as precipitation, inﬁltration, evaporation, redistribution, capillary rise, root-water uptake, water accumulation at the ground surface, surface runoﬀ, andsoil moisture storage, MODFLOW considers all of thegroundwater ﬂow processes. The combined software
thus dramatically expands capabilities of the two separate modules for application to a wide range of spatial
and temporal scales.
Figure 3 schematically shows the invoked couplingof the HYDRUS and MODFLOW models. The ground-water modeling domain for MODFLOW was discretized
into grids/blocks using ﬁnite-diﬀerences as described by
Harbaugh et al. (2000). Based on similarities in soilhydrology and topographical characteristics, the discretized MODFLOW domain can be divided into zones
that comprise one or more grids of the MODFLOW
model. A separate HYDRUS soil proﬁle is then assigned to each of these zones. The approach assumes
that the HYDRUS soil proﬁle adequately represents vadose zone ﬂow for the entire zone. The HYDRUS vertical soil proﬁles are subsequently discretized verticallyinto ﬁnite elements. Since the numerical solution of the
Richards equation requires much ﬁner temporal (often
as low as few seconds) as well as spatial discretizations
than traditional groundwater models, the two models(HYDRUS and MODFLOW) interact (i.e., exchange
information about the groundwater recharge rate andthe groundwater level) only at the end of each MODFLOW time step, with HYDRUS generally requiring
multiple time steps during each MODFLOW time step
(Seo et al. 2007). MODFLOW receives the recharge ﬂux
from HYDRUS and calculates a new water table depth
for a particular time step. A new water table depth isthen assigned as the pressure head bottom boundary
condition in the HYDRUS package for the next MODFLOW time step (Fig. 3). Twarakavi et al. (2008) recently compared the HYDRUS package to other contemporary modeling approaches and evaluated its performance for three case studies of increasing complexity
and spatial and temporal scales. They concluded thatthe HYDRUS package provided a much more eﬃcient
alternative and could account better for vadose zone
processes than other existing MODFLOW packages.
For large-scale ground water ﬂow problems, the HYDRUS package provides an optimal trade-oﬀ between
computational eﬀort and accuracy of model simulationsof coupled vadose zone – groundwater problems.

HYDRUS modules for modeling subsurface ﬂow and contaminant transport
References

Parkhurst D.L. & Appelo C.A.J. 1999. User’s guide to PHREEQC (Version 2) – A computer program for speciation, batchreaction, one-dimensional transport and inverse geochemical calculations. Water-Resources Investigations, Report 99–
4259, Denver, Co, USA, 312 pp.
Reynolds K.A. 2004. Groundwater Vulnerability to Microbial
Contamination. Water Conditioning & Puriﬁcation 56: 28–
30.
Rockhold M.L., Yarwood R.R., Niemet M.R., Bottomley P.J.
& Selker J.S. 2005. Experimental observations and numerical modeling of coupled microbial and transport processes invariably saturated sand. Vadose Zone J. 4: 407–417.
Sansoulet J., Cabidoche Y.-M., Cattan P., Ruy S. & Šimůnek J.
2008. Spatially distributed water ﬂuxes in an Andisol underbanana plants: experiments and 3D modelling. Vadose Zone
J. 7(Special Issue “Vadose Zone Modeling”): 819–829.
Seo H.S., Šimůnek J. & Poeter E.P. 2007. Documentation of the
HYDRUS Package for MODFLOW-2000, the U.S. Geological
Survey Modular Ground-Water Model, GWMI 2007–01,. Int.
Ground Water Modeling Center, Colorado School of Mines,
Golden, CO, 96 p.
Schijven J. & Šimůnek J. 2002. Kinetic modeling of virus transport at ﬁeld scale. J. Contam. Hydrol. 55: 113135.
Seuntjens P., Nowack B. & Schulin R. 2005. Root-zone modelingof heavy metal uptake and leaching in the presence of organic
ligands. Plant Soil. 265: 61–73.
Šimůnek J. & Suarez D.L. 1993. Modeling of carbon dioxidetransport and production in soil: 1. Model development. Water Resour. Res. 29: 487–497.
Šimůnek J. & Suarez D.L. 1994. Two-dimensional transportmodel for variably saturated porous media with major ion
chemistry. Water Resour. Res. 30: 1115–1133.
Šimůnek J., He C., Pang J.L. & Bradford S.A. 2006a. Colloidfacilitated transport in variably-saturated porous media: Numerical model and experimental veriﬁcation. Vadose Zone J.
5: 1035–1047.
Šimůnek J., van Genuchten M.Th. & Šejna M. 2008. Development
and applications of the HYDRUS and STANMOD software
packages, and related codes. Vadose Zone J. 7 (Special Issue“Vadose Zone Modeling”): 587–600.
Šimůnek J. & van Genuchten M.Th. 2008. Modeling nonequilibrium ﬂow and transport with HYDRUS. Vadose Zone J. 7
(Special Issue “Vadose Zone Modeling”): 782–797.
Suarez D.L. & Šimůnek J. 1993. Modeling of carbon dioxidetransport and production in soil: 2. Parameter selection, sensitivity analysis and comparison of model predictions to ﬁeld
data. Water Resour. Res. 29: 499–513.
Tebes-Stevens C. & Valocchi A.J. 1998. Numerical Solution Techniques for Reaction Parameter Sensitivity Coeﬃcients in Multicomponent Subsurface Transport Models. Hydraulic Engineering Series No. 59. UILI-ENG-98–2009. Department ofCivil and Environmental Engineering, University of Illinois
at Urbana-Champaign.
Twarakavi N.K.C., Šimůnek J. & Seo H. S. 2008. Evaluatinginteractions between groundwater and vadose zone using
HYDRUS-based ﬂow package for MODFLOW. Vadose Zone
J. 7 (Special Issue “Vadose Zone Modeling”): 757–768.

A
ut
h
co o
py r'sBauer J., Herbst M., Huisman J.A., Weihermüller L. & Vereecken
H. 2008. Sensitivity of simulated soil heterotrophic respiration
to choice of temperature and moisture reduction functions.
Geoderma 145: 17–27.
Bradford S.A., Šimůnek J. & Walker S.L. 2006. Transport andstraining of E. coli O157:H7 in saturated porous media. Water Resour. Res. 42, W12S12, doi:10.1029/2005WR004805,
12 pp.
Buchner J.S., Šimůnek J., Lee J., Rolston D.E., Hopmans J.W.,
King A.P. & Six J. 2008. Evaluation of CO2 ﬂuxes from anagricultural ﬁeld using a process-based numerical model. J.
Hydrology 361: 131–143.
Gargiulo G., Bradford S.A., Šimůnek J., Ustohal P., Vereecken H.
& Klumpp E. 2007a. Transport and deposition of metabolically active and stationary phase Deinococcus radiodurans in
unsaturated porous media. Environ. Sci. Technol. 41: 1265–
1271.
Gargiulo G., Bradford S. A., Šimůnek J., Ustohal P., Vereecken
H. & Klumpp E. 2007b. Bacteria transport and deposition
under unsaturated conditions: the role of the matrix grain
size and the bacteria surface protein. J. Contam. Hydrol. 92:
255–273.
Gargiulo G., Bradford S.A., Šimůnek J., Ustohal P., Vereecken H.
& Klumpp E. 2008. Bacteria transport and deposition underunsaturated conditions: the role of bacteria surface hydrophobicity. Vadose Zone J. 7: 406–419.
Ginn T.R., Wood B.D., Nelson K.E., Scheibe T.D., Murphy E.M.
& Clement T.P. 2002. Processes in microbial transport in thenatural subsurface. Adv. Water Resour. 25: 1017–1042.
Harbaugh A.W., Banta E.R., Hill M.C. & McDonald M.G.
2000. MODFLOW-2000, the U.S. Geological Survey modular
ground-water model user guide to modularization concepts
and the ground-water ﬂow process. Denver, CO, Reston, VA:
U.S. Geological Survey.
Herbst M., Hellebrand H.J., Bauer J., Huisman J.A., Šimůnek J.,
Weihermüller L., Graf A., Vanderborght J. & Vereecken H.
2008. Multiyear heterotrophic soil respiration: evaluation ofa coupled CO2 transport and carbon turnover model. Ecological Modelling 214: 271–283.
Jacques D. & Šimůnek J. 2005. User Manual of the Multicomponent Variably-Saturated Flow and Transport Model HP1, Description, Veriﬁcation and Examples, Version 1.0, SCKCENBLG-998, Waste and Disposal, SCKCEN, Mol, Belgium, 79
pp.
Jacques D., Šimůnek J., Mallants D. & van Genuchten M.Th.
2006. Operator-splitting errors in coupled reactive transport
codes for transient variably saturated ﬂow and contaminanttransport in layered soil proﬁles. J. Contam. Hydrol. 88: 197–
218.
Jacques D., Šimůnek J., Mallants D. & van Genuchten M.Th.
2008a. Modeling coupled hydrological and chemical processes:
Long-term uranium transport following mineral phosphorus
fertilization. Vadose Zone J. 7(Special Issue “Vadose Zone
Modeling”): 698–711.
Jacques D., Šimůnek J., Mallants D. & van Genuchten M.Th.
2008b. Detailled modeling of coupled water ﬂow, solute transport and geochemical reactions: Migration of heavy metals ina podzol soil proﬁle. Geoderma 145: 449–461.
Pang L. & Šimůnek J. 2006. Evaluation of bacteria-facilitatedcadmium transport in gravel columns using the HYDRUS
colloid-facilitated solute transport model. Water Resour. Res.
42, W12S10, doi:10.1029/2006WR004896.

469

Received October 22, 2008
Accepted January 22, 2009