Numerical Modelling of Hydrodynamics for Water Resources – Garcia-Navarro & Playán (eds)© 2008 Taylor & Francis Group, London, ISBN 978-0-415-44056-1

Modeling Subsurface Water Flow and Solute Transport with
HYDRUS and Related Numerical Software Packages
J. Šimůnek
Department of Environmental Sciences, University of California Riverside, Riverside, CA 92521, USA

M. Th. van Genuchten
US Salinity Laboratory, USDA-ARS, 450 West Big Springs Road, Riverside, CA 92507, USA

M. Šejna
PC-Progress, Ltd., Anglická 28, 120 00 Prague, Czech Republic

ABSTRACT: HYDRUS and related software packages represent a system of numerical models
for simulating variably-saturated water, heat and solute movement in the subsurface. The software packages can simulate flow and transport in one-dimensional (e.g., HYDRUS-1D, HP1,
SOILCO2, and UNSATCHEM), two-dimensional (SWMS_2D, CHAIN_2D, or HYDRUS-2D),and three-dimensional (SWMS_3D or HYDRUS (2D/3D)) transport domains. While water flowin all of the models is simulated in a relatively similar manner, different solute transport modules exist in the different codes. Some models can simulate only transport of solutes independent of one another, whereas other models include biogeochemical modules simulating complex
interactions between the various solutes. For example, HYDRUS-1D and UNSATCHEM models can also consider the transport of major ions and their mutual reactions, such as dissolutionprecipitation, cation exchange, or aqueous complexation. The wetland module in HYDRUS(2D/3D) furthermore considers biochemical processes associated with the flow of wastewater
through constructed or natural wetlands. Finally, HP1, a model resulting from coupling HYDRUS-1D and PHREEQC, accounts for a broad range of instantaneous and kinetic chemicaland biological reactions, including aqueous complexation, precipitation-dissolution, cation exchange, surface complexation, and redox reactions. In this manuscript we review the various
models and provide references to many of their applications.

1 INTRODUCTION
We describe in this paper a number of numerical models for simulating water flow and movement of heat and various solutes in the subsurface. The models were developed during the past
15 years or more in close collaboration between the University of California, Riverside, the U.S.
Salinity Laboratory in Riverside, PC-Progress in Prague, Czech Republic, and SCK•CEN, Mol,
Belgium. This manuscript serves as a companion paper to an article that will appear in the upcoming “Vadose Zone Modeling” Special Section of Vadose Zone Journal (Šimůnek et al., inpress). While that article qualitatively describes the various conceptual processes forming thebasis of the different programs, here we focus on the implemented mathematical models. We
provide the governing equations for the various water flow, heat transport, and solute transportmodules, and discuss into which HYDRUS programs they were implemented. We additionallydiscuss the numerical techniques used to discretize the governing flow and transport equations,and the importance of the developed graphical user interfaces to ensure wide acceptance of thesoftware packages. We further discuss selected applications of the HYDRUS programs to demonstrate different features of the codes. Our primary focus is on the programs listed in Table 1,
i.e., HYDRUS-1D, HP1, HYDRUS-2D, and HYDRUS (2D/3D).

95

Numerical Modelling of Hydrodynamics for Water Resources – Garcia-Navarro & Playán (eds)© 2008 Taylor & Francis Group, London, ISBN 978-0-415-44056-1

Table
1. Synopsis of the HYDRUS suite of windows-based software packages.
________________________________________________________________________________________________________Model
Ver Dim Brief Description
Reference________________________________________________________________________________________________________HYDRUS-1D

3.0 1D

Water flow, heat and solute transport, carbon dioxide Šimůnek et al.
transport, transport of major ions, major ion chemistry, (2005)
transport of colloids and bacteria.
HP1
1.0 1D Water flow, transport of heat and multiple components, Jacques &mixed equilibrium/kinetic biogeochemical reactions. Šimůnek (2005)
HYDRUS-2D
2.0 2D Water flow, heat and solute transport (discontinued,Šimůnek et al.
replaced with HYDRUS (2D/3D)).
(1999)HYDRUS (2D/3D) 1.0 2/3D Water flow, heat and solute transport, transport ofŠimůnek et al.
colloids and bacteria wetland module.
(2006b), Šejna &Šimůnek (2007)________________________________________________________________________________________________________

Ver – version, OS – operating system, Dim – dimensions, DOS – Disc Operating System, 1D – onedimensional, 2D - two-dimensional and axisymetrical three-dimensional, 3D – three-dimensional.

2 WATER FLOW
2.1 Uniform flow
Uniform variably-saturated water flow in all of these models is described using the Richards
equation:⎛⎞⎤∂θ (h)∂ ⎡∂h=+ K izA ⎟ ⎥ - S (h)⎢ K (h) ⎜⎜ K ijA⎟⎥∂t∂xi ⎣⎢
∂x j
⎝
⎠⎦

(1)

where θ = volumetric water content [L3L-3]; h = pressure head [L]; K = unsaturated hydraulicconductivity [LT-1]; KijA = components of a dimensionless anisotropy tensor KA (which reduces tothe unit matrix when the medium is isotropic); S = general sink/source term [L3L-3T-1] accounting
for root water uptake; t = time [T]; and xi = spatial coordinate [L].
Because of strongly nonlinear makeup of the Richard equation, only a relatively few simplified analytical solutions can be derived. Most practical applications of (1) require a numericalsolution, which can be obtained using a variety of methods such as finite differences or finite
elements. Equation (1) is generally referred to as the mixed form of the Richards equation since
it contains two dependent variables (i.e., the water content and the pressure head). Solutions of(1) require knowledge of two soil hydraulic functions describing the soil water retention andhydraulic conductivity properties of the porous material. The soil water retention curve (or also
called the soil water characteristic curve, the capillary pressure-saturation relationship, or the pF
curve) describes the relationship between the water content and the pressure head. The HYDRUS packages permit the use of analytical functions suggested by Brooks & Corey (1964),
van Genuchten (1980), Kosugi (1996), and Durner (1994):
λ⎧⎛ he ⎞
⎪
⎜ ⎟
⎝h⎠
⎪
⎪
1
⎪
[1+ |α h|n ]mθ - θr ⎪Se ( h ) ==⎨θs - θr ⎪ 1⎧ ln ( h / hm ) ⎫erfc ⎨
⎬
⎪2
2σ ⎭
⎩
⎪
wi
⎪
ni mi
⎪ ∑
⎩ i [1+ |αi h| ](2)

respectively, where Se = effective saturation [-];θr and θs = residual and saturated water contents[L3L-3], respectively; he = air-entry value [L]; λ = pore-size distribution index that characterizes

96

Numerical Modelling of Hydrodynamics for Water Resources – Garcia-Navarro & Playán (eds)© 2008 Taylor & Francis Group, London, ISBN 978-0-415-44056-1

the width of the pore-size distribution (or the steepness of the retention function); α [L-1], n [-]and m=1/1/n [-] = empirical shape parameters; hm = median pressure head [L] for whichSe(hm)=0.5 and corresponding to the median pore radius; σ = standard deviation of the logtransformed soil pore radius [-] characterizing the width of the pore-size distribution; and wi =weighting factors for the overlapping regions.
Equation (1) contains the hydraulic conductivity function which characterizes the ability of asoil to transmit water, and as such is inversely related to the resistance to water flow. Similarly
as for the soil water retention curve, HYDRUS allows analytical models for the hydraulic conductivity function as suggested by Brooks & Corey (1964), van Genuchten (1980), Kosugi(1996), and Durner (1994):⎧
Se2 / λ + l + 1⎪
2
Sel [1- (1- Se1/ m ) m ]
⎪
⎪
2
⎪⎪
⎧
⎡ ln ( h / hm )
⎤ ⎫⎪
l ⎪1
Se ⎨ erfc ⎢+ σ ⎥⎬K ( Se ) = K s ⎨
2
2
σ
⎪
⎣
⎦ ⎪⎭
⎩
⎪
2
⎪
l
1 /m m1
1 /m m2⎪ w1Se1 + w2 Se2 w1α1[1 - (1 - S e1 1 ) ] + w2α 2 [1 - (1 - S e2 2 ) ]⎪2( w1α1 + w2α 2 )
⎪⎩

(

)(

(3)

)

respectively, where Ks = saturated hydraulic conductivity [LT-1]; and l = pore-connectivity parameter [-].
2.2 Nonequilibrium flow
The Richards equation (1) predicts a uniform flow process in the vadose zone, although possibly
modified macroscopically by allowing the use of spatially variable soil hydraulic properties(e.g., as dictated by different soil horizons, or varying laterally). Unfortunately, the vadose zone
can be extremely heterogeneous at a range of scales, from the microscopic (e.g., pore scale) tothe macroscopic (e.g., field scale). Some of these heterogeneities can lead to a preferential flow
process that is very difficult to capture macroscopically with the standard Richards equation.
One obvious example of preferential flow is the rapid movement of water and dissolved solutes
through macropores (e.g., between soil aggregates, or caused by earthworm pathways or decayed root channels) or within rock fractures, with much of the water bypassing (shortcircuiting) the soil or rock matrix. However, many other causes of preferential flow exist, suchas flow instabilities caused by soil textural changes or water repellency (Hendrickx & Flury
2001; Šimůnek et al. 2003), as well as lateral funneling of water due to inclined or other textural
boundaries (e.g., Kung 1990).
Preferential flow in macroporous soils and fractured rocks can be described using a variety of
dual-porosity or dual-permeability models (Gerke & van Genuchten 1993; Šimůnek et al. 2003).
Dual-porosity and dual-permeability models both assume that the porous medium consists of
two interacting pore regions, one associated with the inter-aggregate, macropore, or fracture system, and one comprising the micropores (or intra-aggregate pores) inside soil aggregates or therock matrix. Dual-porosity models assume that water in the matrix is stagnant, while dualpermeability models allow also for water flow within the soil or rock matrix.
2.2.1 Dual-porosity model
The dual-porosity models implemented in the HYDRUS-1D and HYDRUS (2D/3D) packages
assume that water flow is restricted to the macropores (or the inter-aggregate pores and fractures), and that water in the matrix (the intra-aggregate pores or the rock matrix) does not moveat all. This conceptualization partitions the liquid phase into mobile (flowing, inter-aggregate), θmo,and immobile (stagnant, intra-aggregate), θim, regions [L3L-3]:

θ = θ mo + θim

(4)

97

Numerical Modelling of Hydrodynamics for Water Resources – Garcia-Navarro & Playán (eds)© 2008 Taylor & Francis Group, London, ISBN 978-0-415-44056-1

The dual-porosity formulation for water flow is then based on a mixed formulation of the
Richards equation (1) to describe water flow in the macropores (the preferential flow pathways)
augmented with a mass balance equation to describe the moisture dynamics in the matrix as follows (Šimůnek et al. 2003):⎛⎞⎤∂θ mo (hmo )∂ ⎡∂h=⎢ K ( hmo ) ⎜⎜ K ijA mo + K izA ⎟⎟ ⎥ − Smo ( hmo ) − Γ w∂t∂xi ⎣⎢∂x j⎝⎠ ⎦⎥∂θim ( him )
=
∂t

(5)

− Sim ( him ) + Γ w

where him and hmo = pressure heads for both regions [T-1]; Sim and Smo = sink terms for both regions [T-1]; and Γw = transfer rate for water from the inter- to the intra-aggregate pores [T-1].
Examples of applications of dual-porosity models to a range of laboratory and field data involving transient flow and solute transport are given by Šimůnek et al. (2001), Castiglione et al.
(2003), Zhang et al. (2004), Köhne et al. (2004a, 2005), Kodešová et al. (2005) and Haws et al.
(2005).
2.2.2 Dual-permeability modelResearch versions of HYDRUS-1D and the two-dimensional computational module of HYDRUS (2D/3D) implement the dual-permeability approach suggested by Gerke & van Genuchten (1993). This approach applies Richards equations to each of two pore regions. The flow
equations for the macropore or fracture (subscript f) and matrix (subscript m) pore systems inthis approach are given by:∂θ f ( h f )
∂t

=

⎛
⎞⎤
∂h f
∂ ⎡
+ K izA ⎟ ⎥ - S f ( h f ) - Γ w⎢ K f ( h f ) ⎜⎜ K ijA⎟w∂xi ⎢⎣
∂x j
⎝
⎠ ⎥⎦

⎛
⎞⎤
∂θ m ( hm )
∂ ⎡
∂h
Γw
=
⎢ K m ( hm ) ⎜⎜ K ijA m + K izA ⎟⎟ ⎥ - Sm ( hm ) +1- w∂t∂xi ⎣⎢
∂x j
⎝
⎠ ⎦⎥

(6)

respectively, where w = ratio of the volumes of the macropore or fracture domain and the totalsoil system [-]. Note that the water contents θf and θm in (1) have different meanings than in (5)where they represented water contents of the total pore space (i.e., θ = θmo + θim), while in (6)they refer to water contents of the two separate (fracture or matrix) pore domains such that θ =wθf + (1-w)θm.
Examples of applications of dual-permeability models to a range of laboratory and field data
involving transient flow and solute transport are given by Pot et al. (2005), Kodešová et al.
(2005), and Köhne et al. (2006).
2.3 Water flow and vapor transport
The Richards equation (1) considers only water flow in the liquid phase and ignores the effectsof the vapor phase on the overall water mass balance. While this assumption is justified for the
majority of applications, a number of problems exist in which the effect of vapor flow can not
be neglected. Vapor movement is often an important part of the total water flux when the soil
moisture becomes relatively low. Scanlon et al. (2003) showed that water fluxes in deep vadosezone profiles of the arid and semiarid regions of the western U.S. are often dominated by thermal vapor fluxes. Nonisothermal liquid and vapor flow in an unpublished research version of
HYDRUS is described using the equation (e.g., Saito et al. 2006):⎛ A ∂h⎞∂θ ( h )∂ ⎡∂T ⎤=+ KizA ⎟ + ( K LT + K vT ) KijA⎢( K + K vh ) ⎜⎜ Kij⎥ - S (h )⎟∂t∂xi ⎢⎣∂x j∂x j ⎥⎦
⎝
⎠

98

(7)

Numerical Modelling of Hydrodynamics for Water Resources – Garcia-Navarro & Playán (eds)© 2008 Taylor & Francis Group, London, ISBN 978-0-415-44056-1

where θ = total volumetric water content, being the sum (θ=θl +θv) of the volumetric liquid water content, θl, and the volumetric water vapor content, θv (both expressed as an equivalent watercontent), T = temperature [K]; K = isothermal hydraulic conductivity for the liquid phase [LT-1];
KLT = thermal hydraulic conductivity for the liquid phase [L2K-1T-1]; Kvh = isothermal vapor hydraulic conductivity [LT-1]; KvT = thermal vapor hydraulic conductivity [L2K-1T-1]. Overall waterflow in (7) is given as the sum of isothermal liquid flow, isothermal vapor flow, gravitationalliquid flow, thermal liquid flow, and thermal vapor flow. Since several terms of (7) are a function of temperature, this equation should be solved simultaneously with the heat transport equation (11) to properly account for temporal and spatial changes in soil temperature.
Examples of HYDRUS-1D applications that consider the simultaneous movement of water inthe liquid and vapor phases are given by Scanlon et al. (2003) and Saito et al. (2006). Scanlon etal. (2003) evaluated response of several semiarid (High Plains, Texas) and arid (Chihuahuan
Desert, Texas; Amargosa Desert, Nevada) sites to paleoclimatic forcing using water-potential
and Cl profiles, and by modeling nonisothermal liquid and vapor flow and Cl transport. They
used HYDRUS-1D to show that upward nonisothermal liquid and vapor flow for at least 1 to 2
kyr in the High Plains and for 12 to 16 kyr in the Chihuahuan and Amargosa desert sites is required to reproduce measured upward water potentials and Cl profiles. Saito et al. (2006) additionally combined a numerical solution of the water flow (7) and heat transport (11) equations
with water and energy balances at the soil surface. Their analyses of the distributions of the liquid and vapor fluxes versus depth showed that the soil water dynamics close to the soil surface
is strongly associated with the soil temperature regime.
2.4 Water flow with freezing and thawing
Although considering both liquid and vapor flow, equation (7) still neglects one phenomenon
that is important in cold climatic regions, i.e., freezing and thawing. To include freezing andthawing processes, (7) needs an additional term that accounts for frozen water (e.g., Hansson etal. 2004) as follows:∂θ ( h ) ρ i ∂θi ( T )+=∂tρ w ∂t⎛ A ∂h⎞∂ ⎡∂T ⎤=+ KizA ⎟ + ( K LT + K vT ) KijA⎢( K + K vh ) ⎜⎜ Kij⎥ - S (h)⎟∂xi ⎢⎣
∂
∂
x
x
j
j
⎥⎦
⎝
⎠

(8)

where θi = volumetric ice content [L3L-3]; ρi = density of ice [ML-3] (931 kg m-3); and ρw = density of liquid water [ML-3] (~1000 kg m-3). This equation is even more tightly coupled with theheat transport equation (12) through the second term on the left side, since the change in thevolumetric ice content depends on temperature and can be expressed using the generalized
Clapeyron equation.
2.5 Overland flow
Traditionally, surface and subsurface flow processes have been considered mostly separately,with some of the processes being severely simplified when coupled surface and subsurface flow
is simulated. For example, existing surface irrigation models typically still use empirical infiltration functions such as the Philip, Kostiakov, modified Kostiakov, and Branch infiltration equations, rather than rigorously simulating subsurface water flow and solute transport using the
Richards equation for variably-saturated flow, and advection-dispersion type equations for solute transport. To overcome this deficiency, the following overland flow equation has been incorporated in HYDRUS-2D for use in combination with the atmospheric boundary condition at
the surface boundary (Šimůnek 2003):

99

Numerical Modelling of Hydrodynamics for Water Resources – Garcia-Navarro & Playán (eds)© 2008 Taylor & Francis Group, London, ISBN 978-0-415-44056-1∂h ∂Q+= q( x , t )∂t ∂x
Q = αhm

S 1/ 2
α=
n

where

(9)

m = 5/3

andwhere h = unit storage of water (or mean depth for smooth surfaces) [L]; Q = discharge per unit
width [L2T-1], x = distance coordinate over the soil surface [L]; q(x,t) = rate of local input, or lateral inflow (i.e., local precipitation minus local infiltration) [LT-1]; α [L2-mT-1] and m [-] = parameters related to slope, surface roughness, and the flow conditions (laminar or turbulent flow),usually evaluated using Manning’s hydraulic resistance law; S = slope [-]; and n = Manning’s
roughness coefficient for overland flow [-].
Figure 1 shows a simple application of the overland flow module in HYDRUS-2D to simulate
flow in a hillside transect in which the soil material in the middle third of the transect (between
33 and 66 m) has a 100 times higher saturated hydraulic conductivity than the remaining soil (25versus 0.25 m d-1). Because of this, the middle section can accommodate the infiltration of waterfrom both the rainfall itself and from runoff from the upgradient part of the hillslope. The soiltransect had a slope of 0.01, while the roughness coefficient n was assumed to be equal to 0.01.

Overland flow
Infiltration

0.5
100 m
0.24

0.26

0.28

0.30

0.32

0.34

0.36

0.38

0.40

0.42

0.44

Figure 1. Calculated contours of the water content in a hillsidetransect following rainfall and overland flow (Šimůnek 2003).

2.6 HYDRUS package for MODFLOW
Although computer power has increased tremendously during the last few decades, large scale
three-dimensional applications evaluating water flow in the vadose zone are often still prohibitively expensive in terms of computational resources. To overcome this problem, Seo et al.
(2006, 2007) developed a computationally efficient one-dimensional unsaturated flow HYDRUS package and linked it to the three-dimensional modular finite-difference ground watermodel MODFLOW-2000 (Harbaugh et al. 2000). The HYDRUS unsaturated flow package used
HYDRUS-1D to simulate one-dimensional vertical variably-saturated flow. MODLOW zone arrays were used to define the cells to which the HYDRUS package was applied. MODFLOW
used the time-averaged flux from the bottom of the unsaturated zone as recharge, and calculated
a water table depth which was then used as a pressure head bottom boundary for HYDRUS.
Twarakavi et al. (in press) provided a comparison of the HYDRUS package to other MODFLOW packages that evaluate processes in the vadose zone and presented a field application
demonstrating the functionality of the package.

100

Numerical Modelling of Hydrodynamics for Water Resources – Garcia-Navarro & Playán (eds)© 2008 Taylor & Francis Group, London, ISBN 978-0-415-44056-1

3 HEAT TRANSPORT
3.1 Heat transport without vapor transport
All HYDRUS models consider heat transport due to conduction and convection with flowing water. Neglecting the effects of water vapor diffusion, the heat transport can be described as:

C (θ )

∂T
∂ ⎛
∂T ⎞
∂T
=
⎜⎜ λij (θ )⎟⎟ - Cw qi∂t ∂xi ⎝∂x j ⎠∂xi(10)

where λij(θ) = apparent thermal conductivity of the soil [MLT-3K-1] (e.g. Wm-1K-1); C(θ) and Cw =volumetric heat capacities [ML-1T-2K-1] (e.g. Jm-3K-1) of the porous medium and the liquid phase, respectively; and qi = fluid flux density [LT-1]. The HYDRUS models at each time step first numerically solve the Richards equation (1) and then the heat transport equation (10) using calculated values of the water content and the fluid flux.
3.2 Heat transport with vapor transport
When the effects of water vapor diffusion can not be neglected, the heat transport must be expanded
to the form (e.g., Saito et al. 2006):

C (θ )

∂T
∂θ
∂ ⎛
∂T ⎞
∂T
∂q T
∂q
+ L0 v =− Cv vi − L0 vi⎜⎜ λij (θ )⎟⎟ - Cw qi
∂t
∂t
∂xi ⎝
∂x j ⎠
∂xi
∂xi
∂xi
10

10
T=0
t=0.25 d
t=1
t=5
t=25

Depth [cm]

8

6

T=0
t=0.25 d
t=1
t=5
t=25

8

6

4

4

2

2

a

b

0

0
0

0.05

0.1

0.15

0.2

0

0.02

Water Content [-]
10

0.04

0.06

Total Flux [cm/d]
10

c

d

T=0

8

t=0.25 d8

t=1

Depth [cm]

t=5
6

6

4

4

T=0
t=0.25 d
t=1
t=5
t=25

2

t=25

2

0

0
0

10

20

Temperature [C]

30

0

2

4

6

8

10

Concentration [-]Figure 2. Water content (a), total flux (b), temperature (c), and solute concentration(d) distributions in a 10-cm long vertical soil sample with zero water fluxes acrossthe top and bottom boundaries, and with temperature increasing from top to bottom.

101

(11)

Numerical Modelling of Hydrodynamics for Water Resources – Garcia-Navarro & Playán (eds)© 2008 Taylor & Francis Group, London, ISBN 978-0-415-44056-1

where L0 = volumetric latent heat of vaporization of liquid water [ML-1T-2] (e.g., Jm-3); and qvi =vapor flux density [LT-1]. In equation (11), the total heat flux density is defined as the sum of the
conduction of sensible heat as described by Fourier’s law (the first term on the right side), sensible heat by convection of liquid water (the second term) and water vapor (the third term), andof latent heat by vapor flow (the forth term) (de Vries 1958).
As an example, Figure 2 shows calculated water content, total flux, temperature and concentration profiles for a 10-cm long soil sample with zero water fluxes at both top and bottom
boundaries, and with a temperature gradient along the sample. Increasing temperatures (Fig. 2c)from the top to the bottom of the sample cause vapor flow (Fig. 2b) from the warmer bottom endof the sample toward the colder end. Water evaporates at the warmer end, flows upward as vapor and condensates at the colder end. Water contents correspondingly decrease at the warmer
end and increase at the colder bottom (Fig. 2a). As a consequence of changing water contents, apressure head gradient develops in the sample, leading to water flow in a direction opposite tovapor flow. A steady-state is eventually reached when the upward vapor flow fully balances thedownward liquid flow (Fig. 2b). Since water evaporates at the bottom of the sample and condensates at the top, solute becomes more concentrated near the bottom and more diluted near the top(Fig. 2d). Also, the concentration profile should eventually reach steady-state, although at amuch larger time, when the advective downward solute flux balances the upward diffusive flux.
3.3 Heat transport with freezing and thawing
Similar to water flow, considering freezing and thawing in the heat transport equation requires
an additional term to account for the energy stored in the frozen water (e.g., Hansson et al.
2004):C (θ )∂T∂θ∂θ
− L f ρ i i + L0 v =
∂t
∂t
∂t
∂ ⎛
∂T ⎞
∂T
∂q T
∂q
=
− Cv vi − L0 vi⎜⎜ λij (θ )⎟⎟ - Cw qi∂xi ⎝∂x j ⎠
∂xi
∂xi
∂xi

(12)

where Lf = latent heat of freezing [L2T-2] (~3.34e+5 J kg-1). The first term on the left-hand side
represents changes in the energy content, whereas the second and third terms represent changes
in the latent heat of the frozen and vapor phases, respectively. Terms on the right-hand side represent respectively, soil heat flow by conduction, convection of sensible heat with flowing water,transfer of sensible heat by diffusion of water vapor, and transfer of latent heat by diffusion of water
vapor. Hansson et al. (2004) discussed the extreme nonlinearity of (12) caused by the dependency of the apparent volumetric heat capacity ( = C − L f ρ i dθi / dT ) on temperature. While theapparent heat capacity increases by about two orders of magnitude for silty clay when the freezing point is reached, the increase for sand is almost five orders of magnitude (Fig. 3). For sands
the increase in the apparent heat capacity due to freezing becomes negligible below about -0.05
oC, which corresponds to a pressure head of -6.2 m. At that temperature almost all soil water isfrozen, except for the residual water content. For fine-textured soils the increase in the apparent
heat capacity extends to much lower temperatures, which reflects the fact that for these soils asignificant amount of water remains unfrozen at slightly subzero temperatures.
Hansson et al. (2004) used HYDRUS-1D to evaluate data from laboratory column freezing
experiments. The experiments involved 20-cm long soil columns with an internal diameter of 8
cm that were exposed at the top to a circulating fluid with a temperature of –6˚C. Water and soilin the columns froze from the top down during the experiment, with the freezing process inducing significant water redistribution within the soil. Values of the total water content calculatedwith HYDRUS-1D compared well with measured values.
Figure 4 shows temperatures, pressure heads, liquid water contents, and concentrations at different depths in a loamy soil profile subject to fluctuating surface temperatures. Surface temperatures followed the sinusoidal curve with a mean value of -1oC and a daily amplitude of 5oC.
Changes in temperature can be observed down to a depth of about 10 cm (Fig. 4a). The Clapeyron equation predicts that a temperature decrease of 1oC below zero corresponds to a pressure
head decrease of about -120 m. This shows that negative temperatures can cause extremely102Numerical Modelling of Hydrodynamics for Water Resources – Garcia-Navarro & Playán (eds)© 2008 Taylor & Francis Group, London, ISBN 978-0-415-44056-1

negative pressure heads in the soil profile (Fig. 4b). These pressure heads have their corresponding water contents (Fig. 4c), which follow immediately from the soil water retention curve. The
remaining water within the soil profile stays in a frozen state. From Figure 4 it is obvious thatthe process of freezing is very similar to drying, with frozen water remaining in the larger pores
and liquid water in the smaller pores. Concentrations fluctuate in a similar fashion as the othervariables since solutes are excluded from frozen water and remain in the liquid phase (Fig. 4d).

Figure 3. Apparent volumetric heat capacity (Jm-3K-1) for three soil texturalclasses (sand, loam, and silty clay).
4

0

3

Pressure Head [cm]Temperature [C]

2
1
0
-1
-2
-3

0 cm
1
3.5
10

-4
-5

a

-6

-0.5

0

0.5

1

1.5

-20000-40000
0 cm
1
3.5
10

-60000b-80000
-0.5

2

0

0.5

Time [d]
0.38

2

1.5

2

0 cm
1
3.5
10

1.5

Concentration [-]0.34

Water Content [-]

1.5

1.6

0.36

0.32
0.3
0.28
0.26
0 cm
1
3.5
10

0.24
0.22
0.2

1

Time [d]

c

-0.5

1.4
1.3
1.2
1.1
1

0

0.5

1

1.5

2

d

-0.5

Time [d]

0

0.5

1

Time [d]Figure 4. Temperatures (a), pressure heads (b), water contents (c), and concentrations (d) at different
depths in a loamy soil profile subject to fluctuating surface temperatures.

103

Numerical Modelling of Hydrodynamics for Water Resources – Garcia-Navarro & Playán (eds)© 2008 Taylor & Francis Group, London, ISBN 978-0-415-44056-1

4 SOLUTE TRANSPORT
While water flow in the various HYDRUS software packages is simulated in a relatively similarmanner (i.e., using the Richards equation to describe uniform flow, or using dual-porosity or
dual-permeability models for nonequilibrium flow), different solute transport modules exist in
the different codes. Some modules simulate single-component transport, while others containbiogeochemical modules for simulating potentially very complex interactions between the various solutes. In this section we provide brief descriptions of the various solute transport modules.
4.1 Single component solute transport
Most HYDRUS predecessors (e.g., SWMS_2D and SWMS_3D) allowed solutes to exist only inthe liquid and solid phases. The newer windows-based HYDRUS models (i.e., HYDRUS-1D, HYDRUS-2D and HYDRUS (2D/3D)) assume that solutes can exist in all three phases (liquid, solid,
and gaseous) and that the decay and production processes can be different in each phase. While in
the HYDRUS predecessors solutes could be transported only by advection and dispersion in the liquid phase, the more recent versions assume that they can be transported also by diffusion in the gas
phase. In the various solute transport formulations below we will use the following general conservation equation for solutes in a variably saturated rigid porous medium∂M∂J= − Si − φ
∂t
∂xi

(13)

where M = mass of solute [ML-3]; JSi = solute flux density [ML-2T-1]; and φ = coefficient representing various sinks and sources [ML-3T-1].
4.1.1 Equilibrium transport
The SWMS_2D and SWMS_3D models (Šimůnek et al. 1994, 1995) assume that solutes reside
only in the liquid and solid phases. The various terms of (13) are then given by:

M = θc + ρs
J Si = −θ Dijw∂c+ qi c∂x j(14)

φ = − μwθ c - μs ρ s + γ wθ + γ s ρ - Scr
where c, s = solute concentrations in the liquid [ML-3] and solid [MM-1] phases, respectively; μw, μs= first-order rate constants for solutes in the liquid and solid phases [T-1], respectively; γw, γs = zeroorder rate constants for the liquid [ML-3T-1] and solid [T-1] phases, respectively; ρ = soil bulk density[ML-3]; S = sink term in the water flow equation (1); cr = concentration of the sink term [ML-3]; Dijw
= dispersion coefficient tensor [L2T-1] for the liquid phase.
The HYDRUS models allow solutes to exist also in the gaseous phase, as well as permit them
to be involved in sequential first-order decay chain reactions. The different terms of (13) then become:M k = θ ck + ρ sk + av g k∂c∂g
J Si ,k = −θ Dijw,k k − av Dijg,k k + qi ck
∂x j
∂x j

k = 1, N s(15)

φk = − ( μw,k + μw' ,k ) θ ck - ( μs ,k + μs' ,k ) ρ sk − ( μ g ,k + μ g' ,k ) av g k - Scr ,k ++ ϑ ( μw' ,k −1θ ck −1 + μ s' ,k −1ρ sk −1 + μ g' ,k −1 a vg k −1 ) + γ w,kθ + γ s ,k ρ + γ g ,k av
where g = solute concentration in the gaseous phase [ML-3]; μg = first-order rate constant forsolutes in the gas phase [T-1]; μw’, μs’, and μg’ = similar first-order rate constants providing con-104

Numerical Modelling of Hydrodynamics for Water Resources – Garcia-Navarro & Playán (eds)© 2008 Taylor & Francis Group, London, ISBN 978-0-415-44056-1

nections between individual species in the decay chain; γg = zero-order rate constant for the gas
phase [ML-3T-1]; av = air content [L3L-3]; Dijg = diffusion coefficient tensor [L2T-1] for the gas
phase; and Ns = number of solutes involved in the chain reaction. The subscripts w, s, and g correspond with the liquid, solid and gas phases, respectively; while the subscript k represents thekth chain number. The parameter ϑ in (15) is zero for k=1. The indicial notation used above assumes summations over indices i and j (i,j = 1,2,3), but not over index k.
Interactions between the solid and liquid phases in the HYDRUS models may be describedusing nonlinear equations, while interactions between the liquid and gaseous phases are assumed
to be linear and instantaneous, i.e.,

s=

Kd cβ
1 + ηc β(16)

g = kgc(17)

respectively, where Kd, β, η, kg = empirical coefficients [L3M-1], [-], [L3M-1], [-], respectively. The
Freundlich, Langmuir, and linear adsorption equations are special cases of equation (16). When β=1,equation (16) reduced to the Langmuir equation, when η=0 the equation becomes a Freundlich isotherm, and when β=1 and η=0, equation (16) defines a linear adsorption isotherm. Solute transport
without adsorption is described with Kd=0. We note here that the SWMS_2D and SWMS_3D models considered only linear sorption (i.e., β=1 and η=0).
Typical examples of HYDRUS applications involving sequential first-order decay chains are
the transport of nitrogen species (e.g., Hanson et al. 2006), pesticides, chlorinated aliphatic hydrocarbons (Schaerlaekens et al. 1999; Casey & Šimůnek 2001), hormones (Casey et al. 2003),
and explosives (Dontsova et al. 2006). The ability of HYDRUS to account for diffusion in thegas phase was further used by Wang et al. (1997, 2000) to simulated the transport of volatilecontaminants (fumigants) such as methyl bromine and 1,3-dichloropropene.
4.1.2 Nonequilibrium transport
HYDRUS-1D and the two-dimensional module of HYDRUS (2D/3D) offer a large number of
options to deal with nonequilibrium flow and transport. A complete list of these options is givenby Šimůnek et al. (in press) who divided them into three groups: a) physical nonequilibriumtransport models, b) chemical nonequilibrium transport models, and c) physical and chemical
nonequilibrium transport models. Physical nonequilibrium models include 1) the MobileImmobile Water Model, 2) the Dual-Porosity Model, 3) the Dual-Permeability Model, and 4)
the Dual-Permeability Models with Immobile Water. Chemical nonequilibrium models include
1) the One Kinetic Site Model, 2) the Two-Site Model, and 3) the Two-Kinetic Sites Model. Finally, physical and chemical nonequilibrium transport models include 1) the Dual-Porosity
Model with one Kinetic Site, and 2) the Dual-Permeability Model with Two-Site Sorption.
An example of physical nonequilibrium transport as implemented in HYDRUS is the Mobile-Immobile Water Model that uses the concept of two-region, dual-porosity type solute transport (van Genuchten & Wierenga 1976). Solute exchange between the two liquid regions is thenmodeled as a first-order process:∂θim cim∂s+ (1 − f mo ) ρ im = Γ s - φim∂t∂tΓ s = α ( cmo − cim ) + Γ w c *(18)where cmo and cim = concentrations of the mobile and immobile regions [ML-3], respectively; sim= sorbed concentration of the the immobile region [-]; φim = sink/source term that accounts forvarious zero- and first-order or other reactions in immobile region [ML-3T-1]; fmo = fraction ofsorption sites in contact with the mobile water content [-]; α = mass transfer coefficient [T-1]; andΓs = mass transfer term for solutes between the mobile and immobile regions [ML-3T-1]. The
first equation of (18) is a mass balance for the immobile (micropore) domain, while the second
equation (Γs ) describes the rate of mass transfer between the mobile and immobile domains.
The second (advective) term of Γs in Eq. (18) is equal to zero for the Mobile-Immobile Model
since that model does not consider water flow between the two regions. In the Dual-Porosity

105

Numerical Modelling of Hydrodynamics for Water Resources – Garcia-Navarro & Playán (eds)© 2008 Taylor & Francis Group, London, ISBN 978-0-415-44056-1

Model, c* is equal to cmo for Γw>0 and cim for Γw<0. Note that although the HYDRUS models do
consider diffusive transport in the gas phase also for nonequilibrium applications, for simplicity
we did not include this process in (18) and subsequent equations.
Still much more sophisticated models for physical nonequilibrium transport may be formulated. For example, Pot et al. (2005), Köhne et al. (2006), and Šimůnek et al. (in press) considered not only a dual-permeability model with flow and transport in both regions, but also a dualpermeability model that additionally divides the matrix domain into mobile and immobile
subregions. Pot et al. (2005) and Köhne et al. (2006) used the HYDRUS-1D model with thesedual-permeability submodels to successfully simulate bromide transport in laboratory soil columns at different flow rates or for transient flow conditions, respectively.
An example of chemical nonequilibrium transport in HYDRUS is the Two-Site Model thatassumes that the sorption sites can be divided into two fractions:s = se + sk
e

(19)
-1

where sorption, s [MM ], on one fraction of the sites (often referred to as type-1 sites) is assumed
to be instantaneous, while sorption, sk [MM-1], on the remaining (type-2) sites is considered to betime-dependent. Sorption on the type-2 nonequilibrium sites is assumed to be a first-order kinetic
rate process:∂s k= α k ρ ⎡⎣ sek - s k ⎤⎦ - φk
ρ
∂t

(20)

where αk = first-order rate constant describing the kinetics of the sorption process [T-1]; sek =sorbed concentration that would be reached at equilibrium with the liquid phase concentration[MM-1]; φk = sink/source term that accounts for various zero- and first-order or other reactions at
the kinetic sorption sites [ML-3T-1].
Although we listed the Two-Kinetic Sites Model above among the chemical nonequilibrium
models, this model could equally well be viewed as a physical and chemical nonequilibriumtransport model. This because both physical and chemical nonequilibrium processes could be responsible for solutes associated with each of the two sites. For example, while the chemical
process of attachment/detachment could be associated with one fraction of sorption sites, thephysical process of straining could be associated with the second fraction of sorption sites. In
the Two-Kinetic Sites Model the sorption sites are again divided into two fractions:s = s1k + s2k(21)

but now both sorption sites are assumed to be kinetic. In (21), s1k and s2k = concentrations of thefirst and second fraction of kinetic sorption sites [MM-1], respectively. The kinetic processes arethen described as follows:∂s1k= ka1θ c − kd 1ρ s1k - φk 1∂t∂s k
ρ 2 = ka 2θ c − kd 2 ρ s2k - φk 2
∂t

ρ

(22)

where ka1, ka2 = attachment coefficients for the first and second fraction of kinetic sorption sites[T-1], respectively; kd1, kd2 = detachment coefficients for the first and second fraction of kinetic
sorption sites [T-1], respectively; and φk1, φk2 = sink/source terms for the first and second fractionof kinetic sorption sites [T-1], respectively. As mentioned above, the two kinetic sites may be
used to describe different processes. While the first kinetic process could be used for chemical
attachment, the second kinetic process could represent physical straining of colloids or pathogenic microorganisms (e.g., Bradford et al. 2004, 2007; Gargiulo et al. 2007). Note that in (22)
we do not give the nonlinear blocking coefficients to account for, for example, Langmuirian
blocking to attachment or depth-dependent straining that are considered in the HYDRUS models(e.g., Bradford et al. 2004).
Details about the remaining nonequilibrium models, as well as of various applications, can befound in the HYDRUS manuals as well as in Šimůnek et al. (in press).

106

Numerical Modelling of Hydrodynamics for Water Resources – Garcia-Navarro & Playán (eds)© 2008 Taylor & Francis Group, London, ISBN 978-0-415-44056-1

4.2 Multicomponent Solute Transport
The transport of reactive contaminants in the subsurface is generally affected by a large numberof nonlinear and often interactive physical, chemical and biological processes. Proper simulatingthese processes requires a reactive transport code that couples the physical processes of waterflow and advective-dispersive transport with a range of biogeochemical processes. The soil solution is always a mixture of many ions which may be involved in mutually dependent chemical
processes, such as complexation reactions, cation exchange, precipitation-dissolution, sorptiondesorption, volatilization, redox reactions, and degradation (Šimůnek & Valocchi 2002; Appelo& Postma 2005). The transport and transformation of many contaminants is further mediated bysubsurface aerobic or anaerobic bacteria. Bacteria for example can catalyze redox reactions inwhich organic compounds (e.g., chlorinated hydrocarbons) act as the electron donor and inorganic substances (oxygen, nitrate, sulfate, or metal oxides) as the electron acceptor.
To be able to evaluate these type of systems, several different multicomponent modules weredeveloped for the HYDRUS models. In addition to the single component solute transport module described above, HYDRUS-1D also includes the UNSATCHEM module (Šimůnek &Suarez 1993, 1994) that considers the transport of major ions, and their mutual reactions, suchas dissolution-precipitation, cation exchange, or aqueous complexation. HYDRUS-1D has also
been coupled with PHREEQC, resulting in the HP1 model (Jacques & Šimůnek 2005; Jacqueset al. 2006) that accounts for a wide range of instantaneous and kinetic chemical and biologicalreactions, including aqueous complexation, precipitation-dissolution, cation exchange, surface
complexation, and redox reactions. Finally, the CW2D wetland module (Langergraber & Šimůnek 2005, 2006) in HYDRUS (2D/3D) considers biochemical processes associated with the flowof wastewater through constructed or natural wetlands. Below we briefly review these models
and indicate some of their applications.
Šimůnek & Valocchi (2002) divided geochemical transport models into two major groups:
those with specific chemistry and more general models. Models with specific chemistry are generally restricted to certain prescribed chemical systems and thus are usually constrained to specificapplications. The major ion chemistry UNSATCHEM module in the HYDRUS-1D software package and the CW2D wetland module in HYDRUS (2D/3D) are typical examples of models withspecific chemistry. On the other hand, models with generalized chemistry provide users with much
more freedom in simulating chemical systems, and thus permit a much broader range of applications. Users then either can select species and reactions from large geochemical databases, or areable to define their own species with particular chemical properties and reactions. The HP1 model
represents a typical example of a model with generalized chemistry.
4.2.1 Major ion chemistry
As an alternative to the single-component approach, a major ion chemistry module based on the
UNSATCHEM model (Šimůnek & Suarez 1993, 1994) was incorporated in HYDRUS-1D. This
module considers the transport of major ions and carbon dioxide in soils. CO2 transport is assumed to be governed by diffusion in both the liquid and gas phases, by advection in the liquid
phase, and by respiration by both soil microorganisms and plant roots (Šimůnek & Suarez1993):∂ ( ca av + cwθ ) ∂ ⎛∂c ⎞ ∂ ⎛∂c ⎞ ∂q c ∂qcw
= ⎜ av Da a ⎟ + ⎜ θ Dw w ⎟ - a a - Scw + P∂t∂z ⎝∂z ⎠ ∂z ⎝∂z ⎠ ∂z∂z
P = γ s +γ p(23)where cw, ca = volumetric concentrations of CO2 in the dissolved and gas phases [L3L-3], respectively; Da = effective soil matrix diffusion coefficient of CO2 in the gas phase [L2T-1]; Dw = effective soil matrix dispersion coefficient of CO2 in the dissolved phase [L2T-1]; qa = soil air flux
[LT-1]; q = soil water flux [LT-1]; Scw = dissolved CO2 removed from the soil by root water uptake [L3L-3T-1]; P = CO2 production/sink term [L3L-3T-1]; and γs, γp = CO2 production rates bysoil microorganisms or plant roots [L3L-2T-1], respectively (Šimůnek and Suarez, 1993).
The major variables of the UNSATCHEM chemical system are Ca, Mg, Na, K, SO4, Cl, NO3,
H4SiO4, alkalinity, and CO2. The model accounts for equilibrium chemical reactions between
these components such as complexation, cation exchange and precipitation-dissolution, e.g.,107

Numerical Modelling of Hydrodynamics for Water Resources – Garcia-Navarro & Playán (eds)© 2008 Taylor & Francis Group, London, ISBN 978-0-415-44056-1( 2+ ) (CO32- )K CaCO3 = Ca
(CaCO3o)K Mg,Ca =2+(Ca 2+ )1/ 2
2+
2+ 1/ 2Ca (Mg )

Mg

(24)

C
K SP
= (Ca 2+ ) (CO32− )respectively, where the parentheses denote ion activities, KCaCO3 = equilibrium constant of the
CaCO30 complexed species [-]; KMg,Ca = Gapon selectivity coefficient for Mg2+ and Ca2+ [-]; andKSPC = calcite solubility product [-]. Similar algebraic equations hold for all species of thechemical system.
For the precipitation-dissolution of calcite (CaCO3) and the dissolution of dolomite(MgCO3⋅3H2O), either equilibrium or multicomponent kinetic expressions can be used, including both forward and backward reactions. Other precipitation-dissolution reactions considered inhydromagnesite
the
UNSATTCHEM
modelsinvolve
gypsum(CaSO4⋅2H2O),(Mg5(CO3)4(OH)2⋅4H2O), nesquehonite (MgCO3⋅3H2O), and sepiolite (Mg2Si3O7.5(OH)⋅3H2O).
Since the ionic strength of soil solutions can vary considerably in time and space and often reachhigh values, both the modified Debye-Hückel (Truesdell & Jones 1974) and Pitzer (Pitzer 1979)
expressions were incorporated into the model, thus providing options for calculating single-ion
activities. The module also considers the effects of solution composition on the unsaturated soil
hydraulic properties. The new UNSATCHEM module of HYDRUS-1D enables quantitative
predictions of processes involving major ions, such as simulations of the effects of salinity on
root water uptake and plant growth, evaluation of alternative irrigation, salinity and crop management practices, evaluation of water suitability for irrigation, and estimation of the amount ofwater and amendment required to reclaim soil profiles to desired levels of salinity and ESP (exchangeable sodium percentage).
A recent application of HYDRUS-1D and its major ion chemistry module is given by Gonçalves et al. (2005) who simulated solute transport in three lysimeters irrigated with different
quality waters over a time period of three years. HYDRUS-1D successfully described fieldmeasurements of not only the overall salinity, but also of individual soluble cations as well as
SAR (Sodium Adsorption Ratio) (Fig. 5) and ESP (Gonçalves et al. 2005).

Figure 5. Measured and simulated sodium adsorption ratios (SAR) for lysimeters A, B, and C(Gonçalves et al. 2005).

4.2.2 Chemical and biological reactions in HP1HYDRUS-1D was recently coupled also with the PHREEQC geochemical code (Parkhurst &Appelo 1999) to create a new comprehensive simulation tool, HP1 (acronym for HYDRUS1DPHREEQC) (Jacques & Šimůnek 2005; Jacques et al. 2006; Šimůnek et al. 2006a). This newcode contains modules simulating (1) transient water flow in variably-saturated media, (2) thetransport of multiple components, (3) mixed equilibrium/kinetic biogeochemical reactions, and(4) heat transport. HP1 is a significant expansion of the individual HYDRUS-1D and PHRE-

108

Numerical Modelling of Hydrodynamics for Water Resources – Garcia-Navarro & Playán (eds)© 2008 Taylor & Francis Group, London, ISBN 978-0-415-44056-1

EQC programs by preserving most of their original features. The code still uses the Richards
equation for simulating variably-saturated water flow and advection-dispersion type equations
for heat and solute transport. However, the loosely coupled program can simulate also a broad
range of low-temperature biogeochemical reactions in water, the vadose zone and in ground water systems, including interactions with minerals, gases, exchangers and sorption surfaces basedon thermodynamic equilibrium, kinetic, or mixed equilibrium-kinetic reactions. HP1 uses theoperator-splitting approach with no iterations during one time step (a non-iterative sequential
modeling approach). Jacques et al. (2006) evaluated the accuracy of the operator-splitting approach for a kinetic reaction network (i.e., sequential and parallel kinetic degradation reactions)by comparing HP1 with an analytical solution for TCE-degradation, as well as for mixed equilibrium and kinetic reactions involving different flow conditions (steady-state and transient).
Jacques & Šimůnek (2005), and Šimůnek et al. (2006a) and Jacques et al. (2006; in press),demonstrated the versatility of HP1 on 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 andtransport of major cations (Na+, K+, Ca2+, and Mg2+) and heavy metals (Cd2+, Zn2+, and Pb2+) ina 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.
4.2.3 Wetland module
A multi-component reactive transport model CW2D (Constructed Wetlands 2D) (Langergraber& Šimůnek 2005, 2006) was developed to model the biochemical transformation and degradation processes in subsurface-flow constructed wetlands. The model was incorporated into the
HYDRUS (2D/3D) variably-saturated water flow and solute transport software package. Constructed wetlands have become increasingly popular for removing organic matter, nutrients,
trace elements, pathogens, or other pollutants from wastewater and/or runoff water. Such wetlands involve a complex mixture of water, substrate, plants, litter, and a variety of microorganisms to provide optimal conditions for improving water quality. The water flow regime in subsurface-flow constructed wetlands can be highly dynamic and requires the use of transientvariably-saturated flow model. The biochemical components defined in CW2D include dissolved oxygen, three fractions of organic matter (readily- and slowly-biodegradable, and inert),
four nitrogen compounds (ammonium, nitrite, nitrate, and dinitrogen), inorganic phosphorus,
and heterotrophic and autotrophic micro-organisms. Organic nitrogen and organic phosphorus
were modeled as part of the organic matter. The biochemical degradation and transformation
processes were based on Monod-type rate expressions, such as for NO3-based growth of heterotrophs on readily biodegradable COD (denitrification):r = μDN

K DN,O2

cNO3

K DN,NO2

K DN,O2 + cO2 K DN,NO3 + cNO3 K DN,NO2 + cNO2
cCR
f N,DN cXH
K DN,CR + cCR(25)

We refer to Langergraber & Šimůnek (2005, 2006) for a detailed discussion of the terms in (25).
All process rates and diffusion coefficients were assumed to be temperature dependent. Heterotrophic bacteria were assumed to be responsible for hydrolysis, mineralization of organic matter(aerobic growth) and denitrification (anoxic growth), while autotrophic bacteria were assumed
to be responsible for nitrification, which was modeled as a two-step process. Lysis was considered to be the sum of all decay and sink processes. Langergraber & Šimůnek (2005, 2006) demonstrated the model for one- and two-stage subsurface vertical flow constructed wetlands.
Model simulations of water flow, tracer transport, and selected biochemical compounds werecompared against experimental observations.

109

Numerical Modelling of Hydrodynamics for Water Resources – Garcia-Navarro & Playán (eds)© 2008 Taylor & Francis Group, London, ISBN 978-0-415-44056-1

5 NUMERICAL SOLUTIONS
All models discussed in this manuscript are numerical models that use various techniques tosolve the governing partial differential equations. While spatial derivatives are invariably approximated using the Galerkin-type linear finite element schemes, temporal derivatives in themodels are approximated using finite differences. We always used a fully implicit (backward) finite difference scheme for water flow, and a Crank-Nicholson time-centered scheme for heat andsolute transport. The mass-conservative method proposed by Celia et al. (1990) was used in all HYDRUS models. This method has been shown to provide excellent results in terms of minimizing themass balance error. Time and space discretization of the Richards equation generally leads to anonlinear system of algebraic equations. These equations are then linearized and solved using
the Picard iteration method. Higher-order approximations for the time derivative as derived by van
Genuchten (1976) are used for discretization of the transport equations. Depending upon the size ofthe problem, the matrix equations resulting from discretization of the governing equations are solved
using either Gaussian elimination for banded matrices, or the conjugate gradient method for symmetric matrices and the ORTHOMIN method for asymmetric matrices (Mendoza et al. 1991).
The operator-splitting approach with no iterations during one time step (non-iterative sequential approach) is used in all multicomponent models (i.e., for the major ion chemistry module of
HYDRUS-1D, in HP1, and in the wetland CW2D module of HYDRUS (2D/3D)).
Three-dimensional applications in general require a large number of finite elements to discretize realistically large transport domains. Even with the fast personal computers currently available, it is still very difficult to solve within a reasonable time problems having more than about
half million nodes. To decrease the required computational time, Hardelauf et al. (2007) parallelized SWMS_3D to develop PARSWMS that distributes problems with a large number of elements over multiple processors working in parallel. The PARSWMS code was developed for the
LINUX or UNIX workstations with installed free-wares MPI, PETSc and PARMETIS. Hardelauf et al. (2007) demonstrated that doubling the number of processors may lead to a decrease in
the computational time of up to 48 %.
6 GRAPHICAL USER INTERFACES
One major problem often preventing the use of numerical codes is the extensive work required fordata preparation, finite element grid design, and graphical presentation of the output results. Hence,a more widespread use of numerical models requires techniques which make it easier to create, manipulate and display large data files, and which facilitate interactive data management (Šimůnek et
al. 1999). Introducing such techniques frees users from cumbersome manual data processing, andenhances the efficiency in which programs are being implemented for a particular example. To
avoid or simplify the preparation and management of relatively complex input data files and tographically display final simulation results, we were very keen in developing interactive graphicsbased user-friendly interfaces (GUIs) for the MS Windows environments. From 1995 on this has resulted into several of the software packages described above. While earlier versions (until 1998)
were still 16-bit applications, all software packages from 1998 on are 32-bit applications fully compatible with all Windows operating systems released after Windows 95.
7 CONCLUSSIONSOver the last 15 years the close collaboration between the University of California Riverside,and the U.S. Salinity Laboratory, and more recently with PC-Progress in Prague, Czech Republic, and SCK•CEN, Mol, Belgium, resulted in the development of a large number of computer
tools that are currently being used worldwide for a variety of applications involving the vadosezone. The numerical models may be used to simulate not only one- and multi-dimensional water
flow, heat transport and solute transport in variably saturated media, but also the transport ofcarbon dioxide, major ions, trace elements, colloids, viruses, bacteria, and wastewater. The
popularity of the modeling tools is in large part due to their ease of use because of the availabil-

110

Numerical Modelling of Hydrodynamics for Water Resources – Garcia-Navarro & Playán (eds)© 2008 Taylor & Francis Group, London, ISBN 978-0-415-44056-1

ity of interactive graphical user interfaces. While the earlier models have proved to be useful for
analyzing subsurface flow and processes, several more complex models and modules are being
developed at present. One example is the HP1 code that couples the HYDRUS-1D water flow
and solute transport model with the PHREEQC geochemical model. Since the resulting combined program brings together the state-of-the-art in both flow and transport modeling and modeling of geochemical reactions, we believe that HP1 represents a very flexible and promising
tool for evaluating complex interactions of various processes in the vadose zone.
ACKNOWLEDGEMENT
This paper is based on work supported in part by the Terrestrial Sciences Program of the Army
Research Office (Terrestrial Processes and Landscape Dynamics and Terrestrial System Modeling and Model Integration), by the National Science Foundation Biocomplexity programs #0410055 and NSF DEB 04-21530, and by SAHRA (Sustainability of semi-Arid Hydrology andRiparian Areas) under the STC Program of the National Science Foundation, Agreement No.
EAR-9876800.
REFERENCES
Appelo, C. A. J. & Postma, D. 2005. Geochemistry, Groundwater and Pollution, 2nd ed., A. A. Balkema
Publ., Leiden, The Netherlands, 649 p.
Bradford, S. A., Bettehar, M., Šimůnek, J. & van Genuchten, M. Th. 2004. Straining and attachment of colloids in physically heterogeneous porous media, Vadose Zone J. 3(2), 384-394.
Bradford, S. A., Šimůnek, J., Bettahar, M., van Genuchten, M. Th. & Yates, S. R. 2006. Significance ofstraining in colloid deposition: evidence and implications. Water Resour. Res., 42, W12S15,
doi:10.1029/2005WR004791, 16 pp.
Brooks, R. H. & Corey, A. T. 1964. Hydraulic properties of porous media. Hydrol. Paper No. 3, Colorado State Univ., Fort Collins, CO.
Casey, F. X. M. & Šimůnek, J. 2001. Inverse analyses of the transport of chlorinated hydrocarbons subject
to sequential transformation reactions. J. of Environ. Quality, 30(4), 1354-1360.
Casey, F. X. M., Larsen, G. L., Hakk, H. & Šimůnek, J. 2003. Fate and transport of 17β-Estradiol in soilwater systems. Environ. Sci. Technol., 37(11), 2400-2409.
Castiglione, P., Mohanty, B. P., Shouse, P. J., Šimůnek, J., van Genuchten, M. Th., & Santini, A. 2003. Lateralwater diffusion in an artificial macroporous system: Theory and experimental evidence. Vadose Zone J., 2,
212-221.
Celia, M. A., Bououtas, E. T. & Zarba, R. L. 1990. A general mass-conservative numerical solution for the unsaturated flow equation. Water Resour. Res., 26, 1483-1496.
de Vries, D. A. 1958. Simultaneous transfer of heat and moisture in porous media. Transactions, American Geophysical Union, 39(5), 909-916.
Dontsova, K. M., Yost, S. L., Šimůnek, J., Pennington, J. C. & Williford, C. 2006. Dissolution and transport of TNT, RDX, and Composition B in saturated soil columns. J. of Environ. Quality, 35, 20432054.
Durner, W. 1994. Hydraulic conductivity estimation for soils with heterogeneous pore structure. Water Resour. Res., 32(9), 211-223.
Gargiulo, G., Bradford, S. A., Šimůnek, J., Ustohal, P., Vereecken, H. & Klumpp, E. 2007. Transport anddeposition of metabolically active and stationary phase Deinococcus Radiodurans in unsaturated porous media. Environ. Sci. and Technol., 41(4), 1265-1271.
Gerke, H. H., & van Genuchten, M. Th., 1993. A dual-porosity model for simulating the preferential
movement of water and solutes in structured porous media. Water Resour. Res. 29, 305-319.
Gonçalves, M. C., Šimůnek, J., Ramos, T. B., Martins, J. C., Neves, M. J. & Pires, F. P. 2006. Multicomponent solute transport in soil lysimeters irrigated with waters of different quality. Water Resour. Res.,
42, W08401, doi:10.1029/2006WR004802, 17 pp.
Hanson, B. R., Šimůnek, J. & Hopmans, J. W. 2006. Numerical modeling of urea-ammonium-nitrate fertigation under microirrigation. Agric. Water Management, 86, 102-113.
Hansson, K., Šimůnek, J., Mizoguchi, M., & Lundin, L.-Ch. 2004. Water flow and heat transport in frozen soil: Numerical solution and freeze/thaw applications. Vadose Zone Journal, 3(2), 693-704.

111

Numerical Modelling of Hydrodynamics for Water Resources – Garcia-Navarro & Playán (eds)© 2008 Taylor & Francis Group, London, ISBN 978-0-415-44056-1

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 groundwater flow process. Denver, CO, Reston, VA, U.S. Geological Survey.
Hardelauf, H., Javaux, M., Herbst, M., Gottschalk, S., Kasteel, R., Vanderborght, J. & Vereecken, H.
2007. PARSWMS: a parallelized model for simulating 3-D water flow and solute transport in variably
saturated soils. Vadose Zone Journal, 6(2), 255-259.
Haws, N. W., Rao, P. S. C., Šimůnek, J. & Poyer I. C. 2005. Single-porosity and dual-porosity modelingof water flow and solute transport in subsurface-drained fields using effective field-scale parameters. J.
of Hydrology , 313(3-4), 257-273.
Hendrickx, J. M. H., Dekker, L. W., & Boersma, O. H. 1993. Unstable wetting fronts in water repellent fieldsoils. J. Environ. Qual. 22, 109-118.
Jacques, D. & Šimůnek, J. 2005. User Manual of the Multicomponent Variably-Saturated Flow and Transport Model HP1, Description, Verification and Examples. Version 1.0, SCK•CEN-BLG-998, Waste
and Disposal, SCK•CEN, 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 flow and contaminant transport in layeredsoil profiles. J. Contam. Hydrology, 88, 197-218.
Jacques, D., Šimůnek, J., Mallants, D. & van Genuchten, M. Th. 2008. Modeling coupled hydrologicaland chemical processes in the vadose zone: a case study on long term uranium migration following
mineral p-fertilization, Vadose Zone J., Special Issue “Vadose Zone Modeling”, (in press).
Kodešová, R., Kozák, J., Šimůnek, J. & Vacek, O. 2005. Field and numerical study of chlorotoluron transport in the soil profile: Comparison of single and dual-permeability model. Plant, Soil and Environment, 51(7), 310-315.
Köhne, J. M., Köhne, S., Mohanty, B. P. & Šimůnek, J. 2004. Inverse mobile-immobile modeling oftransport during transient flow: Effect of between-domain transfer and initial water content. Vadose
Zone J., 3(4), 1309-1321
Köhne, J. M. & Mohanty, B. P. 2005. Water flow processes in a soil column with a cylindrical macropore: Experiment and hierarchical modeling. Water Resour. Res., 41, W03010,
doi:10,1029/2004WR003303.
Köhne, J. M., Mohanty, B. P. & Šimůnek, J. 2006. Inverse dual-permeability modeling of preferential water flow in a soil column and implications on field-scale solute transport prediction. Vadose Zone Journal, 5, 59-76.
Kosugi K. 1996. Lognormal distribution model for unsaturated soil hydraulic properties. Water Resour. Res.,32(9), 2697-2703.
Kung, K-J. S. 1990. Preferential flow in a sandy vadose zone. 2. Mechanism and implications. Geoderma,
46, 59-71.
Langergraber, G. & Šimůnek, J. 2005. Modeling Variably-Saturated Water Flow and Multi-Component
Reactive Transport in Constructed Wetlands. Vadose Zone Journal, 4, 924-938.
Langergraber, G. & Šimůnek, J. 2006. The Multi-component Reactive Transport Module CW2D for Constructed Wetlands for the HYDRUS Software Package, Manual – Version 1.0, HYDRUS Software Series 2, Department of Environmental Sciences, University of California Riverside, Riverside, CA, 72 pp.
Mendoza, C. A., Therrien, R. & Sudicky, E. A. 1991. ORTHOFEM User`s Guide, Version 1.02. Waterloo
Centre for Groundwater Research, Univ. of Waterloo, Waterloo, Ontario, Canada.
Parkhurst, D. L. & Appelo, C. A. J. 1999. User’s guide to PHREEQC (version 2) – A computer programfor speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations. Water
Resources Investigation, Report 99-4259, Denver, Co, USA, 312 pp.
Pitzer, K. S. 1979. Activity Coefficients in Electrolyte Solutions. Chap. 7, CRC Press, Boca Raton, Fl.
Pot, V., Šimůnek, J., Benoit, P., Coquet, Y., Yra, A. & Martínez-Cordón, M.-J. 2005. Impact of rainfall
intensity on the transport of two herbicides in undisturbed grassed filter strip soil cores. J. of Contaminant Hydrology, 81, 63-88.
Saito, H., J., Šimůnek, J. & Mohanty, B. 2006. Numerical analyses of coupled water, vapor and heat transport in the vadose zone. Vadose Zone Journal, 5, 784–800.
Scanlon, B., Keese, K., Reedy, R. C., Šimůnek, J. & Andraski, B. 2003. Variations in flow and transport
in thick desert vadose zones in response to paleoclimatic forcing (0 – 90 kyr): Monitoring, modeling,
and uncertainties. Water Resour. Res. 39(7), 1179, doi:10.1029/2002WR001604, 13.1-13.7.
Schaerlaekens, J., Mallants, D., Šimůnek, J., van Genuchten, M. Th. & Feyen, J. 1999. Numerical simulationof transport and sequential biodegradation of chlorinated aliphatic hydrocarbons using CHAIN_2D. J. of
Hydrological Processes, 13(17), 2847-2859.
Šejna, M. & Šimůnek, J. 2007. HYDRUS (2D/3D): Graphical User Interface for the HYDRUS Software
Package Simulating Two- and Three-Dimensional Movement of Water, Heat, and Multiple Solutes in
Variably-Saturated Media, published online at www.-pc-progress.cz, PC-Progress, Prague, Czech Republic.

112

Numerical Modelling of Hydrodynamics for Water Resources – Garcia-Navarro & Playán (eds)© 2008 Taylor & Francis Group, London, ISBN 978-0-415-44056-1

Seo, H. S., Šimůnek, J. & Poeter, E. P. 2006. Application of the Unsaturated Flow (UNSF) Package,
Modflow Conference Proceedings, ed. E. Poeter, Ch. Zheng, M. Hill, and J. Doherty, International
Ground Water Modeling Center, Colorado School of Mines, May 22-24, Golden, CO.
Seo, H. S., Šimůnek, J. & Poeter, E. P. 2007. Documentation of the HYDRUS package for MODFLOW2000, the U.S. Geological Survey Modular Ground-Water Model. GWMI 2007-01. International
Ground Water Modeling Center, Colorado School of Mines, Golden, Colorado, 96 pp.
Šimůnek, J. & Suarez, D. L. 1993. Modeling of carbon dioxide transport and production in soil: 1. Model development. Water Resour. Res., 29(2), 487-497.
Šimůnek, J. & Suarez, D. L. 1994. Two-dimensional transport model for variably saturated porous media withmajor ion chemistry. Water Resour. Res., 30(4), 1115-1133.
Šimůnek, J., Vogel, T. & van Genuchten, M. Th. 1994. SWMS_2D simulating water flow and solute transport
in two-dimensional variably saturated media. Version 1.2, IGWMC - FOS 61, International Ground Water
Modeling Center, Colorado School of Mines, Golden, Colorado, 196 pp.
Šimůnek, J., Huang, K. & van Genuchten, M. Th. 1995. The SWMS_3D code for simulating water flow andsolute transport in three-dimensional variably saturated media. Version 1.0, Research Report No. 139, U.S.
Salinity Laboratory, USDA, ARS, Riverside, California, 155 pp.
Šimůnek, J., Šejna, M. & van Genuchten, M. Th. 1999. The HYDRUS-2D software package for simulating two-dimensional movement of water, heat, and multiple solutes in variably saturated media. Version 2.0, IGWMC - TPS - 53, International Ground Water Modeling Center, Colorado School ofMines, Golden, Colorado, 251 pp.
Šimůnek, J., Wendroth, O., Wypler, N., & van Genuchten, M. Th. 2001. Nonequilibrium water flow characterized from an upward infiltration experiment. European J. of Soil Sci., 52(1), 13-24.
Šimůnek, J. & Valocchi, A. J. 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, 15111536.
Šimůnek, J. 2003. HYDRUS-2D Code Modification: Modeling Overland Flow and Dynamic Interactions
Between Plants and Water Flow in a Hillslope Transect. Idaho National and Environmental Laboratory, Final report for project under contract number 00021013, 55 pp.
Šimůnek, J., Jarvis, N. J., van Genuchten, M. Th. & Gärdenäs, A. 2003. Nonequilibrium and preferential
flow and transport in the vadose zone: review and case study. Journal of Hydrology, 272, 14-35.
Šimůnek, J., van Genuchten, M. Th. & Šejna, M. 2005. The HYDRUS-1D software package for simulating
the one-dimensional movement of water, heat, and multiple solutes in variably-saturated media. Version
3.0, HYDRUS Software Series 1, Department of Environmental Sciences, University of California Riverside, Riverside, CA, 270 pp.
Šimůnek, J., Jacques, D., van Genuchten, M. Th. & Mallants, D. 2006a. Multicomponent geochemicaltransport modeling using the HYDRUS computer software packages. J. Am. Water Resour. Assoc.,
42(6), 1537-1547.
Šimůnek, J., van Genuchten, M. Th. & Šejna, M. 2006b. The HYDRUS Software Package for Simulating Twoand Three-Dimensional Movement of Water, Heat, and Multiple Solutes in Variably-Saturated Media.
Technical Manual, Version 1.0, PC Progress, Prague, Czech Republic, pp. 241.
Šimůnek, J. & van Genuchten, M. Th. 2008. Modeling nonequilibrium flow and transport with HYDRUS. Vadose Zone Journal, Special Issue “Vadose Zone Modeling”, (in press).
Šimůnek, J., van Genuchten, M. Th. & Šejna, M. 2008. Development and applications of HYDRUS andrelated software packages. Vadose Zone Journal, Special Issue “Vadose Zone Modeling” ”, (in press).
Truesdell, A. H. & Jones, B. F. 1974. WATEQ, a computer program for calculating chemical equilibria
of natural waters. J. Res. U. S. Geol. Surv., 2(2), 233-248.
Twarakavi, N. K. C., Šimůnek, J. & Seo, H. S. 2008. Evaluating interactions between groundwater andvadose zone using HYDRUS-based flow package for MODFLOW, Vadose Zone J., Special Issue“Vadose Zone Modeling”, (in press).
van Genuchten, M. Th. 1976. On the accuracy and efficiency of several numerical schemes for solving theconvective-dispersive equation, in Finite Elements in Water Resources, edited by W. G. Gray et al., Pentech Press, London, pp. 1.71-1.90.
van Genuchten, M. Th. & Wierenga. P. J. 1976. Mass transfer studies in sorbing porous media, I. Analyticalsolutions. Soil Sci. Soc. Am. J. 40, 473-481.
van Genuchten, M. Th. 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J., 44, 892-898.
Wang, D., Yates, S. R. & Gan, J. 1997. Temperature effect on methyl bromide volatilization in soil fumigation. Journal of Environmental Quality, 26(4), 1072-1079.
Wang, D., Knuteson, J. A. & Yates, S. R. 2000. Two-dimensional model simulation of 1,3dichloropropene volatilization and transport in a field soil. Journal of Environmental Quality, 29(2),639-644.

113

Numerical Modelling of Hydrodynamics for Water Resources – Garcia-Navarro & Playán (eds)© 2008 Taylor & Francis Group, London, ISBN 978-0-415-44056-1

Zhang, P., Šimůnek, J., & Bowman, R. S. 2004. Nonideal transport of solute and colloidal tracers throughreactive zeolite/iron pellets. Water Resour. Res., 40, doi:10.1029/2003WR002445.

114

Numerical Modelling of Hydrodynamics for Water Resources – Garcia-Navarro & Playán (eds)© 2008 Taylor & Francis Group, London, ISBN 978-0-415-44056-1

OUTLINE:
1. Introduction
2. Water Flow
2.1. Uniform Flow
2.2. Nonequilibrium Flow
2.3. Water Flow and Vapor Transport
2.4. Water Flow with Freezing and Thawing
2.5. Overland Flow
2.6. HYDRUS package for MODFLOW
3. Heat Transport
3.1. Heat Transport without Vapor Transport
3.2. Heat Transport with Vapor Transport
3.3. Heat Transport with Freezing and Thawing
4. Solute Transport
4.1. Single Component Solute Transport
4.1.1. Equilibrium Transport
4.1.2. Nonequilibrium Transport
4.2. Multicomponent Solute Transport
4.2.1. Major Ion Chemistry
4.2.2. Chemical and Biological Reactions in HP1
4.2.3. Wetland Module
5. Numerical Solutions
6. Graphical Interfaces
7. Summary and Conclusions

115