J. Hydrol. Hydromech., 66, 2018, 2, 133–142
DOI: 10.1515/johh-2017-0050

New features of version 3 of the HYDRUS (2D/3D) computer software
packageJiří Šimůnek1*, Miroslav Šejna2, Martinus Th. van Genuchten3, 4
1

Department of Environmental Sciences, University of California, Riverside, CA 92521, USA.
PC-Progress, Ltd., Prague, Czech Republic.
3
Center for Environmental Studies, CEA, São Paulo State University, UNESP, Rio Claro, SP, Brazil.
4
Department of Earth Sciences, Utrecht University, Utrecht, Netherlands.
2*

Corresponding author. E-mail: Jiri.Simunek@ucr.edu
Abstract: The capabilities of the HYDRUS-1D and HYDRUS (2D/3D) software packages continuously expanded during the last two decades. Various new capabilities were added recently to both software packages, mostly by developing
new standard add-on modules such as HPx, C-Ride, UnsatChem, Wetland, Fumigant, DualPerm, and Slope Stability.
The new modules may be used to simulate flow and transport processes in one- and two-dimensional transport domains
and are fully supported by the HYDRUS graphical user interface (GUI). Several nonstandard add-on modules, such as
Overland, Isotope, and Centrifuge, have also been developed, but are not fully supported by the HYDRUS GUI. The objective of this manuscript is to describe several additional features of the upcoming Version 3 of HYDRUS (2D/3D),
which was unveiled at a recent (March 2017) HYDRUS conference and workshop in Prague. The new features include a
flexible reservoir boundary condition, expanded root growth features, and new graphical capabilities of the GUI. Mathematical descriptions of the new features are provided, as well as two examples illustrating applications of the reservoir
boundary condition.
Keywords: HYDRUS; Reservoir boundary condition; Pumping well; Root growth; Graphical user interface.

INTRODUCTION
The HYDRUS (2D/3D) software package and its various
predecessors (e.g., UNSAT, SWMS-2D, CHAIN-2D, and HYDRUS-2D) have a long history that goes back to the early
1970s as documented in detail by Šimůnek et al. (2008, 2016a).
The HYDRUS conferences often serve as a forum where newversions of HYDRUS (2D/3D) are unveiled. For example,
Version 2 of HYDRUS (2D/3D) and several of its add-on modules (UnsatChem, DualPerm, C-Ride, HP2) were first
introduced in 2013 at the Third International Conference on“HYDRUS Software Applications to Subsurface Flow andContaminant Transport Problems,” held in Prague, Czech Republic. Similarly, Version 3 of HYDRUS (2D/3D) was
unveiled in 2017 at the Fifth HYDRUS conference in Prague.
Similarly, as Version 2, Version 3 is supported by a number ofstandard and non-standard specialized HYDRUS modules thatwere summarized by Šimůnek et al. (2016a).
The objective of the present paper is to describe in detail two
main computational modules that have been included in Version 3 of HYDRUS (2D/3D): a new "Reservoir" boundary
condition” (BC) and a dynamic root system. We provide here
not only mathematical descriptions of the new standard computational modules but also two examples illustrating the use of
the reservoir BC: i.e., for evaluating the performance of a
pumping well as well as of a falling head experiment. Additionally, selected new graphical capabilities of the software
package are described.
THE RESERVOIR BOUNDARY CONDITION
Previous versions of HYDRUS (2D/3D) required that pressure heads at boundaries representing external water bodies be
specified as external input. No feedback was possible to relate
temporal changes of the water level in these external water

bodies in response to interactions with the subsurface, such as
infiltration or exfiltration. HYDRUS users typically assumed
that the water level in an external water body (e.g., a furrow orborehole) was constant (e.g., Hinnell et al., 2009; Lazarovitch
et al. 2009; Warrick et al., 2007), or they needed to prescribe
the dynamics of the water level themselves.
Version 3 of the standard two-dimensional computational
module of HYDRUS (2D/3D) offers a new system-dependent
boundary condition, further referred to as the reservoir boundary condition. This option allows users to consider a reservoir
that is external to the HYDRUS transport domain, while water
can be added (injected) to or removed (pumped) from thereservoir. Flow into or out of the reservoir through its interface
with the subsurface depends on the prevailing conditions in theflow domain (e.g., the position of the groundwater table) andexternal fluxes. Since mass balances of water and solute in theexternal reservoir are constantly being updated (based on all
incoming and outgoing fluxes), the boundary conditions along
the transport domain are dynamically adjusted depending upon
the water level in the reservoir.
Reservoir boundary conditions potentially have a large number of applications such as estimating dynamically the waterlevel in wells, in furrows during irrigation, and in wetlands.
Another application concerns a relatively new approach to landdevelopment known as Low-Impact Development (LID), which
is a “green” approach to storm-water management that seeks tomimic the natural hydrology of a site using decentralized microscale control measures (Coffman, 2002). Low-impact
development practices range from the use of bioretention cells,
infiltration (dry) wells or trenches, stormwater wetlands, wetponds, level spreaders, permeable pavements, swales, green
roofs, vegetated filter and buffer strips, sand filters, smaller
culverts, to water harvesting systems (Brunetti et al., 2016,2017). We provide later two applications where the reservoir
BC is used to evaluate fluxes into or out of wells, as well as to133Jiří Šimůnek, Miroslav Šejna, Martinus Th. van Genuchten

Fig. 1. Three different types of the reservoir boundary condition: a well, a furrow, and a wetland. In the figure hw is the water level in thereservoir [L], S is the volume of water in the reservoir ([L2] or [L3] for two-dimensional and axisymmetric systems, respectively), Qp is thepumping rate (positive for removal of water, negative for adding water) ([L2T–1] or [L3T–1] for two-dimensional and axisymmetric systems,
respectively), c is the solute concentration in reservoir water [ML–3], cp is the solute concentration in injected water [ML–3], P and E refer tothe precipitation and evaporation rates [LT–1], rw is the radius of the well [L], a is a half-width of the furrow [L], α is the slope of the furrowside [–], rmax is the maximum width (radius) of the wetland [L], zmax is the maximum depth of the wetland [L], and b represents the width of
the water surface ([L] or [L2] for two-dimensional and axisymmetric systems, respectively).

simulate a falling head experiment.
The reservoir BC is implemented in HYDRUS (2D/3D) for
three different geometries (Fig. 1): a well, a furrow, and awetland. While the furrow BC is implemented only for twodimensional transport domains, the well and wetland BCs are
implemented for both two-dimensional planar and axisymmetric domains.
The Well Reservoir Boundary Condition
The position of the water level in the well is obtained bysolving the following mass balance equation:dVw
dh= πrw2 w = Qin ( t ) − Q p ( t )
dt
dt

(1)

which in HYDRUS is implemented in terms of the finite difference discretization:

hwj +1 − hwj= Qin − Q pΔtΔt
hwj +1 = hwj + 2 ( Qin − Q p )
πrw
πrw2

(2)

where Vw is the volume of water in the well [L3], rw is the wellradius [L], hw is the water level in the well [L], Qin is the waterflow rate into the well from the soil profile across the well wall(or its screened part) [L3T–1], Qp is the pumping rate [L3T–1], Δt
is the time step [T], and hwj and hwj +1 are water levels in thewell at the previous and current time levels [L], respectively(Fig. 1a). The fluxes Qin and Qp can have negative values, inwhich case they represent outflow from the well (infiltration
into the soil profile) and the water flow rate entering the well,respectively.
Equation (1) can be solved providing that the initial positionof the water level in the well, hw,init, and the pumping rate are

134

known. The parts of the boundary below and above the waterlevel in the well are then assigned internally to be (timevariable) pressure head (Dirichlet) and seepage face boundary
conditions, respectively. During execution, HYDRUS calculates which part of the seepage face boundary is active (with aprescribed zero pressure head) and which is inactive (with aprescribed zero flux). HYDRUS also calculates and reports
fluxes separately across these two parts of the boundary along awell, and uses these fluxes in Eq. (1). Note that this option canalso be used for 2D geometries. The surface area πrw2 is thenreplaced with the width, rw, of the well bottom.
A similar mass balance equation as for water flow is usedalso for solute transport:

dVw c
dV
dc
= Vw
+ c w = Qin cin − Q p c p
dt
dt
dt
dc
Vw
= Qin cin − Q p c p − c ( Qin − Q p )dt(3)where c is the solute concentration in the reservoir [ML–3], cin issolute concentration associated with mass transfer between thesoil profile and water reservoir [ML–3] (cin is equal to c whenwater infiltrates into the soil profile and equal to the soluteconcentration in the soil profile when water exfiltrates into thereservoir), and cp is the solute concentration associated withpumping or injection (equal to c for pumping or to theconcentration of water being injected into the reservoir) [ML–3].
An energy balance equation similarly as Eq. (3) was implemented also to account for heat transport and time-variabletemperatures in the reservoir.
The Furrow Reservoir Boundary Condition
Similarly as for wells, the reservoir BC for furrows within
HYDRUS is based on mass balance considerations to determine
the position of the water level in the furrow, hw [L]. The volume
of water, S [L2], in the furrow, and its change in time depending

New features of version 3 of the HYDRUS (2D/3D) computer software package

upon the inflow and outflow rates, are described using theequations (Šimůnek et al., 2016b):

h 
hw2
1
1 
S = h( a + b) = hw  a + a + w  = hw a +22 tan α 2 tan αdS
= Qin ( t ) − Q p ( t ) + ( P − E )bdtdh
h dhw h  dha w+ w= a + w  w =dt tan α dt tan α  dt

Similarly, as for the well and furrow boundary conditions,users can specify a reservoir boundary condition along theboundary of the wetland (Fig. 1c). The bottom coordinates of awetland are described using the following equations:
(4)

 r 
z = zmax  rmax 

p

(7)

1/ p

 z 
r = rmax  zmax = Qin ( t ) − Q p ( t ) + ( P − E )bor in terms of their finite difference discretization:
hwj  hwj +1 − hwj= Qin − Q p + ( P − E )ba +
tan α  Δt
Δt tan α
hwj +1 = hwj + j( Qin − Qp + Pb − Eb )hw + a tan α

The Wetlands Reservoir Boundary Condition(5)

where a is the half-width of the bottom of the furrow [L], b is
the half-width of the water level [L], hw is the water level in the
furrow [L], Qin the water infiltration rate from the furrow to thesoil profile across the furrow walls [L2T–1], Qp is the pumping
rate (positive for removal of water, negative for adding water)[L2T–1], Δt is the time step [T], P and E are precipitation andevaporation rates [LT–1], respectively, and hwj and hwj +1 arewater levels in the furrow at the previous and current timelevels [L], respectively (Fig. 1b). Similarly, as for the well BC,
parts of the boundary below and above the water level in the
furrow are then assigned time-variable pressure head (Dirichlet)
and seepage face boundary conditions, respectively.
A similar solute mass balance approach as for wells is usedalso for the furrow BC (Šimůnek et al., 2016b):d ( Sc )

dc
dS
=S
+c
= Q p ( t ) c p − Qin ( t ) c
dt
dt
dt
dc
S
= Q p ( t ) c p − Qin ( t ) c − c Q p ( t ) − Qin ( t ) + ( P − E )b  (6)dtdc
S= Q p ( t ) ( c p − c ) − c [ ( P − E )b ]dt
where cp is the solute concentration of irrigation water (fertigation) [ML–3; or dimensionless], and c is the average soluteconcentration of the furrow water and hence the concentrationof the infiltrating water [ML–3]. Note that we assume that precipitation and evaporation fluxes are devoid of solutes and that
there is instantaneous and complete mixing of solute within the
furrow. Precipitation will, therefore, lead to dilution of solute in
the furrow water, while evaporation will lead to increasing
concentrations. The furrow BC module requires as input the soluteconcentration in the irrigation water (cp) and the timing of fertigation (i.e., the beginning and end time of fertigation). The module
then calculates the solute concentration in the furrow water, c,which is subsequently used in a Cauchy (concentration flux)
boundary condition together with the local infiltration flux calculated with HYDRUS. Similarly, as the standard computational
module of HYDRUS, the Reservoir boundary condition cannot be
used to account for precipitation/dissolution of mineral phases in
the reservoir or its boundaries.

where zmax is the depth of a wetland [L], rmax is the radius of awetland [L], and p is a shape parameter (Fig. 1c). The volume
of water in the wetland, Vw ([L2] or [L3] for two- and threedimensional problems, respectively), and the water depth, hw[L], for two-dimensional problems are calculated using:
1/ p

 z 
dV = rdz = rmax  zmax V=rmax( zmax

dz
(8)

rmax

z dz =
) 
(z )
1/ p

1/ p

1/ p

max

 p  1/ p +1z
 1+ p p 1 + p  ( z max )1/ p  1+ p
hw = 
V p  rmax
(9)and for axisymmetric problems using: z 
dV = πr dz = π ( rmax )  zmax 

2/ p

2

2

dz
(10)
2

( rmax ) z 2/ p dz = π  rmax   p  z 2/ p +1
2/ p
1/ p  ( zmax )  zmax   2 + p 
2

V =π

p

 1  2 + p   z 1/ p 2  2+ p
max
hw =   V π  p  rmax  

(11)

A similar solute mass balance equation as for the well and furrow reservoir BCs is also used for the wetlands reservoir BC.
Before continuing with the two applications, we note thatwater flow and solute transport can occur in both directions
between a reservoir and the transport domain depending upon
the prescribed initial and boundary conditions. For example, areservoir can be initially empty while groundwater table ispresent in the transport domain. Water will then flow from thetransport domain into the reservoir, with HYDRUS dynamically adjusting the boundary conditions at the reservoir walls. The
code will assign a time-variable pressure head boundary condition to the reservoir walls that are below the water level in thereservoir, and a seepage face boundary condition to the walls
above the water level. At the same time, water may be pumped
from or added to the reservoir using other boundary conditions.
Alternatively, water could be initially present in or be added to
the reservoir, while the transport domain is initially still dry.
Water will then flow from the reservoir into the transport domain, in which case the water level in the reservoir and corresponding boundary conditions will be adjusted dynamically
based on the fluxes into and out of the reservoir.

135

Jiří Šimůnek, Miroslav Šejna, Martinus Th. van Genuchten

Additionally, the interface between the reservoir and thetransport domain is also selected by the software users. For
example, in the case of a well reservoir, users can assume that
only the bottom or the walls of the well are permeable, or thatsome part of the well is lined (impermeable to water) or incontact only with certain soil horizons. Since the initial andboundary conditions for the reservoir and the transport domain
are specified independently, many different combinations ofthese various scenarios can hence be accommodated with the
new reservoir boundary condition module is HYDRUS. Two
examples of these features are given in the next section.
RESERVOIR BOUNDARY CONDITIONAPPLICATIONS
We now provide two examples illustrating the use of thewell reservoir boundary condition. The first example pertains toa pumping well, with water flowing from the transport domain
into the well from which water is being pumped. The waterlevel in the well depends on the balance between water inflow
and pumping. The second example considers falling head infiltration, with water flowing from the reservoir into the transportdomain, while the water level in the reservoir is correspondingly being adjusted. This example represents a classic falling head
infiltration experiment. Additional examples involving furrow
reservoir boundary conditions for both water flow and solutetransport problems are given by Siyal et al. (2012) and Šimůnek
et al. (2016b).

Fig. 2. The transport domain and its discretization for the pumping
well example. Values in the figure are in meters.

A pumping wellThis example considers conditions typical of a pumping
well. The transport domain in the example is assumed to be 5 mwide and 3 m deep. The well has a diameter of 0.25 m, is 2 m
deep, and has an initial water level of 0.50 m above the wellbottom. The soil profile is considered to be homogeneous andassumed to consist of loam, with the default van Genuchten(1980) parameters taken from the HYDRUS soil catalog basedon textural class-averaged pedotransfer functions derived by
Carsel and Parrish (1988). Initial conditions in the profile wereassumed to be in static equilibrium with a groundwater tablelocated 50 cm below the soil surface. Inflow into the well occurred through both the well bottom and its sides. Timevariable pressure head and seepage face boundary conditions
were applied at the well boundary below and above the waterlevel in the well, respectively. A constant pressure head boundary condition, corresponding to the initial pressure head conditions was specified along the right side of the domain, while noflow boundary conditions were imposed at all other boundaries.
The pumping rate from the well was set equal to 0.03 m3/h.
The transport domain was discretized into an unstructured
triangular finite element mesh, with a targeted element size of
0.075 m. The size of finite elements was adjusted to 0.03 m at
the well bottom and 0.3 m at the right side of the domain. The
transport domain and its discretization are shown in Figure 2.
Figure 3 shows the water level and water volume in the well,as well cumulative water fluxes, as a function of time. Initially,
the pumping rate was slightly smaller than the inflow rate ofwater from the soil into the well through the well walls below(Exfiltration) and above (Seepage) the water level in the well.
As a result, the water level in the well and the well water volume slightly increase. However, the inflow of water into thewell after about 10 hours was about equal to, and then became
smaller than, the pumping rate, leading to a gradual but slow
decrease in the water level in the well. Figure 4 shows calcu-

136

Fig. 3. Calculated water levels and water volumes in the well, andcumulative water fluxes as a function of time for the pumping well
example. Note that the terms "Exfiltration" and "Seepage" represent cumulative flow from the transport domain into the well,while "Pumping" refers to the removal of water from the well.

Fig. 4. Calculated pressure head (m) profile and position of the watertable after 48 hours of pumping for the pumping well example.

lated pressure head profiles and the position of the water table
after 48 hours of pumping.
While relatively simple, the example demonstrates the use of
the new reservoir boundary condition. The problem in actuality
can be made much more complex, for example by assuming a

New features of version 3 of the HYDRUS (2D/3D) computer software package

heterogeneous soil profile, making the pumping rate time variable, assuming different parts of the well to be screened, andincluding contaminant transport among essentially all otherfeatures available in the standard HYDRUS (2D/3D) software
package.
Falling head infiltration
In this example, the well reservoir boundary condition isused to simulate a falling head experiment when the water
reservoir is placed on top of the soil. The transport domain inthis example was 1.25 m wide and 1.30 m deep, and discretized
into a structured finite element mesh consisting of 380 nodes
and 684 elements. The soil profile was again considered to behomogeneous and to consist of loam, with the default van
Genuchten (1980) parameters taken from the HYDRUS soil
catalog. The radius of the reservoir placed on top of the soil
surface was 0.20 m, and the initial water level in the reservoir
0.15 m above the soil surface.
Figure 5 shows that the water level and the water volume in
the reservoir decreased as a result of water infiltrating into thesoil profile. At the end of the experiment, all water initially inthe reservoir infiltrated into the soil profile. Note that in Figure
5 the “Volume” at time zero is equal to cumulative “Infiltration” at the end of the experiment. Figure 6 shows the calculated pressure head profiles at 0.2 d when water is still infiltrating
into the profile, and at 0.5 d during redistribution when thereservoir is empty. The plots hence represent situations during
and after the falling head infiltration experiment.
ROOT GROWTH
Previous Version 2 of HYDRUS(2D/3D) included a relatively comprehensive macroscopic root water and solute uptake
model (Šimůnek and Hopmans, 2009) to account for the effects
of both water and salinity stress on root water uptake, while
additionally accounting for possible active and passive root
contaminant or nutrient uptake. Root water and solute uptake
both could furthermore be treated as being either noncompensated or compensated, while users could select thedegree of compensation (Šimůnek and Hopmans, 2009).
HYDRUS-1D additionally allowed users to externally prescribe a time-variable rooting depth, either using the logistic
growth function or a tabulated form. This feature was thus far
not available in HYDRUS (2D/3D), which has been noted by

Fig. 5. Calculated water levels and water volumes in the reservoir
on top of the soil, and cumulative infiltration volumes into the soil,
as a function of time for the falling head experiment.

several HYDRUS users (e.g., Karandish and Šimůnek, 2016;
Ramos et al., 2012; Roberts et al., 2009; among others). The
spatial distribution of roots in the root zone was assumed toremain constant during the simulations, and only the intensityof potential uptake (as a result of time-variable potential transpiration rates) could be specified as input. The HYDRUS
models also did not allow the spatial extent of the rooting zone
to change actively as a result of certain environmental stresses(e.g., Hartman et al., 2018).
To overcome these deficiencies, several studies either further modified the HYDRUS models (or their predecessors suchas CHAIN-2D or SWMS-3D), or coupled the models withexisting crop growth or root growth models. For example,
Javaux et al. (2008, 2013) developed R-SWMS, a threedimensional root growth model that couples the model ofSomma et al. (1998) (based on SWMS-3D) with the root architecture model of Doussan et al. (1998). Zhou et al. (2012) coupled HYDRUS-1D with the WOFOST crop growth model(Boogaard et al., 1998), and used the resulting code to simulate
the growth and yield of irrigated wheat and maize (Li et al.,2012, 2014). Han et al. (2015) similarly coupled HYDRUS-1D
with a simplified crop growth version used in SWAT, based on
the EPIC root growth model (Williams et al., 1989). They usedthe coupled model to simulate the contribution and impact ofgroundwater on cotton growth and root zone water balance.

Fig. 6. Calculated pressure head (cm) profiles at 0.2 (left) and 0.5 d (right) for the falling head experiment.

137Jiří Šimůnek, Miroslav Šejna, Martinus Th. van Genuchten

Wang et al. (2014, 2015) further coupled EPIC with CHAIN-2Dand HYDRUS-1D to assess the effects of furrow and sprinkler
irrigation, respectively, on crop growth. Finally, Hartmann et al.
(2018) implemented into both HYDRUS-1D and HYDRUS(2/3D) a root growth model developed by Jones et al. (1991).
Their coupled model assumes that various environmental growthstress factors can influence root development under suboptimal
conditions. Of the different software, only the models of Hartmann et al. (2018) are directly available from the HYDRUS
website as a non-standard HYDRUS module (i.e., not fully
supported by the HYDRUS GUI).
To extend the capabilities of the standard module of HYDRUS (2D/3D), a simple root growth model with similar capabilities as those in HYDRUS-1D was implemented into Version
3. The rooting depth, LR, can now be either constant (the standard approach) or variable during the simulations. For annual
vegetation, a growth model is required to simulate changes inrooting depth with time. Time-variable rooting depth values canbe provided either using a table on input, or calculated with theprogram assuming that the actual rooting depth is the product ofthe maximum rooting depth, Lm [L], and a root growth coefficient, fr(t) [–], given by:

LR ( t ) = Lm f r ( t )(12)

For the root growth coefficient, fr(t), we use the classicalVerhulst-Pearl logistic growth function

f r (t ) =L0
L0 + ( Lm − L0 ) e − rt(13)

where L0 is the initial value of the rooting depth at the beginning of the growing season [L], and r the growth rate [T–1]. The
growth rate can be calculated either from the assumption that
50% of the rooting depth will be reached after 50% of thegrowing season has elapsed or from given data of the rooting
depth at a specified time. The same approach is also used forthe horizontal extent of the rooting zone, except that the maximum extent of the rooting zone in the horizontal direction isused instead of the maximum rooting depth.
When a variable rooting depth is considered, the spatial distribution of roots must be described using either the Vrugt(Vrugt et al., 2001a, 2001b) or Hoffman and van Genuchten(Hoffman and van Genuchten, 1983) functions. The Vrugt twodimensional root distribution function is implemented in HYDRUS as follows (Vrugt et al., 2001a, 2001b):
 pz

*

pr

*




z 
x  − Zm z − z + X m x − x b ( x, z ) =  1 − 1 −eXm 
 Z m (14)

where Xm and Zm are the maximum rooting lengths in the x- andz- directions [L], respectively; x and z are distances from theorigin of the plant (or tree) in the x- and z- directions [L],respectively; px [–], pz [–], x* [L], and z* [L] are empiricalparameters (Vrugt et al., 2001b), and b(x,z) denotes the twodimensional spatial distribution of the potential root water
uptake rate [–]. The parameters x* and z* indicate the location in
the profile having the maximum rooting density, while px and pz
are assumed to be zero for x > x* and z > z*, respectively. We
refer to Vrugt et al. (2001a, 2001b) for different configurationsof the normalized spatial distribution of the potential root water
uptake rate.
Alternatively, one can use the following function (Hoffman
and van Genuchten, 1983):

1.667


Zm

 2.0833  L − z b( z ) = 
1 −

Zm 
 Zm 
0
z > L − 0.2Z m
z ∈ ( L − Z m ; L − 0.2 Z m ) (15)z < L − Zm

where L is the x-coordinate of the soil surface [L] and Zm is therooting depth [L].
Figure 7 shows an example of the two-dimensional spatial
distribution of root water uptake at 10, 20, 30, and 40 d in atransport domain that is 100 cm wide and 150 cm deep. The
plots are for the following root growth parameters: L0x = L0z = 5
cm at t = 5 d, Zm = Lmz = 100 cm at t = 45 d, Xm = Lmx = 75 cm,
x* = 0, z* = 20 cm, px = pz = 1. The root growth rate r for thesevariables was equal to –0.0309 d –1.

Fig. 7. Dynamic two-dimensional spatial distribution of root water uptake at 10, 20, 30, and 40 d for root growth parameters given in the
text. The transport domain is 100 cm wide and 150 cm deep. The scale varies from zero (blue) to one (red) and is dimensionless.

138

New features of version 3 of the HYDRUS (2D/3D) computer software package

THE HYDRUS (2D/3D) GRAPHICAL USER INTERFACE
The popularity of the HYDRUS software packages is inlarge part due to its inclusion of a large number of subsurface
and near-surface water flow and solute transport processes(Šimůnek et al., 2016a). However, equally important has been
the development of advanced graphical user interfaces (GUIs).
The GUI of Version 3 of HYDRUS (2D/3D) now includes asignificant number of major new features, as well as a large
number of smaller corrections and improvements. Although theexternal look has been preserved and is similar to Version 2,
Version 3 internally uses the latest software development tools
and libraries, which is critical to ensure compatibility with new
Windows operating systems and development of more efficient
software in the future. This section provides a brief overview ofthe most important improvements and extensions.
One important objective with new version was to increase
the performance and capacity of the GUI. The limit of Version
2 of HYDRUS (about 1 million finite elements) was extended
by almost an order to 10 million finite elements. The extension
was achieved by overall optimization of the code and developing a 64-bit version of HYDRUS, which can now use all physical memory available on modern PCs. This feature permitted aconsiderable increase in the size of HYDRUS simulations, withthe main limiting factor now being the speed of the calculation
module(s).
We further replaced many old HYDRUS (2D/3D) software
components for charts and tables in Version 3 by new components, which also included fixing some known errors. Now all
dialogs containing charts or tables are resizable, thus making
entering and viewing data more comfortable. The GUI also
supports the Unicode and can now correctly display special
characters such as Greek symbols. Although manipulating
graphical objects in HYDRUS, in general, was relatively easy
thanks to drag-and-drop features, some operations (such as
rotations) had to be defined numerically in a dialog box. Version 3 offers a new graphical tool called “Manipulator” formore user-friendly transformations of selected objects in two or
three dimensions.
We also improved several features related to Mesh-sections,
3D mesh clipping and slicing. The optimization of a GUI for
work with relatively large FE meshes required certain changes
in using mesh-sections. While mesh-sections are still fully
supported, the program by default does not generate as manymesh-sections as in Version 2. However, we created a newgraphical tool called “Clipper” that can now be used to cut orslice a 3D mesh. Figures 8 and 9 show examples of these two
applications. Another useful function related to mesh-sections
is the possibility to select mesh nodes or elements by selecting
geometric objects in a data-explorer tree.
The graphical display of velocity vectors for millions ofmesh nodes is not only slow but usually also too complex forviewing. We improved the new GUI to allow displays of velocity vectors either at mesh nodes or raster points, while the rasterparameters (such as the density of points) are fully adjustable.
Figure 10 shows an example.
Previous Version 2 of HYDRUS (2D/3D) allowed one toshow flowing particles, but only in two dimensions. Flowing
particles are hypothetical objects that can be defined at any
location of the transport domain by users, with the program
then calculating trajectories and positions of these particles withtime considering unretarded convective transport. The “flowing
particles” feature is now available for both 2D and 3D projects.
Another new feature of Version 3 is the calculation and display of streamlines for a given steady-state flow velocity field.

Streamlines are one of the most commonly used graphicalrepresentations of CFD (computational flow dynamics) resultsin that they can display very clearly the flow direction, such as
shown in Figure 11. The software now also has an option to run
flow animations, i.e., the movement of particles along streamlines (Fig. 12). Users can save this animation as a video file anduse in presentations of HYDRUS results.
The ability of charts displaying the distribution of various
quantities along boundary-lines and cross-section lines has also
been extended in Version 3. The GUI now allows the display ofresults at multiple time-layers simultaneously in one chart. Thisfeature allows one to compare very easily values of a given
quantity at different times. Note that cross-section lines andcharts now can be also applied to mesh slices created by the
Clipper tool as described earlier and shown in Figure 8.
We further improved the 3D graphics of Version 3 of
HYDRUS, including allowing the rendering to be smoother andfaster. A transparent mode (object translucency) is now available also for the graphical display of results. An example of thisis shown in Figure 13. We included at the same time new options for automatic numbering of isolines and displaying positions of minimum/maximum values of a current quantity.

Fig. 8. Example showing 3D Clipping. When working with 3Dmodels, the visible mesh (or geometric objects) can be clipped toview values inside the domain. The clipper is controlled by agraphical tool (a manipulator).

Fig. 9. Example showing 3D Slicing. The slicer is another tool
for better visualization of 3D objects. The number of slices,
their position and rotation can be set either in a dialog orgraphically.

139Jiří Šimůnek, Miroslav Šejna, Martinus Th. van Genuchten

Fig. 10. Rasters for Vector Fields. Vector fields can be displayed in mesh-independent rasters, i.e., regular 2D/3D grids with an arbitrary
density of points. Rasters can also be used with the Clipper and Slicer.

Fig. 11. Example showing 3D Streamlines. Streamlines can becalculated for a given time layer (i.e., steady-state flow) from a
set of seed points. Seed point sources of several types can be set
graphically.

Fig. 12. Example of Particle Animation. The movement ofparticles can be animated as flowing-lines or flowing-particles.

Fig. 13. Example of the Transparent mode of the HYDRUS GUI. The Transparent mode is available for all types of graphical displays of
scalar and vector fields. Isolines are automatically numbered, and min/max values of a particular variable in the transport domain areautomatically identified.

140

New features of version 3 of the HYDRUS (2D/3D) computer software package

Fig. 14. Example showing Export to ParaView. ParaView is an open-source, multi-platform data analysis and visualization application.
New visualization options include custom post-processing based on scripting or plugin modules. Details are at www.paraview.org.

Although we believe that HYDRUS post-processing has allof the most important features needed for project execution,some users may require additional functions specific to theirown needs. The new version of HYDRUS offers a robust solution by exporting all results to VTK files. VTK is a well-known
open source library intended for the visualization of scientificdata. Exported HYDRUS results can then be opened in ParaView, such as illustrated in Figure 14, which is a free program
based on VTK. Since both software products (VTK and ParaView) are open source, users have full control over the exported data and can implement any special post-processing as
judged optimal for their application.
CONCLUSIONS
While HYDRUS (2D/3D) and its various predecessors, suchas UNSAT, SWMS-2D, and CHAIN-2D, have been available
to the public for over 20 years, the software package is underconstant development, with new features being continuously
added. This manuscript describes two main new capabilities (ageneral reservoir boundary condition and a dynamic root
growth module) that recently were incorporated in the standard
computational module. The capability to simulate the root
growth in HYDRUS (2D/3D) has regularly been requested forquite some time by many users, and thus may well become apopular new feature. Similarly, the “Reservoir” boundary condition offers new capabilities for application to a large numberof field problems, such as well pumping or the drying of wells,
furrow irrigation design and analysis, falling head infiltration
experiments, and a range of new approaches to land development known as Low-Impact Development (LID). We also
summarized in this paper several new capabilities of the HY-

DRUS (2D/3D) graphical user interface. Continuous updatingof the HYDRUS GUI not only guarantees that the software will
be fully operational with new Window operating systems, butalso will offer more convenient tools, options, and ease of use
to its users.
REFERENCES
Boogaard, H., van Diepen, C., Rotter, R., Cabrera, J., van Laar,
H., 1998. WOFOST 7.1: User's guide for the WOFOST 7.1
crop growth simulation model and WOFOST control center
1.5. Tech. Rep. DLO Winand Staring Ctr., Wageningen, theNetherlands.
Brunetti, G., Šimůnek, J., Piro, P., 2016. A comprehensive
analysis of the variably-saturated hydraulic behavior of agreen roof in Mediterranean climate. Vadose Zone Journal,
15, 9, 15 p. DOI: 10.2136/vzj2016.04.0032.
Brunetti, G., Šimůnek, J., Turco, M., Piro, P., 2017. On the useof surrogate-based modeling for the numerical analysis of
Low Impact Development techniques. Journal of Hydrology,
548, 263–277. DOI: 10.1016/j.jhydrol.2017.03.013.
Carsel, R.F., Parrish, R.S., 1988. Developing joint probability
distributions of soil water retention characteristics. Water Resources Research, 24, 755–769.
Coffman, L.S., 2002. Low-impact development: An alternativestormwater management technology. In: France, R.L. (Ed.):
Handbook of Water Sensitive Planning and Design. Lewis
Publ., Boca Raton, FL, pp. 97–123.
Doussan, C., Pages, L., Vercambre, G., 1998. Modelling of thehydraulic architecture of root systems: An integrated approach to water absorption - Model description. Ann. Bot.,81, 213–223.

141Jiří Šimůnek, Miroslav Šejna, Martinus Th. van Genuchten

Han, M., Zhao, C., Šimůnek, J., Feng, G., 2015. Evaluating theimpact of groundwater on cotton growth and root zone waterbalance using Hydrus-1D coupled with a crop growth model.
Agricultural Water Management, 160, 64–75. DOI:
10.1016/j.agwat.2015.06.028.
Hartmann, A., Šimůnek, J., Aidoo, M.K., Seidel, S.J., Lazarovitch, N., 2018. Modeling root growth as a function of different environmental stresses using HYDRUS. Vadose Zone
Journal, 16 p. DOI: 10.2136/vzj2017.02.0040.
Hinnell, A.C., Lazarovitch, N., Warrick, A.W., 2009. Explicit
infiltration function for boreholes under constant head conditions. Water Resour. Res., 45, W10429. DOI:
10.1029/2008WR007685.
Hoffman, G.J., van Genuchten, M.T., 1983. Soil properties andefficient water use: Water management for salinity control.
In: Taylor, H.M., Jordan, W.R., Sinclair, T.R. (Eds.): Limitations and Efficient Water Use in Crop Production. Am.
Soc. of Agron., Madison, WI, pp. 73–85.
Javaux, M., Schröder, T., Vanderborght, J., Vereecken, H.,
2008. Use of a three-dimensional detailed modeling approach for predicting root water uptake. Vadose Zone Journal, 7, 1079–1088. DOI: 10.2136/vzj2007.0115.
Javaux, M., Couvreur, V., Vanderborght, J., Vereecken, H.,
2013. Root water uptake: From three-dimensional biophysical processes to macroscopic modeling approaches. Vadose
Zone Journal, 12, 4, 16 p. DOI: 10.2136/vzj2013.02.0042.
Jones, C., Bland, W., Ritchie, J., Williams, J.R., 1991.
Simulation of root growth. Modeling Plant and Soil Systems-Agronomy Monograph, Segoe Rd., Madison, WI
53711, USA, ASA-CSSA-SSSA, 31, pp. 91–123.
Karandish, F., Šimůnek, J., 2016. A comparison of numericaland machine-learning modeling of soil water content withlimited input data. Journal of Hydrology, 543, 892–909.
DOI: 10.1016/j.jhydrol. 2016.11.007.
Lazarovitch, N., Warrick, A.W., Furman, A., Zerihun, D.,2009. Subsurface water distribution from furrows describedby moment analyses. J. Irrig. Drain. Eng., 135, 1, 7–12.
Li, Y., Kinzelbach, W., Zhou, J., Cheng, G.D., Li, X., 2012.
Modelling irrigated maize with a combination of coupledmodel simulation and uncertainty analysis, in the northwestof China. Hydrol. Earth Syst. Sci., 16, 1465–1480.
Li, Y., Zou, Q., Zhou, J., Zhang, G., Chen, C., Wang, J., 2014.
Assimilating remote sensing information into a coupled hydrology-crop growth model to estimate regional maize yieldin arid regions. Ecological Modelling, 291, 15–27.
Ramos, T.B., Šimůnek, J., Gonçalves, M.C., Martins,
J.C., Prazeres, A., Pereira, L.S., 2012. Two-dimensionalmodeling of water and nitrogen fate from sweet sorghum irrigated with fresh and blended saline waters. Agricultural
Water
Management,
111,
87–104.
DOI:
10.1016/j.agwat.2012.05.007.
Roberts, T.L., Lazarovitch, N., Warrick, A.W, Thompson, T.L.,
2009. Modeling salt accumulation with subsurface drip irrigation using HYDRUS-2D. Soil Sci. Soc. Am. J., 73, 1,
233–240. DOI: 10.2136/sssaj2008.0033.

142

Siyal, A.A., Bristow, K.L., Šimůnek, J., 2012. Minimizing
nitrogen leaching from furrow irrigation through novel
fertilizer placement and soil management strategies. Agricultural Water Management, 115, 242–251. DOI:
10.1016/j.jconhyd.2012.03.008.
Šimůnek, J., Hopmans, J.W., 2009. Modeling compensated root
water and nutrient uptake. Ecological Modeling, 220, 4,
505–521. DOI: 10.1016/j.ecolmodel.2008.11.004.
Šimůnek, J., van Genuchten, M.T., Šejna, M., 2008. Development and applications of the HYDRUS and STANMOD
software packages and related codes. Vadose Zone Journal,
7, 2, 587–600. DOI: 10.2136/VZJ2007.0077.
Šimůnek, J., van Genuchten, M. T., Šejna, M., 2016a. Recent
developments and applications of the HYDRUS computer
software packages. Vadose Zone Journal, 15, 7, 25 p. DOI:
10.2136/vzj2016.04.0033.
Šimůnek, J., Bristow, K.L., Helalia, S.A., Siyal, A.A., 2016b.
The effect of different fertigation strategies and furrow
surface treatments on plant water and nitrogen use. Irrigation
Science, 34, 1, 53–69. DOI: 10.1007/s00271-015-0487-z.
Somma, F., Hopmans, J.W., Clausnitzer, V., 1998. Transient
three-dimensional modeling of soil water and solutetransport with simultaneous root growth, root water and nutrient uptake. Plant Soil, 202, 281–293.
van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil
Sci. Soc. Am. J., 44, 987–996.
Vrugt, J.A., Hopmans, J.W., Šimůnek, J., 2001a. Calibration ofa two-dimensional root water uptake model. Soil Sci. Soc.
Amer. J., 65, 4, 1027–1037.
Vrugt, J.A., van Wijk, M.T., Hopmans, J.W., Šimůnek, J.,
2001b. One-, two-, and three-dimensional root water uptake
functions for transient modeling. Water Resources Research,
37, 10, 2457–2470.
Wang, J., Huang, G., Zhan, H., Mohanty, B.P., Zheng, J.,
Huang, Q., Xu, X., 2014. Evaluation of soil water dynamics
and crop yield under furrow irrigation with a twodimensional flow and crop growth coupled model. Agricultural Water Management, 141, 10–22.
Wang, X., Huang, G., Yang, J., Huang, Q., Liu, H., Yu, L.,2015. An assessment of irrigation and fertilization practices:
Sprinkler irrigation of winter wheat in the North China
Plain. Agricultural Water Management, 159, 197–208.
Warrick, A.W., Lazarovitch, N., Furman, A., Zerihun, D.,2007. Explicit infiltration function for furrows. J. Irrig.
Drain. Eng., 133, 4, 307–313.
Williams, J.R., Jones, C.A., Kiniry, J.R., Spanel, D.A., 1989.
The EPIC crop growth model. Trans. ASAE, 32, 2, 497–
511.
Zhou, J., Cheng, G., Li, X., Hu, B.X., Wang, G., 2012. Numerical modeling of wheat irrigation using coupled HYDRUS
and WOFOST models. Soil Sci. Soc. Amer. J., 76, 2, 648–
662. DOI: 10.2136/sssaj2010.0467.
Received 14 August 2017
Accepted 4 October 2017