Developments and applications of the HYDRUS computer software packages since 2016

The HYDRUS codes have become standard tools for addressing many soil, agricultural, environmental, and hydrological problems requiring the evaluation of various subsurface physical, chemical, and biological processes. There are now many thousands of HYDRUS users worldwide, with thousands of applications of the HYDRUS models appearing in the peer-reviewed literature. In this manuscript, we provide an overview of the capabilities of the most recent Version 5 of HYDRUS, focusing primarily on features implemented since 2016. We briefly describe the standard HYDRUS model and its standard and nonstandard specialized add-on modules that significantly expand the capabilities of the software packages. The standard add-on modules include HPx, UNSATCHEM, Wetland, Furrow, PFAS, COSMIC, DPU, SLOPE Cube, and Particle Tracking. Recent developments of the HYDRUS Package for MODFLOW are also described, along with additional capabilities incorporated into the graphical user interface supporting HYDRUS. Also discussed are new or improved options to simulate the fate and transport of environmental isotopes,


INTRODUCTION
During the past several decades, the HYDRUS codes have become standard tools for addressing many soil, agricultural, environmental, and hydrological problems involving a range of physical, chemical, and biological processes in the subsurface, particularly in the vadose zone.The HYDRUS models remain one of the most widely used numerical tools for simulated processes in the subsurface.There are many thousands of HYDRUS users (and institutions using HYDRUS) around Vadose Zone Journal the world, with many applications (also several thousand) of the models being published in the peer-reviewed literature (see, e.g., Google Scholar).It is beyond the scope of this paper to describe all applications for which the HYDRUS models have been used.A comprehensive list of publications summarizing many one-dimensional (1D) applications of HYDRUS (over 2000) can be found at https://www.pc-progress.com/en/Default.aspx?h1d-references.A similar list (over 1000) of two-dimensional (2D) and three-dimensional (3D) applications of HYDRUS and its predecessors can be found at http:// www.pc-progress.com/en/Default.aspx?h3d-references.
The early history of the HYDRUS models and their predecessors, which can be traced back to the early 1970s, was reviewed by Šimůnek et al. (2008).More recent accounts of the codes are described by Šimůnek et al. (2012, 2013, 2016, 2018).Our objective here is not to provide a detailed history of the HYDRUS models.Instead, we will focus mainly on developments after those previous reviews appeared, especially on the new capabilities of the current HYDRUS software package and its specialized add-on modules.We thus aim to describe the unique capabilities of the current HYDRUS software package relative to those described in the earlier reviews by Šimůnek et al. (2016, 2018).
The manuscript is organized as follows.After this introduction, we provide in Section 2 a brief history of the Windows-based HYDRUS models.Section 3 summarizes the current capability of the standard HYDRUS model, after which we provide information about the standard (Section 4) and nonstandard (Section 5) HYDRUS modules currently available.In Section 6, we next discuss several new features we plan to include in upcoming releases of the HYDRUS software package, as well as the most common applications of the models.Section 7 discusses the HYDRUS use in undergraduate-and graduate-level courses and available textbooks.Finally, in Section 8, we will provide some summaries and conclusions, along with a personal statement by the senior author in the Appendix.

A BRIEF HISTORY OF THE WINDOWS-BASED HYDRUS MODELS
As documented in our previous reports (Šimůnek et al., 2008, 2016), the development of the HYDRUS software packages during the first two decades of their existence followed two separate paths.On the one hand, the HYDRUS-1D software package (its four versions were released in 1998, 1998, 2005, and 2008) focused on water flow and contaminant transport processes in 1D porous media.While the capabilities and a list of processes that these sequential packages could consider increased dramatically in its 20-year-long history (Šimůnek et al., 2008, 2016), the HYDRUS-1D package during its entire existence remained in the public domain and

Core Ideas
• An overview of the capabilities of the most recent Version 5 of HYDRUS.• The description of the standard and nonstandard specialized HYDRUS add-on modules.• A review of selected popular applications of the HYDRUS models.• Description of new or improved upcoming HYDRUS add-on modules.
was freely downloadable.HYDRUS-2D, on the other hand, was a 2D commercial software package.Its two versions were released in 1996 and 1998.In 2007, the HYDRUS-2D package was extended to consider three dimensions and renamed HYDRUS (2D/3D).That software package was released in three versions in 2007, 2011, and 2017, with each version adding new computational and graphical capabilities and specialized add-on modules (Šimůnek et al., 2016, 2018).Interestingly, based on the number of citations, the manuscript of Šimůnek et al. (2016), which was published in Vadose Zone Journal, was ranked "as the most dominant article in the CTM (contaminant hydrology modeling) research field" by Guleria and Chakma (2022), who carried out a bibliometric analysis of research articles published on this topic since 2010.

THE CURRENT HYDRUS MODEL
The last significant development in the family of the HYDRUS models occurred in April 2022 when the two previously independent software packages, HYDRUS-1D (Version 4, for 1D applications) and HYDRUS (2D/3D) (Version 3, for 2D and 3D applications), were merged into a single software package called HYDRUS.This new HYDRUS package is listed as Version 5 to avoid confusion about the independent Version 4 of HYDRUS-1D and Version 3 of HYDRUS (2D/3D).The HYDRUS software package (Version 5) is a Windows-based application for simulating water, heat, and solute movement in 1D, 2D, or 3D variably saturated porous media.The program was additionally extended with many special (old and new) standard add-on modules (see Table 1).

Dimensionality and geometrical shapes
The HYDRUS program is a finite element model that can simulate the movement of water, heat, and multiple solutes in variably saturated media.It can analyze the flow and transport of water and solutes in unsaturated, partially T A B L E 1 Standard HYDRUS add-on modules.

Wetland 2D
The Wetland module was developed to model biochemical transformation and degradation processes in subsurface-flow-constructed wetlands.Two biokinetic model formulations can be chosen in the Wetland module: (1) CW2D (constructed wetland two-dimensional) (Langergraber & Šimůnek, 2005) and ( 2) CWM1 (constructed wetland model 1) (Langergraber et al., 2009).Aerobic and anoxic transformation and degradation processes for organic matter, nitrogen, and phosphorus are described in CW2D, while aerobic, anoxic, and anaerobic processes for organic matter, nitrogen, and sulfur are considered in CWM1.

C-Ride 1D, 2D
The C-Ride module simulates variably saturated water flow, colloid transport, and colloid-facilitated solute transport in porous media.The module accounts for transient variably saturated water flow, both colloid and solute movement due to advection, diffusion, and dispersion, as well as for solute movement facilitated by colloid transport.

Furrow 2D
The Furrow module is a hybrid finite volume/finite element module describing coupled surface-subsurface flow and transport processes during furrow irrigation and fertigation (Brunetti, Šimůnek, & Bautista, 2018).The numerical approach combines a 1D description of water flow and solute transport in an open channel with a 2D description of water flow and solute transport in the subsurface.

PFAS 1D, 2D, 3D
The PFAS module includes options to consider sorption to the air-water interface and solute concentration effects on surface tension and viscosity (Silva et al., 2020).
Particle Tracking 1D Zhou et al. (2021) implemented the Particle Tracking algorithm into HYDRUS.This module can calculate soil water travel times and ages for different locations in the soil profile.

COSMIC 1D
The COSMIC module developed by Brunetti, Šimůnek, et al. (2019) calculates aboveground neutron fluxes using the physically based COsmic-ray Soil Moisture Interaction Code (COSMIC) of Shuttleworth et al. (2013).

Fumigant 1D, 2D, 3D
The Fumigant module includes options to simulate the fate and transport of fumigants in soils (e.g., removal of a tarp, temperature-dependent tarp properties, additional fumigant injections).

HyPar 2D, 3D
HyPar is a parallelized version of the standard 2D and 3D HYDRUS computational modules.

SLOPE Cube 2D, 3D
The SLOPE Cube (Slope Stress and Stability) module was 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 initiations and conduct slope stability analyses under variably saturated soil conditions.

SLOPE Classic 2D
The SLOPE Classic module is intended mainly for stability checks of embankments, dams, earth cuts, and anchored sheeting structures.
saturated, or fully saturated porous media.The program can model 1D vertical or horizontal flow, 2D vertical or horizontal flow, 3D radial symmetry about the vertical axis, or fully 3D domains.Invoked 3D domains may range from simple hexahedral domains to more complex layered systems and almost entirely general 3D domains (with very few constraints on the system's complexity).The program can handle flow regions delineated by irregular boundaries.The flow region, furthermore, may be composed of nonuniform soils with an arbitrary degree of local anisotropy.
Vadose Zone Journal

Considered processes in HYDRUS
The HYDRUS program numerically solves the Richards equation for variably saturated water flow and advectionconduction/dispersion equations for heat and solute transport in one, two, or three dimensions.The flow equation incorporates a sink term to account for water uptake by plant roots.The heat transport equation considers transport due to conduction and convection with flowing water.The transport equations also include provisions for nonlinear nonequilibrium reactions (sorption) between the solid and liquid phases, linear equilibrium reactions between the liquid and gaseous phases (allowing simulations of volatile chemicals), zeroorder production, and two first-order degradation reactions, one of which can be used to represent parent-daughter degradation reactions (e.g., degradation of organic substances, or nitrification of nitrogen species).Chemical nonequilibrium solute transport can be accounted for using a two-site sorption model (e.g., van Genuchten & Wagenet, 1989), which divides the sorption sites into two groups.Sorption on the first fraction of sorption sites is instantaneous, while sorption on the remaining sites is kinetic.In addition, 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 (e.g., van Genuchten & Wierenga, 1976).Attachment/detachment theory, including filtration theory, is additionally included to enable simulations of the transport of viruses, colloids, nanoparticles, and bacteria (e.g., Bradford et al., 2004).
The dispersion tensor includes terms reflecting the contribution of molecular diffusion, tortuosity, and hydrodynamic dispersion.

Boundary conditions
The water flow part of the model can deal with (constant or time-varying) prescribed pressure head and flux boundaries, as well as boundaries controlled by atmospheric conditions.Soil surface (atmospheric) boundary conditions may change during the simulation from prescribed flux to prescribed head-type conditions (and vice versa), accounting for infiltration, ponding, surface runoff, or actual evaporation processes.
The code can also handle a seepage face boundary, through which water leaves the saturated part of the flow domain, and free drainage boundary conditions.Nodal drains can be represented by a simple relationship derived from analog experiments.Additional boundary conditions can deal with various surface and subsurface drip irrigation, triggered irrigation, and other unique irrigation scenarios (e.g., Gärdenäs et al., 2005;Lazarovitch et al., 2023).Special interactions between different types of boundary conditions, such as an atmospheric boundary condition on a nonactive part of a seepage face or dynamic changes in water levels of external water bodies, can also be considered.Finally, interactions between surface water reservoirs (e.g., wells [e.g., Sasidharan et al., 2018Sasidharan et al., , 2019Sasidharan et al., , 2021]], furrows [e.g., Bristow et al., 2020], and wetlands) and the subsurface can also be accounted for in HYDRUS (Šimůnek et al., 2018).For solute transport, the code supports both (constant and varying) prescribed concentration (Dirichlet or firsttype) and concentration flux (Cauchy or third-type) boundary conditions.Additionally, a special boundary condition permitting gaseous diffusion through a stagnant boundary layer above the soil surface is implemented for solute transport, thereby allowing analyses of the fate and transport of volatile chemicals such as fumigants (e.g., Spurlock et al., 2013).

Soil hydraulic property models
The soil hydraulic properties (SHPs) in HYDRUS can be described using the analytical functions of Brooks and Corey (1964), Durner (1994), van Genuchten (1980), or Kosugi (1996), as well as modified van Genuchten-type analytical functions that improve the description of the hydraulic properties near saturation (Vogel & Císlerová, 1988).Alternatively, SHPs can be defined using lookup tables, thus expanding HYDRUS's ability to consider other analytical models (e.g., Iden et al., 2021;Peters et al., 2021) or other shapes.HYDRUS also incorporates hysteresis, assuming that drying scanning curves are scaled from the main drying curve and wetting scanning curves from the main wetting curve (e.g., Kool & Parker, 1987;Lenhard et al., 1991).The code further implements a scaling procedure to approximate variability in the hydraulic properties of a given soil profile using a set of linear scaling transformations that relate the individual soil hydraulic characteristics to those of a reference soil (Vogel et al., 1991).Additionally, HYDRUS can generate stochastic spatial distributions of scaling factors for the pressure head, water content, and hydraulic conductivity, which can be either completely independent or constrained by Miller-Miller similitude (Miller & Miller, 1956).The scaling factors can be either normally or lognormally distributed.

Numerical methods
The governing equations are solved numerically in HYDRUS using a Galerkin-type linear finite element method applied to a triangular or tetrahedral elements network.Integration in time is achieved using an implicit (backward) finite difference scheme for water flow under both saturated and unsaturated conditions, and central or implicit finite element schemes for solute transport.The resulting equations are solved iteratively by linearization and subsequent Gaussian elimination for banded matrices, a conjugate gradient method for symmetric matrices, or the ORTHOMIN method for asymmetric matrices.In order to improve solution efficiency in transient problems, certain measures are taken, including automatic time step adjustment and ensuring that the Courant and Peclet numbers are kept within preset levels.The water content term in the Richards equation for transient flow is evaluated using the mass-conservative method proposed by Celia et al. (1990).Upstream weighing is included as an option for solving the transport equation to minimize numerical oscillations.

Inverse problems
The standard modules of HYDRUS (not the specialized add-on modules except C-Ride and DualPerm in 1D) also implement a Marquardt-Levenberg-type parameter optimization algorithm (Marquardt, 1963) for inverse estimation of soil hydraulic and/or solute transport and reaction parameters from measured transient or steady-state flow and/or transport data.The procedure allows the estimation of several unknown parameters from observed water contents, pressure heads, concentrations, and instantaneous or cumulative boundary fluxes (e.g., infiltration or outflow data).The parameter estimation procedure can also include additional retention or hydraulic conductivity data.A penalty function can be used to constrain optimized parameters to remain in some feasible region (Bayesian estimation).
The capability of the standard HYDRUS module to directly calibrate and optimize unknown parameters is very popular with HYDRUS users, with nearly all HYDRUS applications involving some calibration or fitting process.Most applications use the implemented Marquardt-Levenberg algorithm, a gradient-based optimizer that provides a statistical basis for assessing parameter uncertainty through confidence intervals.However, this local search strategy cannot guarantee convergence to the global minimum.To address this limitation, multiple studies have coupled HYDRUS externally with global search algorithms, such as Particle Swarm Optimizations (e.g., Brunetti et al. 2016a), Genetic Algorithms (e.g., Vrugt et al., 2001), and Shuffled Complex Evolution approaches (e.g., Vrugt et al., 2004;Wöhling et al., 2008).In an attempt to balance the exploration capabilities of global search strategies and the exploitation features of gradient-based optimizers, Brunetti, Stumpp, et al. (2022) combined the Marquardt-Levenberg algorithm embedded in HYDRUS with the Comprehensive Learning Particle Swarm Optimization strategy.The resulting algorithm outperformed both the Marquardt-Levenberg local search strategy and the global optimizer, Shuffled Complex Evolution, on multiple benchmark problems.
Optimization, sensitivity analyses, and uncertainty assessments are increasingly being implemented in conjunction with HYDRUS.Global sensitivity analysis methods, such as Sobol, Morris, E-FAST, and Random Balance Design-FAST, have been used to investigate the influence of model parameters (e.g., soil hydraulic and solute reactive transport parameters) on the fitting accuracy (e.g., Brunetti et al., 2016b) or on specific variables of interest (Brunetti, Kodešová, et al., 2021;Brunetti, Papagrigoriou, et al., 2021).Similarly, HYDRUS has been coupled externally with Bayesian inference algorithms for calibration and uncertainty assessment purposes (e.g., Brunetti, Šimůnek, et al., 2020;Brunetti, et al., 2019;Guo et al., 2020;Schübl et al., 2022;Vrugt et al., 2008;Wöhling & Vrugt, 2011).In this perspective, Brunetti et al. (2023) recently compared Markov chain Monte Carlo ensemble samplers for vadose zone inverse modeling, highlighting the advantages and limitations of using these algorithms in conjunction with HYDRUS.Samplers have been embedded into the HYDRUS source code, thus making it possible to incorporate them in the HYDRUS software package.

STANDARD HYDRUS ADD-ON MODULES
The capabilities of the HYDRUS software packages have been expanded significantly by including various standard and nonstandard specialized add-on modules described in the following two sections.The standard add-on modules are fully incorporated and supported by the HYDRUS graphical user interfaces (GUIs) and documented in detail in the technical and user manuals and via online help.This is not the case for several nonstandard add-on modules (described below in Section 5), which require additional work outside the GUI.
Common to all HYDRUS add-on modules is that they simulate variably saturated water flow and heat and solute transport in porous media.The specialized add-on modules provide additional capabilities, such as general reactive transport (using the HP1 and HP2 models) or reactive transport with specific chemistry (notably the Wetland and UNSATCHEM modules).Other modules provide additional modeling capabilities, such as to account for preferential flow (the DualPerm module), colloid-facilitated solute transport (the C-Ride module), or the transport of polyfluoroalkyl substances (PFAS) (the PFAS module) or fumigant (the Fumigant module) compounds.Some add-on modules are available for all three dimensions, while others are available only for a particular dimension (e.g., 1D or 2D; see Table 1).In the subsections below, we will discuss only the modules that have been added since 2016, were expanded meaningfully (e.g., having some new capabilities or new dimensions), or Vadose Zone Journal have been used in interesting applications.We will not discuss modules already available in 2016 (and thus described in Šimůnek et al., 2016).

HPx
The current HPx modules (Jacques et al., 2018; HP1 for 1D and HP2 for 2D) couple the PHREEQC-3 geochemical code (Parkhurst & Appelo, 2013) with the HYDRUS computational modules.The 1D version HP1 (the first HP1 versions were coupled with PHREEQC-2; Parkhurst & Appelo, 1999) was first released in 2005 (Jacques & Šimůnek, 2005;Jacques et al., 2008aJacques et al., , 2008b) ) and used successfully in many applications.The 2D equivalent of HP1, that is, HP2, was released in 2012.HPx is a comprehensive simulation module that can be used to simulate (1) transient water flow, (2) the transport of multiple components, including diffusion in the gas phase and exchange with the atmosphere, (3) mixed equilibrium/kinetic biogeochemical reactions, and (4) heat transport in 1D or 2D variably saturated porous media.The HPx modules can simulate a broad range of low-temperature biogeochemical reactions in water, the vadose zone, and/or groundwater systems, including interactions with minerals, gases, exchangers, and sorption surfaces based on thermodynamic equilibrium, kinetic, or mixed equilibrium/kinetic reactions.Additionally, the framework allows the use of flow and transport parameters based on the geochemical state at a given time and location in the domain (Jacques et al., 2018); example applications can be found in Xie et al. (2015) and Valdes-Abellan et al. (2017).
The most significant changes implemented recently in the HPx software are (1) more flexibility, reuse, and reproducibility of model structures and input, (2) the use of alternative geochemical solvers, and (3) better integration of graphical user elements (see Section 2.4 in the Supporting Information).Elements contributing to the first point are (1) full implementation of an alternative scripting language (MyBasic, https://paladin-t.github.io/my_basic),which is a more modern and structured BASIC-like scripting language with object-oriented programming possibilities, (2) options to run Python (Version 3.8) scripts during run time, (3) dynamic definition and use of inline variables in the input files, and (4) a graphical component to define global variables (more details are in Section 2.4 in the Supporting Information).Using the graphical component to define global variables, relatively complex parameter input tables can be created for a given project.Examples (available on request) have been developed to assess the acidification of pyrite-rich soils (Mosley et al., unpublished data, 2024), the implementation of the Tippings model for sorption of heavy metals on organics (Lofts & Tipping, 2011), and cement hydration and aging.
Recently (and currently only available upon request), the PHREEQC geochemical calculations can be replaced by ( 1) ORCHESTRA (Meeussen, 2003), which solves similar geochemical systems as PHREEQC but also allows for the implementation of alternative geochemical models because the model equations are not in the source code but defined in separate text files, and (2) script functions defined by users in the input file.An application of the latter is replacing a (time-consuming) full geochemical solver with a (faster) statistical metamodel derived using machine learning algorithms.Laloy and Jacques (2022) used Python scripts and libraries to replace geochemical calculations with a deep neural network or a k-nearest neighbors model for 1D or 2D cement leaching problems.Depending on the complexity of the model, speedups compared with a single-thread run of the full geochemical model ranged between 10 and 30 times, with results very close to theoretical values calculated by assuming a zero cost for the geochemical calculations.Another example of the U transport application is given by Prasianakis et al. (unpublished data, 2024).Some interesting applications have been implemented in the HPx module since 2016.In addition to the carbon cycling study of Jia et al. (2021), Thaysen, Jessen, et al. (2014) and Thaysen, Jacques, et al. (2014) demonstrated that HPx and other generic coupled reactive transport models provide a flexible framework for implementing soil organic matter dynamics such as those of multiple carbon pools, the effects of biomass, and the fate of gaseous carbon dioxide.Jacques et al. (2018) demonstrated the effects of soil environmental variables (e.g., water content and temperature) on the degradation rate constant and the migration of solid organic matter by bioturbation.Another set of applications concerns the fate of major elements and contaminants in soil profiles and saturated media (e.g., Bruneel et al., 2021;Jiang et al., 2021;Latrille et al., 2021).Colloids and colloid-facilitated solute transport were simulated using HPx by Zhou et al. (2016) and Makselon et al. (2017).The latter used the DVLO theory to calculate the attachment efficiency during variable geochemical conditions.Applications of pyrite oxidation in acid sulfate soils are given by Knabe et al. (2018), who implemented pyrite oxidation under nitratereducing conditions in aquifers, and by Salo et al. (2023) and Mosley et al. (unpublished data, 2024), who simulated geochemical changes in acid sulfate soils following prolonged drought events due to climate change.Tremosa et al. (2020) assessed the weathering of shale tailings by combining transient water flow, mobile-immobile regions, and gas diffusion with geochemical reactions.Besides environmental applications, coupled reactive transport tools are also useful for more engineering-oriented problems, such as the durability of building materials and waste repositories.Although not many applications have been published on these topics using HPx (e.g., Santiago et al., 2023), examples of other similar generic coupled reactive transport models illustrate the feasibility and interest of modeling complex, large-scale engineering Vadose Zone Journal applications (see Claret et al. [2022] for deep radioactive waste repositories).
The HPx modules are general reactive transport models in which users need to define the geochemical system themselves.As a result, these modules provide users with a lot of freedom in designing their particular geochemical system, which facilitates their application to a very wide range of subsurface contaminant transport problems.On the other hand, the following two modules (i.e., UNSATCHEM and Wetland) are modules with specific chemistry in which the geochemical system is predefined, which restricts their application to problems involving the prescribed chemical system.However, the two modules with specific chemistry are much easier to use and are computationally far more efficient than the more general HPx models.

UNSATCHEM
The UNSATCHEM geochemical module (Šimůnek & Suarez, 1994(Šimůnek & Suarez, , 1997) ) is the oldest HYDRUS module.This module was initially a standalone software package simulating carbon dioxide's 1D transport and production and the transport of major ions in variably saturated porous media, including major ion equilibrium and kinetic nonequilibrium chemistry.UNSATCHEM simulates the transport of major ions (i.e., Ca 2+ , Mg 2+ , Na + , K 2+ , SO 4 2− , CO 3 2− , and Cl − ) and their equilibrium and kinetic geochemical interactions, such as complexation, cation exchange, and precipitation-dissolution (e.g., of calcite, gypsum, and dolomite).UNSATCHEM was first incorporated in 2005 as a standard add-on module into Version 3 of HYDRUS-1D, but then also in 2012 as a 2D module into HYDRUS (2D/3D) (Šimůnek et al., 2016), and in 2022 as a 3D module into HYDRUS (2D/3D).The current Version 5 of HYDRUS thus supports UNSATCHEM in all three dimensions.

Wetland
The Wetland module simulates aerobic, anoxic, and anaerobic transformation and degradation processes for organic matter, nitrogen, phosphorus, and sulfur during the treatment of polluted wastewater in subsurface constructed wetlands (Langergraber & Šimůnek, 2012).The Wetland module uses two biokinetic model formulations (constructed wetland twodimensional [CW2D] of Langergraber & Šimůnek [2005] and constructed wetland model 1 [CWM1] of Langergraber et al. [2009]) to account for complex conditions that may occur in various types of wetlands.CW2D considers aerobic and anoxic transformation and degradation processes for organic matter, nitrogen, and phosphorus, whereas CWM1 considers aerobic, anoxic, and anaerobic processes for organic matter, nitrogen, and sulfur.Both biokinetic model formulations were developed for constructed wetlands treating municipal wastewater; they can be applied only to 2D problems.No significant changes have been implemented into the Wetland module since 2016.However, there have been many interesting applications of this module.For example, Pucher et al. (2017) used the HYDRUS Wetland module to optimize the depth of a zeolite layer in a one-stage vertical flow wetland.Similarly, Pucher and Langergraber (2018) used the module to evaluate the functioning of vertical flow wetlands with filter media of different grain sizes.Yasinta et al. (2019) and John et al. (2020) simulated aeration intensity in a saturated vertical upflow constructed wetland.Jayswal and Rodríguez (2021) evaluated the effects of temperature and loading frequency on the performance of vertical subsurfaceflow constructed wetlands.Khan et al. (2022) furthermore assessed the capability of the HYDRUS wetland module to simulate flow and nitrogen removal processes in pilot-scale treatment peatlands under frost and no-frost conditions.A list of other applications of the Wetland module can be found on the HYDRUS website (https://www.pc-progress.com/en/Default.aspx?h3d2-wetland).

Furrow
The Furrow module simulates coupled surface-subsurface flow and transport processes during furrow irrigation and fertigation (Brunetti, Šimůnek, & Bautista, 2018).To reduce the complexity and computational cost, the solver uses a pseudo-3D numerical approach (in 1D and 2D) (Figure 1).The 3D furrow domain is horizontally divided into N elements, associated with a cell-centered 1D finite volume scheme used to resolve overland flow processes spatially.Each element j is characterized by a 2D finite element discretization of the soil domain, which is used to simulate subsurface processes.Zero-inertia approximations of the hydrodynamic wave and advection equation are used to describe water flow and solute transport in the open channel.This reduces model complexity and computational cost and enables one to simulate open-and blocked-ended channels.The standard 2D HYDRUS computational module calculates infiltration, describes water flow and solute transport, and may include root water and nutrient uptake in each element j.The overland and subsurface modules are coupled sequentially in time.In particular, at time t, water levels in the channel are used in HYDRUS as Dirichlet boundary conditions to calculate infiltration in all elements.The resulting infiltration vector is then used to solve overland flow using a fourth-order Runge-Kutta scheme, thereby calculating the water level at the next time step.Soil water contents, pressure heads, and concentrations in each element are stored and updated during the simulation, thus enabling one to thoroughly appraise surface and subsurface water and solute movement along the furrow.
The robustness of the proposed Furrow module was examined and confirmed by mesh and time-step sensitivity analyses (Brunetti, Šimůnek, & Bautista, 2018).The model was theoretically validated by comparison with simulations conducted using the well-established WinSRFR model and experimentally validated by comparison with field-measured data from a furrow fertigation experiment conducted in the United States (Brunetti, Šimůnek, & Bautista, 2018).

PFAS
The PFAS module was developed to simulate, under dynamic vadose zone conditions, the fate and transport of PFAS, often called "forever chemicals."When evaluating the transport of PFAS compounds in the subsurface, one needs to consider, in addition to standard transport, sorption, and degradation processes (as considered in the standard HYDRUS modules), sorption onto the air-water interface (AWI) (Silva et al., 2020;Silva, Šimůnek, Guelfo, et al., 2022).In that case, the term in HYDRUS that accounts for the distribution of a solute between the liquid and gaseous phases (i.e., for volatile chemicals such as ammonia, pesticides, or fumigants) is replaced by a new term accounting for solute sorption onto the AWI as follows: where a v is the air content [L 3 L −3 ], c g is the gas-phase concentration [ML −3 ], t is time [T], Γ is the interfacial adsorbed concentration (i.e., mass or moles per unit area of the interface, or ML −2 ), and A aw is the air-water interfacial area, which has dimensions of [L 2 L −3 ] or [L −1 ].Under transient flow conditions, the air-water interfacial area is a dynamic variable that increases during drying and decreases during wetting, with a corresponding exchange of PFAS compounds between the two phases across the AWI.The air-water interfacial area A aw in the PFAS module is calculated directly from the pressure head-saturation relationship (e.g., Bradford et al., 2015): where σ aw is the air-water surface tension [MT −2 ], P aw is the capillary pressure [ML −1 T −2 ], θ s is the saturated water content [L 3 L −3 ], ρ w is the density of water [ML −3 ], and g is the gravitational acceleration [LT −2 ].Since there is still an ongoing discussion about the suitability of Equation ( 2) (e.g., Silva, Šimůnek, & McCray, 2022) for simulating the fate and transport of PFAS, we additionally implemented a user-specified constant k m , which allows one to scale the interfacial area A aw .The interfacial adsorbed concentration, Γ, can be calculated directly using the Freundlich-Langmuir adsorption isotherm as: where Γ max [ML −2 ] is the maximum surface concentration (i.e., at the solubility limit), and K H (for linear sorption, The two-site sorption concept (similar to sorption onto the solid phase) can be used for solute sorption on the AWI: where sorption, Γ e [ML −2 ], on one fraction of the AWI sites (the Type-1 sites) is assumed to be instantaneous, while sorption, Γ k [ML −2 ], on the remaining (Type-2) AWI sites is considered to be time dependent.Sorption on the Type-2 nonequilibrium AWI sites is assumed to be a first-order kinetic rate process.The mass balance equation for the Type-2 AWI sites in the presence of degradation is given by where ω AWI is the first-order rate constant [T −1 ].
A scaling technique is additionally used in the PFAS module to express concentration effects on the soil hydraulic functions.The primary impact of surface-active organic solutes on unsaturated flow is through the dependence of the soil water pressure head, h [L], on surface tension, σ [MT −2 ], and the contact angle, ω [-], and of the hydraulic conductivity, K [LT −1 ], on kinematic viscosity, ν [L 2 T −1 ].The capillary pressure head of the target liquid can then be scaled with surface tension σ as follows: where α h is the pressure head scaling factor [-], h 0 is the reference soil water pressure head [L], and σ 0 is the reference surface tension [MT −2 ].The dependence of surface tension on solute concentration can be described as (Adamson & Gast, 1997;Henry et al., 2001): where a [L 3 M −1 ] and b [-] are constants for the compound of interest, σ(c) is the surface tension at concentration c, and σ 0 is the surface tension at a reference concentration c 0 .Similarly, the concentration-dependent unsaturated hydraulic conductivity, K(θ, c), of the target liquid can then be scaled with the kinematic viscosity ν as follows: where α K is the hydraulic conductivity scaling factor [-] and K 0 is the hydraulic conductivity at reference concentration c 0 .
The dependence of viscosity on the solute concentration is (Smith & Gillham, 1999) where d [L 3 M −1 ] and e [-] are constants for the compound of interest, ν(c) is the kinematic viscosity at concentration c, and ν 0 is the kinematic viscosity at reference concentration c 0 .
The options to consider equilibrium sorption to the AWI and the concentration effects on viscosity and surface tension were first implemented in Version 5.01 of HYDRUS only for 1D and 2D, but later extended to 3D in Version 5.03.Two-site sorption to the AWI has become available only in Version 5.03.Since the PFAS module has been released only recently, there have not yet been any published applications of this module other than from our team (Silva et al., 2020;Silva, Šimůnek, & McCray, 2022).This is despite immense interest in this module by the research community because of public and regulatory agencies' increased awareness about the environmental dangers of polyfluoroalkyl compounds.

Particle Tracking
The Particle Tracking algorithm from Šimůnek (1991) was recently implemented into the Particle Tracking module of HYDRUS-1D (Zhou et al., 2021) and then in Version 5 of HYDRUS as a 1D add-on module.The module tracks the position of hypothetical particles, either initially defined in the soil profile or released at its boundaries (e.g., with rainfall or irrigation).New particles may be created or leave the soil profile through the soil surface and the bottom, depending on the actual development of the moisture profile.The results of this module can be used to calculate soil water travel times and water ages at different locations in the soil profile.
The Particle Tracking algorithm is based on water balance calculations using soil water profiles simulated by numerically solving the Richards equation.The first monitored particle below the soil surface is at depth z = z 0 and time t = t 0 .Between this particle and the soil surface is an amount of water W 0 [L], which can be calculated by integrating the water content between positions z 0 and L: During time interval (t 0 , t 1 ), an amount of water N [L] passes through the soil surface depending on surface fluxes Vadose Zone Journal such as rainfall, irrigation, and evaporation: where E(t) [LT −1 ] is actual evaporation and I(t) [LT −1 ] is actual infiltration from precipitation or irrigation.During the same interval, the layer between the soil surface and the monitored particle is depleted by water uptake from the root zone by a transpired amount S [L]: where S(z,t) is the sink (extraction) term [T −1 ].At time t 1 , one hence has the following quantity of water W 1 [L] between the soil surface and the monitored particle: Equation ( 13) shows that W 1 is the initial amount of water enriched by infiltration and reduced by evaporation and root water uptake (RWU).The monitored particle is then located at a depth of z = z 1 , which can be calculated from the mass balance equation: (, ).
(14) By repeating the solution of this equation for the time sequence (t 0 , t 1 , . . ., t n ), one obtains a sequence of depths (z 0 , z 1 , . . ., z n ) defining the trajectory of the observed particle.
The locations of the other particles further down in the soil profile can be calculated analogously, except that the amount of water is now balanced between the next two neighboring particles.The amount of water W 0 (W 1 ) between two neighboring particles at time t 0 (t 1 ) is as follows: During time interval (t 0 , t 1 ), the amount of water between the two neighboring particles is depleted by the transpiration amount S: The resulting mass balance equation has then the following form: (, ). (17) Again, by repeating the solution of this equation for the time sequence (t 0 , t 1 , . . ., t n ), one obtains a trajectory (z 0 , z 1 , . . ., z n ) of the particle in question.
The Particle Tracking algorithm itself proceeds as follows: From the known position of the particles at the beginning of the time interval, the pre-solved development of the moisture profile, and the actual values of infiltration, evaporation, and transpiration, the new position of the first monitored particle is calculated using Equation ( 14).The new positions of all other particles are then calculated using Equation ( 17).New particles may be created or leave the profile at the soil surface and the bottom of the soil profile, depending upon the transient development of the moisture profile.Those particles can be released after a certain amount of water passes through the soil surface or with each rainfall.In the latter case, the age of each rainfall can be tracked throughout the soil profile.
By calculating the trajectories of flowing particles, the movement of inert substances not subject to dispersion can be modeled.The average concentration within particular regions between the two neighboring particles can be determined easily from the initial concentration (defined for individual regions) and from a decrease in the water content of these layers due to the extraction term S. The initial position of particles can be defined geometrically (at specified depths) or based on mass balance calculations (by water storage).Similarly, the release of new particles at the boundary can be defined chronologically (at specified times) or meteorologically (rainfall events or depths).
We demonstrate here some typical results of the Particle Tracking module using a dataset reported by Stumpp et al. (2012), who provided all necessary details about the data (which can also be downloaded from the HYDRUS website).Figure 2 shows the trajectories of hypothetical particles initially in the soil profile and generated at the soil surface.The Particle Tracking module additionally provides various particle statistics, such as the initial, final, and total number of particles, the number of particles released at the top and bottom of the profile, the beginning (when released) and end (when leaving the profile) days for each particle, and the particle residence/travel time (Figure 3).Users can use this information to obtain a frequency distribution of particle travel times.
The Particle Tracking module additionally provides information about the age of individual particles since they are released, from which one can interpolate the age of water at  different locations in the soil profile (Figure 4) and the average age of water taken up by roots at any time (Figure 5).This information may be useful for designing various irrigation schemes or for selected ecological applications.
The average age of water taken up by roots (the entire root system), RWUAge, at any particular time is evaluated as follows: where S is RWU [T −1 ] and Age is the water age at a particular location (interpolated from the age of particles evaluated by the Particle Tracking module).Winter wheat was grown in the lysimeters reported by Stumpp et al. (2012).The root zone increased from 0 to 100 cm during one growing cycle of about 90 days.Since the roots were relatively shallow early in the growth cycle, they took up mainly relatively young water from the surface layer.Once the roots reached deeper layers, they started taking up older water from deeper soil layers.However, since the root density was assumed to decrease linearly from the top to the bottom of the root zone, most water taken up by roots was relatively young water throughout the growing cycle (Figure 5).

COSMIC
The COsmic-ray Soil Moisture Interaction Code (COSMIC) module (Brunetti, Šimůnek, et al., 2019) was mainly developed to enable the assimilation of cosmic-ray neutron probe (CRNP) data by HYDRUS.The CRNP measures fast neutrons produced by nuclear interactions between incoming cosmic rays and elements of the Earth's atmosphere (Zreda et al., 2008).When cosmic rays reach the soil surface, they penetrate

DPU
Plant roots generally take up water and solutes moving in the rhizosphere due to the transpiration demand.While water sustains plant growth and its hydraulic functioning, chemicals are bioaccumulated in different plant tissues, where they are eventually metabolized into daughter compounds.The uptake and accumulation of environmental contaminants in edible parts of a plant can pose a major risk to human health.Hence, it is paramount to develop numerical tools to track the fate of chemicals in the soil-plant system, which can also be used for phytoremediation design purposes.
The DPU (i.e., Dynamic Plant Uptake) module addresses this need and extends the capability of HYDRUS to simulate the translocation and transformation of neutral compounds in the soil-plant domain (Brunetti, Kodešová, et al., 2019).As shown schematically in Figure 7, this is accomplished by coupling HYDRUS with a modified version of the multicompartment DPU model of Trapp (2007), which discretizes the plant into four compartments (i.e., roots, stem, leaves, and fruits).The following processes are considered in the DPU module: (a) translocation of compounds with the transpiration stream from roots to the stem and from the stem to leaves and fruits, (b) volatilization of compounds in the stem, leaves, and fruits, (c) gaseous and particle deposition from air to the stem, leaves, and fruits, (d) dilution of compounds by plant growth in all compartments, and (e) sorption and metabolization of compounds in all compartments.At each time step, HYDRUS computes and passes the actual root water and solute uptake flux to the DPU part, which calculates the mass balance and concentrations of the compound(s) in each compartment.The resulting numerical framework, which accounts for multiple metabolization pathways in plant tissues, was successfully used to predict the translocation and transformation of carbamazepine in green pea plants (Brunetti, Kodešová, et al., 2021).The currently available DPU module is fully incorporated into the HYDRUS software package and currently exists for 1D and 2D applications while still being limited to neutral compounds.Brunetti, Kodešová, et al. (2022) extended and validated the DPU module for ionizable chemicals, but their approach is not yet included in the HYDRUS software package.

SLOPE Cube
The SLOPE Cube (or SLOPE 3 ; Slope Stress and Stability) module was developed in cooperation with Dr. Ning Lu from the Colorado School of Mines as a supplemental 2D module of HYDRUS (2D/3D) to simulate 2D transient fields of soil moisture, soil suction, suction stress, total and effective stresses, slope displacement, and a local factor of safety.The module was originally available as a 2D add-on in 2015 and extended to 3D in 2019.The module, which uses a unified effective stress approach for saturated and unsaturated conditions (Lu et al., 2010), is intended to predict spatially and temporally infiltration-induced landslide initiation and carry out slope stability analyses for variably saturated soil conditions.Transient moisture and pressure head fields are obtained directly with HYDRUS and then used to compute effective stress and displacement fields in hillslopes (Lu & Godt, 2013).
Instead of using a one-factor safety methodology for the entire slope, common in conventional slope stability analyses, the SLOPE Cube module computes a field of the safety factors for the whole domain within the hillslope (Lu et al., 2012), thus allowing identification of potential failure surface zones or surfaces.Three new advancements were implemented in SLOPE Cube far beyond traditional physically based landslide models.One is the utilization of pressure head-based effective stress to unify the description of effective stress distribution in hillslopes under variably saturated conditions.Another advancement is the computation of the transient hillslope displacement field caused by variations in effective stress.A third advancement involves implementing the recently established concept of a local factor of safety to capture the evolution of stress paths toward a failure state in hillslopes.The implemented hydromechanical framework has been extensively validated through decades of field monitoring, case studies, and laboratory and field testing, as documented thoroughly by Lu and Godt (2013).

HYDRUS Package for MODFLOW
There have also been further improvements in the "HYDRUS Package for MODFLOW" (HPM) (Beegum et al., 2018(Beegum et al., , 2019(Beegum et al., , 2020)).The HPM, originally developed by Twarakavi et al. (2008), accounts for the main processes and factors affecting fluxes in the vadose zone as incorporated in HYDRUS-1D.These include precipitation, infiltration, evaporation, redistribution, capillary rise, plant water uptake, water ponding on the soil surface, surface runoff, and soil moisture storage.The resulting package provided information about groundwater recharge to the 3D modular finite-difference groundwater model MODFLOW (Harbaugh et al., 2000).The coupled model (HPM+MODFLOW) proved to be effective in addressing spatially variable saturated-unsaturated hydrologic processes at the regional scale by allowing for spatially and temporarily variable water fluxes at the soil surface and in the root zone, considering complex layering in the unsaturated zone, and permitting alternating recharge and discharge fluxes.
The original coupled algorithm assumed that the groundwater table in the HPM profile remained constant throughout an entire MODFLOW time step.This often resulted in unrealistic sudden inflow and/or outflow fluxes at the bottom of the HPM profile when the groundwater table was updated at the next time step.Beegum et al. (2018) developed a methodology to eliminate these sudden bottom fluxes by updating/modifying the pressure head profile at the bottom of the HPM profile after every MODFLOW time step.They showed that the new coupling algorithm effectively eliminated unrealistic sudden variations in the bottom flux of the HPM profiles.
A major limitation of HPM was that it could not simulate solute transport along with water flow and thus could not be used for evaluating groundwater contamination scenarios.This limitation was overcome by Beegum et al. (2019), who integrated the solute transport part of HYDRUS-1D into HPM for the unsaturated zone and MT3DMS for the saturated zone.The resulting update of the coupled model can now simulate solute transport involving many biogeochemical processes and reactions in the vadose zone, including firstorder degradation, volatilization, linear or nonlinear sorption, one-site kinetic sorption, two-site sorption, and two-kineticsite sorption.Beegum et al. (2019) verified the developed model using HYDRUS (2D/3D) and showed that the model effectively simulates the movement of water and contaminants in saturated-unsaturated flow domains.Beegum et al. (2020) further demonstrated the use of MODFLOW with the HYDRUS-1D package and MT3DMS by investigating atrazine contamination in the Zwischenscholle aquifer in Germany.Additional applications of HPM have been reported by Dandekar et al. (2018), Samuel et al. (2020), Szymkiewicz et al. (2018), Zeng et al. (2019), Dai and Guo (2020), and Sarma and Singh (2021), among others.

NONSTANDARD HYDRUS ADD-ON MODULES
Several nonstandard computational modules significantly expand the capabilities of the HYDRUS-1D, HYDRUS (2D/3D), and HYDRUS software packages.Although they still can be run from the standard HYDRUS GUIs, users are usually required to provide manually an additional input file with supplementary information needed for a particular module or to interpret selected input and/or output variables differently from the standard versions.Users may also need to prepare their own graphical output from the output text files.Most of the nonstandard modules can be downloaded freely from the HYDRUS website, together with many examples demonstrating their use and brief descriptions of the theories behind the modules and their implementation (for details, see Table 2 or the HYDRUS website).
A few of the standard HYDRUS add-on modules were initially posted only as nonstandard modules, but they were later included as standard packages with full graphical support when the scientific community showed considerable interest in their use.Examples are C-Ride, Fumigant, and PFAS.Nonstandard HYDRUS modules currently listed in Table 2 likely will be included in the standard HYDRUS software packages in the future.
We are currently developing several additional nonstandard modules, which we intend to convert into standard packages.They involve the modules discussed below in Section 6, such as the Isotope module (Section 6.1), the extended RWU module (Section 6.2), and the Brunswick module (Section 6.3).

DISCUSSION OF RECENT DEVELOPMENTS AND APPLICATIONS
In this section, we will first discuss some of the most recent developments, including processes and options we intend to include in upcoming releases of HYDRUS as either standard or nonstandard add-on modules, depending upon the interest of the research community.New features include options to consider the fate and transport of environmental isotopes (Section 6.1), double vegetation as required in multicropping systems (Section 6.2), more physically based RWU models compared to the currently available Feddes model (Section 6.3), and more complex descriptions of the SHPs, such as the capillary and noncapillary parts of the water retention and hydraulic conductivity functions (Section 6.4).
The second part of this section provides an overview of some of the most common applications of the HYDRUS models or those we see as obtaining broader prominence.They include alternative irrigation applications (Section 6.5), evaluations of the hydraulic functioning of various low-impact developments (LIDs) (Section 6.6), and applications involving different managed aquifer recharge (MAR) technologies (Section 6.7).

Environmental isotopes
The standard third-type (Cauchy type) boundary condition defining the solute flux at the upper boundary is as follows: T A B L E 2 Nonstandard HYDRUS add-on modules.

Isotopes 1D
The Isotopes module for HYDRUS-1D modifies the upper solute transport boundary condition by allowing isotopes not to accumulate at the upper boundary due to evaporation (in contrast to the standard treatment of solutes during evaporation, where solutes stay behind in the soil while water is removed) without considering fractionation processes (Stumpp et al. 2012).This option (with constant fractionation) is now available in Version 5 of HYDRUS.

ColThetaTr 1D
The ColThetaTr module for HYDRUS-1D (both direct and inverse) simulates transient water content effects on colloid transport and the attachment/detachment of colloids from the air-water interface (Bradford et al., 2015).

RootGrowth 1D, 2D
The RootGrowth module was developed by Hartmann et al. (2018) as a supplemental module of the HYDRUS software package (both HYDRUS-1D and HYDRUS [2D/3D]) to model root growth as a function of different environmental stresses.

Freeze 1D
The Freeze module considers freezing/thawing processes and their impact on flow and transport (Hansson et al., 2004;Zheng et al., 2021).

Tillage
The Tillage module performs tillage mixing (averaging or homogenizing water contents and solute concentrations due to tillage over a specified depth) at specified times (Mallman et al., 2014).

Overland 2D
The Overland module considers overland flow (described using a kinematic wave equation) on an atmospheric boundary (van Genuchten & Šimůnek, 2004).

Carbon 2D
The nonstandard Carbon module extends the standard UNSATCHEM module to also consider the production and transport of carbon dioxide in soils, as in the 1D version of UNSATCHEM (Šimůnek & Suarez, 1993).
where q 0 represents the water flux [LT −1 ] and c 0 is the concentration associated with the flux [ML −3 ].When the water flux is directed into the soil profile (e.g., rainfall or irrigation), c 0 is the concentration of the incoming fluid.No solutes are assumed to leave the soil with evaporation at the atmospheric boundary, and hence c 0 in Equation ( 19) is then zero.This standard treatment of the upper boundary condition for solute transport was modified by Stumpp et al. (2012) to account for isotope transport.They assumed that isotope fractionation processes could be neglected and that the relative concentration of isotopes (δ content) will not increase at the upper boundary due to evaporation, and thus, that water and solutes will move at similar rates during evaporation.This contrasts with the standard treatment of solutes during evaporation in HYDRUS-1D, where solutes stay behind in the soil while water is removed.This modification resulted in a nonstandard Isotope module, which has been successfully used by Huang et al. (2015), Sprenger et al. (2015, 2016), Brinkmann et al. (2018), and Nasta et al. (2023), among others.For example, Sprenger et al. (2016) used the modified model to infer soil water residence times at different depths, while Brinkmann et al. (2018) used the model to estimate the residence time distribution of soil water and to identify the temporal origin of water taken up by Picea abies and Fagus sylvatica.
Ignoring evaporative enrichment, as done in the modified HYDRUS-1D, still leads to underestimations of the 2 H and 18 O concentrations in the topsoil, which may be more sig-nificant in regions with higher evaporative losses (Sprenger et al., 2018).For this reason, we modified the Isotope module of Stumpp et al. (2012) further by considering constant evaporative fractionation when c 0 in ( 19) is equal to the surface concentration, c(L,t), multiplied by a fractionation factor r Fract .This factor equals 0.0 when all solutes are left behind in the soil (the standard treatment of solutes), equals 1.0 when fractionation is completely neglected and hence when water and solutes (or isotopes) move at similar rates, and falls between 0 and 1 (close to 1) when evaporation fractionation is considered to be constant.This modified Isotope module is now a standard option in HYDRUS (Version 5).
The Isotope module was further modified by Zhou et al. (2021), who considered the dependence of evaporation fractionation on environmental factors, such as atmosphere temperature and relative humidity, using the Craig and Gordon (1965) and Gonfiantini (1986) models.The resulting formulation has been used successfully in multiple applications (e.g., Post et al., 2022;Zhou et al., 2022Zhou et al., , 2023)).For example, Post et al. (2022) used the code to estimate groundwater recharge rates using soil water isotope profiles in two contrasting dune types on Langeoog Island in Germany.

Double vegetation
Intercropping systems have been widely used worldwide due to their high economic benefits, field productivity, and land-, water-, and nitrogen-use efficiency.Intercropping is an agricultural practice where two crops (e.g., tomato and corn) are grown simultaneously in one field.The practice is often coupled with drip irrigation and mulching to increase crop yields and land-and water-use efficiencies.However, evaluating interspecies nutrition competition and the effects of different spatial arrangements of intercropping species on the spatial distributions of soil water content, temperature, and nitrogen concentration is much more complex than when a single crop species is present.The standard HYDRUS models consider only single vegetation with a single root zone and a single stress response function.The sink term in the Richards equation representing RWU for such a system is described using the macroscopic approach of Feddes et al. (1978): where S is RWU [T −1 ], α is the stress response function [-], b represents the spatial distribution of roots in the root zone [L −1 ], and T p is the potential transpiration [LT −1 ].The HYDRUS code was modified to consider the effects of multiple species and their competition for various resources (e.g., water and nutrients) as follows (Chen et al., 2020(Chen et al., , 2022)): where subscripts c and t refer to two vegetations (e.g., corn and tomato).Each crop has its own root distribution and potential transpiration rate, thus reflecting different growth cycles with different root growth and leaf area parameters.The model additionally allows the specification of different root water and solute uptake parameters for the different crops, which affects the spatial distribution of water and nutrients between irrigation cycles.While only the 1D formulation is given above, the same formulation was also implemented in the 2D code.
The modified HYDRUS code (its 2D version) with double vegetation was used successfully to evaluate soil nitrate dynamics and interspecies water and nitrogen competition in a tomato-corn intercropping system with different spatial arrangements (i.e., different number of rows of different species) by Chen et al. (2020Chen et al. ( , 2022)).

RWU models
The spatiotemporal pattern of RWU depends on multiple factors, including soil moisture availability and plant root biomass distribution.The nonuniformity of soil moisture and its impact on RWU is often mitigated by various coexisting processes, such as compensated root water uptake (CRWU) and hydraulic redistribution (HR).CRWU represents the abil-ity of plants to preferentially extract water from the wet zone to offset a reduction in RWU from the water-stressed zone (e.g., Jarvis, 1989;Šimůnek & Hopmans, 2009).On the other hand, HR refers to the passive transfer of soil moisture through the root system from a relatively wet soil zone to drier soil layers under reduced transpiration demand.RWU in most or all previous HYDRUS versions is described using the macroscopic approach of Feddes et al. (1978).The approach is now generalized in HYDRUS such that the model can simulate root water and solute uptake, account for the effects of both water and salinity stress on RWU, and additionally consider possible active and passive root contaminant or nutrient uptake (Šimůnek & Hopmans, 2009).Root water and solute uptake processes can be treated as being noncompensated or compensated (Jarvis, 1989), while users can select the degree of compensation (Šimůnek & Hopmans, 2009).Environmental stresses and compensation in this approach are considered purely empirical, with the model unable to simulate additional processes such as hydraulic lift.
To overcome the limitation of HYDRUS not being able to consider these processes and to demonstrate the importance of simultaneously considering CRWU and HR in estimating daily transpiration and RWU distributions under nonuniform soil moisture conditions, Thomas et al. (2024) implemented multiple mechanistic RWU models, including those of de Jong van Lier et al. (2006Lier et al. ( , 2008Lier et al. ( , 2009Lier et al. ( , 2013)), Couvreur et al. (2012), Couvreur, Vanderborght, Beff, et al. (2014), Couvreur, Vanderborght, Draye, et al. (2014), and Nimah andHanks (1973a, 1973b), in addition to the existing models of Jarvis (1989) andFeddes et al. (1978).The applications demonstrated much flexibility to account for the various processes involved.Implementing the models into the HYDRUS platform to simultaneously consider CRWU and HR and making them available to HYDRUS users should improve the accuracy of RWU predictions, and it will undoubtedly facilitate and encourage future research in this area.

Soil hydraulic properties
As noted earlier, SHPs can be described in HYDRUS using multiple analytical models for the water retention and hydraulic conductivity functions, including those by van Genuchten (1980), Brooks and Corey (1964), Durner (1994), and Kosugi (1996), as well as modified van Genuchten functions (Vogel & Císlerová, 1988).Alternatively, SHPs can be defined using lookup tables, thereby expanding HYDRUS's ability to consider other analytical models or more general shapes.
Much research has been carried out recently to derive analytical models that improve the description of SHPs in the very dry range (e.g., Iden et al., 2021;Peters, 2013;Peters et al., 2021;Seki et al., 2023;Weber et al., 2019Weber et al., , 2020)).While these models can be accommodated in HYDRUS using the input of lookup tables, we nevertheless included one such formulation directly in HYDRUS.In collaboration with Stathis Diamantopoulos and Tobias Weber, we implemented their generalized modular framework (Weber et al., 2019), which partitions the SHP functions into capillary and noncapillary (e.g., water in corners and/or films) parts.The water retention curve is, for this purpose, considered to be a weighted sum of a parametric capillary retention function, S c (h) (one of the analytical functions referred to above), and a new general model for the noncapillary saturation function, S nc (h) (Weber et al., 2019): where θ cs [-] and θ ncs [-] are the saturated water contents of the capillary and the noncapillary parts, respectively.The hydraulic conductivity function, K(h) [LT −1 ], is then obtained similarly.The capillary part of K(h) is obtained by adopting Mualem's integral, while the noncapillary part of K(h) is estimated directly from the new noncapillary saturation function (Weber et al., 2019), leading to where K c (h) is the hydraulic conductivity of liquid water in completely filled capillaries and K nc is the hydraulic conductivity of noncapillary (e.g., corner flow, film flow) liquid water [LT −1 ].Additionally, K rc (S c ) [-] is the relative hydraulic conductivity associated with completely filled capillaries, K rnc (S nc ) [-] is the saturation-dependent relative conductivity of the noncapillary part, and K s is the traditional total saturated hydraulic conductivity [LT −1 ], which is partitioned into the saturated capillary conductivity, K sc , and the noncapillary conductivity, K snc .
Since the generalized modular framework of Weber et al. ( 2019) can be combined with all soil hydraulic functions currently available in HYDRUS (i.e., Brooks & Corey, 1964;Kosugi, 1996;van Genuchten, 1980;Vogel & Císlerová, 1988), it provides significant additional flexibility to how SHPs can be described in HYDRUS.The resulting module currently has a working title, the Brunswick module.Lastly, a pedotransfer function is also provided to convert the parameters of van Genuchten (1980) to the parameters of the Weber et al. (2019) model (Weber et al., 2020).

Irrigation applications
The number and breadth of HYDRUS applications in the literature have increased dramatically in recent years and are much larger than expected when we first started developing the models.The type of applications is very broad, ranging from environmental applications simulating the transport of different solutes and particle-like substances to agricultural problems evaluating alternative irrigation schemes and the effects of plant covers on the soil water balance and groundwater recharge (Šimůnek et al., 2016).Irrigation applications have always been one of the most common HYDRUS applications.They have been the focus of hundreds of studies evaluating different irrigation methods (e.g., basin, sprinkler, furrow, and surface and subsurface drip), their scheduling (e.g., the timing of irrigation and its amount), and solute-related factors (e.g., fertigation, chemigation, salinization, and sodification) (Lazarovitch et al., 2023).For example, HYDRUS has been used in over 200 peerreviewed journal articles to assess various issues associated with drip irrigation.While a complete understanding of complex irrigation systems requires laborious, time-consuming, and expensive field investigations, which invariably involve only a limited number of treatments, fully calibrated processbased models such as HYDRUS can quickly evaluate different irrigation management strategies without the need for laborintensive fieldwork.The HYDRUS models for these reasons have become invaluable tools for studying complex and interactive water flow and solute transport processes in various irrigation systems.Lazarovitch et al. (2023) provided a recent review of the current modeling capabilities of HYDRUS to evaluate different irrigation methods.They provided not only a description of how HYDRUS can be used to model different irrigation techniques and solute-related topics such as fertigation, chemigation, and salinization/sodification, but also a long list of references of manuscripts in which HYDRUS was used for these applications.Other topics discussed in Lazarovitch et al. (2023) include special boundary conditions implemented in HYDRUS for various irrigation methods, optimization of irrigation systems, special irrigation methods, and the effects of spatial variability.The manuscript emphasizes the advantages and opportunities of using HYDRUS to describe basic processes in the root zone of irrigated plants that support sustainable irrigated agriculture.Especially valuable to HYDRUS users, the manuscript provides a link to a webpage from which all the HYDRUS project files of the discussed examples can be downloaded.
Another more narrowly focused review was written by Morianou et al. (2023).They looked at HYDRUS applications for simulations of water dynamics, RWU, and solute transport under drip irrigation of the four most common categories of tree crops (citrus, olive, avocado, and deciduous fruit/nuts).Morianou et al. (2023) identified 22 papers dealing with this topic, most of which involved surface drip irrigation, while three also focused on subsurface drip irrigation.

Low-impact developments
There has been an explosion of applications of the HYDRUS models to simulate flow and transport in various LIDs that have proven to be valuable alternatives for stormwater management and hydrological restoration of urban areas.LIDs include many different systems such as bioretention cells, infiltration wells/trenches, stormwater wetlands, wet ponds, level spreaders, permeable pavements, swales, green roofs, vegetated filter/buffer strips, sand and gravel filters, smaller culverts, or water harvesting systems (Brunetti et al., 2017).One of the limiting factors in the wide use of LIDs is a lack of adequate modeling tools that can be used to design LIDs that function properly for particular climate conditions.LID modeling requires an accurate description of multiple and interacting hydrological processes.Although several stormwater models can be applied to LID analyses, most of them do not incorporate accurate descriptions of the hydrological processes involved.The HYDRUS models may overcome this limitation.
Most LID applications with HYDRUS involved analyses of the overall hydrological functioning and effects of substrates used in green roofs, bioretention cells, and permeable pavements.Green roofs are complex systems consisting of a growing medium (usually a soil, which may be layered) and a waterproofing membrane, which are fully or partially covered with vegetation and may include additional features such as a drainage layer and an irrigation system.By retaining, delaying, and evapotranspirating stormwater, green roofs increase the resilience of urban areas against flooding.Their hydrological efficiency is strongly influenced by the hydraulic properties of the soil substrate, which is designed to be lightweight, to support vegetation, and to avoid surface water ponding.These features are generally achieved by mixing carbon-enriched soils with lightweight aggregates (e.g., volcanic lapillus or pumice), which confer a multimodal pore size distribution to the substrate.The resulting hydraulic behavior is strongly unsaturated and can be adequately described by the Richards equation.Many evaluated the effectiveness of green roofs using HYDRUS, including Brunetti, Papagrigoriou et al. (2020), Peczkowski et al. (2018), Sandoval et al. (2017), Soulis et al. (2017), andZhang, Lin et al. (2021).Several studied various factors affecting the hydrological functioning of green roofs, such as different substrates (Zhang, Lin et al., 2021), a storage layer (Li et al., 2019), storm events (Peng et al., 2019), compaction (Sandoval & Suárez, 2019), different substrate depths and vegetation covers (Soulis et al., 2017), having no vegetation (Brunetti, Porti, et al., 2018), and layering (Brunetti et al., 2016b;Wang et al., 2021).HYDRUS can be used to simulate not only water flow in green roofs, but also solute transport.Recently, Brunetti, Papagrigoriou, et al. (2021) used HYDRUS to reproduce observations from an experimental campaign about nitrogen turnover in green roofs irrigated with domestic wastewater.HYDRUS predicted high nitrogen turnover rates supported by the unsaturated conditions in the soil.
Strong unsaturated hydraulic behavior is also a common trait of permeable pavements, which is another LID consisting of a surface layer (e.g., porous concrete or turf blocks), a bedding layer (e.g., coarse sand), and a base layer (e.g., stones).The multilayered structure increases the stormwater retention capacity, delays subsurface runoff peaks, and avoids surface runoff.The hydrological efficiency of permeable pavements can be increased further by infiltrating subsurface runoff into the underlining soil, thereby reducing the pressure on urban drainage systems.Analyses of the hydrological functioning of permeable pavements with HYDRUS were carried out by Brunetti et al. (2016a), Brunetti, Šimůnek, Turco, et al. (2018), da Silva Ribas et al. (2021), and Turco et al. (2017), among others.
Finally, a bioretention cell is a stormwater device designed to capture and treat the first portion of runoff from impervious surfaces.The first flush into the cell often contains a large portion of the pollutants that leave an impermeable area.Examples of HYDRUS applications to evaluate the functioning of bioretention cells are given by Stewart et al. (2017), Li et al. (2021), Lisenbee et al. (2021), andHečková et al. (2022).

Managed aquifer recharge
Evaluating problems associated with MAR is another area with a dramatically increasing number of HYDRUS applications.Uncontrolled groundwater use worldwide has resulted in overexploitation, pollution, and other unintended consequences such as land subsidence, destruction of groundwaterdependent ecosystems, and alterations in the hydrological cycle.MAR, which refers to the intentional recharge of groundwater for future recovery or environmental benefit, has emerged as an important tool for achieving groundwater sustainability.MAR uses various technologies for intentional diversion, including the transport, storage, infiltration, and recharge of excess surface water (Sasidharan et al., 2021;Sprenger et al., 2017).
As an example, MAR can be accomplished using technologies such as aquifer storage and recovery, flooding land (Flood-MAR), flooding agricultural land (Ag-MAR), infiltration basins (IB), and vadose zone recharge devices such as stormwater disposal wells, vadose zone trenches, or drywells (Sasidharan et al., 2021).HYDRUS has been used to evaluate many of these technologies.The most common application has been studies of the hydrological functioning of drywells (e.g., Glass et al., 2020;Helles & Mogheir, 2023;Moreno et al., 2023;Sasidharan et al., 2018Sasidharan et al., , 2019Sasidharan et al., , 2020Sasidharan et al., , 2021) ) and infiltration basins (e.g., Cannavo et al., 2018;Glass et al., 2020;Lu et al. 2021;Nieec et al., 2021;Qi et al., 2021).The reservoir boundary condition in HYDRUS, for which the water level in a well or infiltration basin depends on the rate of water injection or pumping along with water infiltration or exfiltration, has most often been used in these applications (e.g., Glass et al., 2020;Sasidharan et al., 2018Sasidharan et al., , 2019)).
Various factors influencing infiltration from drywells and infiltration basins, such as spatial soil heterogeneity (e.g., Sasidharan et al., 2019Sasidharan et al., , 2020Sasidharan et al., , 2021) ) and clogging (Glass et al., 2020), have been evaluated in the above studies.Several have compared the hydrological functioning of drywells with infiltration basins (e.g., Glass et al., 2020;Sasidharan et al., 2021), the surface area required by these technologies (Sasidharan et al., 2021), or how infiltration from an infiltration basin can be enhanced substantially using drywell techniques (Helles & Mogheir, 2023;Sasidharan et al., 2021).While Sasidharan et al. ( 2021) considered a single drywell in the center of a circular infiltration basin, which only requires a 2D axisymmetrical formulation, Helles and Mogheir (2023) considered a more complicated 3D transport domain with an infiltration basin containing 18 drywells.
The references in this section have all used a multidimensional version of HYDRUS.A much higher number of MAR applications (not further given here) used HYDRUS-1D.The generic geochemical module PHREEQC in HPx was recently used to represent redox conditions and site geochemical evolutions during MAR, which made it possible also to account for alterations in the hydraulic conductivity due to clogging or field operations (Rodríguez-Escales et al., 2020).Given the increasing shortage of water resources worldwide, HYDRUS undoubtedly will be used in future studies evaluating current and new MAR technologies.

HYDRUS USE IN CLASSROOMS
Because of their ease of use and the advanced GUIs, the HYDRUS models have become popular tools for teaching the principles of water, solute, and heat movement in soils and groundwater in undergraduate-and graduate-level courses.
Even users with little direct knowledge of soil physics and related disciplines and limited mathematical expertise can benefit from the models.As a result, the HYDRUS software packages have been discussed extensively also in several soil physics and hydrology-related textbooks (e.g., Jury & Horton, 2004;Lazarovitch & Warrick, 2013;Radcliffe & Šimůnek, 2010;Rassam et al., 2003Rassam et al., , 2018;;Shukla, 2014) that can be readily used in university classrooms.
In 2003, Rassam et al. published a book that provided detailed guidance on using HYDRUS-2D to simulate variably saturated water flow and solute transport in soils.The book mainly focused on the technical aspects of defining transport domains, describing their discretization, and interpreting simulated results rather than explaining the underlying flow and transport processes in soils.Thus, the book was aimed at users familiar with the involved physical processes.The book was later translated into Japanese by Nobuo Toride and Mitsuhiro Inoue (Rassam et al., 2004).On the other hand, Soil Physics by Jury and Horton (2004) was likely the first textbook that attempted to explain basic flow and transport processes to students using the HYDRUS-1D software package.Radcliffe and Šimůnek (2010) later wrote a textbook called Soil Physics with Hydrus that covers different soil physics topics and is widely used in university classrooms worldwide.They used the HYDRUS family of programs to help students understand the basic topics better.While the RETC software was used to evaluate the unsaturated SHPs, HYDRUS-1D was used to demonstrate the processes of water infiltration, evaporation, and percolation in soils with varying textures and layers.They also showed how the software can be used for heat and solute transport in these systems, including how physical and chemical nonequilibrium conditions can affect water flow and contaminant transport.The HYDRUS (2D/3D) software was used further to describe 2D flow within field soils, hillslopes, boreholes, and capillary fringes.The software can also evaluate how different transport and reaction parameters affect solute transport.Using the information in this book, students can run HYDRUS and related models for different scenarios and parameters, thereby gaining more insight into the physics of water flow and contaminant transport.The book by Radcliff and Šimůnek (2010) is also suitable for self-studies on how to use the HYDRUS models since all examined problems can be downloaded from the HYDRUS website (https://www.pc-progress.com/en/Default.aspx?h3d-book).
A book called Exercises in Soil Physics was published by Lazarovitch and Warrick (2013) to complement existing soil physics and vadose zone hydrology textbooks by providing practical exercises.The text explores various topics related to soil physics, including solid-phase characteristics, soil water relations, saturated water flow, unsaturated flow, field water flow processes, chemical fate and transport, heat and energy transport, soil gases and their transport, and soil variability.Using HYDRUS-1D, the book further provides solutions to several problems related to variably saturated water flow and RWU.Additionally, STANMOD is used to solve solute transport problems involving sorbing, nonsorbing, degrading, nondegrading, and volatile solutes with different degrees of dispersion.Finally, the book uses ROSETTA and RETC to calculate the soil water retention curve and to determine the unsaturated soil hydraulic functions of van Genuchten (1980) and other soil hydraulic models.
In his textbook Soil Physics: An Introduction, Shukla (2014) published an extensive HYDRUS-1D tutorial focused on the unsaturated zone of a sandy loam furrow-irrigated onion field (Deb et al., 2011).He covered coupled liquid water, water vapor, and heat transport.Readers are provided a detailed description of the HYDRUS-1D input and output windows, including instructions on obtaining the necessary input parameters and interpreting the output.
In 2018, Rassam et al. published an introductory eBook, The HYDRUS-1D Software Package for Simulating the Movement of Water, Heat, and Multiple Solutes in Variably Saturated Media: Tutorial, containing a series of example problems for version 4.16 of HYDRUS-1D.The eBook aimed to provide HYDRUS-1D users with hands-on software package experience and to help them familiarize themselves with the overall organization of the GUI, including the main input and output dialog windows.The tutorial eBook discusses many examples of various complexities that can be easily downloaded from the HYDRUS website (https://www.pcprogress.com/en/Default.aspx?h1d-tut-TutorialBook).
A very recent contribution to the HYDRUS textbooks was provided in Chinese by Hu (2023).Entitled Theory and Applications of Environmental Soil Physics with the HYDRUS Model, the book includes theoretical developments and detailed instructions on using the HYDRUS-1D and HYDRUS (2D/3D) software packages.Although HYDRUS is already widely used in China, this textbook will make it even more accessible to the Chinese audience, and it will likely lead to its broader use in classrooms and for various applications in the country.

SUMMARY AND CONCLUSIONS
We believe that the HYDRUS suite of codes is well-suited to address the many environmental challenges humanity and our planet face, including soil and groundwater pollution problems, a lack of fresh water, soil salinization, climate change, and many other issues.The HYDRUS codes may be used for relatively simple 1D water flow and solute transport problems as well as for far more extensive multidimensional flow and transport applications at the field scale, including relatively complex problems involving a range of biogeochemical reactions.Examples of the latter include the HPx modules that couple the HYDRUS software packages with the PHREEQC geochemical code, which can consider a wide range of biogeochemical reactions, and the Wetland, UNSATCHEM, and DPU modules, which consider various more specific geochemical systems.Much attention has recently been paid especially to the PFAS module for evaluating the fate and transport of PFAS in the subsurface and the Particle Tracking module for determining recharge travel times or the average age of water taken up by plant roots.An important aspect of the HYDRUS models is that they are supported by relatively sophisticated interactive GUIs, making them easy to use and likely a major reason for their wide use.The GUI is continuously updated to make the codes as attractive as possible for nonexpert users and to facilitate classroom learning.Feedback from users is constantly being used to improve the codes, to identify particular strengths and weaknesses of the models, and to define additional processes or features that could be included in the models.We believe the codes and associated manuals are essential in transferring vadose zone research and development technologies to both the scientific community and practitioners.
For these reasons, we believe that the HYDRUS models have provided and will continue to provide an important service to the soil, environmental, and hydrological communities evaluating various fluid flow and biogeochemical processes in especially variably saturated subsurface systems.How to cite this article: Šimůnek, J., Brunetti, G., Jacques, D., van Genuchten, M. Th., & Šejna, M. (2024).Developments and applications of the HYDRUS computer software packages since 2016.Vadose Zone Journal, e20310.https://doi.org/10.1002/vzj2.20310APPENDIX A: PERSONAL NOTE BY THE SENIOR AUTHOR This VZJ special issue is a tribute to Rien van Genuchten, who, while a coauthor of this manuscript, has not been directly involved in writing this appendix.Allow me (a senior author of this article) to make some personal remarks on how Rien has influenced my career and those of many others.

AU T H O R C O N T R I B U T I O N S
My journey in numerical modeling started in the mid-1980s when I became part of a small group in Prague, at that time far behind the Iron Curtain, that was trying to develop one-and multidimensional, variably saturated water flow and solute transport models on a personal computer (an Atari ST).I am using "personal" in its literal meaning instead of mainframes and computers owned by institutions.The Atari was the personal property of Prof. Milena Císlerová, who enabled a few of us (including Tomáš Vogel and Radka Kodešová) to work on her computer to test various numerical schemes and techniques.This was when most governmental research institutions in Czechoslovakia had no or only minimal computer access.Our work was greatly helped by a visit from Rien in 1985 (I did not meet Rien at that time as I was in a conscription service in the army), who left with us many of his reprints and technical reports documenting his early models, as well as a box of punch cards with the source code.These reports documented in detail alternative numerical schemes (from finite differences to linear and Hermitian cubic finite element schemes with and without mass lumping) implemented in Rien's codes.I was an eager learner and tester of these various numerical schemes.
It is well known that Rien has been a mentor and promoter of many people in our research discipline, especially his numerous visitors, coworkers, and collaborators, relentlessly facilitating their continued career development in various ways.Especially remarkable was his support of scientists from non-Western countries that were not free or had just started opening their borders.He invited young scientists from Hungary, Czechoslovakia, East Germany, the Soviet Union, and Iran, providing them with new opportunities that were otherwise unavailable to them at the time.I wrote a letter to Rien 15391663, 0, Downloaded from https://acsess.onlinelibrary.wiley.com/doi/10.1002/vzj2.20310,Wiley Online Library on [15/02/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License in the late 80s, informing him about my work, and soon, I received an invitation to join his group in Riverside.Unfortunately, I did not obtain an exit visa (a special visa that the governments in Eastern Europe used at the time to control their population), and I had to wait until the communist government was overthrown in 1989 before I could embark on my trip to California in 1990.
Rien has been a great mentor to me since I arrived in the United States as a graduate student.He became my mentor and later a friend and a partner after I graduated.I continued to work with him as a postdoc, a research assistant, and eventually a professor.Since then, we have developed (with significant contributions from other coauthors of this manuscript, as well as some others) some of the most widely used numerical models in subsurface hydrology.Before join-ing the University of California, Riverside's faculty (in 2002), I traveled with Rien to all six continents to teach students and researchers worldwide how to use our numerical models and developed many friendships and collaborations.Our teaching job was tremendously helped by the work of Miroslav Šejna, my old high school friend from Czechia, who in the mid-90s developed the first versions of the sophisticated GUIs that made HYDRUS codes much easier and intuitive to use.The software then made its way very quickly into university classrooms and the offices of environmental consulting companies.
As documented by this manuscript, our close collaboration with Rien and Mirek (and later with Diederik and Giuseppe) continues until today.I will be forever grateful for Rien's mentorship and friendship in the past and future decades.

F
I G U R E 1 Schematic of the modeled transport domain in the Furrow module consisting of a trapezoidal open channel and a soil profile (Q inj is the solute injection rate [L 3 T −1 ], C inj is the solute concentration [ML −3 ], Q ir is the irrigation inflow rate [L 3 T −1 ], ET is evapotranspiration, P is precipitation, h is the water depth in the furrow [L], L is the length of the furrow [L], dx is the discretization length [L], h max is the maximum allowed water depth in the furrow [L], b w is the furrow bottom width [L], T w is the top furrow width [L], and SS is the side slope).

F
Particle trajectories calculated by the Particle Tracking module.F I G U R E 3 Particle statistics and travel times reported by the Particle Tracking module.15391663, 0, Downloaded from https://acsess.onlinelibrary.wiley.com/doi/10.1002/vzj2.20310,Wiley Online Library on [15/02/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Vadose Zone Journal F I G U R E 4 Water age as a function of time and depth calculated with the Particle Tracking module.F I G U R E 5 Precipitation and the average age of water taken up by plant roots, as calculated with the Particle Tracking module.RWU, root water uptake.

F
Schematic of the COSMIC module, showing the coupling of HYDRUS with the COSMIC operator.The soil profile consists of two soil horizons with variable water contents (the left bottom figure).The above ground neutron fluxes are given in N h −1 .to a certain depth and are scattered back into the atmosphere.The scattering process is regulated mainly by hydrogen, which implies that the fast neutron intensity at the near-ground level is directly related to the near-surface soil moisture content.While the penetration depth of CRNP measurements varies between 15 cm for wet soils and 55 cm for dry soils, their measurement footprint can reach a radius of up to 240 m, significantly larger than the typical spatial correlation length of soil moisture patterns.Due to this large sample volume, CRNP data can be used to understand the hydraulic behavior of the near-surface soil layer and infer effective SHPs for land surface models.The COSMIC module addresses this scientific need by coupling HYDRUS with the physically based COSMIC ofShuttleworth et al. (2013) (Figure6).The two models are coupled sequentially in time.At each time step, the soil moisture profile simulated by HYDRUS is passed to the COSMIC operator, which calculates the aboveground neutron intensity.This enables the COSMIC module to simulate the soil's transient hydrological behavior and a time series of aboveground neutron intensities.In its direct mode, the COSMIC module creates standard HYDRUS outputs and stores the time series of neutron intensities in an independent output file.In its inverse mode, COSMIC can be combined with transient CRNP data and other relevant information (e.g., soil moisture) to estimate soil hydraulic parameters inversely(Barbosa et al., 2021; Brunetti, Šimůnek et al., 2019).

F
Schematic of the modeled system in the dynamic plant uptake (DPU) module.