Recent Developments and Applications of the HYDRUS Computer Software Packages

The HYDRUS-1D and HYDRUS (2D/3D) software packages (Šimůnek et al., 2008b) are finite-element models for simulating the one- and two- or three-dimensional movement of water, heat, and multiple solutes in variably saturated media, respectively. The standard versions, as well as various specialized add-on modules, of the HYDRUS programs numerically solve the Richards equation for saturated–unsaturated water flow and convection–dispersion type equations for heat and solute transport. The flow equation incorporates a sink term to account for water uptake by plant roots as a function of water and salinity stress. Both compensated and uncompensated water uptake by roots can be considered. The heat transport equation considers movement by both conduction and convection with flowing water. The governing convection–dispersion solute transport equations are written in a relatively general form by including provisions for nonlinear, nonequilibrium reactions between the solid and liquid phases and linear equilibrium reactions between the liquid and gaseous phases. The transport models also account for convection and dispersion in the liquid phase as well as diffusion in the gas phase, thus permitting the models to simulate solute transport simultaneously in both the liquid and gaseous phases. Hence, both adsorbed and volatile solutes, such as pesticides and fumigants, can be considered.

The solute transport equations further incorporate the effects of zero-order production, first-order degradation independent of other solutes, and first-order decay and production reactions that provide coupling between solutes involved in sequential first-order chain reactions. Typical examples of such first-order degradation chains involve radionuclides, various N species, pesticides, and many organic pollutants. Physical nonequilibrium solute transport can be accounted for by assuming a two-region, dual-porosity type formulation that partitions the liquid phase into mobile and immobile regions. Attachment and detachment processes and related filtration provisions are further included to simulate the transport of viruses, colloids, bacteria, nanoparticles, and nanotubes. Many specialized modules, to be described below, have been developed for both HYDRUS-1D and HYDRUS (2D/3D) to account for processes that cannot be handled by the standard computational modules.

In 2008, Šimůnek et al. (2008b) reviewed the early history of the HYDRUS and STANMOD software packages and related programs and modeling tools such as RETC, ROSETTA, UNSODA, UNSATCHEM, and HP1. Since then, several other HYDRUS-related reviews appeared mostly focusing on a particular version or type of application. For example, Šimůnek et al. (2012b) and van Genuchten et al. (2012) reviewed the issues of calibration and validation of the HYDRUS and STANMOD software packages, respectively, while Šimůnek et al. (2013a) reviewed various specialized add-on modules (e.g., C-Ride, DualPerm, Fumigant, and UnsatChem) developed for HYDRUS (2D/3D). The main objective of this paper is to review various new capabilities of the HYDRUS programs that have been implemented since 2008. We believe that such a review would be beneficial for the HYDRUS community, which has grown dramatically during the past several years. An additional objective is to review major types of applications of the different HYDRUS models and applicable add-ons and to briefly discuss future plans and directions.

In the text below, we use various terms such as software package, code, model, module, and program. Although at times we will use these terms interchangeably, we attempt to use them as follows. Under the terms model and module, we understand both the conceptual and mathematical description of the problem as well as its numerical implementation into a computer program. The term model is a broader term in that it includes not only the main module (e.g., HYDRUS) but also multiple standard and nonstandard modules (e.g., UnsatChem or C-Ride). Under the term program, we understand the numerical implementation of the mathematical model into a computer language. And finally, under the term software package, we understand a collection of individual files and resources (such as a GUI, help files, manuals, computational modules, and test examples) that are put together to provide certain functionality.

A major development with respect to HYDRUS-1D occurred in 2008 when Version 4 (Šimůnek et al., 2008a) was released (Table 1). Version 4 substantially enhanced the capabilities of the model compared with Version 3. Version 4.01 additionally considered vapor flow and the fully coupled transport of water, vapor, and energy (Saito et al., 2006); an option to evaluate potential evapotranspiration using the Penman–Monteith combination equation (FAO, 1990) or the Hargreaves equation (Hargreaves, 1994); an option to generate intraday variations in the evaporation and transpiration rates from their daily values; and full graphical support for the HP1 program (Jacques et al., 2008a,b).

A detailed description of additional modifications and the different new options available in various HYDRUS-1D subversions (from 4.04 to 4.17) are given in Table 1. They include options to (i) evaluate tortuosity using the models of Moldrup et al. (1997, 2000) as an alternative to the Millington and Quirk (1961) model (Version 4.06), (ii) calculate soil surface temperatures and actual evaporation fluxes for bare soils using the surface energy balance (Saito et al., 2006) (Version 4.07), (iii) provide support to the HYDRUS package for MODFLOW (Twarakavi et al., 2008) (Version 4.07), (iv) consider both uncompensated and compensated root water uptake as well as both passive and active solute uptake (Šimůnek and Hopmans, 2009) (Version 4.08), (v) use field capacity as a possible initial condition using an equation suggested by Twarakavi et al. (2009) (Version 4.16), (vi) trigger surface irrigation when a prescribed pressure head is reached at a specified depth (Dabach et al., 2013) (Version 4.16), and (vii) allow drainage fluxes to horizontal drains to occur either through the bottom of the soil profile or through a vertically distributed region in the saturated zone (Version 4.17).

Version 4 of HYDRUS-1D—similarly as Version 3—supports two add-on modules simulating (i) carbon dioxide transport and production (Šimůnek and Suarez, 1993) and major ion reactions and transport (the UnsatChem module) (e.g., Šimůnek and Suarez, 1994), and (ii) the transport and general biogeochemical reactions between many different ions (the HP1 module) (Jacques et al., 2008a,b). More details about the HP1 module are given below. Additionally, Version 4 of HYDRUS-1D supports an add-on module simulating water flow and solute transport in dual-permeability porous media (Gerke and van Genuchten, 1993). This module, contrary to two other add-on modules (i.e., UnsatChem and HP1), can be run in both direct and inverse (calibration) mode. However, external optimization tools are required to run UnsatChem and HP1 in the inverse mode (e.g., Jacques et al., 2012). While several applications of the UnsatChem and HP1 modules are described below, an overview of applications of the DualPerm module were given by Köhne et al. (2009a,b)..

In addition to the standard HYDRUS-1D add-on modules, which are fully supported by the HYDRUS-1D GUI and documented in detail in the HYDRUS-1D manuals and via online help, several additional nonstandard modules exist that can be freely downloaded from the HYDRUS website (http://www.pc-progress.com/en/Default.aspx?h1d-library) together with many examples demonstrating their use as well as brief descriptions of the theories behind the modules and their implementation. The nonstandard computational modules significantly expand the capabilities of the HYDRUS-1D software. Although they can still be run from the standard HYDRUS-1D GUI, users are usually required to manually provide an additional input file with supplementary information needed for a particular module or to interpret selected input and output variables differently from the standard versions. Users may also need to prepare their own graphical output from the output text files. Six nonstandard computational modules have been developed so far. They pertain to centrifugal forces, freeze–thaw processes, colloid-facilitated transport, colloid transport with transient water contents, isotope transport, and root growth. The nonstandard modules were developed mostly by ourselves, as well as by various colleague as part of their research, and may become standard HYDRUS modules in the future if sufficient interest exists. They are briefly described below.

1. Centrifugal Forces: This nonstandard computational module considers centrifugal forces in addition to gravitation and capillarity. Since this module can simulate, in both direct and inverse modes, water flow and solute transport in a transient centrifugal field (Šimůnek and Nimmo, 2005), it can be used to analyze data collected using high-speed centrifuges. Note that high-speed centrifugal methods during the last few decades have become relatively standard in many fields (such as in soil physics, the petroleum industry, and environmental and geotechnical engineering) for measuring saturated and unsaturated hydraulic conductivities or for studying various flow and transport processes. Example applications of this module are given by Nakajima and Stadler (2006) and van den Berg et al. (2009).

2. Freezing and Thawing: In addition to fully coupled transport of water, vapor, and energy, this nonstandard module considers the effects of freezing and thawing on water flow and solute and heat transport processes (Hansson et al., 2004). The module is not a standard in HYDRUS-1D since it runs only for unsaturated soils and becomes unstable when the medium reaches full saturation. The freezing and thawing module has been used in studies by Watanabe et al. (2007) and Kurylyk and Watanabe (2013).

3. Colloid-Facilitated Transport: This nonstandard computational module is essentially a one-dimensional version of the C-Ride add-on module of HYDRUS (2D/3D). The module considers particle transport and particle-facilitated solute transport (Šimůnek et al., 2006). The term particle is a general term used here for many substances having a relatively small but finite size (such as viruses, bacteria, pathogens, nanoparticles, and nanotubes), the transport of which is usually described using a convection–dispersion type equation with separate attachment, detachment, and straining terms. Particle-facilitated solute transport (often referred to also as colloid-facilitated transport when solutes are transported sorbed to colloids) is often observed for many strongly sorbing contaminants such as heavy metals and radionuclides. This computational module has been used by Pang and Šimůnek (2006) among others.

4. Colloid Transport with Changing Water Contents: This module can simulate particle transport, similarly as the standard HYDRUS-1D computational modules, while additionally considering the effects of changes in the water content on colloid and bacteria transport and attachment and detachment to or from solid–water and air–water interfaces (e.g., Bradford et al., 2015). For example, when the air–water interface disappears during imbibition, particles residing on this interface are released into the liquid phase. Similarly, during drainage, particles residing at the solid–water interface may be detached from this interface by capillary forces and released into the liquid phase or become attached to the air–water interface.

5. Isotope Transport: This nonstandard module is a modified version of the standard solute transport formulation to account for isotope transport (Stumpp et al., 2012). The module assumes that fractionation processes can be neglected and that the relative concentration of isotopes (their δ content) does not increase at the upper boundary because of evaporation. This is in contrast to the standard formulation during evaporation in HYDRUS-1D, where solutes concentrate at and near the soil surface when water evaporates. Water and solutes in the modified module will move at similar rates. The isotope content taken up by roots during transpiration is then equal to the soil solute concentration without having a fractionation effect (Stumpp et al., 2012). This module has been successfully used also by Stumpp and Hendry (2012), Huang et al. (2015), and Sprenger et al. (2015).

6. Root Growth: This nonstandard computational module can simulate root growth and its dependence on various environmental factors (Hartmann and Šimůnek, 2015). The root growth module is based on approaches developed by Jones et al. (1991). The model assumes that various environmental factors, characterized by growth stress factors, can influence root development under suboptimal conditions. Root growth and the development of root-length density then depend on these environmental factors when a stress factor approach is used. A similar approach was implemented in the 2D part of HYDRUS (2D/3D) (Hartmann and Šimůnek, 2015).

A detailed list of recent developments, additional modifications, and new options in various versions (1.07 to 2.05) of HYDRUS (2D/3D) is given in Table 2. A major new release occurred in 2011 when Version 2 of HYDRUS (2D/3D) with its 3D-Professional Level was made available. This version not only supports complex general three-dimensional geometries that can be designed using three-dimensional objects of general shapes but also includes multiple specialized add-on modules that significantly expand the number of processes that HYDRUS (2D/3D) can consider and which were not available with the main standard module. The add-on specialized modules (i.e., Fumigant, UnsatChem, Wetland, DualPerm, C-Ride, Slope, and Slope Cube) are described below.

A number of special boundary conditions were implemented into Version 2 of HYDRUS (2D/3D). These boundary conditions include (i) a gradient bottom boundary condition (in addition to the unit [free drainage] gradient boundary condition), (ii) a subsurface-drip boundary condition involving a drip characteristic function that reduces irrigation fluxes based on the back pressure as described by Lazarovitch et al. (2005), (iii) a surface-drip boundary condition with a dynamic wetting radius (Gärdenäs et al., 2005), (iv) a seepage-face boundary condition with a specified pressure head (to accommodate a particular suction applied at the bottom of lysimeters), and (v) a triggered-irrigation boundary condition to allow irrigation to be triggered at a specified boundary of the domain when the pressure head at a particular observation node within the domain drops below a certain value (Dabach et al., 2013).

Two- and three-dimensional applications often require a large number of finite elements to discretize large transport domains. Even with powerful personal computers currently available, it is virtually impossible to solve problems having more than about half a million nodes within a reasonable computational time. To decrease the required computational time, Hardelauf et al. (2007) parallelized an earlier three-dimensional computational module of HYDRUS (2D/3D), called SWMS_3D (Šimůnek et al., 2008b), to obtain the ParSWMS code, which distributes problems with a large number of elements over multiple processors working in parallel. While ParSWMS simulates water flow and solute transport in 3D domains, it does not consider some of the advanced features of HYDRUS such as dual-porosity systems, hysteresis, and nonlinear and nonequilibrium solute transport. The ParSWMS program was developed for the LINUX or UNIX workstations using the installed freeware MPI, PETSc, and PARMETIS. Hardelauf et al. (2007) demonstrated that an increase in the number of processors produces a proportional decrease in computational time. Although the parallelized ParSWMS program cannot be run on Windows-based PCs, since it requires LINUX or UNIX, its input and output are fully supported by the HYDRUS GUI (Version 2).

An alternative to ParSWMS is the HYPAR module (acronym for HYdrus PARallelized), which is a parallelized version of the standard 2D and 3D HYDRUS computational modules. HYPAR uses parallel programming tools to take advantage of new multicore and multiprocessor computers to significantly speed up time-consuming simulations, especially those requiring a large number of finite elements. HYPAR currently supports only calculations in a direct (forward) mode but not inverse (parameter estimation) computations. HYPAR similarly does not support any specialized add-on modules such as HP2, UnsatChem, Wetland, and C-Ride.

Several completely new specialized add-on modules have been developed gradually for Version 2 of HYDRUS (2D/3D) to account for various processes not available in the standard software package. These new modules include the HP2, C-Ride, DualPerm, UnsatChem, Wetland, and Fumigant modules. All of these modules simulate water flow and various solute transport processes in two-dimensional, variably saturated transport domains and are fully supported by the HYDRUS GUI. Many of the processes included in these specialized modules of HYDRUS (2D/3D) are currently also available as part of HYDRUS-1D (described above).

1. The C-Ride Module: The C-Ride module simulates the transport of particle-like substances (e.g., colloids, viruses, bacteria, and nanoparticles) as well as considers particle-facilitated solute transport (Šimůnek et al., 2006). Particle-facilitated transport is often observed for many strongly sorbing contaminants such as heavy metals, radionuclides, pharmaceuticals, pesticides, and explosives (see references in Šimůnek et al., 2006). These contaminants are predominantly associated with the solid phase, which is commonly assumed to be stationary. However, they may also sorb or attach to mobile and deposited (colloidal) particles such as microbes, humic substances, suspended clay particles, and metal oxides, which then can act as pollutant carriers and hence provide a rapid-transport pathway for the pollutants. The C-Ride module fully accounts for the dynamics of particles themselves (e.g., attachment and straining) as well as for solute transfer between different phases such as kinetic or equilibrium sorption to the soil phase and kinetic sorption to mobile or deposited colloids. A schematic of the particle-facilitated solute transport model as implemented into C-Ride is shown in Fig. 1.

2. The DualPerm Module: The DualPerm module simulates preferential and nonequilibrium water flow and solute transport in dual-permeability media using the approach suggested by Gerke and van Genuchten (1993). The module assumes that the porous medium consists of two interacting and overlapping regions: one associated with the interaggregate, macropore, or fracture system and one consisting of micropores (or intra-aggregate pores) inside soil aggregates or within the soil or rock matrix. Water flow can occur in both regions, albeit at different rates. We note that this module cannot be applied to systems involving discrete fracture and macropore networks. Modeling details are provided by Šimůnek and van Genuchten (2008). Many applications of this HYDRUS (2D/3D) module, as well as of the corresponding 1D module, are given by Köhne et al. (2009a,b). Figure 2 shows an example for the infiltration of water from a tension disc infiltrometer (having a disc radius of 10 cm) into a 50-cm-wide and 150-cm-deep soil domain. Shown are calculated pressure head profiles in the matrix and fracture domains for different ratios of the anisotropy hydraulic conductivity coefficients (i.e., K_x^A/K_z^A = 1, 10, and 0.1).

3. The UnsatChem Module: The geochemical UnsatChem module has been implemented into all 1D, 2D, and 3D HYDRUS versions. UnsatChem considers the transport of major ions (i.e., Ca²⁺, Mg²⁺, Na⁺, K⁺, SO₄²⁻, CO₃²⁻, and Cl⁻) in conjunction with most or all relevant equilibrium and kinetic geochemical reactions such as complexation, cation exchange, and precipitation and dissolution (e.g., of calcite, gypsum, and dolomite). Table 3 lists the various chemical species considered in UnsatChem. Possible applications of this module include studies evaluating the sustainability of alternative irrigation systems, salinization and reclamation of agricultural soils, and the disposal of brine waters from mining operations (e.g., oil and gas production, shale fracking, or coal seam fracking). Ever since its introduction some two decades ago (Šimůnek and Suarez, 1994), the Unsatchem module (especially its 1-D version) has been used widely in many applications (as described in the Salinization and Sodification section).

4. The Wetland Module: The Wetland module simulates aerobic, anoxic, and anaerobic transformation and degradation processes for organic matter, N, P, and S during treatment of polluted wastewater in subsurface constructed wetlands (Langergraber and Šimůnek, 2012). Constructed wetlands are engineered water treatment systems that optimize the treatment processes taking place in natural environments. They have become popular since they can be very efficient in treating different types of polluted water using sustainable, environmentally friendly approaches. A large number of physical, chemical, and biological processes are simultaneously active and may mutually influence each other in constructed wetlands. The Wetland module uses two biokinetic model formulations to account for complex conditions that may occur in various types of wetlands: CW2D of Langergraber and Šimůnek (2005) for aerobic and anoxic conditions and CWM1 of Langergraber et al. (2009), which also considers anaerobic conditions. The two Wetland modules were tested by Pálfy and Langergraber (2014) and Pálfy et al. (2015). Additional references of Wetland module applications can be found at http://www.pc-progress.com/en/Default.aspx?h3d2-wetland.

5. The Fumigant Module: The Fumigant module implements multiple additional options for simulating processes related to the application and subsurface transport of fumigants, which are not available in the standard HYDRUS models. This module allows users to specify additional injections of fumigants into the transport domain at a specific location at a specific time as well as to consider the presence or absence of a surface tarp, the temperature dependence of tarp properties, and the removal of tarp at a certain time. The Fumigant module has been used recently to investigate the effects of different application scenarios (such as tarped broadcast, tarped bedded shank injection, or tarped drip line-source application) and various factors (such as initial water content or tarp permeability) on fumigant volatilization (Nelson et al., 2013; Spurlock et al., 2013a,b). Figure 3 summarizes one example.

6. The Classic Slope Module: One frequent application of HYDRUS has been to obtain subsurface flow conditions (i.e., relative saturations and water fluxes) for subsequent slope-stability analyses using other programs. This motivated us to develop the Slope Stability (SLOPE) add-on module intended mainly for stability tests of embankments, dams, earth cuts, and anchored-sheeting structures. The influence of water is modeled using the distribution of pore pressure, which is imported automatically from HYDRUS runs into the SLOPE module at specified times, each of which can be analyzed separately. The slip surface in the SLOPE module is considered to be circular and is evaluated using the Bishop, Fellenius–Petterson, Morgenstern–Price, or Spencer method (Lu and Godt, 2013). More details can be found in the user manual of this module.

7. The Slope Cube Module: While the SLOPE module is based on classical engineering soil mechanics theories and uses the effective stress approach only for saturated conditions, a new add-on module, SLOPE Cube (slope, stress, and stability) was recently developed to provide a unified effective stress approach for both saturated and unsaturated conditions (Lu et al., 2010). The module is intended to predict spatially and temporally infiltration-induced landslide initiation and to carry out slope stability analyses under variably saturated soil conditions. Transient moisture and pressure-head fields are directly obtained from the HYDRUS-2D model and subsequently used to compute the effective stress field of hillslopes (Lu and Godt, 2013). Furthermore, instead of the methodology of one-slope-for-one-factor safety in the classical slope stability analysis, the SLOPE Cube module computes fields of the factor of safety in the entire domain within hillslopes (Lu et al., 2012), thus allowing identification of the development of potential failure surface zones or surfaces.

As with HYDRUS-1D, several additional nonstandard computational HYDRUS (2D/3D) add-on modules were developed that are not fully supported by HYDRUS (2D/3D) nor have they been fully documented. These nonstandard add-on modules can again be downloaded from the HYDRUS website (http://www.pc-progress.com/en/Default.aspx?h3d-applications) together with many examples demonstrating their use and a brief description of the theory behind the modules and their implementation. As with HYDRUS-1D, the nonstandard computational modules can still be run from the standard HYDRUS (2D/3D) GUI, but users are usually required to provide an additional input file with supplementary information needed by a particular module or to interpret various input and output variables differently. Three such nonstandard computational add-on modules have been developed thus far:

1. Centrifugal Forces: This nonstandard computational module deals with centrifugal forces in addition to gravitational and capillary forces. In addition to considering processes in 2D transport domains, this module has similar capabilities as the corresponding HYDRUS-1D module (Šimůnek and Nimmo, 2005) as explained above.

2. Overland Flow: The Overland Flow nonstandard module can consider, in addition to subsurface flow and transport, overland flow and transport processes. While the standard HYDRUS modules assume that once the infiltration capacity is exceeded, any excess water is instantaneously removed by surface runoff, this module considers flow of this excess water along the soil surface. The module can account for overland flow (runoff) once the soil infiltration capacity has been reached, can redistribute water on the land surface by moving it to lower parts of a hillside where the water could infiltrate if the local soil infiltration capacity has not been reached, or where it can remain as runoff. While subsurface flow is still described using the Richards equation, overland flow is simulated using the kinematic wave equation (Köhne et al., 2011).

3. Carbon Dioxide Transport and Production: This nonstandard module extends the capabilities of the 2D UnsatChem module discussed earlier. While the standard version of UnsatChem assumes that the spatial distribution of CO₂ concentrations is constant in time (contrary to the 1D UnsatChem model, which considers transient CO₂ transport), this specialized nonstandard module can also simulate CO₂ transport and production (Šimůnek and Suarez, 1993). The module accounts for diffusion of CO₂ in both liquid and gas phases, CO₂ production, and uptake of CO₂ by plant roots. The CO₂ production model considers both microbial and root respiration, which are dependent on water content, temperature, and plant and soil characteristics. The new module was developed so that it can be run using the HYDRUS (2D/3D) graphical user interface, similarly as all other standard and nonstandard add-on modules.

While the 3D-Layered Level of HYDRUS can support only layered geometries that are built above a two-dimensional base, the 3D-Professional Level supports complex general three-dimensional geometries that can be formed from three-dimensional objects (solids) having very general shapes. Three-dimensional objects are formed by boundary surfaces that can be both planar surfaces and curved surfaces (quadrangles, rotaries, pipes, or B-splines). Figure 4 shows examples of various curved surfaces, while Fig. 5 shows how these individual objects can be combined to form complex 3D geometries.

Various spatially variable properties (such as materials, initial conditions, boundary conditions, and domain properties) can be specified in Version 2.0 of HYDRUS either directly on the FEM, as done previously also in Version 1.0, or independently of the FEM on geometric objects (e.g., boundary curves, rectangles, circles, surfaces, solids) as shown in Fig. 6. The main advantage of the latter approach is that when the FEM is changed (e.g., when convergence is not achieved for a given FEM), these properties are not automatically lost but can be reassigned immediately to the new FEM from their initial definition on geometric objects.

Many other improvements were implemented into Version 2.0 of HYDRUS (2D/3D) to make the program easier to use. Particularly useful are options to (i) import domain properties and initial and boundary conditions from existing HYDRUS projects even from projects with a (slightly) different geometry or FEM, (ii) import geometric objects using a variety of formats (e.g., .TXT, .DXF, .SHP), and (iii) display results using isosurfaces. Table 2 lists several other options.

The one-dimensional program HP1 (Jacques et al., 2008a,b), which couples the PHREEQC geochemical program (Parkhurst and Appelo, 1999) with HYDRUS-1D, has been used successfully in many applications since its release in 2005 (see HP1 and HP2 Applications section below). The two-dimensional extension, HP2, was released in 2013 as an add-on module to HYDRUS (2D/3D) (Šimůnek et al., 2012a). HPx, which is an acronym for HYDRUS-PHREEQC-xD (1D or 2D), is a relatively comprehensive simulation module that can be used to simulate (i) transient water flow, (ii) the transport of multiple components, (iii) mixed equilibrium and kinetic biogeochemical reactions, and (iv) heat transport in one- and two-dimensional variably saturated porous media. The HP1 and HP2 modules are suitable for a broad range of low-temperature biogeochemical reactions in water, the vadose zone, and ground water systems including interactions with minerals, gases, exchangers, and sorption surfaces based on thermodynamic equilibrium, kinetic, or mixed equilibrium–kinetic reactions.

HP1 and HP2 both allow thermodynamic equilibrium calculations for multiple chemical reactions and other features such as (i) aqueous speciation with different activity correction models (Davies, extended Truesdell–Jones, B-Dot, Pitzer, and SIT-Specific Ion Interaction Theory), (ii) multisite ion exchange sites with exchange described using different models (Gaines–Thomas, Vanselow, or Gapon), (iii) multisite surface complexation sites with a nonelectrostatic, the Dzombak and Morel or CD_MUSIC models and different options to calculate compositions of the diffuse double layer, (iv) mineralogical assemblages, (v) solid–solutions, and (vi) gas exchange. Kinetic calculations can be used to describe mineral dissolution and precipitation, nonequilibrium sorption processes, biogeochemical reactions including first-order degradation networks, Monod kinetics, and Michaelis–Menten kinetics.

Recent additions to the capabilities of HP1 are (i) diffusion of components (e.g., O₂ or CO₂) in the gas phase and (ii) an option to change the hydraulic and solute transport properties as a function of evolving geochemical state variables. For example, precipitation and dissolution may lead to changes in porosity and corresponding changes in the soil water retention and hydraulic conductivity functions. Similarly, bacterial growth and clogging can affect porosity and corresponding physical properties. HP1 makes it possible to account for changes in (i) the porosity (and hence the saturated water content), (ii) the hydraulic conductivity, (iii) a scaling factor for the pressure head, (iv) aqueous- and gas-phase pore-geometry factors for calculating pore diffusion coefficients, (v) the dispersivity, (vi) the thermal capacity, (vii) the thermal conductivity, and (viii) the thermal dispersivity. HP1 does not require any predefined conceptual or mathematical model to update the flow and transport parameters but rather uses the flexibility of the embedded BASIC interpreter for this purpose. This permits software users to define any user-specific relationship between the geochemical state variables and the transport properties (Jacques et al., 2013).

The HYDRUS Package for MODFLOW was developed by Twarakavi et al. (2008) to account for water fluxes into and through the vadose zone in conjunction with the three-dimensional modular finite-difference ground water model MODFLOW (Harbaugh et al., 2000). The package for MODFLOW consists of two submodels that interact in space and time (Fig. 7): (i) the HYDRUS submodel for flow in the vadose zone and (ii) the MODFLOW submodel for ground water flow. The HYDRUS package considers all of the main processes and factors affecting fluxes in the vadose zone as incorporated in HYDRUS-1D such as precipitation, infiltration, evaporation, redistribution, capillary rise, plant water uptake, water accumulation on the soil surface, surface runoff, and soil moisture storage. Being fully incorporated into the MODFLOW program, the HYDRUS package provides MODFLOW with recharge fluxes into groundwater, while MODFLOW provides HYDRUS with the position of the groundwater table that is used as the bottom boundary condition. The performance of the HYDRUS package was analyzed by Twarakavi et al. (2008) for various case studies involving different spatial and temporal scales. The package has been used in several studies including Deme (2011) and Leterme et al. (2013).

The different versions of HYDRUS models have been used over the years for a large number of applications. We refer to the HYDRUS web site (http://www.pc-progress.com/en/Default.aspx) for an extensive list of various examples. The list currently contains over 850 and 550 references of HYDRUS-1D and HYDRUS (2D/3D) applications, respectively. The types of applications are very broad, ranging between agricultural problems evaluating different irrigation schemes, the effects of plants on the soil water balance and groundwater recharge (see Agricultural Applications section below), to many environmental applications simulating the transport of different solutes and particle-like substances (see Transport of Particle-Like Substances section) as well as evaluating the effects of land use and environmental changes. While many early applications focused mostly on subsurface flow processes, the relatively general formulation of the transport and reaction terms in the HYDRUS models makes it possible to simulate the fate and transport of many different solutes including nonadsorbing tracers, radionuclides (e.g., Pontedeiro et al., 2010; Matisoff et al., 2011; Merk, 2012; Xie et al., 2013), mineral N species (e.g., Li et al., 2015), pesticides (Pot et al., 2005; Dousset et al., 2007; Köhne et al., 2009b), chlorinated aliphatic hydrocarbons (e.g., Kasaraneni et al., 2014; Ngo et al., 2014), hormones (e.g., Casey et al., 2005; Arnon et al., 2008; Chen et al., 2013), antibiotics (e.g., Wehrhan et al., 2007; Unold et al., 2009; Chu et al., 2013; Engelhardt et al., 2015), explosives and propellants (e.g., Dontsova et al., 2006, 2009; Alavi et al., 2011), as well as many particle-like substances such as viruses, colloids, bacteria, nanoparticles, and carbon nanotubes (see Transport of Particle-Like Substances section).

An important advantage of the HYDRUS models is that they are not limited to any particular spatial or temporal scale. HYDRUS-1D has been applied to scales involving very short laboratory soil columns, soil profiles of one to several meters deep (e.g., Ramos et al., 2011, 2012; Neto et al., 2016), as well as to soil profiles several hundred meters deep (Scanlon et al., 2003). HYDRUS (2D/3D) has been used similarly for transport domains ranging from <1 m wide to transects of several tens or hundreds of meters wide and for both laboratory (e.g., Rühle et al., 2013, 2015) and field-scale applications (e.g., Yakirevitch et al., 2010; Pachepsky et al., 2014). Still, we do not recommend HYDRUS for very large 3D domains such as entire catchments (Šimůnek et al., 2012b). Solutions of the Richards equation require relatively fine spatial discretizations especially where and when large pressure gradients may occur such as at and near the soil surface where variable climatological conditions may cause steep gradients in the pressure head. Spatial discretizations of even a relatively small catchment can quickly lead to a finite element mesh (FEM) containing millions of nodes, thus impacting available computational resources. By comparison, no inherent limitations exist for the temporal scale, which can be very short for small-scale laboratory flow studies to hundreds of thousands of years for studies evaluating the effects of the past and current climate (e.g., Scanlon et al., 2003; Leterme et al., 2012) or for long-term environmental risk analyses of radioactive contaminants (Pontedeiro et al., 2010) provided that the material properties, such as soil hydraulic and transport properties, remain constant during the simulation.

A very common use of the HYDRUS models is for inverse estimation of soil hydraulic, solute transport, and heat transport parameters from measured steady-state or transient data. Both HYDRUS models implement a Marquardt–Levenberg-type parameter estimation technique (Marquardt, 1963; Šimůnek and Hopmans, 2002) in such a way that almost any application that can be run in a direct mode (i.e., when all parameters and initial and boundary conditions are specified and predictions are made) can be run equally well in the inverse mode. Hence, the models are effective for various model calibration and parameter estimation applications (Šimůnek et al., 2012b). Because of its generality, the inverse option in HYDRUS has proved to be very popular with many users, leading to a large number of applications. Model calibration and inverse parameter estimation can be performed with the HYDRUS software packages using either a relatively simple, gradient-based, local-optimization approach based on the Marquardt–Levenberg method, which is directly implemented into the HYDRUS models, or more complex global optimization methods (e.g., Vrugt, 2016), which need to be run separately of HYDRUS. We refer readers to a recent review of various HYDRUS applications for model calibration and parameter estimation by Šimůnek et al. (2012b).

It is beyond the scope of this paper to list all possible applications of the HYDRUS models. The breadth of applications is much larger than we expected when we initially started developing the models some 25 yr ago. We briefly note here several different types of applications that are not further discussed below. One is the hydrologic performance of green-roof systems using HYDRUS-1D (Hilten et al., 2008) or HYDRUS (2D/3D) (Palla et al., 2009; Li and Babcock, 2015; Charpentier, 2015; Brunetti et al., 2016). Water flow in highly heterogeneous waste rock piles was evaluated by Fala et al. (2005), Buczko and Gerke (2006), Dawood and Aubertin (2014), and Namaghi et al. (2014). Abramson et al. (2014a,b) further used HYDRUS in a decision support system to investigate the costs and benefits of groundwater access and abstraction for non-networked rural supplies. In yet other studies, Hassan et al. (2008), Finch et al. (2008), Sinclair et al. (2014), and Morrissey et al. (2015) modeled effluent distributions and possible groundwater pollution problems from on-site wastewater treatment systems. We refer to the HYDRUS website for a more complete list of applications.

Agricultural applications of the HYDRUS modules often involve evaluations of various irrigation schemes (e.g., Cote et al., 2003; Ben-Gal et al., 2004; Gärdenäs et al., 2005; Dabach et al., 2013), studies of root water uptake and groundwater recharge (e.g., Turkeltaub et al., 2014; Neto et al., 2016), and the transport of agricultural contaminants (Wehrhan et al., 2007; Unold et al., 2009; Engelhardt et al., 2015). For example, Gärdenäs et al. (2005) used HYDRUS (2D/3D) to evaluate water and N leaching scenarios for three different microirrigation systems (surface and subsurface drip and sprinkler irrigation) and five different fertigation strategies. Siyal et al. (2012) and Šimůnek et al. (2016) similarly used HYDRUS (2D/3D) to evaluate the effect of alternative fertigation strategies and furrow surface treatments on plant water and N use. Li et al. (2014a, 2015) and Dash et al. (2015) used HYDRUS-1D to assess water flow processes and the N balance of a rice paddy field. Others used the HYDRUS models to evaluate the effects of various irrigation practices on soil salinization and sodification (e.g., Corwin et al., 2007; Hanson et al., 2008; Ramos et al., 2011). In the section below, we briefly review applications of the HYDRUS models to drip and furrow irrigation practices, irrigation and soil salinization problems, and groundwater recharge. The examples are included here to show the wide spectrum of applications that are possible with the HYDRUS models. We again refer to the HYDRUS website for many other applications.

Modeling surface or subsurface drip irrigation has been a popular application of HYDRUS (2D/3D). Ever since Skaggs et al. (2004) successfully compared HYDRUS-2D simulations of drip irrigation with experimental observations, the model has been found helpful for evaluating soil water content patterns around drip emitters. Using ISI’s Web of Science, we identified more than 80 manuscripts (listed at http://www.pc-progress.com/en/Default.aspx?h3d-references) in which HYDRUS (2D/3D) was used to simulate drip and trickle irrigation. While the emitters in some studies were simulated as equivalent line sources (Skaggs et al., 2004), other studies considered the emitters to be a point source (Lazarovitch et al., 2009a; Kandelous and Šimůnek, 2010). Kandelous et al. (2011) discussed under what conditions drip emitters can be represented as a point source in an axisymmetrical 2D domain, a line source in a planar 2D domain, or a point source in a fully 3D domain (Fig. 8). They concluded that an axisymmetric 2D representation could be used only before wetting patterns start to overlap and a planar 2D model only after the wetting fronts from neighboring emitters fully merged. Only a 3D model could describe subsurface drip irrigation in its entirety.

HYDRUS has also been used to verify various analytical and empirical models for estimating the position of a wetting front with time, which is useful for designing or operating drip irrigation systems (Cook et al., 2006; Warrick and Lazarovitch, 2007; Lazarovitch et al., 2009a; Hinnell et al., 2010; Kandelous and Šimůnek, 2010). The effects of emitter rate, pulsing, and antecedent water content on water distribution patterns were studied by Skaggs et al. (2010). Dabach et al. (2015) evaluated optimal tensiometer placement for high-frequency subsurface drip irrigation management in heterogeneous soils. The effects of high-frequency pulsing of drip irrigation in heterogeneous soils were also studied by Assouline et al. (2006) and Mubarak et al. (2009).

Soil water and salinity distributions under different treatments of drip irrigation were simulated by Hanson et al. (2008, 2009), Roberts et al. (2008, 2009), Shan and Wang (2012), Selim et al. (2012, 2013), and Phogat et al. (2014), among others. Still others used the HYDRUS models to evaluate N leaching for different fertigation strategies using drip irrigation (Li et al., 2004, 2005; Gärdenäs et al., 2005; Hanson et al., 2006; Ajdary et al., 2007).

The HYDRUS (2D/3D) software, and its predecessors such as SWMS-2D and HYDRUS-2D, also has been used widely to simulate water flow and solute transport in furrow irrigation systems. We identified more than 25 papers addressing these topics (e.g., Benjamin et al., 1994; Abbasi et al., 2003a,b, 2004; Rocha et al., 2006; Wöhling et al., 2004a,b, 2006; Mailhol et al., 2007; Warrick et al., 2007; Wöhling and Schmitz, 2007; Wöhling and Mailhol, 2007; Crevoisier et al., 2008; Lazarovitch et al., 2009b; Ebrahimian et al., 2012, 2013a,b; Siyal et al., 2012; Zerihun et al., 2014; Šimůnek et al., 2016). A more complete list is given at http://www.pc-progress.com/en/Default.aspx?h3d-references. Still, we note that the HYDRUS models as such only consider processes in the subsurface and not overland flow. Hence, when a two-dimensional soil profile perpendicular to the actual furrow is considered, they cannot fully account for flow in the third dimension such as the advance and recession of water in a furrow. A full three-dimensional model that accounts for surface fluxes in the furrow and subsurface flow processes is required to fully describe complex three-dimensional furrow-irrigated systems (e.g., Wöhling et al., 2004b, 2006; Wöhling and Schmitz, 2007; Wöhling and Mailhol, 2007; Zerihun et al., 2014).

A typical early application of HYDRUS-2D to furrow irrigation is given by Benjamin et al. (1994) who simulated fertilizer distributions in the soil profile following broadcast fertilization using conventional- and alternate-furrow irrigation. Abbasi et al. (2003a,b, 2004) in later studies obtained close agreement between measured and predicted soil water contents and solute concentrations along a blocked-end furrow cross-section using HYDRUS-2D. Mailhol et al. (2007) and Crevoisier et al. (2008) similarly found good results with HYDRUS-2D when simulating pressure heads, nitrate concentrations, and N leaching in seasonal studies of conventional- and alternate -furrow irrigated systems while including both root water and nutrient uptake. Rocha et al. (2006) further performed a sensitivity analysis to investigate the effects of different soil hydraulic properties on flow processes below furrows.

In related work, Wöhling and Schmitz (2007) developed a numerical program that coupled HYDRUS-2D with a 1D surface flow and a crop growth model. Their code was used to predict advance and recession times, soil water contents, and crop yield (Wöhling and Mailhol 2007). Ebrahimian et al. (2012) subsequently used the HYDRUS-1D and HYDRUS-2D models to simulate water flow and nitrate transport processes following conventional-furrow irrigation, fixed alternate-furrow irrigation, and variable alternate-furrow irrigation using different fertigation strategies. Ebrahimian et al. (2013a,b) similarly used the 1D surface and 2D subsurface models to study scenarios that could minimize nitrate losses in two different alternate-furrow fertigation systems.

In a more recent study, Šimůnek et al. (2016) developed a furrow-irrigation submodule for HYDRUS (2D/3D) (Fig. 9) to evaluate the effects of furrow soil surface treatment and fertigation timing on root water and solute uptake, deep drainage, and solute leaching in a loamy soil. Simulations showed that although more water was lost by evaporation in treatments with plastic placed along the furrow bottom than in the control treatments, more water was available for transpiration and less water was drained from the soil profile for these treatments. While some of the above studies involved only simulations (e.g., Rocha et al., 2006; Warrick et al., 2007; Lazarovitch et al., 2009b), several used HYDRUS (2D/3D) to calibrate and test predictions against experimental data (e.g., Abbasi et al., 2003a,b, 2004; Wöhling and Mailhol 2007; Crevoisier et al., 2008; Zerihun et al., 2014), thus providing confidence that the model can adequately describe these complex systems.

Saline waters are used often for irrigating agricultural crops in regions having limited water resources, thus potentially causing salinization and sodification of irrigated agricultural lands. Efficient irrigation and leaching management practices are critical in these regions to prevent or limit soil salinization when rainfall is not sufficient to leach accumulated salts during or following irrigation. The HYDRUS models have been used in several studies to evaluate the sustainability of various irrigation schemes with respect to salinization and sodification processes, to assess reclamation of saline or sodic soils, and to evaluate the movement of salts after the accidental release (or possible beneficial application) of saline waters resulting from mining operations (e.g., Jakubowski et al., 2014). Such problems can be addressed with HYDRUS using two approaches. One would be to use the standard HYDRUS models by assuming that salinity behaves more or less like an inert tracer and hence is now subject to chemical reactions (e.g., Hanson et al., 2008; Dudley et al., 2008; Roberts et al., 2009; Groenveld et al., 2013). An alternative is to use the UnsatChem module, which considers the transport and reactions between major ions (e.g., Gonçalves et al., 2006; Ramos et al., 2011). While the former approach does not permit such processes as cation exchange, dissolution of mineral amendments (e.g., gypsum or calcite), or precipitation of these minerals when the soil solution becomes oversaturated; the latter approach allows one to consider those geochemical processes and the effects of salts and soil water quality on soil properties.

The UnsatChem modules (especially its 1D version) has been used in many applications as exemplified by Kaledhonkar and Keshari (2006), Kaledhonkar et al. (2006, 2012), Schoups et al. (2006), Skaggs et al. (2004), Corwin et al. (2007), and Rasouli et al. (2013), among others. Gonçalves et al. (2006) and Ramos et al. (2011, 2012) demonstrated the applicability of these modules to simulating multicomponent major ion transport in soil lysimeters irrigated with waters of different quality. While Gonçalves et al. (2006) used the UnsatChem module of HYDRUS-1D (Fig. 10), Ramos et al. (2011) used both the standard HYDRUS-1D and UnsatChem modules. Ramos et al. (2011) compared results obtained with the two modules and discussed their respective advantages and disadvantages. For example, the UnsatChem module requires much more input information (e.g., the solution composition of irrigation waters and Gapon exchange constants for all soil horizons) and runs much slower (∼20 times) than the standard HYDRUS-1D model. While both HYDRUS-1D modules were used by Ramos et al. (2011) to describe field data of the water content and overall salinity as expressed in terms of electrical conductivity (EC), the UnsatChem module was additionally used to describe the concentrations of individual soluble cations as well as of the Na adsorption ratio and the exchangeable Na percentage. Whereas EC values were calculated using different methodologies (treated as a nonadsorbing tracer in the standard module and calculated from concentrations of individual ions in UnsatChem), the two modules produced very similar results during the irrigation seasons. The main differences were found when soil water contents decreased significantly below field capacity, in which case the standard HYDRUS transport module simply increased EC linearly as the soil dried out, while the UNSATCHEM module produced a nonlinear increase in EC as a result of cation exchange (Ramos et al., 2011). Larger differences in EC values predicted with the two modules would have been observed if the soil solution had become oversaturated with respect to calcite and gypsum.

Important conclusions about the practical implications of salinity management were obtained in several studies such as by Corwin et al. (2007) and Hanson et al. (2008). Corwin et al. (2007) used the UnsatChem module to demonstrate that leaching requirements would be lower when estimated with a transient modeling approach than when using a more standard steady-state approach. Adapting leaching requirements based on the transient approach would produce significant savings in terms of irrigation water volumes. Hanson et al. (2008) showed that while the conventional or water balance approach for estimating leaching fractions predicts little or no leaching when applied water levels are less than potential evapotranspiration, field data and HYDRUS modeling showed considerable leaching around the drip lines. The spatially varying soil wetting patterns that occur during drip irrigation causes localized leaching near the drip lines (Hanson et al., 2008), thus allowing for more profitable production of various crops (e.g., processing tomato [Solanum lycopersicum L.]) than with other irrigation methods.

The HYDRUS models now include a relatively comprehensive macroscopic root water and solute uptake module (Šimůnek and Hopmans, 2009) to account for both water and salinity stress effects on water uptake while also accounting for possible active and passive root solute uptake. Root water and solute uptake furthermore can be treated as being either noncompensated or compensated.

HYDRUS-1D allows users to externally prescribe a time-variable rooting depth either using the logistic growth function or in a tabulated form. Such a feature is currently not available in HYDRUS (2D/3D), which forces the spatial distribution of roots in the root zone to remain constant during the simulations. Both models also do not allow the spatial extent of the rooting zone to change actively as a result of environmental stresses. To overcome these deficiencies, several studies either further modified the HYDRUS models (or their predecessors such as CHAIN-2D or SWMS-3D) or coupled the models with various crop-growth or root-growth models. For example, Javaux et al. (2008, 2013) developed R-SWMS, a three-dimensional root-growth model that couples the model of Somma et al. (1998) (based on SWMS-3D) with the model of Doussan et al. (1998).

For these same reasons, Zhou et al. (2012) coupled HYDRUS-1D with the WOFOST (Boogaard et al., 1998) crop-growth model and used the resulting model to simulate the growth and yield of irrigated wheat (Triticum aestivum L.) and maize (Zea mays L.) (Li et al., 2012, 2014b). Han et al. (2015) similarly coupled HYDRUS-1D with a simplified crop-growth version used in the Soil and Water Assessment Tool (SWAT) to simulate the contribution and impact of groundwater on cotton (Gossypium hirsutum L.) growth and root zone water balance. Wang et al. (2014a, 2015b) coupled the crop-growth model EPIC (Williams et al., 1989) with CHAIN-2D and HYDRUS-1D to assess the effects of furrow and sprinkler irrigation, respectively, on crop growth. Hartmann and Šimůnek (2015) furthermore implemented into both HYDRUS-1D and HYDRUS (2D/3D) the root-growth model developed by Jones et al. (1991). Their model assumes that various environmental factors as characterized by growth stress factors can influence root development under suboptimal conditions.

The governing convection–dispersion solute transport equations, as solved numerically in the HYDRUS models, allow consideration of kinetic attachment and detachment processes of particle-like substances to the solid phase. The term particle-like substance is used to represent colloids, viruses, pathogens, bacteria, nanoparticles, nanotubes, and related constituents, whose subsurface transport is often modeled using the convection–dispersion equation with certain attachment, detachment, and straining terms. This approach is used widely even though the various constituents can have dramatically different shapes and sizes with sizes varying from nanometers to micrometers. Modeling their transport represents one of the most popular applications of the HYDRUS models. We identified more than 80 manuscripts in which HYDRUS was used for simulating the transport of particle-like substances (see references at http://www.pc-progress.com/en/Default.aspx?h1d-lib-bacteria).

The particle-transport option was used first by Schijven and Šimůnek (2002) to simulate the transport of viruses at the field scale using both HYDRUS-1D and HYDRUS-2D. They modified the models by including two kinetic attachment–detachment processes involving two different sorption sites and then used the programs to simulate the removal of bacteriophages MS2 and PRD1 by dune recharge and deep-well injection. Many others since then have used the HYDRUS models to simulate virus transport in laboratory columns (e.g., Torkzaban et al., 2006a,b; Zhang et al., 2012; Frohnert et al., 2014) as well as field systems (e.g., Schijven et al., 2013).

The HYDRUS models have been used similarly as tools to understand and predict various complexities of colloid and microbial transport in the subsurface under different conditions. For example, Bradford et al. (2002, 2003, 2004) evaluated the effects of attachment, straining, and exclusion on the fate and transport of colloids in saturated porous media. Gargiulo et al. (2007a,b, 2008) evaluated the effects of such factors as matrix grain size, water content, metabolic activity, and surface proteins on bacterial transport and deposition in saturated and unsaturated media. Torkzaban et al. (2010) and Bradford et al. (2012) additionally evaluated the effects of dynamic changes in the solution ionic strength on the transport and release of colloids and microorganisms in soils. Bradford et al. (2015) further considered the effects of changing water contents on E. coli D21g transport and attachment and detachment to or from solid–water and air–water interfaces. We emphasize that at present a specialized nonstandard HYDRUS-1D module must be used to consider the effects of changes in solution chemistry and water contents on the transport and release of colloids (see Nonstandard Modules section)

The HYDRUS models are increasingly being used also to simulate the fate and transport of various nanoparticles and nanotubes in the environment. For example, Liang et al. (2013a,b), Ren and Smith (2013), Cornelis et al. (2013), Neukum et al. (2014), and Wang et al. (2015a) evaluated the sensitivity of the transport and retention of stabilized silver nanoparticles to various physicochemical factors in column studies and undisturbed soil. Kasel et al. (2013a,b) and Mekonen et al. (2014) evaluated the effects of input concentration, grain size, and saturation on the transport of multiwalled carbon nanotubes. Such studies are important for providing new knowledge about the processes affecting the environmental fate of particle like substances, which in turn allows us to continuously update the HYDRUS models.

As discussed below (HYDRUS Applications Published in VZJ), the HYDRUS models are often used (we identified 34 papers) in studies involving the use of various geophysical methods including electrical resistivity tomography (ERT), ground-penetrating radar (GPR), cosmic-ray sensing (CRS), and electric magnetic resonance. For example, electrical resistivity surveys and HYDRUS modeling were used by Batlle-Aguilar et al. (2009) to investigate axisymmetrical infiltration patterns and by Lehmann et al. (2013) to observe the evolution of soil wetting patterns preceding a hydrologically induced landslide. A large number of studies involved the complementary use of HYDRUS modeling and GPR data (e.g., Laloy et al., 2012; Jadoon et al., 2012; Scholer et al., 2013; Moghadas et al., 2013; Busch et al., 2013; Léger et al., 2014; Tran et al., 2014) or cosmic-ray neutron probes (e.g., Franz et al., 2012; Bogena et al., 2013; Lv et al., 2014; Villarreyes et al., 2014). While the depth of penetration for GPR may be up to 10 to 15 m, its spatial extent is quite limited. By comparison, CRS monitors water contents mainly near the soil surface but over much larger areas. Although CRS methods do not provide a horizontal or vertical resolution for soil moisture, it averages water contents over tens of hectares and thus can provide very useful data for agriculture and hydrological models at the hectometer scale. Other studies using magnetic resonance imaging and time-lapse electromagnetic induction are given by Pohlmeier et al. (2009) and Robinson et al. (2012), respectively. Additional applications of HYDRUS, in conjunction with geophysical methods, are discussed below.

Historically, one of the most common applications (∼40 papers) of the HYDRUS models have been to estimate subsurface water fluxes and groundwater recharge and how these processes are affected by soil surface and root zone conditions such as precipitation, evaporation, and the presence or absence of plants (e.g., Scanlon et al., 2002, 2003; Garcia et al., 2011; Kurtzman and Scanlon, 2011; Kodešová et al., 2014; Fan et al., 2015; Turkeltaub et al., 2014; Dafny and Šimůnek, 2016; Neto et al., 2016). For example, Dafny and Šimůnek (2016) showed that for the coastal plain of Israel, groundwater recharge dramatically decreases as a percentage of precipitation from ∼30% to ∼10 and 1% for conditions with bare sandy loess and sandy loess with vegetative covers of 26 and 50%, respectively.

The impact of changing land use on groundwater recharge was investigated in several other studies (e.g., Le Coz et al., 2013; Ibrahim et al., 2014; Turkeltaub et al., 2014, 2015). Of these, Le Coz et al. (2013) found an increase in groundwater recharge because of changes from rainfed to irrigated cropping conditions in a semiarid region. Turkeltaub et al. (2015) evaluated the impact of switching crop type on water and solute fluxes in deep vadose zones. Similarly, changes in groundwater recharge in response to the expansion of rainfed cultivation in the Sahel, West Africa, were evaluated by Ibrahim et al. (2014). Another related application is to anticipate the sensitivity of groundwater recharge to changes in climate in response to greenhouse effects (e.g., Leterme et al., 2012; Newcomer et al., 2014; Pfletschinger et al., 2014; Wine et al., 2015). Additional applications of HYDRUS for evaluating groundwater recharge are given below in the Groundwater Recharge Applications section and on the HYDRUS website.

The versatility of HP1 was demonstrated by Jacques et al. (2008a,b) by means of several examples including (i) the transport of heavy metals (Zn²⁺, Pb²⁺, and Cd²⁺) subject to multiple cation exchange reactions, (ii) transport with mineral dissolution of amorphous SiO₂ and gibbsite [Al(OH)₃], (iii) heavy-metal transport in a porous medium having a pH-dependent cation exchange complex, (iv) infiltration of a hyperalkaline solution in a clay sample (this example considered kinetic precipitation or dissolution of kaolinite, illite, quartz, calcite, dolomite, gypsum, hydrotalcite, and sepiolite), (v) long-term transient flow and transport of major cations (Na⁺, K⁺, Ca²⁺, and Mg²⁺) and heavy metals (Cd²⁺, Zn²⁺, and Pb²⁺) in a soil profile, (vi) cadmium leaching in acid sandy soils, (vii) radionuclide transport, and (viii) long-term uranium migration in agricultural field soils following mineral P fertilization.

More recent HP1 applications include evaluations of (i) laboratory and field experiments involving the treatment of Hg-contaminated soils with activated C (Bessinger and Marks, 2010; Leterme et al., 2014), (ii) CO₂ production and transport in bare and planted mesocosms (Thaysen et al., 2014a), (iii) the effects of lime and concrete waste on vadose zone C cycling (Thaysen et al., 2014b), (iv) chemical degradation of concrete during leaching with rain and different types of water (Jacques et al., 2010), and (v) the effects of chemical degradation on the hydraulic properties of concrete such as porosity, tortuosity, and the hydraulic conductivity (Jacques et al., 2013). Jacques et al. (2012) additionally combined HP1 with the general optimization UCODE program (Poeter et al., 2005) to inversely optimize hydraulic, solute transport, and cation exchange parameters pertaining to column experiments subject to transient water flow and solute transport with cation exchange.

HP1 has recently been used also to solve a number of benchmark problems that were developed for model developers to demonstrate model conformance with norms established by the subsurface science and engineering community (Steefel et al., 2015). These benchmarks involved (i) multirate surface complexation and 1D dual-domain multicomponent reactive transport of U(VI) (Greskowiak et al., 2015), (ii) generation of acidity as a result of sulfide oxidation and its subsequent effect on metal mobility above and below the water table (Mayer et al., 2015), and (iii) implementation and evaluation of permeability–porosity and tortuosity–porosity relationships associated with mineral precipitation and dissolution processes (Xie et al., 2015).

The versatility of the two-dimensional HP2 was demonstrated recently by Šimůnek et al. (2012b) on several examples: (i) sodic soil reclamation using furrow irrigation to demonstrate the cation exchange features of HP2 and (ii) the release and migration of uranium from a simplified uranium mill tailings pile toward a river. These examples included the processes of water flow, solute transport, precipitation and dissolution of the solid phase, cation exchange, complexation, and many other reactions.

Vadose Zone Journal (VZJ) has been a frequent outlet for manuscripts documenting various HYDRUS applications. HYDRUS-1D and HYDRUS (2D/3D) were used in over 100 and 50 VZJ papers, respectively. This means that almost 20% of peer-reviewed journal articles using the HYDRUS models have been published in VZJ. This trend continued in recent years: 18, 14, and eight papers using HYDRUS appeared in VZJ in 2013, 2014, and 2015, respectively. In the sections below, we provide an overview of HYDRUS applications that have appeared in VZJ in recent years and which partly mirror the main types of applications discussed above.

The largest number of HYDRUS papers in VZJ simulated subsurface water fluxes and groundwater recharge (e.g., Dickinson et al., 2014; Pfletschinger et al., 2014; Rieckh et al., 2014; Turkeltaub et al., 2014; Fan et al., 2015; and Guber et al., 2015). Of these, Guber et al. (2015) used HYDRUS-2D to evaluate a new subsurface water retention technology consisting of subsurface polyethylene membranes installed within the soil profile to improve root-zone water storage and to limit downward recharge fluxes. Fan et al. (2015) used both HYDRUS-1D and HYDRUS (2D/3D) to model the effects of plant canopy and roots on soil moisture and deep drainage in forested ecosystems. Dickinson et al. (2014) used HYDRUS-1D to verify the appropriateness of a proposed screening tool for delineating areas with constant groundwater recharge. Turkeltaub et al. (2014) used data collected with a deep vadose zone monitoring system to calibrate HYDRUS-1D and subsequently used the software to investigate the temporal characteristics of groundwater recharge and how recharge may be affected by climate change. Similarly, Pfletschinger et al. (2014) used HYDRUS-1D to evaluate the effects of climate shifts in arid areas on groundwater recharge. Rieckh et al. (2014) further used HYDRUS-1D to evaluate water and dissolved C fluxes in an eroding soil landscape and their dependence on terrain position, while Le Coz et al. (2013) used HYDRUS-1D to evaluate how a change from rainfed to irrigated cropping in a semiarid region will affect groundwater recharge.

The second largest group of HYDRUS applications published in VZJ comprised studies that use data collected with various geophysical methods (e.g., Montzka et al., 2013; Grunat et al., 2013; Moghadas et al., 2013; Ganz et al., 2014; Thoma et al., 2014; Lv et al., 2014; Dimitrov et al., 2014, 2015; and Persson et al., 2015). For example, several issues related to data assimilation, which involved both HYDRUS modeling and ERT or GPR were studied by Grunat et al. (2013), Moghadas et al. (2013), Ganz et al. (2014), Thoma et al. (2014), and Persson et al. (2015). Of these various studies, Persson et al. (2015) used HYDRUS-2D to simulate laboratory experiments involving dye movement in a glass tank. They successfully compared modeled horizontal velocities with those obtained by image analysis and ERT. Experimental and numerical results both showed that horizontal velocities in the capillary fringe are more or less identical to those in the saturated zone. Ganz et al. (2014) used HYDRUS-3D to simulate ponded infiltration into a water-repellent sand and successfully compared their numerical results with ERT observations. They discussed the importance of considering hysteresis for water repellent soils. Lv et al. (2014) calibrated HYDRUS-1D using soil moisture measurements from a network of time-domain transmissometry (TDT) probes and then compared both measured and modeled water content values against cosmic-ray neutron probe estimates. Finally, a series of papers by Dimitrov et al. (2014, 2015) and Montzka et al. (2013) used the HYDRUS-1D model to inversely derive soil hydraulic parameters and surface soil water contents using L-band brightness temperatures. All of these studies demonstrate how numerical modeling of subsurface flow processes can be used to optimize the analysis of geophysical data.

As discussed above, the HYDRUS models are often used to evaluate the transport of particle-like substances such as colloids, bacteria, viruses, or nanoparticles. Two manuscripts addressing these topics appeared in VZJ. Wang et al. (2014b) used HYDRUS to study physical and chemical factors influencing the transport and fate of E. coli in soil affected by preferential flow, while Wang et al. (2015a) evaluated the transport and retention of polyvinylpyrrolidone-coated silver nanoparticles in natural soils.

Several other applications of the HYDRUS software packages models have appeared in VZJ. Two such applications in 2015 focused on the effects of root water uptake on soil moisture dynamics and deep drainage or recharge. Fan et al. (2015) used both HYDRUS-1D and HYDRUS (2D/3D) to model the effects of spatial distributions of the plant canopy, rainfall, and roots on soil moisture and deep drainage in a coastal sand dune forest of subtropical Australia. Périard et al. (2015) used HYDRUS (2D/3D) to simulate root water uptake by romaine lettuce (Lactuca sativa L.) and to evaluate the effect of moisture deficit on tip burn, a physiological disorder that can lead to a complete loss of harvest. A similar HYDRUS-1D study for evaporation was performed later by Huang et al. (2013).

The HYDRUS models were further used in a large number of studies to inversely optimize various soil hydraulic and solute transport parameters (Rühle et al., 2015; Qu et al., 2014; Lv et al., 2014; Caldwell et al., 2013; Schelle et al., 2013). Of these studies, Qu et al. (2014) used HYDRUS-1D to inversely estimate van Genuchten (1980) soil hydraulic parameters from field soil water content measurements at multiple locations to evaluate the spatial variability of the soil water content. Lv et al. (2014) calibrated HYDRUS-1D by optimizing soil hydraulic parameters using soil moisture measurements from a network of TDT probes, while Zhao et al. (2013) used the multistep outflow method to determine the soil hydraulic properties of a frozen soil.

Several specialized HYDRUS modules, as discussed above, have been used also in multiple VZJ publications. For example, Spurlock et al. (2013a,b) used the fumigant module to evaluate soil fumigant transport and volatilization to the atmosphere for different types of fumigant applications. The HP1 module was used further in a study by Thaysen et al. (2014b) to evaluate the effects of lime and concrete waste on carbon cycling in the vadose zone. Skaggs et al. (2014) used the UnsatChem module in a global sensitivity analysis to simulate crop production with degraded waters, whereas Lassabatere et al. (2014) used the dual-permeability flow module to evaluate a new analytical model for calculating cumulative infiltration into dual-permeability soils.

As numerical models, such as the HYDRUS software packages, are becoming increasingly more accurate, comprehensive, and numerically efficient, their application to a large number of theoretical and practical problems is becoming more and more widespread. For these reasons, the Windows-based HYDRUS models are now rapidly becoming useful tools for teaching the principles of water, solute, and heat movement in soils and groundwater, even for users with very little direct knowledge of soil physics and related disciplines and with limited mathematical expertise. As a result, the HYDRUS software packages have been used to advantage in several soil physics and hydrology related textbooks (e.g., Rassam et al., 2003; Radcliffe and Šimůnek, 2010; Lazarovitch and Warrick, 2013; Shukla, 2013). Below, we give a brief account of the more recent books and conference proceedings.

Radcliffe and Šimůnek (2010), in their textbook Soil Physics with HYDRUS: Modeling and Applications, describe a broad range of relatively standard soil physics topics. They used various tools from the HYDRUS family of programs (Šimůnek et al., 2008b) to make the topics more accessible to students. For example, the RETC software is used to describe and quantify the unsaturated soil hydraulic properties, while HYDRUS-1D software was used to demonstrate infiltration, evaporation, and percolation processes of water in soils having different textures and layering. The software is also used to demonstrate various heat and solute transport problems in these systems including the effect of physical and chemical nonequilibrium conditions. The HYDRUS (2D/3D) software is used further to describe two-dimensional flow in field soils, hillslopes, boreholes, and within capillary fringes. The effects of various transport and reaction parameters on solute transport are also evaluated. Using information in this book, users can run HYDRUS and related models for different scenarios and with different parameters, thus obtaining more insight into the physics of water flow and contaminant transport. The book can also be used for self-study on how to use the HYDRUS models.

Another book, Exercises in Soil Physics, was edited by Lazarovitch and Warrick (2013) to complement available soil physics and vadose zone hydrology texts by providing additional practical exercises. The topics of soil physics are explored using nine categories: solid phase, soil water relations, saturated water flow, unsaturated flow, field water flow processes, chemical fate and transport, heat and energy transport, soil gases and transport, and soil variability. Several problems involving variably saturated water flow and root water uptake are solved using HYDRUS-1D. Some of the solute transport problems involved sorbing, nonsorbing, degrading, nondegrading, and volatile solutes with different degrees of dispersion and are solved using STANMOD. Finally, ROSETTA and RETC are used in forward calculations of the soil water retention curve and for inverse calculation of the soil hydraulic properties of the van Genuchten and other soil hydraulic models.

A very extensive HYDRUS-1D tutorial, Soil Physics: An Introduction, was published by Shukla (2013). This textbook focused on coupled liquid water, water vapor, and heat transport in the unsaturated zone of a sandy loam, furrow-irrigated onion field (Deb et al., 2011). Readers are provided with a very detailed description of most HYDRUS-1D input and output windows used in the tutorial including details on how the required input parameters can be obtained and how the output is to be interpreted.

Three special workshops dedicated to various applications of the HYDRUS models have been conducted since 2008. The second, third, and fourth HYDRUS workshops or conferences were organized in Prague, Czech Republic (in 2008), in Tokyo, Japan (in 2008), and again in Prague (in 2013), respectively. A large number of HYDRUS applications presented at these conferences have been published in the conference proceedings (Šimůnek and Kodešová, 2008; Saito et al., 2008; and Šimůnek et al., 2013b), which can be downloaded freely from the HYDRUS website (http://www.pc-progress.com/en/Default.aspx).

As numerical models are becoming much more efficient, comprehensive, and numerically accurate, their application to a large number of theoretical and practical problems is becoming increasingly widespread. This is true not just for the HYDRUS models but also for other models addressing various soil, hydrologic and environmental science, and engineering problems such as the TOUGH models (Finsterle et al., 2008), STOMP (White et al., 2008), SWAP (van Dam et al., 2008), VS2DI (Healy, 2008), and many other models as discussed by Vereecken et al. (2016). As we noted earlier in our 2008 paper (Šimůnek et al., 2008b), we believe that these various models and modeling tools have served, and will continue to serve, an extremely important role in vadose zone research.

In this paper we illustrated a large number of applications of HYDRUS-1D and HYDRUS (2/3D) and its standard and nonstandard specialized add-on modules that significantly expanded the versatility of the models. The popularity of the HYDRUS models and related models (notably, STANMOD, RETC, UNSATCHEM, and HP1) is reflected by their increasing use in a variety of applications and publications. That these models serve a purpose is certainly reflected by the number of downloads from the HYDRUS website (http://www.pc-progress.com/en/Default.aspx). HYDRUS-1D has been downloaded more than 40,000 times since the program was made freely available (7738 times in 2015 alone), STANMOD 6000 times (>1000 times in 2015), and RETC nearly 11,000 times (2260 times in 2015). The website received nearly 140,000 visitors in 2015, while more than 30,000 people are registered users, mostly from the United States, China, Germany, France, Australia, Colombia, Israel, and Turkey (in this order). We hope to continue further development and improvement of these models in the near future as part of a continual cycle of improvement.

While much effort has gone into the development of the HYDRUS models, we also realize that model development and validation and verification never end. In terms of future work, one major priority for us is to formalize most or all of the nonstandard modules that thus far are included only in an approximate manner and without much documentation. In terms of HYDRUS-1D, these nonstandard modules deal with centrifugal forces, freeze–thaw processes, colloid-facilitated transport, colloid transport with changing water contents, isotope transport, and root growth (see HYDRUS 1D: Nonstandard Modules section above). Nonstandard HYDRUS (2D/3D) modules concern centrifugal forces, overland flow, and CO₂ transport and production [see HYDRUS (2D/3D): Nonstandard Modules section above]. Especially important is the coupling of the HYDRUS models with surface runoff processes to produce a more comprehensive surface–vadose zone–groundwater modeling environment. Also needed in the future are further improvements in the accuracy and computational efficiency of the numerical solutions of the governing equations to facilitate larger-scale applications and continual updates of some of the components of the HYDRUS software packages and related models (like RETC and STANMOD) to make them more compatible with the 64-bit Windows 10 operating systems and future Windows versions.

We further realize that models remain a reflection of what is known—or thought to be known—about prevailing subsurface water flow and solute processes and our ability to capture those processes in usable mathematical formulations and related computer software. Many scientific and organizational challenges remain in this respect to advance systematic modeling of all of the physical, chemical, and biotic processes operative in the vadose zone and relevant connections with both groundwater and aboveground surface hydrologic and atmospheric processes. We refer to Vereecken et al. (2016) for a wide-ranging discussion of these aspects within the general context of modeling soil processes.

We greatly acknowledge USDA NIFA (projects CA-R-ENS-5047-RR and CA-R-ENS-7274-H) for the financial support.