LAND AND WATER

The HYDRUS-1D Software Package for Simulating the
One-Dimensional Movement of Water, Heat, andMultiple Solutes in Variably-Saturated Media:
Tutorial
Version 1.00, July 2018

David Rassam, Jirka Šimůnek, Dirk Mallants,
andMartinus Th. van Genuchten

ISBN 978-1-4863-1001-2
CSIRO Land and Water

Citation
The HYDRUS-1D Software Package for Simulating the Movement of Water, Heat, and Multiple
Solutes in Variably Saturated Media: Tutorial, Version 1.00. By Rassam, D., Šimůnek, J., Mallants,
D., and M. Th. van Genuchten (2018). CSIRO Land and Water, Australia, pp. 183. – ONLINE
ISBN: 978-1-4863-1001-2.

Copyright© Commonwealth Scientific and Industrial Research Organisation 2018. To the extent permitted bylaw, all rights are reserved and no part of this publication covered by copyright may be reproduced
or copied in any form or by any means except with the written permission of CSIRO.

Cover image. Courtesy of Dr. Mat Gilfedder

Important disclaimer
CSIRO advises that the information contained in this publication comprises general statements
based on scientific research. The reader is advised and needs to be aware that such information may
be incomplete or unable to be used in any specific situation. No reliance or actions must therefore
be made on that information without seeking prior expert professional, scientific and technicaladvice. To the extent permitted by law, CSIRO (including its employees and consultants) excludes
all liability to any person for any consequences, including but not limited to all losses, damages,
costs, expenses and any other compensation, arising directly or indirectly from using thispublication (in part or in whole) and any information or material contained in it.
CSIRO is committed to providing web accessible content wherever possible. If you are having
difficulties with accessing this document please contact csiroenquiries@csiro.au.
This report documents version 4.16 of HYDRUS-1D, a software package for simulating water, heatand solute movement in one-dimensional variably-saturated media. The software has been verified
against a large number of test cases. However, no warranty is given that the program is completely
error-free. If you do encounter problems with the code, find errors, or have suggestions for
improvement, please contact:
Jirka Šimůnek
Department of Environmental Sciences
University of California Riverside
Riverside, CA 92521
USA
Phone/Fax: (951) 827-7854
Email: Jiri.Simunek@ucr.edu

Preface
This tutorial provides a series of example problems for version 4.16 of HYDRUS-1D, asoftware package for simulating water, heat and solute movement in one-dimensional variably
saturated media. The community of HYDRUS users has been growing steadily over the years
as evidenced in part by the number of downloads (over 10,000 times each of the past several
years) and visits to the HYDRUS web pages at https://www.pc-progress.com/ (on average
several hundred each day). The objective of this introductory HYDRUS-1D tutorial is to give
HYDRUS-1D users a first hands-on experience with the HYDRUS-1D software package andto familiarize themselves with the overall organization of the HYDRUS-1D graphical userinterface, including the main input and output dialog windows. Example problems include
variably-saturated water flow in single- and multi-layered soil profiles with and without root
water uptake. Both short- and long-term climatic boundary conditions are considered, with afocus on calculating soil water balance and groundwater recharge. Subsequent examples deal
with basic advective-dispersive solute transport without solid-liquid interactions, as well as
more advanced contaminant transport problems involving linear or non-linear equilibrium
sorption reactions and chemical non-equilibrium processes described using one-site sorption
kinetics. Physical non-equilibrium (mobile-immobile water type) solute transport is also
discussed, as well as the combined effect of physical and chemical non-equilibrium effects on
solute transport (i.e., a mobile-immobile type model with two-site sorption in the mobile zone).
Two inverse modelling examples consider one-step and multi-step outflow, respectively. The
final example is about coupled heat movement and reactive solute transport. All examples canbedownloaded
fromthe
HYDRUS
website(https://www.pcprogress.com/en/Default.aspx?H1d-tut-TutorialBook).
Additional supporting information is provided in an extended Appendix, which is based on
Frequently Asked Questions (FAQs) previously discussed on the HYDRUS-1D website. Thistutorial provides complementary information to the existing HYDRUS-1D Technical Manual(Šimůnek et al., 2008) and the more theoretical background documentation on soil physicalprocesses implemented in HYDRUS-1D as documented by Radcliffe and Šimůnek (2010). A
comprehensive overview of recent developments and applications of the HYDRUS-1D
software packages is available from Šimůnek et al. (2016).

i

ii

CONTENTS
Preface......................................................................................................................................... i
Introduction ................................................................................................................................ 1
Example 1: Infiltration of Water into a Single-Layered Soil Profile ......................................... 3
Example Description and Objectives ................................................................................. 3
Pre-Processing.................................................................................................................... 3
Post-Processing ................................................................................................................ 17
Example 2: Infiltration of Water into a Two-Layered Soil Profile .......................................... 27
Example Description and Objectives ............................................................................... 27
Pre-Processing.................................................................................................................. 27
Post-Processing ................................................................................................................ 31
Example 3: Root Water Uptake ............................................................................................... 37
Example Description and Objectives ............................................................................... 37
Pre-Processing.................................................................................................................. 37
Post-Processing ................................................................................................................ 42
Example 4: Estimating Groundwater Recharge ....................................................................... 49
Example Description and Objectives ............................................................................... 49
Pre-Processing.................................................................................................................. 49
Post-Processing ................................................................................................................ 53
Example 5a: Basic Solute Transport: Solute Introduced through a Constant Flux Boundary. 59
Example Description and Objectives ............................................................................... 59
Pre-Processing.................................................................................................................. 59
Post-Processing ................................................................................................................ 66
Example 5b: Basic Solute Transport: Solute Introduced through a Constant Pressure Head
Boundary .......................................................................................................................... 71
Example Description and Objectives ............................................................................... 71
Pre-Processing.................................................................................................................. 71
Post-Processing ................................................................................................................ 72
Example 5c: Basic Solute Transport: Profile with Non-Zero Initial Solute Concentration .... 75
Example Description and Objectives ............................................................................... 75
Pre-Processing.................................................................................................................. 75
Post-Processing ................................................................................................................ 76
Example 5d: Basic Solute Transport: The Effects of Solute Dispersion ................................. 79
Example Description and Objectives ............................................................................... 79
Pre-Processing.................................................................................................................. 79
Post-Processing ................................................................................................................ 80
Example 5e: Basic Solute Transport: The Effects of Diffusion .............................................. 83iii

Example Description and Objectives ............................................................................... 83
Pre-Processing.................................................................................................................. 83
Post-Processing ................................................................................................................ 86
Example 6a: Advanced Solute Transport: Equilibrium (instantaneous) Linear Sorption ....... 89
Example Description and Objectives ............................................................................... 89
Pre-Processing.................................................................................................................. 90
Post-Processing ................................................................................................................ 92
Example 6b: Advanced Solute Transport: Non-Linear Equilibrium Non-Linear Sorption ..... 97
Example Description and Objectives ............................................................................... 97
Pre-Processing.................................................................................................................. 98
Post-Processing .............................................................................................................. 100
Example 6c: Advanced Solute Transport: One-Site Chemical Nonequilibrium Sorption .... 103
Example Description and Objectives ............................................................................. 103
Pre-Processing................................................................................................................ 104
Post-Processing .............................................................................................................. 105
Example 6d: Advanced Solute Transport: Physical Nonequilibrium (Dual-Porosity)
Transport ........................................................................................................................ 109
Example Description and Objectives ............................................................................. 109
Pre-Processing................................................................................................................ 110
Post-Processing .............................................................................................................. 111
Example 6e: Advanced Solute Transport: Combined Physical and Chemical
Nonequilibrium .............................................................................................................. 117
Example Description and Objectives ............................................................................. 117
Pre-Processing................................................................................................................ 117
Post-Processing .............................................................................................................. 119
Example 7a: Inverse Modeling: One-Step Outflow Method ................................................. 127
Example Description and Objectives ............................................................................. 127
Pre-Processing................................................................................................................ 127
Post-Processing .............................................................................................................. 133
Example 7b: Inverse Modeling: Multi-Step Outflow Method ............................................... 135
Example Description and Objectives ............................................................................. 135
Pre-Processing................................................................................................................ 135
Post-Processing .............................................................................................................. 138
Example 8: Heat Movement and Reactive Solute Transport ................................................. 139
Example Description and Objectives ............................................................................. 139
Pre-Processing................................................................................................................ 139
Post-Processing .............................................................................................................. 143
Appendix I: Notes on Spatial and Temporal Discretization .................................................. 149
Temporal Discretization................................................................................................. 149
Spatial Discretization ..................................................................................................... 151
iv

Appendix II. Boundary Conditions ........................................................................................ 153
Water Flow Boundary Conditions ................................................................................. 154
Solute Transport Boundary Conditions ......................................................................... 165
Appendix III: Unsaturated Soil Hydraulic Properties ............................................................ 167
Water Retention Properties ............................................................................................ 167
Models Describing the SWRC ....................................................................................... 168
Hydraulic Conductivity Functions ................................................................................. 174
References .............................................................................................................................. 179v

vi

Introduction
This tutorial provides a series of example problems for version 4.16 of HYDRUS-1D,a software package for simulating water, heat and solute movement in one-dimensionalvariably saturated media. The software consists of several computation modules suchas the standard direct module (h1d_calc.exe), the standard inverse module(h1d_clci.exe), a dual-permeability direct module (h1d_dual.exe), a dual-permeability
inverse module (h1d_dlin.exe), the UnsatChem module for major ion chemistry andtransport (h1d_unsc.exe), and the HP1 module for multicomponent transport (hp1.exe),
as well as the HYDRUS-1D interactive graphics-based user interface.
The standard HYDRUS-1D module numerically solves the Richards equation forvariably-saturated water flow and advection-dispersion type equations for heat andsolute transport. The flow equation incorporates a sink term to account for water uptake
by plant roots. The flow equation may also consider dual-porosity type flow in which
one fraction of the water content is assumed to be mobile and another fraction
immobile, or dual-permeability type flow involving two mobile regions, one
representing the matrix and one the macropores or fractures. The heat transport equationconsiders transport due to conduction and convection with flowing water. Coupledwater, vapor, and energy transport can be considered as well. The solute transport
equations consider advective-dispersive transport in the liquid phase, as well as
diffusion in the gaseous phase. The transport equations also include provisions fornonlinear nonequilibrium reactions between the solid and liquid phases, linear
equilibrium reactions between the liquid and gaseous phases, zero-order production,and two first-order degradation reactions: one which is independent of other solutes,
and one which provides the coupling between solutes involved in sequential first-orderdecay reactions. In addition, physical nonequilibrium solute transport can be accounted
for by assuming a two-region, dual-porosity type formulation which partitions theliquid phase into mobile and immobile regions. Alternatively, the transport equations
include provisions for kinetic attachment/detachment of solute to the solid phase, andthus can be used to simulate the transport of viruses, colloids, or bacteria.
The objective of this introductory HYDRUS-1D tutorial is to give HYDRUS-1D users
hands-on experience with the software and to familiarize them with the overall
organization of the HYDRUS-1D graphical user interface, including the main input andoutput dialog windows. Example problems include variably-saturated water flow insingle- and multi-layered soil profiles with and without root water uptake. The
applications consider both short-term (1 and 30 days) and longer-term (1 year) climaticboundary conditions, with a focus on calculating soil water balance and groundwater
recharge. Subsequent examples deal with basic advective-dispersive solute transport
without solid-liquid interactions, as well as more advanced contaminant transportconsidering linear or non-linear equilibrium sorption reactions and chemical nonequilibrium described using one-site sorption kinetics. Physical non-equilibrium(mobile-immobile) solute transport is also discussed, as well as the combined effect ofphysical and chemical non-equilibrium effects on solute transport (i.e., a mobileimmobile type model with two-site sorption in the mobile zone). Two inverse modellingexamples consider one-step and multi-step outflow, respectively. The final example isabout coupled heat movement and reactive solute transport. All examples can be
downloaded
fromthe
HYDRUS-1D
website(https://www.pcprogress.com/en/Default.aspx?H1d-tut-TutorialBook).
1The following specific examples are considered in this tutorial:
1. Infiltration of water into a single-layered soil profile
2. Infiltration of water into a two-layered soil profile
3. Root water uptake
4. Estimating groundwater recharge
5. Basic solute transport
5a. Solute introduced through a constant flux boundary
5b. Solute introduced through a constant pressure head boundary
5c. Profile with non-zero initial concentration
5d. The effects of solute dispersion
5e. The effects of diffusion
6: Advanced solute transport
6a. Equilibrium (instantaneous) linear sorption
6b. Equilibrium (instantaneous) non-linear sorption
6c. Chemical non-equilibrium one-site sorption
6d. Physical non-equilibrium, dual-porosity (mobile-immobile water)
media
6e. Dual-porosity model with two-site sorption in the mobile zone
(physical and chemical non-equilibrium)
7. Inverse modeling
7a. One-step outflow method
7b. Multi-step outflow method
8. Heat and reactive solute transport (with a stagnant boundary condition forvolatile chemicals)

2

Example 1: Infiltration of Water into a Single-Layered Soil
ProfileExample Description and Objectives
This example problem involves infiltration of water from a ponded soil surface into a
1-m deep single-layered loam soil. The soil profile is initially unsaturated, having an
initial pressure head of -100 cm. Water infiltrates from the saturated soil surface, whichis represented by a ‘Constant Pressure Head’ boundary condition by assuming that thepressure head at the soil surface is 1 cm (corresponding to the ponding depth). Water
drains from the bottom of the soil profile by gravity (represented by a ‘Free Drainage’boundary condition) since the groundwater table is at an unspecified point deep in the
profile (e.g., more than 10 m). This introductory example demonstrates the advance of
the wetting front throughout the soil profile, examines surface and bottom fluxes, runtime information (such as time steps and a number of iterations), mass balances, andother input and output information.

Pre-Processing
Project Manager
The first step when working with any HYDRUS-1D project involves opening the
Project Manager (Fig. 1.1) by clicking on the Project Data Manager Tab in the upper
left corner (third button from the left) of the displayed Graphical User Interface (GUI).
The Project Manager is a tool to manage data of existing projects and helps to locate,
open, copy, delete and/or rename desired projects or their input or output data. A project
represents any particular problem to be solved using HYDRUS-1D. The project name,as well as a brief description of the project, helps to locate a particular problem. Input
and output data for each project are placed in a folder with the same name as the project.
Projects are represented by a file project_name.h1d and the project_name folder.
The current version of the Project Manager gives you considerable freedom inorganizing your projects. The projects are grouped into Work Groups, which can beplaced anywhere in accessible memory (such as on local and/or network hard drives).
The Work Groups can be any existing accessible subdirectory (folder). HYDRUS-1D
is installed together with four default workspaces: Direct, Inverse, Unsatch, and HP1,
which are located in the HYDRUS-1D folder. The Direct and Inverse work groups
contain test examples for direct and inverse problems, respectively, while the Unsatchand HP1 work groups contain examples for carbon dioxide/major ion chemistry andHP1 biogeochemical examples, respectively. We suggest that you create your own
Work Groups (e.g., the My Direct and My Inverse Work Groups) and keep the preprogrammed examples intact for future reference. Projects can be copied with the
Project Manager only within a particular Work Group. However, you can copy
projects between Work Groups using standard file managing software, e.g., Windows
Explorer. In that case, you must copy both the folder of a particular project and the
project_name.h1d file. Another way of copying a project between Work Groups isfirst to open the project, and then using the command Save As (File menu) to save thatproject to a new location.

3

The Project Manager has two Tabs (in the upper left corner): Project Groups andProjects. The Project Groups Tab displays currently defined Project Groups, and the
Project Tab displays projects in a selected particular Project Group. The ProjectManager can be obtained either from the File Menu (the command Project Manager)
or the Toolbar (the third button).

Figure 1.1. The Project Manager dialog window.
- Open the Project Manager (from the File Menu or the third button on the
Toolbar,).
- Select the Projects Tab of the Project Manager.
- Click button “New” and the New Project dialog window (Fig. 1.2) will open.
This command allows you to create a new project with a given name, provides abrief project description, and indicates where the project files will be located.
- Type the Name of the project: “Example 1”- Type the Description of the project: “Infiltration of water into a soil profile”.
- Click button “OK”. A default HYDRUS-1D project is now created, presumably
for simulating water flow in a one-dimensional 1-m deep soil profile. Thisdefault project now can be modified for a specific problem.

4

Figure 1.2. The New Project dialog window.
Main Processes
Figure 1.3 shows the Main Processes dialog window, which is the first window you willencounter when starting a new project. This window will show the processes that can beselected for a particular simulation. When selecting Water Flow (see the mathematicaldescription in Chapter 2 of the Technical Manual1 (Šimůnek et al., 2008)), you can also
consider Vapor Flow and accumulation/melting of Snow at the soil surface. Whenselecting Solute Transport, you can either use the Standard Solute Transport module(Chapter 3), the more specialized Major Ion Chemistry module (Chapter 6), or thebiogeochemical module HP1 (which resulted from coupling HYDRUS-1D with the
PHREEQC biogeochemical code (Jacques and Šimůnek, 2005; Jacques et al., 2006,2018)). Because degradation and other transport and reaction processes are temperaturedependent (usually degradation increases with temperature), one can also consider Heat
Transport (Chapter 4) and the effects of temperature on solute reactions (such as
contaminant degradation (EFSA, 2008) and organic matter decomposition (Jacques et al.,2018)). When plants are present, one should also consider Root Water Uptake and Root
Growth (Chapter 2.2 and 3.5). The UnsatChem module can be used to simulate carbon
dioxide transport (considering diffusion in both liquid and gas phases and convection in
the liquid phase) and production (CO2 Transport) (Chapter 5), as well as major ion solutetransport (Chapter 6). Finally, you can select whether HYDRUS-1D should be run in for
a forward prediction, or in the Inverse Mode (Chapter 9), in which case various
parameters can be calibrated against experimental data and/or experimental data can bedisplayed together with the simulation results.
- Type Heading: “Infiltration of water into a soil profile”- Click the Button “Next” (or press n of the keyboard)
From this point on the program will navigate you through the entire process of entering
input information. Users may either select particular commands from a menu (or the preprocessing part of the desktop) or allow the interface to lead them through the process ofentering input data by selecting the Next button. Alternatively, clicking the Previous
button will return users to the previous window (since the Main Processes dialog window
is the first window, the Previous button is not present).
1

Subsequent references to a particular Chapter refer to the HYDRUS-1D Technical Manual (Šimůnek
et al., 2008). Much of the information can be accessed immediately also by simply pressing F1, or byclicking on the “Help” button (and then, among other options, clicking on “Help on …” and dragging“?” to a particular pre- or post-processing part of the window).

5

Figure 1.3. The Main Processes dialog window.
Geometry Information
The second dialog window one encounters is called Geometry Information (Fig. 1.4).
In this window, one selects the Length Unit (mm, cm, or m), the Number of Soil
Materials that are present in the system (i.e., the number of soil horizons), the Numberof Layers over which a Mass Balance will be calculated, whether flow is vertical or
horizontal (i.e., a Decline from Vertical Axis of 1 indicates flow in a vertical direction,
while 0 indicates flow a horizontal direction), and the Depth of the Soil Profile. Note thatwhen the Length Unit is changed later, all input parameters for a particular project will
be converted automatically from the original units into the new units.

Figure 1.4. The Geometry Information dialog window.
- Length Units: cm- Number of Soil Materials (typically corresponding to the soil horizons): 1(indicates a homogenous soil profile). The maximum allowed number is 100.
- Number of Layers for Mass Balance: 1 (you may select up to 10 mass balance
layers; the default option is a single mass balance layer, which corresponds tothe entire soil profile. You may wish to assign individual mass balance layers to6each soil material (horizon) or a group of soil materials, or apply separate massbalances the root zone and the remainder of the soil profile.)- Decline from vertical axis: 1 (vertical flow)- Depth of the Soil Profile (cm): 100 cm- Click Button “Next”
Time Information
The third dialog window users encounter is called Time Information (Fig. 1.5). Here you
select the Time Unit to be used throughout the application (i.e., years, days, hours, min,or sec), when the simulation starts (Initial Time) and ends (Final Time), and how thetime should be discretized (i.e., Initial, Minimum, and Maximum Time Steps). We
recommend to simply use the default parameters for time stepping in most applications.
The boundary conditions will either be constant in time, or one can specify TimeVariable Boundary Conditions, such as daily values of precipitation, evaporation, andtranspiration. Evapotranspiration fluxes additionally can be evaluated from various
Meteorological Data using the Penman-Monteith combination equation
recommended by FAO (Monteith, 1981, 1990) or the Hargreaves Formula (Jensen etal., 1997). Time-variable boundary conditions can be either constant during specified
time intervals or be subject to Daily Variations (of transpiration or meteorologicaldata) to be Generated by HYDRUS-1D using predefined interpolation functions.
Finally, a surface Energy Balance can be evaluated from specified meteorological data.
Further details on temporal discretization can be found in Appendix I.

Figure 1.5. The Time Information dialog window.
7

- Time Units: Days
Time Discretization:- Initial Time (day): 0- Final Time (day): 1- Initial Time Step (day): 0.001- Minimum Time Step (day): 0.00001- Maximum Time Step (day): 5 (you do not need to change this even though it islarger than the final time).
- Click Button “Next”
Print Information
The fourth dialog window HYDRUS-1D is called Print Information (Fig. 1.6). In thiswindow, you simply select how often the output of the model should be printed. Output
can be printed at every nth Time Step (if 1 then after each time step), at a predetermined
Regular Time Interval, or at specified Print Times, which then will be specified in thenext dialog window.

Figure 1.6. The Print Information dialog window.
- Print Options: Leave as Default- Number of Print Times: 12- Click Button “Select Print Times.” This command will open a Print-Times dialog
window (Fig. 1.7), which by default will list uniformly distributed Print Times atwhich detailed information about pressure heads, water contents, concentrations,water and solute fluxes, and soil-water and solute balances will be printed. Thesedefault Print Times can be edited as needed. Clicking on the Default command
button will cause the print times to be distributed evenly between the initial andfinal time. Clicking on the Default (log) command button will cause the print times
to be distributed between the initial and final time evenly on a log scale. This option
is recommended only for larger times. When using this Default (log) option, thefirst Print Time is often extremely small. One should make sure that the Initial and
Minimum Time Steps (in the Time Information dialog window) are smaller than
this first Print Time.
8

- Click Button “Next”

Figure 1.7. The Print-Times dialog window.
The first four windows are general windows (with general information), which you needfor any application of HYDRUS-1D. All of the dialog windows to be discussed next deal
with a) water flow, b) heat transport, c) solute transport, d) root water uptake, and/or otherinformation.
Water Flow – Iteration Criteria
The Iteration Criteria dialog window is used to specify iteration criteria, which arerequired for the numerical solution of the Richards equation (the equation governing waterflow in the unsaturated zone) and various Time Step Control factors, which depend on
the numerical solution (Fig. 1.8). In general, one should not change any of these values,
which control the fully automated self-adjusting time-stepping procedure of the numericalsolution. The default values are based on many years of experience using the HYDRUS
models and their predecessors for a range of applications. More details about iteration
criteria and time step control can be found in Appendix I.
Internal Interpolation Tables of water contents, hydraulic conductivities, and specificwater capacities are generated from the specified set of unsaturated soil hydraulicparameters. Values of the hydraulic properties are then computed for all numericalnodes during the iterative solution process using linear interpolation between the entries
in these tables. This interpolation approach was found to be much fastercomputationally than direct evaluation of the hydraulic functions over the entire rangeof pressure heads. Interpolation tables can be avoided by setting both the Upper andLower Limits of the Tension Interval to zero. To avoid problems with the numericalsolution, the Internal Interpolation Tables should not be used when the dual-porosity
model of Durner (1994) is used to describe the soil hydraulic functions since this model
does not have a monotonous hydraulic capacity function. When the Brooks and Corey(1964) model is used, either the lower limit should be larger (in absolute value) than1/α (i.e., the air-entry value) or the tables should be avoided altogether (e.g., by fixing
the lower and upper limits at zero).

9Figure 1.8. The Iteration Criteria dialog window.
- Leave default values- Click Button “Next”
Water Flow – Soil Hydraulic ModelHYDRUS-1D allows you to consider different analytical models to describe the soil
hydraulic properties (i.e., water retention and hydraulic conductivity functions). The
unsaturated soil hydraulic properties can be described in HYDRUS-1D using Brooks andCorey (1964), van Genuchten-Mualem (van Genuchten, 1976), modified van Genuchten(Vogel and Císlerová, 1988), Durner (1994), and Kosugi (1996) type analytical
functions. The different analytical models can be selected in the Soil Hydraulic Modeldialog window (Fig. 1.8). They are described in detail in Chapter 2.4 of Šimůnek et al.
(2008) and summarized further in Appendix III. Or, alternatively, you can click on F1
for more information.
The Soil Hydraulic Model window allows you to select additionally various models
for simulating nonequilibrium or preferential flow, such as dual-porosity or dualpermeability models as summarized by Šimůnek and van Genuchten (2008), andwhether or not to consider hysteresis. Hysteresis is a property of soils to have different
wetting and drying branches of the soil water retention curve, including scanning curves(Luckner et al., 1989; Lenhard and Parker, 1992; Vogel et al., 1996). One may further
limit hysteresis to only the soil water retention curve, or also consider hysteresis in theunsaturated hydraulic conductivity (which generally shows less effects of hysteresis
when plotted versus water content)

10Figure 1.9. The Soil Hydraulic Model dialog window.
- Leave as Default (Single Porosity, Models, van Genuchten – Mualem)- Click Button “Next”
Water Flow – Soil Hydraulic Parameters
The Soil Hydraulic Parameters for the hydraulic property model selected in theprevious dialog window can be specified in the dialog window shown in Figure 1.11.

Figure 1.10. The Water Flow Parameters dialog window.
11

It is always preferable to use soil hydraulic parameters based on actual measurements,for example by analyzing observed soil water retention and/or hydraulic conductivity
data using the RETC code (van Genuchten et al., 1991), which can be downloaded fromthe HYDRUS website. Unfortunately, this is not possible when data are limited or
lacking. For this reason, we implemented in the Water Flow Parameters window a
Soil Catalog, which provides average soil hydraulic parameters for twelve different
USDA soil textural classes as derived by Carsel and Parrish (1988). One can open the
Soil Catalog and simply select the desired soil texture (which in our current example
is actually loam). The hydraulic parameters of that soil type are then immediatelytransferred to the table in Fig. 1.10. Users should be aware that the average parametersin the Soil Catalog were derived for agricultural soils and may not be appropriate forrelatively deep soil horizons.
Additionally, we also implemented in the Water Flow Parameters window the optionof selecting pedotransfer functions, which relate easily measured soil properties suchas texture and bulk density with the hydraulic parameters. The pedotransfer functions
used in HYDRUS-1D are those developed by Schaap et al. (2001) as incorporated into
their Rosetta module (Neural Network Prediction; Fig. 1.11). While hydraulicparameters obtained with Rosetta may be superior to those from the Soil Catalog, we
decided to keep the soil catalog, mostly for historical reasons.

Figure 1.11. The Rosetta Lite dialog window.
Finally, the temperature dependence of soil hydraulic parameters can be considered by
selecting the Temperature Dependence checkbox in the Water Flow Parameters
window. The influence of temperature on the soil water pressure head is thenquantitatively predicted from the effect of temperature on surface tension (Philip andde Vries 1957), while the influence of temperature on the hydraulic conductivity is
predicted from the effect of temperature on viscosity and the density of water(Constantz, 1982).

12

- θr (residual water content, units depend on the selected length unit; here
cm3/cm3)- θs (saturated water content, units depend on the selected length unit; cm3/cm3)- Alpha (the hydraulic shape parameter α, units depend on the selected length
unit; here 1/cm)- n (the hydraulic parameter n, dimensionless)- Ks (the saturated hydraulic conductivity, units depend on the selected length andtime units; here cm/day)- l (Mualem’s pore connectivity exponent; the default value of 0.5 is used here)- Select Loam from the Soil Catalogue (the default soil texture in HYDRUS-1D)- Click Button “Next”
Water Flow – Boundary Conditions
The Water Flow Boundary Conditions window (Fig. 1.12) allows you to selectboundary conditions for the top and bottom of the soil profile, depending upon theapplication More information on the meaning and application of various boundary
conditions are found in Appendix II (and/or see Chapter 2.7 or press the F1 key).
Various constant or time-variable, pressure heads or water fluxes, as well as
atmospheric boundary conditions, can be used at the soil surface. For the bottom of thesoil profile, you can additionally select free or deep drainage boundary conditions, aseepage face, or horizontal drain conditions. Additionally, in the Initial Condition partof the dialog window, you must specify whether the initial conditions for water flow
are to be given later in terms of the pressure head or the water content.
The water table in the current example is assumed to be located far below the bottomof the soil domain and thus does not affect flow processes in the adopted soil profile(e.g., by capillary rise). This default BC (Free drainage) corresponds to gravity flow
(i.e., flow according to a total unit hydraulic gradient).

Figure 1.12. The Water Flow Boundary Conditions dialog window.
13

- Upper Boundary Condition: Constant Pressure Head (see Appendix II.1.1)- Lower Boundary Condition: Free Drainage (see Appendix II.1.7)- Initial Condition: In Pressure Heads- Click Button “Next”
HYDRUS-1D Guide
The next dialog window (Fig. 1.13) will ask “Do you want to run PROFILE
application?” When selecting OK, the external Profile module, in which users can
specify spatially variable properties (such as material distribution, initial conditions,root water uptake distribution, observation nodes, etc.), will be opened.

Figure 1.13. The HYDRUS-1D Guide dialog window.
- Click Button “OK”
At this point, you will get an additional question "Do you want to save data before
running the Graphical Profile Editor?" Please answer yes since you want the externalProfile module to reflect the information entered thus far in the main module, such as
which processes were selected, what the depth of the soil profile is, and what type ofthe boundary conditions will be used.
Soil Profile – Graphical Editor
HYDRUS-1D discretizes by default 2 the soil profile into 100 finite elements (FE)involving 101 nodes. Since the profile was assumed to be 100 cm, the elements all havea size of 1 cm. This default discretization is appropriate for this simple example.
However, there are other problems for which you may need to refine the FE grid atcertain locations along the modeled soil profile. More details on the spatial
discretization can be found in Appendix I.
- Menu: Conditions->Initial Conditions->Pressure Head, or press the Red Arrow(i.e.,) on the Toolbar- Click the “Edit Condition” Button on the left-side Edit Bar, select with your
Mouse the top node (the selected node should now have a yellow color), thenleft-click the mouse again, and specify a pressure head of 1 cm. This is now the
top boundary condition corresponding to the ponding depth on the soil surface.
Note that when a constant pressure head boundary condition is used, its value isentered into the initial condition and then kept constant during the simulation.
All other nodes will remain at the default initial pressure head of -100 cm.
2

Since the default number of finite elements is always 100, the element size will depend on theassumed total profile length.

14

- Menu: Conditions->Observation Nodes, or pressof the Toolbar
Note that when you move the cursor around the view window, you can monitor
its position, i.e., the node numbers and the depth, at the lower right corner.
- Click the “Insert” Button on the left-side Edit Bar, and specify nodes at 20, 40,
60, 80, and 100 cm. Use a left click to insert an observation node, and a rightclick to end the command (hence you can add several observation nodes after
the initial left click).
- Observation nodes can be deleted using the "Delete" and "Delete All"
commands of the Edit Bar. In addition to serving as locations for which theoutput (pressure head, water content, temperature and concentration at each timestep) is printed, observation nodes also provide necessary output required forspecific applications. For example, observation nodes are used to identify thelocation of experimental data points used in inverse modeling projects, or thelocation of a threshold pressure head for triggered irrigation (see Appendix II).
- Menu: File->Save Data, or use Toolbar- Menu: File->Exit
Soil Profile – Summary
After exiting the "Soil Profile" Editor, the Soil Profile Summary dialog window (Fig.
1.14) appears. This window provides tabular information about the spatial
discretization ("z [cm]"), initial conditions ("h [cm]"), spatial root distribution ("Root[1/cm]"), scaling factors ("Axz", "Bxz", "Dxz"), and material distributions (Mat) withdepth, among other data. This table is adjusted automatically to accommodate
concentration and temperature data if the solute and heat transport options had been
activated in the Main Processes window. In the window, Axz is the nodal value of thedimensionless scaling factor αh [-] associated with the pressure head, Bxz is the nodal
value of the dimensionless scaling factor αΚ [-] associated with the saturated hydraulicconductivity, and Dxz is the nodal value of the dimensionless scaling factor αθ [-]associated with the water content (see Chapter 2.4 or press F1).
The command "Set to Default Values", if used, would set all values in this table to its
default values. The command "Set Initial Conditions equal to Field Capacity" would
set the initial values equal to Field Capacity. Field capacity in HYDRUS-1D is definedas the water content or pressure head corresponding to the hydraulic conductivity of
about 0.01 cm/d (Twarakavi et al., 2009).
- Click Button “Next”
HYDRUS-1D Guide
The next dialog window (Fig. 1.15) will now appear and ask “Do you want to run
HYDRUS-1D application?” When you select OK, another dialog window appears andasks "Do you want to save the input data before executing HYDRUS-1D?". You
should always answer these questions positively by clicking on "Yes", except whenrealizing that some errors may have been made and should be corrected (then click on"No"). Once “Yes” is chosen, the external computational H1D_Calc.exe module will
be executed.

15

Figure 1.14. The Soil Profile Summary dialog window.

Figure 1.15. The HYDRUS-1D Guide dialog window.
- Click Button “OK”
- Alternatively, you can run the computational module using the Menu Command:
Calculation->Execute HYDRUS-1D or a Toolbar button.

Execute HYDRUS-1D
The computations are now being carried out, with parts of the results being displayed
in a DOS window (Fig. 1.16). This window will show selected run-time information,such as time (Time), the number of iterations (ItW), the cumulative number of iteration(ItCum), the surface flux (vTop), the cumulative surface flux (SvTop), cumulative root
water uptake (SvRoot), the cumulative bottom flux (SvBot), the surface pressure head
(hTop), the average pressure head in the root zone (hRoot), and the bottom pressure
head (hBot). When the simulation ends, the Run Time in seconds is displayed, as wellas the note “Calculations have finished successfully”. You are now asked to “Press
Enter to Continue”. Once the keyboard button “Enter” is pressed, the DOS window is
closed and the program returns to the HYDRUS-1D GUI. Detailed output information
16

will not be shown in the DOS window if the user does not check the “Screen Output”
box in the Print Information dialogue window.

Figure 1.16. The DOS window in which the computational module h1d_calc.exe is
executed.

Post-Processing
Once the calculations are finished and the “Enter” button is pressed, the “PostProcessing” section (the right side) of the main window of the HYDRUS-1D GUI will
be populated with commands that can be used to display selected results, which will be
discussed below ((Fig. 1.17)). Program output consists of a number of output files thatare described in detail in Chapter 13 of the technical manual. These output files arelocated in the project folder, the location of which is indicated in the Project Managerof the HYDRUS-1D GUI. The output files are in standard text format and as such canbe opened by any standard text editing programs, such as Notepad or MS Excel.

Figure 1.17. The main window of the HYDRUS-1D GUI with both pre- and postprocessing commands.
17

Observation Points
Pressure heads (Fig. 1.18), water contents, and water fluxes at selected observation
nodes (i.e., at depths of 20, 40, 60, 80, and 100 cm) are displayed under the command“Observation Nodes” at the post-processing section of the HYDRUS-1D GUI (orusing the menu command “Results->Observation Nodes).
Referring to Figure 1.18:• Note how the pressure heads evolve over time at the selected observation points.
This output option is useful when one requires the model to produce predictions
at a specific location to compare with field observations (such as observedpressure heads from a tensiometer or water contents from a TDR probe) and isparticularly relevant when conducting inverse simulations (parameter
identification or model calibration).
• Note that the wetting front passes the selected depths of 20, 60, and 100 cm attimes of approximately 0.08, 0.35, and 0.65 days, respectively. This timing is
dictated by the hydraulic properties of the loam soil.
• When the wetting front approaches a certain depth, the pressure head increases
from the initial value of -100 cm and eventually reaches 1 cm (not shown in Fig.
1.18, where the maximum displayed pressure head is zero), which is equal tothe imposed positive pressure head at the surface.
• The colors in Figure 1.18 show which curve belong to which node (i.e., N1, N2,…). This connection is not immediately shown initially but appears when you
right-click on the figure, select “Legend…”, and then click on “Visible”. Manyother options can be used by right-clicking on the figure, including copying thefigure. You can also alter the labels along the axes (e.g., to use different fonts
or font sizes), or change the scales of the axes by right-clicking on the axis itself.
These changes are possible with all post-processing figures that are generated.

Figure 1.18. Pressure heads versus time at selected observation points (N1 - 20 cm,
N2 - 40 cm, N3 - 60 cm, N4 - 80 cm, and N5 - 100 cm).

18

Profile Information
Pressure heads (Fig. 1.19), water contents (Fig. 1.20), water fluxes, hydraulicconductivities, and hydraulic capacities versus depth at selected Print Times (i.e., every
two hours in this example) are displayed using the command “Profile Information” ofthe post-processing section of the HYDRUS-1D GUI (or using the menu command“Results->Profile Information).
•••

Figure 1.19 shows the advance of the wetting front at selected print times. For
example, after 0.08 days (2 hours), the wetting boundary had reached a depth ofabout 20 cm (similar results are shown in Figure 1.18 for the first observation
point N1).
The wetting front advances in time until it reaches the bottom of the profile atabout 0.75 days. The entire profile then attains a pressure head of 1 cm.
Figure 1.20 portrays the same information for the water content, with the initialwater content equal to 0.24 (as defined by the initial pressure head of -100 cm
and the soil water retention curve of the loam soil). Note that once the wettingboundary passes a certain point, the water content becomes equal to thesaturated water content (θs) of 0.43 (see Water Flow Parameters: Qs=0.43).

Figure 1.19. Pressure head profiles at selected print times (from 0.083 to 0.75 days).

Water Flow - Boundary Fluxes and Heads
Various water fluxes (potential and actual surface fluxes, potential and actual root water
uptake rates, and bottom fluxes, as well as their cumulative values) and pressure heads(at the surface and bottom of the profile, and average root zone pressure heads) versus
time are displayed under the command “Water Flow - Boundary Fluxes and Heads”
of the post-processing section of the HYDRUS-1D GUI (or using the menu command“Results->Water Flow - Boundary Fluxes and Heads") (Fig. 1.21).
19

Figure 1.20. Water content profiles at selected times (from 0.083 to 0.75 days).

Referring to Figure 1.21:• Note that HYDRUS-1D reports downward fluxes (against the direction of thez-axis) as negative and upward fluxes (in the direction of the z-axis) as positive.
This is due to the convention in the numerical solution that the z-axis is positive
upward.
• Water starts infiltrating at the soil surface due to the 1 cm positive pressure head
imposed at the “Constant Head Boundary” and continues to do infiltrate during
the entire simulation.
• The lower “Free Drainage Boundary” allows unit gradient gravity drainage.
However, the initial magnitude of this drainage is so small that it does not show
in this figure at early times (before the arrival of the wetting front).
• As soon as the wetting front reaches the lower boundary (after about 0.70 days;
see also the pressure head profiles in Figure 1.18), the bottom flux increases
significantly. Note that the two lines representing the upper and lower boundary
fluxes become parallel (i.e., inflow = outflow), and the water balance of the soil
profile reaches a steady-state.

20

Figure 1.21. The surface flux (negative and into the soil profile) and the bottom flux
(negative and out of the soil profile) versus time.

Soil Hydraulic Properties
Graphs of the unsaturated soil hydraulic functions (various combinations of pressure
heads, water contents, hydraulic conductivities, and soil water capacities) for soils inthe soil profile are displayed using the command “Soil Hydraulic Properties” of the
post-processing section of the HYDRUS-1D GUI (or using the menu command“Results->Soil Hydraulic Properties") (Fig. 1.22).
•••

Figure 1.22 shows the water retention curve for the selected “Loam” soil. Notethat the horizontal axis is the logarithm of the absolute value of the pressure
head, thus representing matric suctions in terms of pF units (one cannot take thelogarithm of a negative number). Values of 0, 1, 2, 3, and 4 hence represent
pressure heads equal to -1, -10, -100, -1000, and -10,000 cm, respectively). Plots
using pF values often provide a better display of the entire soil water retention
curve, including values near saturation (such as between -1 and -10 pressure
heads in the current example).
Note from Figure 1.22 that the water content at a pressure head of -100 cm (i.e.,
pF=log(|h|) = 2), which is the initial condition, is equal to about 0.24. This valueof 0.24 is the initial water content in Figure 1.20).
The “Soil Hydraulic Properties” command allows users to view the soil water
retention and the hydraulic conductivity functions of the soil or soils being
studies, among other combinations of soil hydraulic properties (the list boxes of
the commend provide many graphing options, including the use of regular of
logarithmic scales).

21

Figure 1.22. Soil water retention curve.

Run Time Information
Selected information about the numerical solution (such as temporal changes in the time
step, the number of iterations necessary to solve the Richards equation at any particulartime level, the cumulative number of iterations, and the dimensionless Peclet and
Courant numbers when solute transport is considered) versus either time or the numberof time-steps are displayed under the command “Run Time Information” of the postprocessing section of the HYDRUS-1D GUI. One could for this purpose use the menu
command “Results->Run Time Information" (Fig. 1.23).
••

Figure 1.23 shows how the time step dynamically changed during thesimulation. The algorithm used to dynamically adjust the time step based on thenumber of iterations needed to solve the Richards equation at any particular time
step is described in Chapter 7.4.2 of the Technical Manual. Note that after about
0.65 days, the time step increases significantly. This reflects the fact that thesolution is then found much easier (less nonlinear) since the profile approaches
steady-state, and HYDRUS-1D hence can use larger time steps.
Similarly, Figure 1.24 shows that after about 0.67 days, the number of iterations
required to achieve a solution drops to the minimum allowed number of 2. Notethat very large spikes in the number of iterations may indicate numericalinstability, thus warranting a thorough check of the solution integrity (e.g.,check mass balance errors).

22

Figure 1.23. Time steps versus time.

Figure 1.24. The number of iterations versus time.

Mass Balance Information
The output file Balance.out containing the mass balance information at selected print
times is displayed under the command “Mass Balance Information” of the postprocessing section of the HYDRUS-1D GUI (or, again, also using the menu command“Results->Mass Balance Information").
Mass balance information is provided for the entire soil profile and each selected“Subregion” at every “Print Time.” Since there is only one subregion in this example,
the data in columns 3 and 4 are identical.

23

Column 1 refers to the following:
Length:The depth of the soil profile
W-volume: The total volume of soil water in the soil profile or selected
subregions
In-flow:
The sum of fluxes flowing into and out of the soil profile (orsubregion)h Mean:The average pressure head along the soil profile (or subregion)
Top Flux: The flux into (at the top) the soil profile
Bot Flux: The flux out of (at the bottom) the soil profile
WatBalT: The total (absolute) water balance error for the entire soil profile
WatBalR: The relative water balance error in % (see Chapter 7.44)
Column 2 gives the units for each variable. The units used are displayed at the
beginning of the Balance.out file. In this particular example, L is in cm and T in days.
•••

Table 1.1 shows that at times up to 0.667 days, the flux into the soil profile ismuch larger than out of the soil profile.
Consequently, “In-flow” is positive, which means that water is being stored in
the soil profile.
Table 1.2 shows that at 0.75 days, “In-flow” is zero, thus indicating that a
steady-state has been attained (inflow is equal to outflow). This statement isconfirmed by the fact that “h Mean” is 1 cm (equal to the imposed upperboundary condition), and the profile is fully saturated (W-volume = Qs x Length= 0.43 x 100 cm = 43 cm). Furthermore, all of the mass balance information for
T > 0.75 is identical.

Table 1.1. Mass balance information at Time=0.6667 days.
Time (T)0.6667

Sub-region number
LengthW-volume
In-flowh MeanTop Flux
Bot Flux
WatBalT
WatBalR

1
[L]
[L]
[L/T]
[L]
[L/T]
[L/T]
[L]
[%]

24

1.00E+02
4.27E+01
2.43E+01-9.23E-01-2.53E+01-4.95E+00
1.91E-05
0.00

1.00E+02
4.27E+01
2.43E+01-9.23E-01

Table 1.2. Mass balance information at Time=0.75 days.
Time (T)0.75

Sub-region number
LengthW-volume
In-flowh MeanTop Flux
Bot Flux
WatBalT
WatBalR

1
[L]
[L]
[L/T]
[L]
[L/T]
[L/T]
[L]
[%]

25

1.00E+02
4.30E+01
0.00E+00
1.00E+00-2.496E+01-2.496E+01
2.86E-05
0.00

1.00E+02
4.30E+01
0.00E+00
1.00E+00

26

Example 2: Infiltration of Water into a Two-Layered Soil
ProfileExample Description and Objectives
The objective of this example is to simulate water flow in a 1-m deep, two-layered soil
profile. The example demonstrates how to specify time-variable atmospheric boundary
conditions (daily values of precipitation and evaporation) and how to interpret
simulated results, such as actual boundary fluxes. The soil profile consists of a 50-cm
soil layer of clay loam underlain by a 50-cm soil layer of sandy loam. The initiallyunsaturated soil profile (having a default initial pressure head = -100 cm) is wetted via
precipitation during the first 7 days of the simulation. Precipitation changes from 6 cm/dduring the first 2 days, to 10 cm/d during the next three days, and finally to 2 cm/dduring the last 2 days of the first week. The ‘Atmospheric Boundary Condition with a
Surface Layer’ is invoked, which allows accumulation or ponding of water on the soil
surface, with surface runoff occurring only if a user-defined threshold (the maximum
ponding depth) has been exceeded. In this example, a different lower boundary
condition is imposed, i.e., a ‘Seepage Face' boundary condition. Water starts draining
from the soil profile through the bottom ' Seepage Face' boundary condition only oncethe bottom of the profile reaches full saturation (the pressure head becomes zero). As
the simulations will show, the soil profile continues to drain until day 10, after which
evaporation (1 cm/d) is activated until the end of the 20-day simulation. Note that thisis a hypothetical example; neglecting evaporation on certain days allows us to focus on
one process at a time (e.g., rainfall, ponding, evaporation).
This example demonstrates the advancement of the wetting front through the twolayered soil profile, calculates prevailing fluxes, and examines the mass balance of eachof the two soil layers. The relevance of the pressure head at the seepage face ishighlighted. The concept of surface ponding due to a potential flux that exceeds thesaturated hydraulic conductivity of the surface soil layer is demonstrated. A suitable
Hydraulic Property Model is chosen to cope with the ponding process.

Pre-Processing
Please create 'Example 2' by copying and modifying 'Example 1' using the ProjectManager.
Project Manager- Select Projects Tab- Select Example 1- Click Button “Copy”- Type New Name: Example 2- Type Description: Infiltration of water into a 2-layer soil profile- Click Button “OK”
- Select Example 2- Click Button “Open”
Main Processes- Heading: Infiltration of water into a 2-layer soil profile27

- Click Button “Next”
Geometry Information- Number of Soil Materials: 2 (indicates a two-layered soil profile)- Number of Layers for Mass Balance: 2 (mass balances will be calculated foreach soil layer)- Click Button “Next”
Time Information
Time Discretization:- Final Time (day): 20- Maximum Time Step (day): 0.1 (We do not recommend to use a large time stepwhen a surface water layer is allowed to develop (i.e., for saturated near-surface
conditions). A larger maximum time step may lead to overshoot in the
maximum allowed water layer depth.)- Check Time-Variable Boundary Conditions (this allows you to specify timevariable precipitation and evaporation fluxes)- Number of Time-Variable Boundary Records: 5- Click Button “Next”
Print Information- Number of Print Times: 6- Button “Select Print Times”- Print Times: 2, 4, 5, 6, 10, 20 days
- Click Button “OK”
- Click Button “Next”
Water Flow – Iteration Criteria- Leave default values- Click Button “Next”
Water Flow – Soil Hydraulic Model- Select the van Genuchten - Mualem model With Air-Entry Value of -2cm (thismodel is recommended for very fine textured materials, such as silt loam or
clayey soils, or for soils having a van Genuchten n parameter smaller than about
1.2; please refer to Appendix III.1.2 and III.2).
- Click Button “Next”
Water Flow – Soil Hydraulic Parameters- Click on the first row: Mat 1- Select from the Soil Catalogue: Clay Loam
- Click on the second row: Mat 2- Select from the Soil Catalogue: Sandy Loam
- Click Button “Next”
Water Flow – Boundary Conditions- Upper Boundary Condition: Atmospheric BC with Surface Layer.
This boundary condition requires you later to specify time-variable surface
fluxes (e.g., rainfall and evaporation rates), while allowing a surface water layer28to build up on the soil surface when some of the rainfall water cannot infiltrate
into the soil profile, before surface runoff starts (see Appendix II.1.5).
- Lower Boundary Condition: Seepage Face; h=0This boundary condition allows water to drain from the soil profile once thepressure head at the bottom of the soil profile reaches a value of 0 (i.e., whenfull saturation is reached) (see Appendix II.1.7).
- Initial Condition: In Pressure Heads- Max h at Soil Surface: 3 cm
Surface runoff starts only once the depth of the surface water layer reaches 3cm.
- Click Button “Next”
Time Variable Boundary Conditions- Fill out the table in the Time-Variable Boundary Conditions dialog window asshown in Figure 2.1.
- The Precip. Column. This column indicates that the rainfall rate is 6 cm/dduring the first two days (between time 0 and 2 d), 10 cm/d during the next 3
days (between time 2 and 5 d), 2 cm/d during the following 2 days (between 5
and 7 d), and no rainfall during the remainder of the simulation (between times 7
and 20 d). If precipitation data are available from a weather station, then theyneed to be converted into rates with units corresponding to those used in
HYDRUS-1D as defined during geometry and time information preprocessing(more details can be found in Appendix II.1.5.
- The Evap. column. This column indicates that there is no evaporation during thefirst 10 days and evaporation of 1 cm/d between day 10 and 20. We assume that
evaporation is negligible during (first seven days) and immediately after (nextthree days) the rainfall event.
- The hCritA column. This parameter indicates that actual evaporation remains
equal to potential evaporation until the pressure head in the top node reaches 50,000 cm. Once the hCritA value is reached, the actual evaporation rate
becomes lower than the potential value since the soil is then too dry to deliver
the potential rate. The hCritA parameter (given here in absolute values) is not
used in any other calculation, except for evaluating evaporation rates. Note thatyou can click on the corner of a cell and copy values downward by dragging the
cursor, similarly as in an Excel spreadsheet; in this case, you can do this for thehCritA column]- Click Button “Next”

29

Figure 2.1. The Time-Variable Boundary Conditions dialog window.
HYDRUS-1D Guide- Do you want to run PROFILE application?
- Click Button “OK”
Soil Profile – Graphical Editor
The Soil Profile editor allows you to specify the initial pressure head condition, thedistribution of soil materials and soil layers versus depth for mass balance calculations,and the location of observation nodes.
- Menu: Conditions->Initial Conditions->Pressure Head or Toolbar- Click the “Edit Condition” Button on the left-side Edit Bar, select with the
Mouse the entire soil profile, left click the mouse again, and specify a pressure
head of -100 cm.
(The entire profile will be at the default initial pressure head of -100 cm)- Menu: Conditions->Material Distribution or Toolbar- Click the “Edit Condition” Button on the left-side Edit Bar, select with the
Mouse node 52 to 101 (left click at node 52, drag the mouse to the bottom node,and left click again).
- In the "Condition Specification" dialog window, set the Material Index equal to
2 and click OK (Material 2 corresponds to the second layer, the sandy loam;
Material Index 1 corresponds to the clay loam layer).
- Menu: Conditions->Subregions or Toolbar- Repeat the same process (as for Material Distribution); set the Subregion Index
equal to 2 for nodes between 52 and 101 and click OK. In this way the Mass
Balance will be calculated separately for Material 1 and 2, together with theMass Balance for the entire profile.
- Menu: Conditions->Observation Nodes or Toolbar- Click the “Insert” Button on the left-side Edit Bar and select nodes at 0, 50, and100 cm (right-click to end the command)- Menu: File->Save Data or Toolbar
30

- Menu: File->Exit
Soil Profile – Summary- Click Button “Next”
Execute HYDRUS-1D- Menu: Calculation->Execute HYDRUS-1D or use the Toolbar button.

Post-Processing
Observation Points
Referring to Figure 2.2:
• Note that the pressure head at the surface node (ON-1 at 0 cm) becomes positive
after two days when the precipitation rate of 10 cm/day exceeds the saturated
hydraulic conductivity of the upper Clay Loam layer (6.24 cm/day). The
pressure head (reflecting the thickness of the surface water layer) continues toincrease up to the preset limit of 3 cm, at which time surface runoff is triggered(day three). At this stage, the boundary condition switches to a “Constant Head"boundary condition and surface runoff is calculated subsequently as thedifference between the precipitation and infiltration rates (both instantaneousand cumulative runoff is stored as output and can be displayed upon request).
Note that during the first time step when the maximum thickness of the waterlayer is reached, the pressure head may slightly rise above 3 cm. HYDRUS-1D
will then force it back to the prescribed 3 cm during the next time step. Thisbreach can be limited by setting a smaller “Maximum Time Step” (say, 0.01
day).
• For times between day 3 and day 5, the profile has attained a steady-state
condition (all three pressure heads are constant; horizontal lines).
• When the precipitation rate drops below the saturated hydraulic conductivity ofthe surface layer (clay loam, 6.24 cm/day) at day 5, all of the precipitation stored
in the surface water layer gradually infiltrates into the soil profile and the soil
surface becomes unsaturated again (at about 5.5 d).
• After the cessation of precipitation (Time > 7 d) and before the activation of
evaporation (at day 11), the pressure heads in the soil profile continue todecrease as a result of drainage due to gravity. During this period, the lower
boundary (ON-3 at 100 cm) remains saturated, and water continues draining outof the bottom of the soil profile through the seepage face.
• At day 11, evaporation is activated, and hence the soil profile starts drying out
at a much faster rate. The pressure head at the soil surface eventually reaches
the minimum allowed pressure head ‘hCritA’ of -50,000 cm (not shown in
Figure 2.2). The drying process progresses down through the soil profile until itreaches the lower Seepage Face boundary at about 13.4 days. Then, the bottom
boundary becomes unsaturated, and hence the seepage flux ceases.

31

Rainfall

6 cm/d10 cm/dEvaporation

2 cm/d1 cm/dFigure 2.2. Pressure heads versus time at selected observation points. The upper partof the figure shows the applied boundary conditions.

Profile Information
Referring to Figure 2.3:
• The pressure head profile at day 2 indicates that the lower boundary has not yet
reached full saturation (pressure head < 0), which means that there is no bottom
flux.
• The pressure head profiles at days 4 and 5 are identical, which confirms that a
steady- state has indeed been attained during this period.
• On the other hand, the bottom pressure head is equal to zero from day 4 to day
10, indicating draining of the soil profile through the bottom boundary.
Referring to Figure 2.4:
• Note the abrupt change in the water content at the interface between the two soil
types. While the pressure head across the material interface is continuous(Figure 2.3), the water content is not so due to the different water retention
properties (and thus retention curves) of the two materials. Also note that while
the pressure head across the material interface is continuous, its gradient aboveand below the interface is discontinuous due to different conductivities of thetwo materials. The pressure head gradients at the material interface need tocompensate the different hydraulic conductivities of the two materials since thewater flux is the same when approaching the interface from above and below.
32

Figure 2.3. Pressure head profiles at selected times. The boundary between layer 1(clay loam) and layer 2 (sandy loam) is at -50 cm.

Figure 2.4. Water content profiles at selected times. The boundary between layer 1(clay loam) and layer 2 (sandy loam) is at -50 cm.

33

Water Flow - Boundary Fluxes and Heads
Referring to Figure 2.5:
• Note that HYDRUS-1D reports downward fluxes (against the upward directionof the z-axis) as negative and upward fluxes (in the direction of the z-axis) as
positive.
• The surface flux between days 0 and 2 is equal to the applied precipitation rateof 6 cm/day since this rate is lower than the saturated hydraulic conductivity ofthe surface soil layer (clay loam, 6.24 cm/day).
• Once the precipitation rate is higher than the saturated hydraulic conductivityof the surface soil layer, the soil surface quickly reaches full saturation, andexcess water (the difference between the rainfall and the infiltration rates) starts
to accumulate at the soil surface.
• Before the surface becomes fully saturated (at time > 2 d), the infiltration rate
can be higher than the saturated hydraulic conductivity of the surface layer
because of the capillary component of the total gradient.
• Once full saturation is reached, and the surface water layer gradually increases
until it reaches its maximum thickness of 3 cm (between 2 and 3 d), theinfiltration rate gradually increases as well because of the additional
contribution to the total gradient by the hydrostatic pressure due to the surface
water layer.
• The infiltration flux stabilizes once the maximum thickness of the surface waterlayer (3 cm) is reached.
• Between 4 and 5 days, the flux into the soil profile is equal to the flux out of thesoil profile, which confirms that the soil profile has attained a steady-state.
• The infiltration rate starts decreasing when the precipitation rate becomes
smaller than the saturated hydraulic conductivity (Ks) but remains
approximately equal to Ks as water from the surface layer infiltrates.
• As soon as the bottom of the profile becomes saturated (between day 2-3),
drainage commences. Water continues flowing through the lower ‘Seepage
Face’ boundary as the profile continues to drain due to gravity.
• The seepage flux ceases at 13.4 days when the pressure head at the lower
boundary becomes negative (desaturation).
• The actual surface evaporative flux (positive) is equal to the applied potential
evaporation rate of 1 cm/day up to about 11.2 days (as long as the surfacepressure head is higher than hCritA), after which it drops as the soil surface
dries out and the hydraulic conductivity can no longer sustain this evaporative
flux.

34

Rainfall

6 cm/d10 cm/dEvaporation

2 cm/d1 cm/dFigure 2.5. The surface (into the soil profile) and bottom (out of the soil profile)
fluxes versus time. The upper part of the figure shows the applied boundary
conditions.
Mass Balance Information•

Table 2.1 shows the mass balance information for the layered soil profile. Note
that, in addition to the mass balance for the entire profile (the first data column),
a separate mass balance is provided for the two subregions: the top 50 cm (thesecond data column) and the bottom 50 cm (the third data column).
Table 2.1. The mass balance information for Time = 10 days.
Time
day
10
Sub-region
number1
2
Lengthcm
100
50
50
W-volume
cm
31.8
18.1
13.7
In-flow
cm/day-0.348-0.194-0.154h Mean
cm
-38.6
-55.1
-22.2
Top Flux
cm/day-0.00225
Bot Flux
cm/day-0.348
WatBalTcm-1.62E-05
WatBalR
[%]
0.00

35

36

Example 3: Root Water Uptake
Example Description and Objectives
The objective of this example is to simulate water flow and root water uptake in a 3-m
deep uniform soil profile and to demonstrate how to simulate root water uptake andassociated crop transpiration processes. The soil profile is initially at hydrostatic
equilibrium, i.e., the pressure head changes uniformly from -300 cm at the top of thesoil profile to 0 cm at the bottom. The profile is wetted via precipitation through an‘Atmospheric Boundary with Surface Runoff’ during the first three days, whileevaporation and transpiration (Evapotranspiration, ET) are active during the entire 30day simulation. The bottom boundary is a ‘Constant Head Boundary’, with a 0-cm
constant pressure head representing the groundwater table. Roots are uniformly
distributed from the surface down to a depth of 200 cm, which is still 100 cm above thegroundwater table.
This example simulates root water uptake and its reduction due to the water stress when
the soil can no longer supply water to the plant roots at an optimal rate. Project results
will show how pressure heads and water fluxes evolve differently at various depths inthe soil profile as a result of the combined processes of infiltration, evaporation, andtranspiration. The default (uniform) spatial discretization is changed so as todemonstrate how to specify a finer grid at some locations. The phenomenon of‘groundwater ET’, whereby the plant sources its water by capillary action from theunderlying groundwater aquifer (despite the fact that roots are well above the watertable), is demonstrated using this example.

Pre-Processing
Project Manager- Open the Project Manager (from the File Menu or the third button on the
Toolbar,).
- Select Projects Tab- Click button “New”- Type the Name of the project: “Example 3”- Type the Description of the project: Root water uptake in a soil profile- Click Button “OK”
Main Processes- Heading: Root water uptake in a soil profile- Check the Root Water Uptake checkbox (in addition to already checked Water
Flow)- Click Button “Next”
Geometry Information- Length Units: cm- Number of Soil Materials: 1 (indicates a homogenous soil profile)- Number of Layers for Mass Balance: 1 (one may select up to 10 mass balance
layers)- Decline from vertical axis: 1 (vertical flow)
37

- Depth of Soil Profile (cm): 300 cm- Click Button “Next”
Time Information- Time Units: Days
Time Discretization:- Initial Time (day): 0- Final Time (day): 30- Initial Time Step (day): 0.001- Minimum Time Step (day): 0.00001- Maximum Time Step (day): 5- Check Time-Variable Boundary Conditions- Number of Time-Variable Records 2- Click Button “Next”
Print Information- Print Options: Leave as Default- Number of Print Times: 11- Click Button “Select Print Times”- Print Times (day): 1, 2, 3, 4, 5, 6, 7, 8, 10, 15, 30- Click Button “OK”
- Click Button “Next”
Water Flow – Iteration Criteria- Leave default values- Click Button “Next”
Water Flow – Soil Hydraulic Model- Check the Single Porosity model: van Genuchten – Mualem- Check "With an Air-Entry Value of -2 cm" (this model is recommended for finetextured materials, such as clayey soils, or for soils with the n parameter smaller
than about 1.2)
- Click Button “Next”
Water Flow – Soil Hydraulic Parameters- Select Clay from the Soil Catalogue- Click Button “Next”
Water Flow – Boundary Conditions- Upper Boundary Condition: Atmospheric BC with Surface Run Off(surface runoff is initiated immediately once the precipitation rate exceed theinfiltration rate. At the same time, the flux boundary condition is replaced withthe pressure head boundary condition with zero pressure head. Under thisboundary, condition excess precipitation does not accumulate at the surface).
- Lower Boundary Condition: Constant Pressure Head(this boundary condition with a zero pressure head will indicate that the
groundwater level is at the bottom of the soil profile and is constant in time).
- Initial Condition: In Pressure Heads- Click Button “Next”

38Root Water and Solute Uptake ModelIn the Root Water and Solute Uptake Model dialog window (Fig. 3.1.) you can selectwhich model to use to account for the effects of water stress on root water uptake,either the Feddes (1977) model or an S-shape function (van Genuchten, 1987) can beused, and whether or not (No Solute Stress) should be considered to account for salinity
stress. Salinity stress can be either additive (Additive Model) or multiplicative(Multiplicative Model), and either the threshold and slope model (Maas, 1990)(Threshold Model) or the S-shaped (S-Shape) function (van Genuchten, 1987) can beused. It is also possible to consider either non-compensated (the Critical Stress Index =1) or compensated (the Critical Stress Index < 1) root water uptake. You canadditionally consider Active Solute Uptake (Šimůnek and Hopmans, 2009) andspecify related parameters. Root water and solute uptake are described in detail in
Chapter 2.2 and 3.5, respectively, of the Technical Manual.
Root Water Uptake Reduction Model- Check Feddes (Default)- Click Button “Next”

Figure 3.1. The Root Water and Solute Uptake Model dialog window.
Root Water Uptake Parameters
Parameters for the water stress response functions are specified in the Root Water
Uptake Parameters dialog window (Fig. 3.2). The Root Water Uptake Parameters
for the water stress response function suggested by Feddes et al. (1978) are describedin detail in the HYDRUS-1D technical manual. Water uptake in this model is assumed
to be zero close to saturation (i.e., wetter than some arbitrary “anaerobiosis point” P0).
Root water uptake is also zero for pressure heads less (more negative) than the wilting
point (P3). Water uptake is considered optimal between pressure heads Popt and P2,
whereas for pressure heads between P2 and P3 (or P0 and Popt) water uptake decreases39(or increases) linearly with pressure head. Root water uptake occurs at its optimal valueeven at full saturation, when both P0 and POpt are set equal to zero.
- Select Pasture (default) from the Database (all parameters will then be
populated automatically with values from the implemented database)- Click Button “Next”

Figure 3.2. The Root Water Uptake Parameters dialog window.
Time Variable Boundary Conditions
In this dialog window, specify a rainfall rate of 2 cm/d during the first two days, anevaporation rate of 0.1 cm/d during all 30 days of simulation, and a transpiration rateof 0.5 cm/d during the first 3 and 1.5 cm/d during the subsequent 27 days.
- Row 1: Time=3 day; Precip=2 cm/day; Evap=0.1 cm/day; hCritA=50,000 cm;
Transp=0.5 cm/day- Row 2: Time=30 day; Precip=0 cm/day; Evap=0.1 cm/day; hCricA=50,000 cm;
Transp=1.5 cm/day- Click Button “Next”
HYDRUS-1D Guide- Do you want to run PROFILE application?
- Click Button “OK”
Soil Profile – Graphical Editor
You can use the Soil Profile editor to adjust the spatial discretization of the soil profile(in this case smaller elements at the top and bottom boundaries, with the element size
increasing by a geometric factor from both top and bottom towards the middle of the
profile), specify the initial pressure head conditions, the root distribution with depth,
and the position of observation nodes.
- Menu: Conditions->Profile Discretization or Toolbar- Click the “Number” button on the left-side Edit Bar: Number of Nodes = 120.
- Click the “Insert Fixed” button: Select the Node at a depth of 150 cm (themidpoint of the profile). These points are used to specify the local nodal density.
40

By default, they are placed at the upper and lower boundaries, but in addition
can be inserted at, or deleted from, any other point in the soil profile.
(Note that when you move the cursor around the view window, you can monitor
its position, i.e., the node numbers and the depth, at the lower right corner)- Click the “Density” button, select the top node, left click the mouse button, in
the Nodal Density dialog window enter the Lower Density of 0.1, and click the
OK button.
- Click the “Density” button again, select the bottom node, left click the mousebutton, in the Nodal Density dialog window enter the Upper Density of 0.3, andclick the OK button.
Note that the FE nodes will get automatically redistributed proportionally basedon the specified densities. The distance between nodes at the top will be (1/0.1=)
10 times smaller than in the middle, while the distance between nodes at thebottom will be (1/0.3=) 3.33 times smaller than in the middle. The nodal density
is a real number between 0.01 and 100, specifying the local density of nodes.
The density can be specified only at fixed points (either at the top or bottom
boundary or at an additional user-defined “Insert Fixed” node). The program
distinguishes between top and bottom density. At the top (bottom) boundary,
only the Lower Density (Upper Density) can be specified. For any inserted fixed
node, the top (bottom) density at a fixed point specifies the relative thickness ofthe elements above (or below) this point. If the top and bottom densities areequal, then the elements above and below the fixed node will have the same
thickness. If the top and bottom density values are different, then the element
thicknesses will be different as well.
- Menu: Conditions->Initial Conditions->Pressure Head or Toolbar- Click the “Edit Condition” button (left side), select with the Mouse the entire
profile, left click the mouse button, in the Condition Specification dialog
window that appears enter Pressure Head Distribution: Top value=-300 cm andBottom value=0 cm, uncheck ‘Use top value for both’, and click the OK button.
- Menu: Conditions->Root Distribution or Toolbar- Click the “Edit Condition” button (left side), select with the Mouse the nodes
down to a depth of 200 cm, left click the mouse again, in the Condition
Specification dialog window that appears enter Roots Linear Distribution: Topvalue=1, check ‘Use top value for both’, and click the OK button. This meansthat roots are distributed uniformly in the top 200 cm of the soil profile.
Typical root distribution functions are shown in Figure 3.3. Note that HYDRUS-1D
can accommodate any arbitrarily measured or prescribed root distribution function
by using the Root Distribution function.
- Menu: Conditions->Observation Nodes or Toolbar- Click the “Insert” Button on the left-side Edit Bar and specify nodes at 0 and
300 cm (right-click to end the command)- Menu: File->Save Data or Toolbar- Menu: File->Exit

41

Figure 3.3. Schematic of possible root distribution functions (Kumar et al., 2014).
Soil Profile – Summary- Click Button “Next”
Execute HYDRUS-1D- Menu: Calculation->Execute HYDRUS-1D, or a Toolbar button.

Post-Processing
Observation Points
Referring to Figure 3.4:
• The pressure head in the surface node (ON-1 at 0 cm), which was initially at apressure head of -300 cm, increases during the first three days of the simulation
as a result of precipitation and infiltration. Subsequently, the pressure headdecreases down to the minimum limit of hCritA (-50,000 cm) as a result of
evaporation (the minimum value shown in Figure 3.4 is -500 cm).
• The pressure head at the bottom node (ON-2 at 300 cm) remains constant at 0
cm (i.e., the Constant Pressure Head boundary condition representing the
groundwater table).
• The full pressure head range (-50,000-0 cm) and total simulation time (0-30 d)
are not shown in Figure 3.4 in order to highlight the changes at early times.

42Figure 3.4. Pressure heads versus time at selected observation points.

Profile Information
Referring to Figure 3.5:
• The pressure head profiles at days 1 and 3 show the typical advance of thewetting front during the precipitation period, up to day 3.
• The pressure head at day 6 displays a drying profile due to the combined effectsof transpiration and evaporation; drying is most pronounced at the soil surface
due to the evaporation boundary. The effect of transpiration on the soil pressure
head and water content distributions are depth-dependent and visualized bymeans of the root water uptake function.
Referring to Figure 3.6:
• As a result of the low transpiration rate (0.5 cm/day), root water uptake at t = 1
day is limited to just above 0.002 day-1; water uptake is maximal from about 25 cm down to the bottom of the roots (-200 cm).
• When transpiration increases at the beginning of day 3 from 0.5 to 1.5 cm/day,root water uptake in the top of the profile is suboptimal owing to the pressure
head being greater (i.e., less negative) than the optimal pressure head (Popt = 25 cm; see Figure 3.2).
• At time = 6 and 10 days, root water uptake has tripled following the increasein transpiration. At t = 6 days water uptake is nearly uniform in the entire rootzone. As the profile starts drying out at the top (see pressure head evolution in
Figure 3.5), water uptake will progressively be from deeper depths.
Referring to Figure 3.7:
• Time = 1 day: Near the surface, the downward (negative) flux reflects
infiltration due to precipitation (the surface flux equals evaporationprecipitation: 0.1-2=-1.9 cm/d) with the wetting front approaching a depth of
100 cm. However, below a depth of 130 cm, there is a small positive upward
flux from the constant head boundary (the water table) to the soil profile driven
by plant uptake. Note the zero-flux point at a depth of about 130 cm.
43

•

•

•

Time = 3 days: The wetting front advances deeper into the soil profile andprovides plants with its water demand (and hence there is no upward positive
flux from the water table at this time). Water infiltrating through the upper"atmospheric" boundary (due to precipitation) drains through the bottom"constant pressure head" boundary.
Time = 5 days: The upper 70 cm of the profile is drying out due to evaporation.
As a result, a net positive upward flux occurs out of the soil profile through theupper "atmospheric" boundary. As the upper part of the profile gets drier, roots
start to source their water via an upwards flux from the bottom "constant head"boundary (note the positive flux below a depth of 150 cm).
Time = 30 days: The top 150 cm of the soil profiles becomes so dry that the
upward flux practically ceases and the pressure head at the soil surface reaches
the limiting hCritA value (-50,000 cm). Transpiration (root water uptake) alsoceases in this part of the profile where the pressure head drops to or below thewilting point (P3) of -8000 cm (note the pressure head distribution in Figure3.8). It is the high pressure head gradient between the lower part of the root zone
(below 130 cm) and the "constant head" boundary that drives the upward flux
from the water table to satisfy part of the plant water demand (in this case, roots
between depths of 130 cm and 200 cm are sourcing their water demand fromgroundwater). Note that the fine-textured clay soil has a high air-entry value anda high capacity to maintain the soil close to saturation, thus promoting thisphenomenon.

Figure 3.5. Pressure head profiles at selected print times. The bottom boundary
condition corresponds to a groundwater table with h = 0 cm.

44

Figure 3.6. Root water uptake at selected print times.

Figure 3.7. Flux profiles at selected print times.

45

Figure 3.8. Pressure head profile at Time = 30 days.

Water Flow - Boundary Fluxes and Heads
Referring to Figure 3.9:
• The upper boundary flux is equal to precipitation minus evaporation (2 cm/day– 0.1 cm/day) during the first three days. After day 3 when precipitation ceases,
the flux becomes positive (the upward flux from the atmospheric boundary) due
to evaporation. The actual flux remains equal to the potential flux up to day 12(the first stage of evaporation), after which it drops below the potential value(the second stage of evaporation), until it becomes very small (about 0.01
cm/day at day 30).
• Actual root water uptake remains equal to potential (0.5 cm/day) for a shortperiod of time, after which it starts to decrease. Potential transpiration increases
to 1.5 cm/day after 3 days. Actual transpiration continues to drop until it isapproximately equal to the upward flux from the "constant head" bottom
boundary. At that time, the soil profile attains a steady-state condition.
• The lower boundary flux is initially zero since the initial condition was equal tothe hydrostatic condition. Due to root water uptake drawing water up into the
profile, the bottom flux starts gradually increasing (representing an upward flux
from the water table). After about two days, the precipitation front reaches thebottom boundary and reverses the bottom flux to be negative. Following the
cessation of precipitation, the trend reverses again due to the root water demand.
The flux eventually reaches a steady-state, satisfying about 23% of the root
water demand (the long-term upward flux = 0.366 cm/day, potential
transpiration = 1.5 cm/day).

46

Rainfall

Transpiration
0.5 cm/d

1.5 cm/d

Evaporation
0.1 cm/d

2 cm/dFigure 3.9. The surface flux (into the soil profile), root water uptake (from within thesoil profile), and the bottom flux (out of the soil profile) versus time. The upper partof the figure shows the applied boundary conditions.

47

48

Example 4: Estimating Groundwater Recharge
Example Description and Objectives
The objective of this example is to simulate water flow and root water uptake in a 100cm deep uniform soil profile and to evaluate flow across the bottom boundary of the
profile (i.e., groundwater recharge). The soil profile is initially at a uniform pressure
head equal to -100 cm. The boundary condition for the soil surface is an ‘Atmospheric
Boundary Condition with Surface Runoff’. The bottom boundary condition is ‘Free
Drainage’ (i.e., the groundwater table is assumed to be deep enough to not influencewater redistributions in the soil profile). Roots are uniformly distributed from the soil
surface down to a depth of 50 cm. A 366-day time series of meteorological data fromthe Tully region in north Queensland, Australia, is provided, which allows you to
calculate daily values of potential evapotranspiration. The HYDRUS-1D model is thus
used to calculate actual plant water uptake (i.e., actual transpiration) and to estimate
deep drainage below the root zone, which eventually recharges the underlying aquifer.

Pre-Processing
Project Manager- Select Projects Tab- Click Button “New”- Type Name: Example 4- Type Description: Root uptake with meteorological data
- Click Button “OK”
Main Processes- Heading: Root uptake with meteorological data
- Check “Root Water Uptake” (in addition to a default process "Water Flow")- Click Button “Next”
Geometry Information- Length Units: cm- Number of Soil Materials: 1 (indicates a homogenous soil profile)- Number of Layers for Mass Balances: 1 (users may select up to 10 mass balance
layers)- Decline from vertical axis: 1 (vertical flow)- Depth of Soil Profile (cm): 100 cm- Button “Next”
Time Information- Time Units: Days
Time Discretization:- Initial Time (day): 0- Final Time (day): 366- Initial Time Step (day): 0.001- Minimum Time Step (day): 0.00001- Maximum Time Step (day): 5- Check "Time-Variable Boundary Conditions"49

- Number of Time-Variable Boundary Records: 366- Check "Meteorological Data"- Number of Meteorological Records: 366- Select "Penman-Monteith Equation" (Default)- Button “Next”
Note that contrary to previous examples with Time-Variable Boundary Conditions
when users specify directly the potential transpiration and evaporation rates, in thisexample these rates (or a combined evapotranspiration rate) will be evaluated from
Meteorological Data (e.g., solar radiation, wind speed, air temperature and humidity)
using the Penman-Monteith combination equation recommended by FAO(Monteith, 1981, 1990).
Print Information- Print Options: Leave as Default- Number of Print Times: 12- Button “Select Print Times”- Button “Default”- Button “OK”
- Button “Next”
Water Flow – Iteration Criteria- Leave default values- Button “Next”
Water Flow – Soil Hydraulic Model- Check Single Porosity Model, van Genuchten – Mualem (default)- Button “Next”
Water Flow – Soil Hydraulic Parameters- Soil Loam from the Soil Catalogue (default)- Button “Next”
Water Flow – Boundary Conditions- Upper Boundary Condition: Atmospheric BC with Surface Run Off- Lower Boundary Condition: Free Drainage- Initial Condition: In Pressure Heads- Button “Next”
Root Water and Solute Uptake Model- Water Uptake Reduction Model- Check Feddes (Default)- Button “Next”
Root Water Uptake Parameters- Select Pasture from the Database (default)- Button “Next”

Time-Variable Boundary Conditions
50

Prepare the precipitation time series in a spreadsheet, select the first cell in the‘Time’ column and then paste values from the spreadsheet using the Ctrl+V keyboard
command. The time series of rainfall data is displayed in Figure 4.1.
20

Rainfall (cm/day)

16

12

8

4

0
0

100

200

300

400

Time (days)

Figure 4.1. Time series with daily rainfall values.
Meteorological Parameters
In the Meteorological Parameters dialog window (Fig. 4.2), users specify parameters
needed to calculate potential evapotranspiration using either the Penman-Monteith orHargreaves equations, using geographical and meteorological parameters, andparameters related to crop growth, LAI, and interception. You must then specify
Latitude (deg) and Altitude (m), Angstrom values as (recommended value is 0.25)
and b (recommended value is 0.5), values a1 (recommended value is 0.9) and b1(recommended value is 0.1) for calculating the effect of a cloudiness factor on longwave radiation, values al (recommended value is 0.34) and bl (recommended value is–0.139) for calculating the effect of emissivity on long wave radiation, and heights
of wind speed and temperature and humidity measurements (usually 200 cm) (seeChapter 2.7 for details).
You must also specify whether or not a crop is present. If a crop is present, then thecrop data may be constant in time or must be entered. If the latter, data can be entered
via a table (the data are then interpolated linearly in time between entered values) or
can be provided in terms of daily values. Relevant crop data are Crop Height, Albedo,
LAI, and RootDepth. When “No Crop” is selected, a (surface) Albedo has to beprovided. A default value of Albedo for simulations with a crop is 0.23.
If the Crop Data Daily option is selected in the Meteorological Parameters window,
then the Crop Height and Root Depth need to be specified at the same times as other51

Meteorological variables. If Crop Data in Tables is selected, then Crop Height andRoot Depth should be specified at selected times in the Growth Data table.
Users must also specify how LAI is calculated. LAI can be calculated from:a) From Crop Height, Clipped Gras: the crop height combined with a formulafor grass (i.e., LAI = 0.24×CropHeight),b) From Crop Height, Alfalfa: the crop height combined with a formula foralfalfa (i.e.,
LAI = 1.5×log(CropHeight)+5.5), orc) From Surface Fraction: based on the surface cover fraction, SCF, (i.e., LAI =-1/ai × ln(1-SCF), where ai is the constant for the radiation extinction by the
canopy [-] = 0.463).
- Enter values as shown below in Figure 4.2. Note that the geographicalparameters refer to the Tully region in north Queensland, Australia.
- Button “Next”

Figure 4.2. The Meteorological Parameters dialog window.
Meteorological Conditions- Prepare the time series for solar radiation (Radiation), maximum (T_max) andminimum (T_min) daily air temperatures, relative humidity or vapor pressure(Humidity), and wind speed (Wind) in a separate spreadsheet, select the first cellin the ‘Time’ column of the Meteorological Conditions dialog window (Fig.
4.3) and then paste values from a spreadsheet using the Ctrl+V keyboard
command.
- Button “Next”

52

Figure 4.3. The Meteorological Conditions dialog window.
HYDRUS-1D guide- Do you want to run PROFILE application?
- Button “OK”
Soil Profile – Graphical Editor- Menu: Conditions->Root Distribution or Toolbar- Click the “Edit Condition” button (left side), select with the Mouse the nodes
down to a depth of 50 cm, left click the mouse again, in the Condition
Specification dialog window that appears enter Roots Linear Distribution: Topvalue=1, check ‘Use top value for both’, and click the OK button. This meansthat roots are distributed uniformly in the top 50 cm of the soil profile(according to the Feddes model in Figure 3.3).
- Menu: Conditions->Observation Nodes or Toolbar- Click the “Insert” Button on the left-side Edit Bar and specify nodes at 50 and
100 cm (right-click to end the command)- Menu: File->Save Data or Toolbar- Menu: File->Exit
Soil Profile – Summary- Button “Next”
Execute HYDRUS-1D- Menu: Calculation->Execute HYDRUS-1D or a Toolbar button

Post-Processing
Profile Information
53

.

•

•

Figure 4.4 shows two pressure head profiles at 152.5 (black line) and 335.5(red line) days, which represent conditions after a wet and a dry period,
respectively. Note that the pressure head in the top 50-cm root zone drops toP3 (the wilting point, -8,000 cm for pasture) after the dry period,
Figure 4.5 shows two corresponding flux profiles at 152.5 and 335.5 days.
Note a much higher flux following the wet period and almost a zero flux along
most of the root zone following the dry period.

Figure 4.4. Pressure head profiles at two selected times.

Figure 4.5. Flux profiles at two selected times.

54

Water Flow - Boundary Fluxes and Heads
•
•

•

•

•

Figure 4.6 shows the estimated potential root water uptake rates derived fromthe meteorological data. Note the seasonal signal, which reflects lower uptake
during winter (which in the southern hemisphere is around July).
Figure 4.7 shows actual root water uptake (compare actual uptake with potential
uptake (Fig. 4.6) and the actual surface flux due to infiltration of precipitation(Fig. 4.8)). Note that actual root water uptake (transpiration) is significantly
smaller during periods without precipitation or immediately after large rainfalls.
Figure 4.9 shows the actual bottom flux, which represents deep drainage that
occurs below the root zone. This flux will continue downwards by gravity until
it eventually recharges the underlying aquifer. Note that most of the recharge
corresponds to large rainfall events when the soil profile is saturated (and highly
conductive).
Figure 4.10 shows the cumulative fluxes; the actual surface flux (infiltration
minus evaporation) is negative (downward) and represents infiltration, theactual bottom flux is negative (outflow) and represents recharge, and the actualroot water uptake (transpiration) is positive (upward). Note that the magnitudeof the surface flux is about equal to that of the root water uptake flux and thebottom flux.
If one is interested in the time history of groundwater recharge, i.e., when theaquifer is actually recharged, one would need to model the entire soil profile
down to the water table. In the current example, only the first 100 cm of a deep
unsaturated zone are considered.

Figure 4.6. Potential root water uptake versus time.

55

Figure 4.7 Actual root water uptake versus time.

Figure 4.8. Actual surface flux (into the soil profile) versus time.

56

Figure 4.9. Actual bottom flux (out of the soil profile) versus time.

Figure 4.10. Cumulative actual (surface, bottom, and transpiration) fluxes versus time.

57

58

Example 5a: Basic Solute Transport: Solute Introduced
through a Constant Flux Boundary
Example Description and Objectives
The objective of Example 5a is to provide an introductory example involving solutetransport under variably-saturated flow conditions. Water flow and solute transport
through a 300-cm deep uniform soil profile with an initial pressure head of -100 cm(the default initial value in HYDRUS-1D) is simulated in this example. The soil profile
is initially solute free. Water infiltrates into the soil profile at a constant flux of 5 cm/dduring the entire 20-day simulation. Infiltrating water has a unit solute concentration
during the first 10 days and is solute free during the second 10 days. Water and solute
drain through the bottom ‘Free Drainage' boundary.
This example demonstrates how solute is introduced into the soil profile using a flux
boundary and how solute breaks through at various locations in the soil profile. The
numerical stability criteria and dimensionless Peclet and Courant numbers (Perrochet
and Berod, 1993) are also discussed.

Pre-Processing
Project Manager- Select the Projects Tab- Button “New”- Name: Example 5a- Description: Basic Solute Transport - Solute through a flux boundary- Button “OK”
Main Processes- Heading: Solute through a flux boundary- Check “Solute Transport” (in addition to "Water Flow", which is checked bydefault)- Button “Next”
Geometry Information- Depth of Soil Profile (cm): 300 cm- Button “Next”
Time Information- Time Discretization:- Final Time (day): 20- Button “Next”
Print Information- Print Options: Leave as default- Number of Print Times: 10- Button “Select Print Times”- Button “OK”
- Print Times: Leave default values
59

- Button “Next”
Water Flow – Iteration Criteria- Leave default values- Button “Next”
Water Flow – Soil Hydraulic Model- Leave as default (Single Porosity, van Genuchten – Mualem)- Button “Next”
Water Flow – Soil Hydraulic Parameters- Leave as default (parameters for Loam)- Button “Next”
Water Flow – Boundary Conditions- Upper Boundary Condition: Constant Flux
- Lower Boundary Condition: Free Drainage- Initial Condition: In Pressure Heads- Button “Next”
Water Flow – Constant Boundary Fluxes- Upper Boundary Flux: -5 (cm/d; note that negative fluxes are downwards; theratio of water flux to the saturated hydraulic conductivity is 0.2, i.e., the soil
remains unsaturated)- Button “Next”
Solute Transport - General Information
The next set of four dialog windows deals with solute transport. In the Solute
Transport dialog window (Fig. 5a.1), you can select a time and space weighting
discretization scheme (although using the default schemes, i.e., the Crank-Nicholson
Scheme with Galerkin Finite Elements, is recommended). For minimizing or eliminating
numerical oscillations, HYDRUS-1D implements the combined Peclet and Courant
number criterion introduced by Perrochet and Berod (1993), with a recommended default
stability criterion of 2. In other words, when the combined effect of Peclet and Courant
number is less than or equal to 2 (Pe×Cr ≤ 2), numerical oscillations in the solutetransport solution should be minimal. Further details about the time and spaceweighting scheme are discussed in Appendix I.

60Figure 5a.1. The Solute Transport dialog window.
The next option available is about whether the solute transport and reaction parameters
are water content and/or temperature dependent (in the "Dependence on Environmental
Factors" section). Indeed, several of the diffusion, zero-order production, first-orderdegradation, and adsorption coefficients may be strongly dependent upon temperature(Chapter 3.4). This could also be relevant when the Major Ion Chemistry module (notdiscussed here) is invoked to calculate carbon dioxide transport and production or
multi-component chemistry related to saline and/or sodic soils (see Chapter 5). The
water content dependence of degradation coefficients is implemented using a modified
equation of Walker (1974); if this option is selected, two parameters need to be provided(Figure 5a.): a solute dependent parameter B (usually 0.7), and a reference pressure head(href) which defines the reference water content (θref) at which the reference degradation
constant (ar) is defined. Water content dependent degradation rates have been reported for
denitrification (Morales et al., 2016) and carbon (Jacques et al., 2018; Sierra et al., 2012)
cycling. In the current example, neither the temperature nor the water content dependence
option is activated.

61

Figure 5a.2. The Water Content Dependent Solute Transport and Reaction Parameters
dialog window.
The next options are about invoking the nonequilibrium solute transport model (in the"Nonequilibrium Solute Transport Models" section). Note that the default option is theequilibrium model. Solute transport in that case occurs in the full water-filled porespace (i.e., only a mobile flow region exists), while solute partitioning between theliquid and solid phase is assumed to instantaneous (i.e., not kinetically controlled).
Details of the different nonequilibrium models are described in Chapter 3 and Šimůnek
and van Genuchten (2008).
When solute transport is nonlinear (e.g., when nonlinear sorption is considered), youneed to specify "Iteration Criteria" for the solution of the solute transport equation.
You also have an option to select whether the "Tortuosity Factor" in the diffusion
coefficient for the water and gas phases should be considered and whether theMillington and Quirk (1961) or Moldrup (Moldrup et al., 2000) model should be used
for this (Table 5a.1). If the "Tortuosity Factor" option is not used, the model uses adefault value of one for the tortuosity factor. This means that the effective diffusion
coefficient for the water phase (Dp,w) is equal to the free water diffusion coefficient(Dw):
𝐷𝐷𝑝𝑝,𝑤𝑤 = 𝐷𝐷𝑤𝑤 𝜏𝜏𝑤𝑤 =𝐷𝐷𝑤𝑤 if 𝜏𝜏𝑤𝑤 = 1(5a.1)

while the pore diffusion coefficient for the gas phase (Dp,g) is defined as
𝐷𝐷𝑝𝑝,𝑔𝑔 =𝐷𝐷𝑔𝑔 𝜏𝜏𝑔𝑔 𝑎𝑎𝑣𝑣
𝜃𝜃

=
62

𝐷𝐷𝑔𝑔 𝑘𝑘𝑔𝑔 𝑎𝑎𝑣𝑣
𝜃𝜃

if 𝜏𝜏𝑔𝑔 = 1(5a.2)

where Dg is the diffusion coefficient in the gas phase, av is volumetric air content(cm3/cm3) and kg is the inverse of the Ostwald coefficient (1/kg = KH Ru TA, where KH is
the Henry's Law constant (mol/(atm×m3)), Ru is the universal gas constant (8.31451 m2 kg
s-2 mol-1 K-1), and TA is absolute temperature (degree K). Values for free water diffusion
coefficient (Dw) are listed in Table 5a.2 for a number of common chemicals (about 1.7
cm2/d for Cl). A good reference publication on solute diffusion coefficients is Flury andGimmi (2002). Typical values for the diffusion coefficient in the gas (air) phase (Dg)are 20,000 cm2/day for CO2 (Jacques et al., 2018), 8,640 cm2/day for radon (SφgaardHansen and Damkjaer, 1987), and 16,940 cm2/day for methane (Marrero and Mason,1972).
Table 5a.1. Tortuosity factors for the water (τw) and gas (τg) phases. Other parameters
are the water content (θ), the saturated water content (θs), and the volumetric air content(av).
Tortuosity model
Millington-Quirk (1961)
Moldrup et al. (2000)

Water phaseτw =Gas phase
a7 / 3τg = v2

θ
θ s2

7/3

θ 
τ w = 0.66   θs 

θs

8/3

av1.5
τg =
θs

Table 5a.2. Diffusion coefficients in free water at 18°C, Dw (source: Li and Gregory,1974).
ChemicalDw (cm2/day)Cations
+
H
7.06
Na+
0.98
+
K
1.44
Mg2+
0.51
2+
Ca
0.58
Anions
OH3.88Cl
1.48
HCO31.02
2SO4
0.77
PO430.53

Users can further specify how many solutes will be considered in the advectivedispersive transport equations (one equation per solute will then be solved). Up to 10
solutes can be defined. Finally, you can specify for how long ("Pulse Duration") soluteis applied at (or released from) the soil surface. Use this option only when constantwater flow boundary conditions are used at the soil surface.

63

- Leave everything as default, except- Pulse Duration [day]: 10 (solute will be applied during the first 10 days, while
solute-free solution will be applied after that)- Button “Next”

Solute Transport - Transport Parameters
The next dialog window is used to specify general solute transport parameters. In the
Solute Transport Parameter dialog window (Fig. 5a.3) you can specify suchparameters as the bulk density, the dispersivity (e.g., Vanderborght and Vereecken,2007), and diffusion coefficients for solutes in both the liquid phase (see Table 5.a.2
for examples) and the gas phase if volatile solutes are present. If the “Use Tortuosity
Factor” option is activated, the liquid (“Diffus. W.”) and gas (“Diffus. G.”) diffusion
coefficients should be the free water (Dw) and air phase (Dg) diffusion coefficients.
Alternatively, if the “Use Tortuosity Factor” option is not activated, the tortuosity isequated to one and the diffusion coefficients entered should be the pore diffusion
coefficients for the water (Dp,w) and gas (Dp,g) phases.
The two parameters "Frac" and "ThIm" are used only with nonequilibrium solutetransport models and should be equal to 1 and 0, respectively, for equilibrium solutetransport (i.e., transport without preferential flow and kinetically controlled sorption).
- Leave as default- Note that Disp. = 30 cm, i.e., one-tenth of the travel distance, which is 300 cm.
- Note that the molecular diffusion coefficient in the liquid phase (Diffus.W.) is
set equal to zero. This process is often neglected since in the majority of
applications it is significantly smaller than hydrodynamic dispersion. On theother hand, it is possible to set this value equal to about 1 cm2d-1, which is anapproximate value for a majority of solutes.
- Button “Next”

Figure 5a.3. The Solute Transport Parameters dialog window.

64

Solute Transport - Reaction Parameters
The Solute Transport and Reaction Parameter dialog window (Fig. 5a.4) is used
next to specify parameters for various chemical and biological processes such as
sorption (Kd, Nu, Beta), Henry's law constant (Henry) for volatile substances, andseveral degradation constants (e.g., SinkWater1, SinkSolid1, SinkGas1, SinkWater0,etc). Degradation rates can be in the form of zero- and first-order processes and need to
be specified for each phase (i.e., liquid, solid, and gaseous) separately. The default
values are all zero. When identifying degradation rates, you need to pay close attention
to how the values (e.g., half-lives) are reported. They may be reported for eachindividual phase (again, either liquid, solid, or gaseous), or for the entire (multi-phase)system. For details see Chapter 3 (or press F1).
- Leave as Default- Note that all parameters are equal to zero except for the Freundlich exponent,
which is equal to 1, i.e., parameters are set for a nonreactive tracer- Button “Next”

Figure 5a.4. The Solute Transport and Reaction Parameters dialog window.

Solute Transport - Boundary Conditions
Next, you need to specify appropriate boundary conditions for solute transport to define
how the soil profile interacts with the outside environment (i.e., at the upper and lower
boundaries). This is done using the Solute Transport Boundary Condition dialog
window (Fig. 5a.5).
- Leave as Default(i.e., Upper Boundary Condition: Concentration Flux BC;
Lower Boundary Condition: Zero Concentration Gradient)
For more details on solute transport boundary conditions, refer to Appendix II.2
- Boundary Condition for Solute 1: 1 (this is the concentration of the infiltratingwater, in mmol/cm3)- Button “Next”

65Figure 5a.5. The Solute Transport Boundary Conditions dialog window.
HYDRUS-1D guide- Do you want to run PROFILE application?
- Button “OK”
Soil Profile – Graphical Editor- Menu: Conditions->Observation Nodes or Toolbar- Button (left side) “Insert”: select nodes at 0, 150, and 300 cm- Menu: File->Save Data or Toolbar- Menu: File->Exit
Soil Profile – Summary- Button “Next”
Execute HYDRUS-1D- Menu: Calculation->Execute HYDRUS-1D or a Toolbar button.

Post-Processing
Observation Points
Referring to Figure 5a.6:• The concentration at the surface (ON-1) increases during the first 10 days while
the solute pulse is active at the upper Constant Flux boundary. Solute (1mmol/cm3) was applied during the first 10 days only (see Figure 5.a.2), but theconcentration never reaches unity. Had the pulse continued indefinitely, theentire profile would have attained a unit concentration. However, owing to the
low infiltration rate, it takes a rather long time for the soil in this example toreach a steady-state water content and to replace the initially solute free soil
solution with 1 mmol/cm3 infiltrating solution). The concentration declined as
the solute-free water infiltrates the profile after day 10.

66

•

•
•

The concentration at a depth of 150 cm (ON-2) reached 0.36 mmol/cm3 after 10
days. Even though solute-free water infiltrated after day 10, the concentration
continued to increase until it peaks at 0.58 mmol/cm3 as flowing water continues
to carry solutes downward from the soil above this observation point. The soil
above this point has a higher concentration.
The solute breaks through at the bottom of the profile (ON-3)) after about 9
days.
Had solute-free water continued to infiltrate the profile indefinitely, all solutewould eventually be leached out of the profile, thus returning the soil column to
its initial state.

Figure 5a.6. Solute concentrations versus time at selected observation points.

Profile Information
Referring to Figure 5a.7:• Note the declining concentration peak along the depth of the profile due tosolute dispersion.
• Note that at 10 d, the concentration at a depth of 300 cm (the lower boundary)is just above zero, indicating that solute has just broken through this boundary(Figure 5a.6).
• If the area under each curve is integrated, we should get the total solute massstored within the soil profile at that particular time. This value is reported as‘ConcVol’ in the ‘Mass Balance Information’.

67

Figure 5a.7. Concentration profiles at selected times.
Solute Fluxes
Referring to Figure 5a.8:• Note that solute fluxes are positive when entering the transport domain andnegative when leaving the transport domain. This is different than those forwater flow, which are always positive upwards and negative downwards.
• The cumulative solute flux entering the soil profile through the upper
atmospheric boundary peaks at 50 mmol/cm2. This because a unit concentration
was present in the precipitation, which entered the profile at a rate of 5 cm/day
for 10 days (i.e., 1 mmol/cm3 × 5 cm/day × 10 days =50 mmol/cm2).
• The sum of the influx and efflux curves at any time defines the amount of solute
stored within the soil profile. This value is reported as ‘ConcVol’ in the ‘Mass
Balance Information’ (e.g., 41.27 mmol/cm2at day 20, which is equal to the
cumulative surface flux + the cumulative bottom flux).

68Figure 5a.8. Cumulative solute fluxes (into and out of soil profile) versus time.
Run Time Information•••

The grid Peclet number (Pe) defines the predominant type of the solute transport(the ratio of the advective and dispersive transport terms) in relation to thecoarseness of the finite element grid. The Peclet number increases when theadvective component of the transport problem dominates the dispersive part,which is exhibited by a steep concentration front. To achieve acceptablenumerical results, the spatial discretization must be kept relatively fine to maintaina low Peclet number (smaller than 2). Figure 5a.9 shows a low Peclet number
throughout the simulation, which indicates a stable simulation.
The second dimensionless number, which characterizes the relative extent of
numerical oscillations, is the Courant number (Cu), which is associated with thetime discretization and should be kept below 1 (which the program automaticallydoes).
A stable HYDRUS-1D simulation maintains the product of these twodimensionless numbers (i.e., the Peclet and Courant numbers) to be under 2. For
more details, refer to the HYDRUS-1D Technical Manual (Chapter 8.4.5).

69

Figure 5a.9. The Peclet and Courant numbers (which are related to stability criteria)
versus time.

70

Example 5b: Basic Solute Transport: Solute Introduced
through a Constant Pressure Head Boundary
Example Description and Objectives
This example is the same as Example 5a, except that the constant flux upper boundary
is replaced by the constant pressure head boundary. The objective of Example 5b isthus to simulate water flow and solute transport through a 300-cm deep uniform soil
profile with an initial pressure head of -100 cm. The soil profile is initially solute free.
Water infiltrates into the soil profile through a constant head boundary (0 cm) during
the entire 20-day simulation period. Infiltrating water has a unit solute concentration
during the first 10 days and is solute free during the second 10 days. Water and solute
drains through the bottom ‘Free Drainage' boundary.
This example demonstrates how solute is introduced into the soil profile using a
constant pressure head boundary and how solute breaks through at various locations inthe soil profile. This could be part of a constant-head infiltration experiment wheretracers are introduced together with water, a furrow irrigation system providing waterand dissolved solutes to a cropped soil (Ebrahimian et al., 2013a,b), water flow
processes and the N balance of a rice paddy field (Dash et al., 2015), or a leaking water
holding pond introducing a mixture of chemicals into the underlying soil (Mallants etal., 2017).

Pre-Processing
The example is created most easily by simply copying and renaming Example 5a andimplementing the required modifications.
Project Manager- Select the Projects Tab- Select Example 5a- Button “Copy”- New Name: Example 5b- Description: Basic Solute Transport - Solute through a flux boundary- Button “OK”
- Select Example 5b- Button “Open”
Since the problem is almost identical as Example 5a (except for the water flowboundary condition), we can skip a large number of input dialog windows and start
directly with the Water Flow Boundary Conditions. Note that the Solute Boundary
Condition is identical to the previous example, i.e., a Concentration Flux Boundary.
Water Flow – Boundary Conditions- Upper Boundary Condition: Constant Pressure Head
- Button “OK”
Soil Profile – Graphical Editor- Menu: Conditions->Initial Conditions->Pressure Head or Toolbar
71.

- Button (left side) “Edit Condition”: Select with the Mouse the first node, left
click the mouse, and specify a pressure head of 0 cm.
- Menu: File->Save Data or Toolbar- Menu: File->Exit
Execute HYDRUS-1D- Menu: Calculation->Execute HYDRUS-1D or a Toolbar button.

Post-Processing
Observation Points
Referring to Figure 5b.1:• The concentration at the surface (ON-1) increases during the first 10 days while
the solute pulse is active at the upper Constant Pressure Head boundary.
However, contrary to Example 5a, the surface concentrations do reach unity.
This because under the constant head flow boundary condition the soil becomessaturated quickly and thus the initially solute-free solution near the soil surface
is displaced quickly. Similarly, the concentration declines at a much faster rate
and all solute is washed out at this location when solute-free water infiltrates
into the profile at about day 15. The same applies to (ON-2).
• The solute breaks through at 300 cm depth after only 2 days (compared to 9days in Example 5a). This is due to the much larger water flux that is now
infiltrating into the soil profile through the Constant Pressure Head boundary.
To be exact, the water flux q under the constant pressure head is five times larger
than under the Constant Flux boundary (q=0.2×Ks for the previous and q=Ks forthe current example). Except at very early times (<0.1 day), the actual surface
flux is approximately equal to the saturated hydraulic conductivity of the soil
(about 25 cm/day), which is five times the constant flux applied in Example 5a.
Profile Information
Referring to Figure 5b.2:• Note that after 10 days (the end of the solute pulse), the minimum soluteconcentration across the entire profile was about 95% of the applied
concentration (the profile almost reached a steady-state corresponding to a unit
concentration).
• At day 20, the applied solute was almost completely washed out from the upper
half of the profile, with little solute remaining in the lower half.

72

Figure 5b.1. Solute concentrations versus time at selected observation points.

Figure 5b.2. Concentration profiles at selected times.

Solute Fluxes
Referring to Figure 5b.3:• The cumulative solute mass that infiltrated into the soil profile through the upper
Constant Pressure Head boundary was about 250 mmol/cm2. This follows(approximate) from having a unit concentration, Ks being about 25 cm/day, ahaving a pulse of 10-day duration (i.e., 1 mmol/cm3 × 25 cm/day × 10 days =250 mmol/cm2).
• The sum of the two curves at any time is the solute stored within the soil profileat that time. At the end of the simulation, the solute mass is equal to about 0.9
73mmol/cm2 (refer to ‘ConcVol’ in the ‘Mass Balance Information’ at day 20),
which shows that almost all of the solute introduced to the soil profile via theupper boundary has been leached out through the lower Free Drainage
boundary.

Figure 5b.3. Cumulative solute fluxes into and out of the soil profile versus time.

74

Example 5c: Basic Solute Transport: Profile with Non-Zero
Initial Solute Concentration
Example Description and Objectives
This example is again based on Example 5a, except that solute is now initially present
in the top 30 cm of the soil profile and solute-free water infiltrates the soil. Example 5cthus simulates water flow and solute transport through a 300-cm deep uniform soil
profile with an initial pressure head of -100 cm. The soil profile has solute mass stored
within the top 30 cm (concentration = 30 mmol/cm3). Solute-free water is introduced
through the upper ‘Constant Flux' boundary’ during the entire 30-day simulation. Water
drains through the bottom ‘Free Drainage' boundary.
This example demonstrates how solute is mobilized and washed out from the soil
profile.

Pre-Processing
You can create this example again by copying and renaming Example 5a andmodifying it.
Project Manager- Select the Projects Tab- Select Example 5a- Button “Copy”- New Name: Example 5c- Description: Basic Solute Transport - Solute initially in the profile- Button “OK”
- Select Example 5c- Button “Open”
Main Processes- Heading: Solute initially in the profile- Button “OK”
Time Information
Time Discretization:- Final Time (day): 30- Button “Next”
Print Information- Print Options: Leave as default- Number of Print Times: 4- Button “Select Print Times”- Print Times (day): 2, 5, 15, 30- Button “OK”
- Button “OK”

75

Solute Transport - General Information- Leave as default, except- Pulse Duration [day]: 30 (a solute-free solution is applied at the boundary)- Button “Next”
Solute Transport - Boundary Conditions- Leave as Default(i.e., Upper Boundary Condition: Concentration Flux BC;
Lower Boundary Condition: Zero Concentration Gradient)- Boundary Condition for Solute 1: 0 (a solute-free solution is applied at theboundary)- Button “Next”
Soil Profile – Graphical Editor- Menu: Conditions->Initial Conditions->Concentration or Toolbar
andselect Concentration 1 from the Quantity listbox on the Edit bar- Button (left side) “Edit Condition”: select with the Mouse the upper 30 cm, left
click the mouse, and specify the initial solute concentration of 30 (mmol/cm3).
- Menu: File->Save Data or Toolbar- Menu: File->Exit
Soil Profile – Summary- Button “Next”
Execute HYDRUS-1D- Menu: Calculation->Execute HYDRUS-1D or a Toolbar button.

Post-Processing
Observation Points
Referring to Figure 5c.1:• The concentration at the surface (ON-1) rapidly drops as solute-free water
infiltrates into the soil profile.
• The initial solute bulge (from 0-30 cm) rapidly disperses and travels down the
profile, resulting in a maximum concentration of 3.13 mmol/cm3 at a depth of
150 cm (ON-2) after 9.6 days.
• Further dispersion results in a maximum concentration of only 2.2 mmol/cm3
after 19.7 days at the lower boundary (ON-3).

76

Figure 5c.1. Solute concentrations versus time at selected observation points.

Profile Information
Figure 5c.2 shows how the initial solute mass travels down the profile with the peakconcentration decreasing down the profile due to dispersion.

Figure 5c.2. Concentration profiles at selected times.

77Solute Fluxes
Figure 5c.3 shows the bottom solute flux, which peaks just after about 20 days, whichalmost coincides with the time when the solute concentration peaks at the same location(see Figure 5c. 5c.1, ON-3 at 300 cm depth). Recall that solute fluxes are positive whenentering the transport domain and negative when leaving the transport domain.

Figure 5c.3. Actual and cumulative bottom solute fluxes versus time.

78

Example 5d: Basic Solute Transport: The Effects of Solute
Dispersion
Example Description and Objectives
The objective of Example 5d is to demonstrate the effect of solute dispersion (throughthe use of different values of the dispersivity) on solute transport. By reducing thedispersivity, advective transport becomes more dominant, leading to sharper
concentration fronts. Example 5d is identical to Example 5c except that it uses a muchsmaller dispersivity parameter (DL = 2 cm) compared to the previous example (DL = 30cm).
Example 5d thus simulates water flow and solute transport through a 300-cm deep
uniform soil profile with an initial pressure head of -100 cm. The soil profile has asolute mass stored within the top 30 cm (concentration = 30 mmol/cm3). Solute-freewater is introduced through the upper ‘Constant Flux' boundary’ during the entire 30day simulation. Water drains through the bottom ‘Free Drainage' boundary.

Pre-Processing
This example will be created by copying and renaming Example 5c and thenmodifying it.
Project Manager- Select the Projects Tab- Select Example 5c- Button “Copy”- New Name: Example 5d- Description: Basic Solute Transport - Showing the effects of dispersivity- Button “OK”
- Select Example 5d- Button “Open”
Main Processes- Heading: Show the effect of dispersivity- Button “OK”
Since the problem is identical to Example 5c (except for the dispersivity value), we canskip all input dialog windows except for the Solute Transport Parameters dialog
window.
Solute Transport - Transport Parameters- Disp. = 2 (cm)- Button “OK”
Execute HYDRUS-1D- Menu: Calculation->Execute HYDRUS-1D, or a Toolbar button.
- The dialog window appears and asks "Do you want to save the input data
before executing HYDRUS-1D?".
79

- Button “Yes” (the input data need to be saved before the computational module
can be run).
- The dialog window appears and warns "All existing output files will bedeleted. Continue?". When output files exist, the program issues a warning thatthey will be deleted if the computational module is run.
- Button “Yes”

Post-Processing
Observation Points
Compare Figure 5d.1 (lower dispersion) to Figure 5c.1 (higher dispersion):• The peak concentrations are higher.
• The peak concentrations appear at a later time.
• Solute advances in a more advective manner (i.e., closer to piston flow).

Figure 5d.1. Solute concentration versus time at selected observation points

Profile Information
Compare Figure 5d.2 (lower dispersion) to Figure 5c.2 (higher dispersion):• Solute moves as a narrower bulge, while having a higher peak concentration.
• Solute breaks through at the lower boundary at a later time.

80Figure 5d.2. Profile concentrations at selected times
Solute Fluxes
Compare Figure 5d.3 (lower dispersion) to Figure 5c.3 (higher dispersion):• The solute flux peak is higher.
• Solute breaks through the lower boundary at a later time.

Figure 5d.3. Bottom solute fluxes versus time81

Run Time Information
Referring to Figure 5d.4:• Note a much higher Peclet number in Example 5d (less dispersion) comparedto Example 5c (high dispersion). Higher Peclet numbers indicate a higher ratioof advective transport to dispersive transport.
• The Peclet number reached a value of 2 during the simulation. We recommendthat when Peclet numbers close to two or higher are obtained, the mass balanceshould be checked (in this case, the relative error in the solute mass balance isonly 0.15%). Note that for much lower dispersivities (and correspondingly
higher Peclet numbers), the current time and space weighting schemes may not
be appropriate and may lead to unacceptable mass balance errors. When thishappens you may need to implement a finer spatial discretization or use various
stabilizing options (e.g., using upstream weighting or a stability criterion). We
refer to the Chapter 8.4.5 of the HYDRUS-1D Technical Manual for details.

Figure 5d.4. Stability Criteria versus time82

Example 5e: Basic Solute Transport: The Effects ofDiffusion
Example Description and Objectives
The objective of Example 5e is to demonstrate the effect of diffusion on solutetransport. Example 5e simulates steady-state water flow and solute transport through avertical saturated 100-cm uniform soil profile. The pressure head is equal to zero
everywhere in the soil column. The same pressure head is also imposed on the two'Constant Pressure Head' boundaries (i.e., the upper and lower boundaries). As a result,flow occurs only due to gravity and the fluid flux is equal to the saturated hydraulicconductivity. Water having a unit solute concentration is introduced through the upper'Constant Pressure Head' during the entire 1,000-day simulation. Because the soil
profile is fully saturated, the soil hydraulic parameters (except for the saturated watercontent and the saturated hydraulic conductivity) are irrelevant. The saturated hydraulicconductivity Ks has a very low value of 0.001 cm/d, under which conditions thecontribution from molecular diffusion to solute transport becomes important. Twosimulations are carried out using diffusion coefficients of 0.10 cm2/day (marked as Low
Diffusion) and 2.0 cm2/day (marked as Typical Diffusion), the latter being a good ruleof-thumb value for the molecular diffusion coefficient in water. Note that we set thedispersivity equal to zero and disabled the tortuosity factor so that the soil diffusion
coefficient is equal to the specified molecular diffusion coefficient. This example
demonstrates the critical role of solute diffusion when flow rates are extremely low.

Pre-Processing
Project Manager- Select Projects Tab- Button “New”- Name: Example 5e1- Description: Basic Solute Transport - Typical Diffusion- Button “OK”
Main Processes- Heading: Show the effect of diffusion- Uncheck "Water Flow"(When the "Water Flow" option is unchecked, the initial pressure head
distribution is kept constant during the simulation. For example, when thepressure head is kept constant and equal to zero, the gravitational flow will
result in a steady-state water flux equal to the saturated hydraulic conductivity.)- Check “Solute Transport”- Button “Next”
Geometry Information- Keep default values- Depth of Soil Profile (cm): 100 cm- Button “Next”
Time Information
83

Time Discretization:- Final Time (day): 1000
- Button “Next”
Print Information- Keep default values- Number of Print Times: 5- Button “Select Print Times”- Button “OK”
- Button “Next”
Water Flow – Iteration Criteria- Keep default values- Click Button “Next”
Water Flow – Soil Hydraulic Model- Keep default values (Single Porosity Models, van Genuchten – Mualem)- Click Button “Next”
Water Flow – Soil Hydraulic Parameters- Click Button “Next”- Qr=0, Qs=0.4, Alpha=0.01 n=1., Ks=0.001 cm/d, l=0.5(note that under saturated conditions, only the saturated water content, Qs, andthe saturated hydraulic conductivity, Ks, are relevant. The values of otherparameters are not important.)- Button “Next”
Water Flow – Boundary Conditions- Upper Boundary Condition: Constant Pressure Head
- Lower Boundary Condition: Constant Pressure Head
- Initial Condition: In Pressure Heads(This setup ensures that the pressure head throughout the soil column and during
the entire simulation will remain constant and equal to zero. As a result, only thegravitational flow will occur.)- Button “Next”
Solute Transport- Keep default values- Uncheck “Use Tortuosity Factor” (to keep the soil diffusion coefficient equal tothe molecular diffusion coefficient in water)- Pulse duration (day): 1000
- Button “Next”
Solute Transport Parameters- Disp. = 2 (cm)
- Diffus. W. = 2 (cm2/day) [or 2.3×10-9 m2/s, a Typical Diffusion transport
simulation]- Button “Next”

84Solute Transport and Reaction Parameters- Keep default values- Button “Next”
Solute Transport Boundary Conditions- Sol. No. 1 Bound. Cond. = 1- Button “Next”
HYDRUS-1D guide- Do you want to run PROFILE application?
- Button “Yes”- Do you want to save data before running Graphical Profile Editor?
- Button “Yes”- All existing output files will be deleted. Continue?
- Button “Yes”
Soil Profile – Graphical Editor- Menu: Conditions->Initial Conditions->Pressure Head or Toolbar- Button (left side) “Edit Condition”: select with the Mouse the entire domain, left
click, specify 'Pressure Head Distribution': Top value = 0, Bottom value = 0, andunchecked "Use top value for both".
- Button “OK”
- Menu: Conditions->Observation Nodes or Toolbar- Button (left side) “Insert”: select nodes at 50 and 100 cm (right-click to end thecommand)- Menu: File->Save Data or Toolbar- Menu: File->Exit
Soil Profile – Summary- Button “Next”
Execute HYDRUS-1D- Menu: Calculation->Execute or a Toolbar button.

Rerun the example with lower diffusion coefficient:
Project Manager- Select Projects Tab- Select Example 5e1- Copy- New Name: Example 5a2- Description: Basic Solute Transport - Low diffusion- Select Example 5e2- Button “OK”
Solute Transport Parameters- Diffus. W. = 0.1 (cm2/day) [or 1.16×10-10 m2/s, [a Low Diffusion simulation]- Button “Next”

85Execute HYDRUS-1D- Menu: Calculation->Execute HYDRUS-1D or a Toolbar button.

Post-Processing
Observation Points
Referring to Figure 5e.1:• The extremely low flow rate (0.001 cm/d) significantly hinders solute transport
dominated by the diffusion process. The low flow rate in this example is due tothe very low Ks and having a zero pressure head gradient.
• For the case of low diffusion (Diffusion coefficient = 0.1 cm2/day),breakthrough curves at selected observation nodes are not shown because thesolute front had not reached the selected nodes within the simulation period.
• Obviously, having a high diffusion coefficient significantly enhances solutetransport since diffusion is then the dominant process at such low flow rates.
• Increasing the dispersivity has a very marginal impact on solute breakthroughat such low flow rates. Conversely, when flow rates are high and solute transport
is dominated by dispersion processes [e.g., Example 5d], changing the diffusion
coefficient hardly has any impact on solute breakthrough.

Figure 5e.1. Solute concentrations versus time at selected observation points.

Profile Information
Referring to Figure 5e.2:• For the case of low diffusion (red lines), the solute concentration profiles showthat the solute front had not yet reached a depth of 50 cm (the location of ON1).
• A higher diffusion coefficient significantly enhanced solute breakthrough,albeit still during very long timeframes.

86Figure 5e.2. Concentration profiles at selected times for low (red lines) and high(black lines) diffusion.
Solute Fluxes
Referring to Figure 5e.3:• The solute mass introduced at the upper constant pressure head boundary is
extremely low due to the very low flow rate.
• A higher diffusion rate leads to more solute mass entering the constant pressure
head boundary where the solute is introduced.
• The cumulative bottom flux is not shown here due to its negligibly small value(only 0.038 mmol/cm2 after 1000 days for the high diffusion case).

Figure 5e.3. Surface solute fluxes versus time for low (a red line) and high (a blackline) diffusion.
87

88

Example 6a: Advanced Solute Transport: Equilibrium(instantaneous) Linear Sorption
Example Description and Objectives
Example 6a simulates flow and transport through an initially solute-free 300-cm deep
uniform, fully saturated soil profile (having an initial pressure head of 0 cm). Water
drains through the bottom ‘Free Drainage' boundary. Water having a unit soluteconcentration is introduced through the upper ‘Constant Pressure Head' boundary’during the entire 20-day simulation. Solute sorption onto the soil solid phase isintroduced by considering linear equilibrium sorption, which is typically represented
mathematically by the distribution coefficient (or the solid-liquid partitioning
coefficient), Kd (Figure 6.a.1). Applying the Kd concept results in solutes being
instantaneously sorbed to the solid phase, and in solute retardation that is described by
the retardation factor R defined as:
R= 1+

ρ Kd
θ

(6a.1)

where Kd (L/kg = cm3/g) is an empirical distribution factor of the linear sorption
isotherm that relates S (the sorbed concentration, i.e., solutes sorbed to soil particles, gsolute/g soil) to C (the concentration in the liquid pore space, g solute/cm3), θ is thewater content (equal to porosity at full saturation, cm3/cm3) and ρ is the bulk density(g/cm3).
S (sorbed concentration)&LangmuirC (dissolved concentration)

Figure 6a.1. General sorption isotherm representation, showing differences between
linear (indicated as ‘Kd’ region) and non-linear sorption (indicated as ‘Freundlich’
region), a transition period left of the precipitation threshold, and solute precipitation(indicated as ‘precipitation’ region). The vertical axis (S) and the horizontal axis (C)are, respectively, the sorbed concentration (on solid phases) and the dissolved
concentration (in the liquid phase). (Source: Wang et al., 2009).

89Substituting Kd = 0.267 (cm3/g), θ =0.4 (cm3/cm3), and ρ = 1.5 (g/cm3) results in an R
value of 2.0 for the present example. The retardation factor (R, dimensionless)represents the relative rate of fluid flow to the transport rate of a reactive sorbing
solute (Freeze and Cherry, 1979) or:

R=

v
vr

(6a.2)

where v is fluid flow rate, νr is the transport rate for the center of mass of the reactive
solute. The retardation factor is equal to 1 for non-reactive (non-sorbing) solutes
travelling in a medium where the entire pore space is effective in transporting the fluid(effective porosity=total porosity).
Example 6a demonstrates the effect of retardation when a solute is sorbed onto the soil
solid phase. This process is modeled here using a simple linear model (termed a linear
isotherm). The effect of retardation is demonstrated by comparing the results from twonearly identical simulations: one having R = 1 (no sorption) and the other having R = 2(with sorption). This project will show the effect of retardation on solute breakthroughfrom the soil profile (or from a laboratory column). Since HYDRUS-1D does not reportsorbed solute masses for the equilibrium model, the mass balance information will be
used to demonstrate how this term can be calculated across the entire profile.

Pre-Processing
We create 'Example 6a1' by copying and modifying 'Example 5b' using the ProjectManager.
Project Manager- Select Projects Tab- Select Example 5b- Click Button “Copy”- Type New Name: Example 6a1- Type Description: Advance Solute Transport - Retardation- Click Button “OK”
- Select Example 6a1- Click Button “Open”
Main Processes- Heading: Simulate the Effects of Retardation- Uncheck "Water Flow" (while in Example 5b transient water flow wasconsidered, in Example 6a1 water flow is steady-state)- Click Button “OK”
Water Flow – Soil Hydraulic Parameters- Soil Catalogue: Loam
- Qs = 0.4 [saturated water content; the value used to get round numbers in themass balance]- Ks = 20 [saturated hydraulic conductivity; the value used to get round numbersin the mass balance]- Button “OK”
90Solute Transport - General Information- Leave default values- Pulse duration (day): 20 [Located in the bottom right corner; to obtain acontinuous solute application throughout the entire simulation period of 20
days]- Button “Next”
Solute Transport - Transport Parameters- Bulk. D. = 1.5 (g/cm3)- Disp. = 2 (cm)- Button “Next”
Solute Transport - Reaction Parameters- Kd = 0.267 (cm3/g, this will invoke retardation)- Button “OK”
Soil Profile – Graphical Editor- Menu: Conditions->Observation Nodes or Toolbar- Button (left side) “Delete”: select the top node at 0 cm to delete this observation
node (right-click to end the command)- Menu: Conditions->Initial Conditions->Pressure Head or Toolbar- Button (left side) “Edit Condition”: select with the Mouse the entire profile andspecify the Pressure Head of 0 cm.
- Menu: File->Save Data or Toolbar- Menu: File->Exit
Soil Profile – Summary- Button “Next”
Execute HYDRUS-1D- Menu: Calculation->Execute HYDRUS-1D or a Toolbar button
Rerun the example without retardation:
Project Manager- Select Projects Tab- Select Example 6a1- Copy- New Name: Example 6a2- Description: Advance Solute Transport - No Retardation- Select Example 6a2- Button “OK”
Solute Transport - Reaction Parameters- Kd = 0 (without retardation)- Button “OK”

91.

Execute HYDRUS-1D- Menu: Calculation->Execute HYDRUS-1D or a Toolbar button.

Post-Processing
Observation Points• Figures 6a.2 and 6a.3 show that retardation delays solute breakthrough by afactor of two (R = 1+0.267×1.5/0.4 = 2) since a fraction of the solute is sorbed
to the soil particles.
• Figures 6a.2 and 6a.3 show that a steady-state concentration of unity (equal tothe concentration applied at the inlet constant pressure head boundary) has
been achieved in the profile.

Figure 6a.2. Liquid phase solute concentrations versus time with (R=2) and without(R=1) sorption at ON-1.

Figure 6a.3. Liquid phase solute concentrations versus time with (R=2) and without(R=1) sorption at ON-2.
•

Figure 6a.4 displays the concept of dimensionless time (Skaggs and Leij, 2002).
The use of dimensionless time is useful since it eliminates the effects of travel92

distance and flow velocity, which are specific to any one experiment.
Dimensionless time as used here is also known as the number of pore volumes
and represents the cumulative outflow volume from a soil profile or laboratory
column, divided by the volume of fluid in the system. Travel time t is the ratioof distance to average pore-water velocity (L/ν). To obtain dimensionless time,T, we multiply time t by this ratio (ν/L) to get:
T =

vt
L

(6a.3)

Multiplying Eq. 6a.3 by A/A, where A is the cross-sectional area of the systembeing considered, gives

T=

vtA QA Veluted= =LA LA Vavailable(6a.4)

where T (number of pore volumes) represents how many times the fluid content
of the column has been displaced by the influent fluid. In this example, T is
calculated as follows:
T=

20t
300 * 0.40

(6a.5)

where the flux is 20 cm/day, the column length is 300 cm, and the water content
is 0.40. Figure 6a.4 shows that the center of solute mass breaks through after
about 1 pore volume (when nonequilibrium processes are present solute does
not break through at one pore volume), and that sorption delays it by anotherpore volume. In fact, one can show that the number of displaced pore volumes
is mathematically equal to R (in reality the area above the curves in Fig. 6a.4).

Figure 6a.4. Liquid phase solute concentrations versus dimensionless time with (R=2)
and without (R=1) retardation at ON-2.
Profile Information
Figure 6a.5 shows solute profile concentrations after 4 and 8 days. Note that identicalprofiles are obtained when doubling the residence time and doubling the retardation
factor.

93

Figure 6a.5. Liquid phase solute concentration profiles with (R=2) and without (R=1)
retardation at different times.
Solute Fluxes•

Figure 6a.6 shows cumulative solute fluxes at the surface (inlet) and bottom(outlet) boundaries. When solute retardation is active (R=2), solutes become
sorbed onto the solid phase and hence less solute mass (since some solute mass
is stored at the sorption sites) is observed at the bottom boundary compared tothe case when retardation is absent. Note that the breakthrough at the end of thecolumn starts at about 5 and 10 days for R equal to 1 and 2, respectively,reflecting the retardation factors ratio of 2.

Figure 6a.6. Cumulative surface and bottom solute fluxes versus time with (R=2) andwithout (R=1) retardation.
94

Mass Balance Information
HYDRUS-1D does not explicitly report sorbed concentrations (S) since they can beeasily calculated from the adsorption isotherm and the liquid phase concentrations (C),i.e., S=Kd×C. The information reported in the mass balance file can be used to derive
the sorbed mass. Note that the calculations below use data from the entire profile (since
only one subregion was used) and hence produce correct results only when a steadystate concentration has been reached [in this case, cMean = 1 mmol/cm3].
The mass balance information refers to the following:
Length:W-volume:
In-flow:
subregion)h Mean:ConcVol:cMean:Top Flux:
Bot Flux:
WatBalT:
WatBalR:

Depth of the soil profile
Total volume of soil water
The sum of fluxes flowing into and out of the soil profile (orThe average pressure head in the soil profile
The total solute mass in the profile (in the liquid phase and sorbed
equilibrium phase)
The mean concentration in the liquid phase
The flux into the soil profile
The flux out of the soil profile
The total water balance error in the soil profile
The relative water balance error in %Table 6a.1 compares mass balance information at the end of the simulation (day 20).
Note that the total combined solute mass in the liquid and sorbed phases (ConcVol) is
240 mmol/cm2 when R=2 compared to 120 mmol/cm2 when R=1.

Table 6a.1. Mass balance information for Time = 20 days with (R=2) andwithout (R=1) retardation.

LengthW-volume
In-flow
H MeanConcVolcMean
Top Flux
Bot Flux
WatBalT
WatBalR
CncBalT
CncBalR[L][L][cm/day][L][mmol/cm2][mmol/cm3][cm/day][cm/day][cm][%][mmol]
[%]

95

R=2

R=1

300
120
0
0
240
1.00

300
120
0
0
120
1.00
-20
-20
0
0

1.20
0.214

0.6
0.088

By definition, Kd is an empirical distribution factor of the linear sorption isotherm that
relates S (sorbed concentration) to C (resident concentration):

Kd =

S
C

(6a.7)

We substitute Equation 6a.7 into Equation 6a.1 and obtain:

R= 1+

Sρ
Cθ

(6a.8)

where Cθ is the mass of solute in the liquid phase and Sρ is the mass of sorbed solute
per unit soil volume.
Referring to Table 6a.1, for the case when R=2:
ConcVol= 240 mmol/cm2
W-volume= 120 cm3
cMean
= 1 mmol/cm3 (this calculation works well here as the profile
has reached a steady state)
Mass in the liquid phase = Cθ = 1 x 120 = 120 mmolSρ [sorbed mass] = ConcVol [total mass] – Cθ [mass in liquid phase] = 240 –
120 = 120 mmol
120
Sρ
=
=
1+
1+
2 [Implemented in this example]
R=
Cθ
120

96

Example 6b: Advanced Solute Transport: Non-Linear
Equilibrium Non-Linear Sorption
Example Description and Objectives
Example 6b simulates water flow and solute transport through a 100-cm horizontal
uniform saturated soil profile having a 10% pressure head gradient imposed across two
Constant Pressure Head Boundaries. Water having a 10 mmol/cm3 solute concentration
is introduced to the initially solute-free soil profile through the upper ‘Constant Pressure
Head’ boundary during the entire 50-day simulation. Sorption is modeled using a nonlinear Freundlich-Langmuir isotherm given by

Kd C βS=
1 + ηC β(6b.1)

where S is the sorbed concentration (solutes sorbed onto the solid phase), C is the liquid
phase resident concentration, and Kd, β and η are empirical or quasi-empiricalconstants. The units for Kd are (L3β M−β) and those for η are (L3M−1). When β=1, Eq.
(6b.1) becomes a Langmuir isotherm, when η=0, the equation becomes a Freundlich
isotherm, and when both β=1 and η=0 the equation reduces to a linear sorption isotherm(Eq. 6a.7). The Freundlich and Langmuir isotherms are shown in Figure 6b.1 fordifferent values of the parameters β and η. For nonzero (positive) values, the Langmuircoefficient η leads to an upper limit (given by Kd/η) that sorbed concentrations cannot
exceed, while the Freundlich exponent β dictates how quickly this saturation limit isachieved. We use this example to demonstrate how solute transport is affected by withnonlinear sorption.

Figure 6b.1. Plots of the Freundlich adsorption isotherm given by Eq. (6b.1, η= 0), with
Kd = 1 and β defined in the caption (left), and the Langmuir adsorption isotherm given
by Eq. (6b.1, β=1), with Kd = 1 and η defined in the caption (right) (source: Mallants etal., 2011)97

Pre-Processing
Project Manager- Select Projects Tab- Button “New”- Name: Example 6b- Description: Advance Solute Transport - Equilibrium nonlinear sorption- Button “OK”
- Select Example 6b- Button “Open”
Main Processes- Heading: Solute transport with non-linear sorption- Uncheck "Water Flow"(When the "Water Flow" option is unchecked, the initial pressure head
distribution is kept constant during the simulation. For example, when thepressure head is kept constant and equal to zero, gravitational flow in the soil
profile will result in a steady-state flux equal to the saturated hydraulicconductivity.)- Check “Solute Transport”- Button “Next”
Geometry Information- Decline from Vertical Axis: 0 (horizontal flow and transport)- Depth of Soil Profile (cm): 100 cm- Button “Next”

Time Information
Time Discretization:- Final Time (day): 50- Button “Next”
Print Information- Print Options: Leave as Default- Number of Print Times: 12- Button “Select Print Times”- Button “Default”- Button “OK”
- Button “Next”
Iteration Criteria- Leave as Default- Button “Next”
Water Flow – Soil Hydraulic Model- Leave as Default- Button “Next”

98Water Flow – Soil Hydraulic Parameters- Soil Catalogue: Loam
- Qs = 0.4 (to get round numbers in the mass balance)- Ks = 20 (to get round numbers in the mass balance)- Button “Next”
Water Flow – Boundary Conditions- Upper Boundary Condition: Constant Pressure Head
- Lower Boundary Condition: Constant Pressure Head
- Button “Next”
Solute Transport- Leave as Default- Pulse duration (day) 50
Iteration Criteria: (Need to be specified for nonlinear solute transport problems.
Nonlinear sorption causes the problem to be nonlinear.)- Absolute Concentration Tolerance = 0.001- Relative Concentration Tolerance = 0.001- Maximum Number of Iterations = 20
- Button “Next”
Solute Transport Parameters- Bulk. D. = 1.5 (g/cm3)- Disp. = 10 (cm)- Button “Next”
Solute Transport and Reaction Parameters – Solute 1- Kd = 0.267 (cm3×2/g2), because Beta = 2)
- Nu = 0.1 (cm3/g)- Beta = 2.0 (-)- Button “Next”
Solute Transport Boundary Conditions- Sol. No. 1 Bound. Cond. = 10 mmol/cm3- Button “Next”
Soil Profile – Graphical Editor- Menu: Conditions->Initial Conditions->Pressure Head or Toolbar- Button (left side) “Edit Condition”: select with the Mouse the entire domain, left
click
- Pressure Head Distribution: Top value = 10 (this in fact a left-side value for thehorizontal soil profile), Bottom value = 0 (a right-side value), uncheck "Use top
value for both"- Button “OK”
- Menu: File->Save Data or Toolbar- Menu: File->Exit
Soil Profile – Summary- Button “Next”
99Execute HYDRUS-1D- Menu: Calculation->Execute HYDRUS-1D or a Toolbar button.

Post-Processing
Observation Points
Figure 6b.2 shows breakthrough curves of the liquid and sorbed concentrations at
depths of 50 and 150 cm in the soil profile. Solute concentrations in the liquid phaseapproach the applied concentration of 10 mmol/cm3 but do not quite reach it. Note thatsorbed concentrations (which are not reported by HYDRUS-1D) were calculated
directly from the liquid phase concentrations using the sorption isotherm (eq. 6b.1).

Figure 6b.2 Solute concentrations in the liquid (solid lines) and sorbed (dashed lines)
phases versus time at depths of 50 and 100 cm.
Profile Information
Figure 6b.3 shows calculated solute concentration profiles at different times. They
confirm that concentrations at the end of the simulation (the blue line) have not yet
reached the steady state concentration of 10 mmol/cm3.

100

Figure 6b.3. Concentration profiles at different times.

101

102

Example 6c: Advanced Solute Transport: One-Site
Chemical Nonequilibrium Sorption
Example Description and Objectives
Example 6c includes three simulations of water flow and solute transport through a
100-cm uniform saturated horizontal soil profile. Water having various soluteconcentrations is introduced into the initially solute-free profile through the upper (leftside) ‘Constant Pressure Head' boundary during the entire simulation period. Chemical
nonequilibrium is modeled using the one-site sorption model (Figure 6c.1), in which afirst-order rate coefficient (Alpha) controls time-dependent (not instantaneous) solute
mass transfer between the liquid and solid phases. For example, the Alpha rate constantof 0.05 d-1 adopted in the three simulations below means that solutes in the liquid phase
are transferred to the sorption sites at a rate of approximately 5% per day. The
Freundlich-Langmuir isotherm is used similarly as in Example 6b, indicating that there
is an upper value for the sorbed phase concentration when an equilibrium state isreached. The one-site chemical nonequilibrium sorption equation is then given by:

 K CβdS= α d β −Sdt 1 + ηC


(6c.1)

where (Alpha in HYDRUS-1D) is the rate coefficient.

Figure 6c.1. Conceptual representation of chemical nonequilibrium one-site sorption.
The parameter Sk is the kinetically sorbed solute concentration on the solid phase.
In the first two simulations ( i.e., 6c1 and 6c2), pressure head gradients of 10% and 5%
are imposed across the two Constant Pressure Head boundaries to demonstrate theeffect of residence time on non-equilibrium sorption. Simulation times are 50 and 100
days in the two simulations, so that the total inlet solute mass is equal in both cases.
The third simulation (6c3) uses a much higher inlet concentration of 1,000 mmol/cm3(compared to 10 mmol/cm3 in the other two simulations), a longer simulation time of
200 days, and a 10% pressure head gradient to ensure that the saturation sorbed
concentration is achieved.
This example demonstrates that time can be an important factor affecting solute masstransfer between the liquid phase and solid phases for the solute transport under103

chemical nonequilibrium conditions. Since mass transfer is time-dependent, theresidence time plays a crucial role since at lower flow rates the mass transfer between
the liquid and solid phases becomes more effective. Since solutes are introduced into
the soil profile at different rates, we will use dimensionless time to enable comparisons
between simulations that have different inlet fluxes.

Pre-Processing
Project Manager- Select Projects Tab- Select Example 6b- Button “Copy”- New Name: Example 6c- Description: Advance Solute Transport - Chemical non-equilibrium- Button “OK”
- Select Example6b- Button “Open”
Main Processes- Heading: Solute transport with chemical non-equilibrium- Button “Next”
Time Information
Time Discretization:- Final Time (day): 50 (Simulation 6c1), 100 (Simulation 6c2),
Final Time (day): 200 (Simulation 6c3)
- Button “OK”
Solute Transport
Non-equilibrium Solute Transport Models- Check One-site sorption model (chemical non-equilibrium)- Pulse duration (day): 50 (Simulation 6c1), 100 (Simulation 6c2), 200(Simulation 6c3)
- Button “Next”
Solute Transport Parameters- Frac. = 0 (all sorption sites are kinetic – as in Figure 6c.1; no equilibrium
sorption)- Button “Next”
Solute Transport and Reaction Parameters – Solute 1- Kd = 0.267 (cm3×2/g2) (keep the same sorption parameters)- Nu = 0.10 (cm3/g)- Beta = 2.0 (-)- Alpha = 0.05/day (need to scroll to the right of the table to find this lastparameter)- Button “Next”
Solute Transport Boundary Conditions- Sol. No. 1 Bound. Cond. = 10 mmol/cm3 (Simulations 6c1 and 6c2)
104

- Sol. No. 1 Bound. Cond. = 1000 mmol/cm3 (Simulation 6c1)
- Button “Next”
Soil Profile – Graphical Editor- Menu: Conditions->Initial Conditions->Pressure Head
- Button (left side) “Edit Condition”: select with the Mouse the entire domain,then left click- Pressure Head Distribution: Top value = 10 (for Simulations 6c1 and 6c3)
- Pressure Head Distribution: Top value = 5 (for Simulation 6c2)
Bottom value = 0; uncheck "Use top value for both" (for all three Simulations)- Button “OK”
- Menu: File->Save Data
- Menu: File->Exit
Soil Profile – Summary- Button “Next”
Execute HYDRUS-1D- Menu: Calculation->Execute HYDRUS-1D or a Toolbar button.

Post-Processing
Solute Fluxes• Dimensionless time is calculated according to Equation 6a.3 (T=vt/L). For
simulations 6c.1 and 6c.3, the pore-water velocity (v) was equal to 5 cm/day(the pressure head gradient (0.10) times Ks (20) divided by the water content(0.4), while for simulation 6c.2, v was equal to 0.05 * 20/0.4 = 2.5 cm/day; L
was 100 cm for all simulations.
• Simulation 6c1 with the higher flux of 2.0 cm/day is marked as ‘Low Residence
Time’, whereas Simulation 6c2 with the lower flux of 1.0 cm/day is marked as‘High Residence Time’.
• Referring to Figures 6c.2: For a longer residence time, more time is available
during which the time-dependent non-equilibrium solute mass transfer process
occurs, resulting in a more cumulative transfer of solutes to the sorption sites.
• Figure 6c.3 shows that for the example with a long residence time, more solute
mass is transferred to the sorption sites, which means that less solute mass
remains in the liquid phase and hence less cumulative solute mass exits thebottom boundary at a given time compared to the simulation with the shorter
residence time.

105

Figure 6c.2 Effect of the residence time on the non-equilibrium solute
mass transfer.

Figure 6c.3. Effect of the residence time on the cumulative bottom solute flux.
Mass Balance Information
Nonequilibrium concentrations are not explicitly reported in the ‘Observation Node’file (Obs_nod.out) of HYDRUS-1D (as opposed to the HYDRUS (2D/3D) software
package). However, they are reported at Print Times in the Mass Balance file
(Balance.out). The same technique as adopted in Example 6b.3 is used to calculate nonequilibrium solute breakthrough at a depth of 50 cm based on the data reported for
Subregion 2 at each Print Time.

106

•
•

•

Consistent with the trend shown in Figure 6c.2, Figure 6c.1 shows a higher nonequilibrium concentration in the case of 'High Resident Time' compared to 'Low
Resident Time'.
Simulations 6c.1 and 6c.2 both did not achieve an equilibrium state, and hence
the non-equilibrium concentration did not reach the saturation concentration of
2.67 mmol/cm3 (Kd/η = 0.267/0.1 = 2.67 mmol/cm3). Increasing the simulation
time and the inlet concentration in simulation 6c.3 did achieve this goal (see theblue time series marked as ‘Steady State’ in Figure 6c.4).
A practical application of modeling sorption processes is for reactive barriers
installed down-gradient of a contaminated groundwater plume. One common
type of reactive barrier is an interception-sorption trench filled with sorptive
material that captures the groundwater contaminants. The sorptive material isthen subsequently disposed of in a land fill when it can longer attract morecontaminants (when it reaches the saturation concentration limit).

Figure 6c.4. Effect of the resident time on the non-equilibrium sorbed
solute concentration.

107

108

Example 6d: Advanced Solute Transport: PhysicalNonequilibrium (Dual-Porosity) Transport
Example Description and Objectives
Next we use Example 6d to simulate water flow and solute transport through a 100-cm
horizontal uniform saturated soil profile having a 20% pressure head gradient imposed
across two Constant Pressure Head boundaries. Solute-free water is introduced to thesoil profile through the ‘Constant Pressure Head' boundary’ on the left during the entire
50-day simulation. The soil structure is such that 90% of the pore space is immobile(similarly as often in fractured rock), while the remaining 10% represents the mobile
fraction. Physical nonequilibrium is considered with a first-order coefficient (Alpha)controlling the time-dependent solute mass transfer process between the mobile andimmobile zones (Figure 6d.1).

Figure 6d.1. Conceptual physical nonequilibrium model for water flow and solutetransport in a dual-porosity medium. In the plots, θ is the water content, θmo and θim
are water contents of the mobile and immobile flow regions; cmo and cim areconcentrations of the corresponding regions, while S is the total solute content of theliquid phase.
Linear sorption is modeled in the mobile phase with a retardation factor Rm = 1+fρKd/θm= 1+1×1.5×0.0267/0.04 = 2.0 with f (dimensionless) being the fraction of sorption sites
in contact with the mobile phase (assumed 1). All other parameters are the same as in
Example 6a, except that the mobile phase is only 10% of the total pore space;additionally, we use Kd=0.0267 to maintain R at 2.0. The soil profile has initially a unit
solute concentration in the mobile phase. Since R=2.0, an equal mass of solute is
initially present in the mobile liquid phase and the sorbed phase. The initial soluteconcentration in the immobile zone is equal to 0.5 mmol/cm3. This variable is referred
to as ‘initial sorbed concentration’ in the ‘Graphical Editor’. Note that no sorption
considered is in the immobile zone since f=1, i.e., all sorption occurs in the mobile
liquid zone. This example is not necessarily typical since most of the sorption usuallyoccurs in the immobile phase where smaller pores are in contact with a much larger
fraction of the surface area than larger pores. The model formulation used in HYDRUS1D is based mostly on early work by van Genuchten and Wierenga (1976).
109

The example illustrates solute transport in a dual-porosity medium subject to physicalnonequilibrium whereby time-dependent solute mass transfer takes place between themobile and immobile zones. The immobile zone represents poorly connected or deadend pores. The model is well suited to simulate flow and transport in systems withpreferential flow paths such as fractured rocks, cracked clays or soils containing
macropores due to decayed root channels or earthworm channels (Mallants et al.,1996a, b)

Pre-Processing
Project Manager- Select Projects Tab- Select Example 6c1- Button “Copy”- New Name: Example 6d- Description: Advance Solute Transport - Physical non-equilibrium- Button “OK”
- Select Example 6d- Button “Open”
Main Processes- Heading: Solute transport with physical non-equilibrium- Keep unchecked "Water Flow"- Keep checked "Solute Transport" and "Standard Solute Transport"- Button “OK”
Print Information- Print Options: Leave as Default- Number of Print Times: 10- Button “Select Print Times”- Print Times: 0.5, 1, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50- Button “OK”
- Button “Next”
Solute Transport
Non-equilibrium Solute Transport Models- Check Dual Porosity (Mobile-Immobile Water) Model (Physical nonequilibrium)
Iteration Criteria - Only for Nonlinear Problems- Absolute Concentration Tolerance: 0 (since this is a linear transport (no nonlinear sorption and no blocking), there is no need to specify iteration criteria forsolute transport)- Relative Concentration Tolerance: 0- Maximum Number of Iteration: 0- Pulse duration (day)" 0- Button “Next”
Solute Transport Parameters- Frac. = 1.00 (all sorption sites are in contact with mobile water; perhaps press
F1 here for a precise definition of the various parameters)110

- Thlmob. = 0.36 (the immobile water content; considering that porosity = 0.4; themobile water content = (0.4-0.36=) 0.04)- Button “Next”
Solute Transport and Reaction Parameters – Solute 1- Kd = 0.0267 (cm3/g)- Nu = 0 (linear sorption isotherm)- Beta = 1 (linear sorption isotherm)- Alpha = 0.02/day (you need to scroll to the right of the table to find this lastparameter)- Button “Next”
Solute Transport Boundary Conditions- Sol. No. 1 Bound. Cond. = 0 mmol/cm3 (no inflow; solute is initially in the soil
profile)- Button “Next”
Soil Profile – Graphical Editor- Menu: Conditions->Initial Conditions->Pressure Head
- Button (left side) “Edit Condition”: select with the Mouse the entire domain andthen left click- Pressure Head Distribution, Top value = 20
Bottom value = 0; uncheck 'Use top value for both'- Button “OK”
- Menu: Conditions->Initial Conditions->Concentration (represents the mobile
phase)- Button (left side) “Edit Condition”: select with the Mouse the entire profile andthen left click the mouse again and specify a concentration 1.0 (mmol/cm3).
- Menu: Conditions->Initial Conditions->Sorbed Concentration (for the dualporosity model, this variable represents the immobile zone)- Button (left side) “Edit Condition”: select with the Mouse the entire profile, thenleft click the mouse again, and specify a concentration of 0.5 (mmol/cm3).
- Menu: File->Save Data
- Menu: File->Exit
Soil Profile – Summary- Button “Next”
Execute HYDRUS-1D- Menu: Calculation->Execute HYDRUS-1D or a Toolbar button.

Post-Processing
Observation Points
Referring to Figure 6d.2:• A sharp drop occurred in the concentration of the mobile phase during the first
3 days of the simulation. This reflects normal dispersive flow in the mobile111

•phase with an open (connected) pore structure and a high flow rate, during
which most of the solutes in the mobile phase are being leached out.
After 3 days, a very long tail developed whereby concentrations continued todecrease at a slower rate. This is due to the slow migration of solutes from theimmobile zone to the mobile phase owing to the concentration gradient between
the immobile and mobile zones (the highest concentration in the immobilezone).

Figure 6d.2 Solute concentration in the mobile phase versus time at a depth of 50 cm.

Profile Information
The concentration profiles in Figures 6d.3 and 6d.4 show sharper fronts in the mobile
phase compared to those in the immobile zone since concentrations in the latterdecrease more uniformly along the profile as a consequence of the slow rate dictatedby the first-order mass transfer parameter.

Figure 6d.3 Concentration profiles for the mobile zone at different times.

112

Figure 6d.4. Concentration profiles for the immobile zone at different times.

Solute Fluxes
Referring to Figure 6d.5:• The nonequilibrium solute flux (i.e., the solute flux between the mobile andimmobile regions) slightly increased at the beginning of the simulationreflecting solute migration from the mobile to the immobile zone (since theinitial concentration of the mobile zone is higher than that of the immobilezone). Following a rapid drop in the concentration of the mobile phase (due to
inflow of solute-free water) as shown in Figure 6d.1, the mass transfer process
reverses and solutes move from the immobile to the mobile phase.
• The bottom flux increased at a relatively fast rate at the beginning of thesimulation, thus reflecting the steep decrease in the mobile concentration (dueto the inlet solute-free flux).
• The two curves become almost parallel at large times since the slow migrationof solutes from the immobile zone becomes almost equal to the solute flux
across the bottom boundary.

113

Figure 6d.5. Cumulative solute fluxes.
Mass Balance Information
The mass balance information at Time=0 (initial conditions) is given in Table 6d.1
where:
W-volume
ConcVol:
ConcVolIm:
cMean
cMeanIm:

Total water volume in the profile
Total Solute Mass in equilibrium region (i.e., in the mobile phase andsorption sites in contact with the mobile phase)
Total Solute Mass in nonequilibrium region (i.e., in the immobile
liquid zone and sorption sites in contact with this zone)
Mean concentration in the mobile liquid zone
Mean concentration in the immobile liquid zone

Table 6d.1. Mass balance information at Time = 0 (initial conditions).
Time[days]

0

Sub-region number
LengthW-volume
In-flow
hMean
ConcVol
cMean
ConcVolIm
cMeanImTop flux
Bot flux

1
[cm]
[cm]
[cm/day][cm][mmol/cm2]
mmol/cm2[mmol/cm2]
mmol/cm2[cm/day][cm/day]

114

100
40
0
10
8
1
18
0.5
-4
-4

100
40
0
10
8
1
18
0.5

Interpretation of initial (at Time=0) solute masses (ConcVol and ConcVolIm) in theMass Balance File
W-volume = Qs × L = 0.40 × 100 =40 cm
ConcVol, Sm
Initial concentration in the mobile zone, Cm = 1 mmol/cm3
Thlmob. (the immobile water content, θim) = 0.36 cm3/cm3
The mobile water content, θm = Qs – ThImob = 0.40 – 0.36 = 0.04 cm3/cm3
Length of the soil profile, L = 100 cm
Bulk Density, ρ = 1.5 g/cm3
The distribution coefficient, Kd = 0.0267 cm3/g
Fraction of sorption sites in contact with the mobile zone, f = 1
Sm = Cm(fρKd+θ m)L=1×(1×1.5×0.0267+0.04)×100=8.005 mmol/cm2
ConcVolIm, Sim
Initial concentration in immobile zone Cim = 0.50 mmol/cm3
Sim=Cim[(1-f)ρKd+θ im)L=0.5×(0×1.5×0.0267+0.36) ×100=18 mmol/cm2

115

116

Example 6e: Advanced Solute Transport: Combined
Physical and Chemical Nonequilibrium
Example Description and Objectives
Example 6e simulates water flow and solute transport through a 100-cm horizontal
uniform saturated soil profile having a 20% pressure head gradient imposed across two
Constant Pressure Head boundaries. Water with a solute concentration of 2 mmol/cm3
is introduced into the soil profile through the upper (left) ‘Constant Pressure Head'boundary’ during the entire 50-day simulation. The soil structure is such that 50% of
the pore space (θs =0.4) is immobile (θim =0.2). Physical nonequilibrium is modeledusing a first-order mass transfer coefficient (Alpha) that controls time-dependent solute
mass transfer between the mobile and immobile zones. All sorption sites are assumed
to be in contact with the mobile phase (f=1) and are split into those subjected toequilibrium sorption (10%; FracM=0.1) and those with kinetic sorption (90%).
Equilibrium sorption is modeled using equilibrium (instantaneous) linear sorption with
R=2.0. Kinetic sorption (chemical nonequilibrium) is modeled using a first-ordercoefficient (Omega) that controls time-dependent solute mass transfer between theliquid phase and the sorption sites of the mobile zone (see Figure 6e.1). The soil profile
initially has a unit solute concentration in the mobile zone, which implicitly implies thatan equal mass of sorbed solutes is initially present on the equilibrium sorption sites (tohonor the linear sorption isotherm with R=2.0). Initial solute concentration in theimmobile zone is equal to 0.5 mmol/cm3. Note that after a sufficiently long run time,
the profile reaches a steady-state situation whereby non-equilibrium sites do reachequilibrium and thus satisfy the condition imposed by the linear sorption isotherm with
R=2.0 (solute masses in the liquid mobile and sorbed mobile phases become equal).
This example demonstrates solute transport with combined physical and chemical nonequilibrium whereby time-dependent solute mass transfer takes place between themobile and immobile zones, in combination with instantaneous and time-dependent
mass transfer between the liquid phase and sorption sites of the mobile zone.

Pre-Processing
Project Manager- Select Projects Tab- Select Example 6d- Button “Copy”- New Name: Example 6e- Description: Solute transport with physical and chemical non-equilibrium- Button “OK”
- Select Example 6e- Button “Open”

117Solute mass transfer withinliquid phase between mobile
and immobile zones

Non-equilibrium
Sorption sites
Pore space
Immobile zone
‘ThImob’

Equilibrium
Sorption sites

Solid soil particles
1 - QsFigure 6e.1. Conceptualization of the dual-porosity model with two-site sorption in
the mobile zone.

Main Processes- Heading: Solute transport with physical and chemical non-equilibrium- Button “Next”
Geometry Information- Number of Layers for Mass Balance: 2 (to obtain the mass balance of each soil
layer)- Button “Next”
Print Information- Print Options: Leave as Default- Number of Print Times: 25- Button “Select Print Times”- Button “Default”- Button “OK”
- Button “Next”
Solute Transport
Non-equilibrium Solute Transport Models- Check Dual Porosity Model with Two-Site Sorption in the Mobile Zone
(Physical and Chemical non-equilibrium)- Pulse duration (day) 50- Button “Next”
118Solute Transport Parameters- Disp = 10 cm- Frac. = 1.00 (all sorption sites are in contact with mobile water)- Thlmob. = 0.20 (immobile water content)- Button “Next”
Solute Transport and Reaction Parameters – Solute 1- Kd = 0.133 (cm3/g)- Nu = 0 (linear sorption)- Beta = 1 (linear sorption)- Frac_M = 0.10 (scroll to the right to find this parameter; fraction of sorption
sites in contact with the mobile region that are subject to instantaneous sorption)- Omega = 0.50 day-1 (first-order mass transfer coefficient representing kinetic
sorption)- α = 0.02 day-1 (first-order mass transfer coefficient representing mass transfer
between mobile and immobile zones)- Button “Next”
Solute Transport Boundary Conditions- Sol. No. 1 Bound. Cond. = 2 mmol/cm3- Button “Next”
Soil Profile – Graphical Editor- Menu: Conditions->Subregions- Button (left side) “Edit Condition”, select with the Mouse node 51, left click themouse
Material Index: 2 (this will select a 1-cm element located at a depth of 50-51
cm)- Menu: File->Save Data
- Menu: File->Exit
Execute HYDRUS-1D- Menu: Calculation->Execute HYDRUS-1D or a Toolbar button.

Post-Processing
Observation Points
Referring to Figure 6e.2:• An initial drop occurs in concentrations at both observation points as solutes
migrate from the mobile to the immobile region, thus reducing the concentrationin the former.
• After about 1 and 3 days, solutes introduced at the upper Constant Pressure
Head boundary reach the depths of 50 cm and 100 cm, respectively, thusincreasing concentrations.
• At 50 days, concentrations at both observation points almost reach a steadystate (approaching 2 mmol/cm3, which is equal to the inlet concentration at theupper boundary).
119

Figure 6e.2. Solute concentrations in the mobile region versus time at depths of 50
and 100 cm.
Profile Information:
Referring to Figure 6e.3, which shows concentration profiles in the mobile zone at 0,2, and 50 days:• The initial profile (Time = 0) has a uniform initial unit concentration of 1
mmol/cm3 in the mobile zone.
• After 2 days, the concentration near the upper (surface) boundary increases dueto the solute (2 mmol/cm3) inflow at the constant pressure head boundary. At
this early time, the concentration increase was restricted to the top 40 cm of the
profile. Concentrations below this point continue decreasing due to themigration of solutes from the mobile to the immobile zone (as noted in Figure6e.1) at early times.
• At 50 days, the concentration profile near the surface (the top 23 cm) has alreadyreached the steady-state concentration of 2 mmol/cm3. The concentration in therest of the soil profile is slightly lower (1.96 mmol/cm3 at the 100-cm depth).

120

Figure 6e.3. Concentration profiles in the mobile zone at different times.
Referring to Figure 6e.4, which shows concentration profiles in the immobile zone at
0, 2, and 50 days:• The initial profile (Time = 0) has a uniform concentration of 0.5 mmol/cm3 in
the immobile zone.
• After 2 days, concentrations in the immobile zone slowly increase across theentire profile. In the lower 40 cm, there is a linear increase driven by the firstorder mass transfer parameter (Omega). Note that the straight line is parallelwith the initial condition (Time=0). The top 50 cm, however, displays anonlinear concentration profile resulting from a higher migration of solutes
from the mobile zone as driven by the concentration gradient caused by the inlet
concentration of 2 mmol/cm3 at the constant pressure head boundary. Had thisgradient been absent, concentrations would have been constant throughout theentire profile (represented by a dotted red line).
• At 50 days, the concentration profile has not yet reached the steady-state
concentration of 2 mmol/cm3 due to slow migration of solutes from the mobile
zone.

121

Figure 6e.4. Concentration profiles in the immobile zone at different times.
Solute Fluxes
Referring to Figure 6e.5:• The cumulative non-equilibrium mass transfer represents the sum of masstransfer between the mobile and immobile zones, and mass transfer between
the liquid phase and chemical non-equilibrium sorption sites (i.e., mass transferfrom the equilibrium phase to the nonequilibrium phase). It seems that nonequilibrium mass transfer has almost approached a steady-state (since the blueline is almost horizontal) (as also indicated by Figs 6e.3 and 6e.4.
• The cumulative surface and bottom fluxes will be used to check against solute
masses reported in the Mass Balance Information.

Figure 6e.5. Cumulative solute fluxes versus time.

122

Mass Balance Information•

Mass balance information at Time=0 (initial conditions) and Time=50 days forthe entire soil profile is given in Table 6e.1 where:

W-volume
ConcVol

cMean
ConcVolIm
cMeanIm
SorbVolIm2

Total water volume in the profile
Total solute mass in equilibrium phases, i.e., the mobile liquid phaseand equilibrium sorption sites of the mobile zone (this excludes themass in the immobile liquid phase and non-equilibrium sites of themobile zone).
Mean concentration of the mobile liquid phase
Total solute mass of the immobile liquid phase
Mean concentration of the immobile liquid phase
Total solute mass of non-equilibrium sorption sites of the mobile zone

Table 6e.1. Mass balance information for the entire profile at various times.
Time (Days)

0

Sub-region number
LengthW-volume
In-flow
hMean
ConcVol
cMean
ConcVolIm
cMeanIm
SorbVolIm2•

50

Steady State

1
[cm]
[cm]
[cm/day][cm][mmol/cm2][mmol/cm3][mmol/cm2][mmol/cm3][mmol/cm2]

100
40
0
10
21.995
1.0
10
0.5
0

0
10
43.63
1.98
38.75
1.93
35.56

0
10
44
2
40
2
36

Interpretation of solute masses (ConcVol, ConcVolIm and SorbVolIm2) in theMass Balance Information

Referring to Table 6e.1:
At Time=0
ConcVol
Initial concentration in the mobile zone, Cm = 1 mmol/cm3
Thlmob. (immobile water content, θim) = 0.20 cm3/cm3
Mobile water content, θm = Qs – ThImob = 0.40 – 0.20 = 0.20 cm3/cm3
Length of the soil profile, L = 100 cm
Bulk Density, ρ =1.5 g/cm3
Distribution coefficient, Kd = 0.133 cm3/g
Retardation factor of the mobile zone:
Rm=1+fρKd/θm=1+1×1.5×0.133/0.20=2.0
123

The fraction of sorption sites in contact with the mobile zone, f=1
The fraction of sorption sites in contact with the mobile zone that is atequilibrium, fm=0.1
Sm=Cm(ffmρKd+θ m)L=1×(1×0.1×1.5×0.133+0.20) ×100=22 mmol/cm2
ConcVolIm
Initial concentration of the immobile zone Cim = 0.50 mmol/cm3
Sim=Cim[(1-f)ρKd+θ im)L=0.5×(0×1.5×0.133+0.20) ×100=10 mmol/cm2
SorbVolIm2
Initial sorbed mass of chemical non-equilibrium sites is equal to zeroAt Time=50 days (compare with flux data)
From Mass Balance data:
Net mass of solutes added to a profile after 50 days = Solute mass at 50 days – initialsolute mass
= (43.63+38.75+35.56) – (22+10) = 85.94 mmol/cm2
Non-equilibrium mass transfer between mobile and immobile zones = 38.75-10 =28.75 mmol/cm2
Non-equilibrium sorption mass transfer = 35.56-0.0 = 35.56 mmol/cm2
Total non-equilibrium mass transfer = 28.75 + 35.56 = 64.31 mmol/cm2
From Solute Flux data:
Cumulative Surface Solute Flux = 400 mmol/cm2
Cumulative Bottom Solute Flux = -314.43 mmol/cm2
Net mass of solutes added to profile after 50 days = 400 – 314.43 = 85.57 mmol/cm2
Cumulative non-equilibrium mass transfer = 64.31 mmol/cm2
Note that the small discrepancy in the total solute mass added is reflected in the error
term (CncBalT) in the Mass Balance information.
At steady state:
ConcVol = 44 mmol/cm2, comprising 40 mmol/cm2 in the liquid mobile phase and 4
mmol/cm2 in the sorbed mobile phase.
SorbVolIm2 is equal to 36 mmol/cm2
Total sorbed mass (equilibrium + non-equilibrium sites) = 4 + 36 = 40 mmol/cm2
Total sorbed mass = total mass in the liquid phase, which satisfies R=2 as anequilibrium state has been reached.
•

Figure 6d.6 shows the solute breakthrough in the immobile zone at a depth of
50 cm using data collated from the Mass Balance Information File (cMeanIm
for Subregion 2).

124

Figure 6e.6. Immobile solute concentrations versus time at a depth of 50 cm.
•

Figure 6d.7 shows the individual mass transfer for the physical (ConcVolIm)and chemical (SorbVolIm2) non-equilibrium regions. If the initial solute massin the immobile zone (10 mmol/cm2) is deducted from the physical nonequilibrium mass values (red line), and then added to the chemical nonequilibrium mass (black line), it will yield the total non-equilibrium masstransfer reported in Figure 6e.4 (blue line).

Figure 6e.7. Non-equilibrium solute masses versus time.

125

126

Example 7a: Inverse Modeling: One-Step Outflow Method
Example Description and Objectives
Example 7a considers the inverse solution for a one-step outflow experiment (Figure7a.1). Data presented by Kool et al. (1985) and used in Example 6 of the HYDRUS-1D
manual (p. 170) are used in the analysis.
An undisturbed core sample, having a length of 3.95 cm and a diameter of 5.4 cm, wasequilibrated at near-saturated conditions in a Tempe pressure cell. The retention curve
was first measured for pressure heads up to -10 m. Water contents at pressure heads of-30 and -150 m were measured on disturbed samples. After resaturation, a pneumaticpressure of 10 m was imposed at the top of the sample, and the cumulative outflow wasrecorded with time. The position of the measuring burette was adjusted manually every
time reading was made to maintain a constant head lower boundary condition. At theend of the experiment, the soil was resaturated, and the saturated hydraulic conductivityof the soil and porous plate were measured with a falling head method.
Three hydraulic parameters (α, n, and θr) are to be estimated by numerical inversion of
the observed cumulative outflow data and the measured water content at a pressure headof -150 m. Since water exits the soil column across a ceramic plate, the flow problem
involves a two-layered system. To be able to simulate flow through the ceramic bottom
plate without having to modify the code, the HYDRUS model must be able to simulate
flow through materials with very high air entry values such that the ceramic plateremains saturated at all times during the outflow experiment. A very high air entry valueof the plate was ensured by specifying parameter α to be 10-20 (1/cm). The soil profile,which consists of a 3.95-cm long soil sample and a 0.57-cm thick ceramic plate, wasdiscretized into 50 nodes with five nodes representing the ceramic plate. Only a few
nodes were used for the ceramic plate since the plate remained saturated during theentire experiment, thus causing the flow process in the plate to be linear.

Figure 7a.1. Schematic of the one-step outflow measurement device (from Green etal., 1998).

Pre-Processing
127Project Manager- Click Button "New"- Type Name: Onestep- Type Description: Onestep Outflow Method- Button "OK"Main Processes- Heading: Onestep Outflow Method- Check "Inverse Solution"- Button “Next”
Inverse Solution
When the option "Inverse Solution" is selected in the Main Processes window, the next
dialog that appears is the "Inverse Solution" dialog window (Fig. 7a.2). In this window,
you must first select which parameters are to be estimated (Estimate …). One canestimate either a) soil hydraulic parameters, b) solute transport parameters, c) heattransport parameters, or any combination of these three sets of parameters. When only
water flow is simulated, then only soil hydraulic parameters can be optimized, and theother sets of parameters (i.e., for solute and heat transport) are disabled.
You must also select the method of Weighting of Inversion Data in the objectivefunction. One can choose between no weighting, weighting by mean ratios (of different
data types in the objective function), or weighting by standard deviations. When noweighting is selected, users need to supply weights for particular data points in the Datafor Inverse Solution dialog window (Fig. 7a.4). When weighting by mean ratio or
standard deviation is selected, the code will calculate either means or standard
deviations for different data sets (e.g., water contents, pressure heads, concentrations,…), and then adjust the weights proportionally. These internal weights can still be
multiplied by weights from the Data for Inverse Solution dialog window.
The objective function for the inverse estimation of solute transport parameters can be
defined using different types of concentrations. Available Concentration Types are:a) the resident concentration in the liquid phase, b) a log-transformation of the resident
concentration in the liquid phase, c) the outflow (flux) concentration, d) the soluteconcentration flux, e) the cumulative concentration (solute) flux, and f) the total
resident concentration. The total resident concentration includes concentrations in thesorbed and nonequilibrium phases.
The maximum number of iterations for the inverse solution is also specified in this
dialog window. If you select a zero for the number of iterations, then only the direct
simulation is carried out. However, you can then still enter the measured data, in whichcase the code compares results of the direct simulation with the measured data.
HYDRUS allows users to optimize up to 15 parameters (NPaD is the maximum allowed
number of optimized parameters). We, however, do not recommend optimizing somany parameters simultaneously. Unsaturated flow problems are inherently ill-posed,
and thus not too many parameters can be successfully optimized simultaneously. When
more parameters are optimized, the problem often becomes nonunique. The problems
128

of unsaturated flow are different than those of saturated flow where it is common to
optimize many more parameters simultaneously. For details see Šimůnek and Hopmans(2002), Hopmans et al. (2002), and Šimůnek et al. (2002).
Finally, the "Number of Data Points in Objective Function" must be specified. This is
the total number of all data points that will be considered in the objective function. This
number needs to be specified in the Data for Inverse Solution dialog window (Fig.
7a.4).

Figure 7a.2. The Inverse Solution dialog window.
- Check "Soil Hydraulic Parameters"- Check "No Internal Weighting" (we want to replicate the example in the same
way as originally solved by Kool et al. (1985)- Max. Number of Iterations:20- Number of Data Points:
10- Button “Next”
Geometry Information- Number of Soil Materials:
2 (soil itself and the ceramic plate)- Number of Layers for Mass Balances: 2 (ditto)- Depth of the Soil Profile (cm):4.52- Button “Next”
Time Information- Time Units:- Final time (hour):- Initial Time Step (hour):- Minimum Time Step (hour):- Maximum Time Step (hour):- Button “Next”

Hours
100
0.001
0.000110

Print Information
129

- Uncheck "Screen Output"(since the direct simulations are carried out many times during the inverse
parameter estimation process, we are not interested in detailed information about
each run, but only in information printed after each inverse solution iteration,such as the value of the objective function and values of optimized parameters)- Number of Print Times:
11- Button "Select Print Times": 0.017, 0.033, 0.05, 0.167, 0.5, 1.33, 2.75, 5.417,
10, 15, 100- Button “Next”
Water Flow - Iteration Criteria- Water Content Tolerance: 0.0001- Pressure Head Tolerance (cm): 0.1- Upper Limit of the Tension Interval: 15000 (to extend the Internal Interpolation
Tables all the way to the lowest recorded value, i.e., the wilting point).
- Button “Next”
Water Flow - Soil Hydraulic Model- Leave as Default- Button “Next”
Water Flow - Soil Hydraulic Parameters
When the soil hydraulic parameters are optimized, a different version of the Water Flow
Parameters dialog window (Fig. 7a.3) appears. In this window, which appears for eachsoil material, you must provide initial estimates of the optimized soil hydraulicparameters, specify which parameters are to be optimized (check the appropriate
boxes), and provide parameter constraints for the optimization. Entering zeros (thedefault values) for the minimum and maximum constraints signifies that the parameters
are unconstrained. Initial values specified for parameters that are not optimized will be
kept constant during the inverse optimization process.

Figure 7a.3. The Water Flow Parameters - Inverse Solution dialog window.
- Initial Estimate:- Fitted:- Second Material:- Button “Next”

θr=0.15, θs=0.388, α=0.025, n=1.5, Ks=5.4, l= 0.5
θr, α, and nθr=0., θs=1., α=1e-20, n=1.001, Ks=0.003, l= 0.5

130Water Flow - Boundary Conditions- Upper Boundary Condition: Constant Flux
- Lower Boundary Condition: Constant Pressure Head
- Button “Next”
Water Flow - Constant BC- Upper Boundary Flux: 0- Button “Next”
Data for Inverse Solution
In the table Data for Inverse Solution (Fig. 7a.4) you must specify the measured datathat define the objective function to be minimized during the parameter optimization
process. Many different types of data can be used to define the objective function. How
the values in the X and Y columns are interpreted depends on the Type and Position
values. Weight is the weight associated with a particular data point. Table 7a.1 lists
selected data types that can be included in the objective function. A complete list ofvarious data types that can be used to define the objective function and how the X andY columns are interpreted is given in the online help and in the manual of HYDRUS1D. You can double-check the correct definition of the data types by clicking on ‘Show
list boxes (not recommended for large data sets’ at the bottom of the window.
Table 7a.1. Data types for the objective function (Inverse Problem).

131

Figure 7a.4. The Data for Inverse Solution dialog window.
The first nine data points in this example are interpreted in the program as (Type=0)
cumulative outflow Q (column Y) at time t (column X) through the bottom boundary(Position=2). Cumulative outflow values are negative because outflow is downward(against the z axis). The last (10th) value is interpreted as (Type=5), i.e., a data point of
the retention curve (water content θ (column Y) at pressure head h (column X) for
material 1 (Position=1)).
- Enter the following data (copy from the Excel file)
0.017-0.0786 02
1
0.033-0.1616 02
1
0.05-0.2097 021
0.167-0.3408 02
1
0.5-0.4456 02
1
1.033-0.49802
1
2.75-0.5614 02
1
5.417-0.5937 02
1
100-0.6824 0
2
1-150000.157
5
1
1- Button “Next”
Soil Profile - Graphical Editor- Menu: Options->Grid or Toolbar: Height: 0.05
- Menu: Condition->Profile Discretization or Toolbar- Button (Edit Bar) "Number": 50- Button (Edit Bar) "Insert Fixed" at 3.95 cm- Button (Edit Bar) "Density": deselect "Use upper", upper density =0.1 at 3.95cm
- Menu: Condition->Initial Condition->Pressure Head or Toolbar- Button (Edit Bar) "Edit Condition"- Select entire profile: Top value=-2, Bottom value=2.52
Deselect "Use top value for both"- Select the lowest node = -1000 cm
132

- Menu: Condition->Material Distribution or Toolbar- Button (Edit Bar) "Edit Condition"- Select the ceramic plate and specify "Material Index"=2- Ditto for "subregions"- Menu: File->Save Data or Toolbar- Menu: File->Exit
Soil Profile - Summary- Button “Next”
Execute HYDRUS- Menu: Calculation->Execute HYDRUS or a Toolbar button.

Post-Processing
Figure 7a.5 shows the measured cumulative outflow curve versus time, as well as thebest fit obtained with HYDRUS-1D. Initial and final parameter estimates are listed inTable 7a.2. Notice the very good fit of the measured data in Figure 7a.5, with R2 being
0.9987. Figure 7a.6 shows a comparison of the predicted and measured retention
curves, as well as a comparison of the diffusivity curve obtained by parameter
estimation and D(θ) values calculated independently using the method of Passioura[1976]. Again, notice the close agreement between predicted and measured values.

Figure 7a.5. Measured and optimized cumulative outflow versus time forthe one-step outflow experiment.

133

Figure 7a.6. Observed and predicted retention curves (left) and calculated andpredicted diffusivities (right).
Table 7a.2. Initial estimates of and optimized parameters for one-step outflow experiment.

Parameter

Initial Estimate

Final Valueθr [cm3/cm3]
θs [cm3/cm3]
α [1/cm]

0.15
0.388+
0.025
1.5
5.4+
0.003+
0.5+

0.166

n [-]
Ks [cm/h]Ks,cer [cm/h]l [-]+Not optimized

134

0.0363
1.42

Example 7b: Inverse Modeling: Multi-Step Outflow Method
Example Description and Objectives
Although initial applications of the inverse approach to one-step outflow data appeared
promising, later studies revealed possible problems in terms of non-uniqueness of the
optimized parameters (e.g., van Dam et al., 1992, 1994). To circumvent uniquenessproblems, van Dam et al. (1994) conducted and analyzed outflow experiments in which
the pneumatic pressure was increased in several smaller steps. Eching and Hopmans(1993) similarly estimated soil hydraulic parameters from simultaneous measurementsof transient cumulative outflow and the soil water pressure head inside of the soil
sample during a multistep outflow experiment.
In this test example, we analyze a multistep outflow experiment with simultaneous
measurement of the pressure head inside the soil sample (Hopmans, personalcommunication). The experimental setup consisted of a 6-cm long soil column in a
Tempe pressure cell modified to accommodate a micro-tensiometer-transducer system.
A tensiometer was installed, with the cup centered 3 cm below the soil surface. The soil
sample was saturated from the bottom and subsequently equilibrated to an initial soil
water pressure head of -25 cm at the soil surface. Pressures heads of -100, -200, -400,
and -700 cm were applied subsequently in consecutive steps at 0, 12.41, 48.12, and105.92 hours, respectively.

Pre-Processing
Project Manager- Button "New"- Name:- Description:- Button "OK"MultistepMultistep Outflow Method

Main Processes- Heading: Multistep Outflow Method- Check "Inverse Solution"- Button “Next”
Inverse Solution- Check "Soil Hydraulic Parameters"- Check "Weighting by Standard Deviations"- Max. Number of Iteration:
50- Number of Data Points:
229- Button “Next”
Geometry Information- Length Units:
cm- Number of Soil Materials:
1- Number of Layers for Mass Balances: 1- Depth of the Soil Profile:
6.0
135

- Button “Next”
Time Information- Time Units:
Hours- Final time (hours):
190.384- Initial Time Step (hours):
0.01- Minimum Time Step (hours):
0.0001- Maximum Time Step (hours):
120- Check Time-Variable Boundary Conditions- Number of Time-Variable Boundary Records: 4- Button “Next”
Print Information- Uncheck "Screen Output"- Number of Print Times:
10- Button "Select Print Times": 0.01, 0.1, 0.2, 19.038, 38.077, 57. 115, 133, 27,
152.31, 171.35, 190.384- Button “Next”
Water Flow - Iteration Criteria- Leave the default values- Button “Next”
Water Flow - Soil Hydraulic Model- Leave the default model- Button “Next”
Water Flow - Soil Hydraulic Parameters- Initial Estimate: Select Loam from the Catalog for Initial Estimate- Fitted: All six parameters- Button “Next”
Water Flow - Boundary Conditions- Upper Boundary Condition: Constant Flux
- Lower Boundary Condition: Variable Pressure Head
- Button “Next”
Water Flow - Constant BC- Upper Boundary Flux: 0- Button “Next”
Variable Boundary Conditions- Enter the following values:
Time
GWL
12.42
-100
48.12
-200
105.92-400
190.384-705- Button “Next”

136Data for Inverse Solution- Copy from the Excel File: Worksheet HYDRUS-1D (One-Step)- Button “Next”
Soil Profile - Graphical Editor- Menu: Condition->Profile Discretization or Toolbar- Button (Edit Bar) "Number": 50- Button (Edit Bar) "Insert Fixed" - at -3 cm- Button (Edit Bar) "Density": Select the top node and specify Density=3- Button (Edit Bar) "Density": Select the middle node (-3 cm) and set Density=2- Menu: Condition->Initial Condition->Pressure Head or Toolbar- Button (Edit Bar) "Edit condition"- Select entire profile: Top value=-25, Bottom value=-19- Deselect "Use top value for both"- Menu: Conditions->Observation Points or Toolbar- Button (Edit Bar) “Insert”: Set Observation Point at -3 cm- Menu: File->Save Data or Toolbar- Menu: File->Exit
Soil Profile - Summary- Button “Next”
Execute HYDRUS- Menu: Calculation->Execute HYDRUS or a Toolbar button.

The computations (using the inverse computational module h1d_clci.exe) are now
being carried out, with parts of the results being displayed in a DOS window (Fig. 7b.1).

Figure 7b.1. The DOS window in which the inverse computational module
h1d_clci.exe is executed.
This window shows selected information about the inverse run, such as the number of
iterations (Iteration), the value of the objective function being minimized (SSQ), andvalues of optimized parameters at each inverse iteration. When the simulation ends, theRun Time in seconds is displayed. You are now asked to “Press Enter to Continue”.
Once the keyboard button “Enter” is pressed, the DOS window is closed, and the137program returns to the HYDRUS-1D GUI. Detailed output information about
individual direct runs are not shown in the DOS window since we unchecked the“Screen Output” box in the Print Information dialogue window.

Post-Processing
Figure 7b.2 below compares the measured and optimized cumulative outflow curves
for the soil sample (left) and measured and optimized pressure heads (right). Excellent
agreement was obtained for both variables. The final fit had an R2 of 0.9995. Table 7b.1
lists initial estimates and final values of the six optimized parameters.

Figure 7b.2. Measured and optimized cumulative bottom flux (left) and pressure
heads inside the soil sample (right) for a multistep outflow experiment.

Table 7b.1 Initial and final parameter values obtained for the multistep
outflow experiment.

138

Example 8: Heat Movement and Reactive Solute Transport
Example Description and Objectives
Example 8 simulates combined heat and solute transport through a 100-cm verticaluniform soil profile that maintains hydrostatic equilibrium conditions throughout theentire simulation. We further consider a stagnant boundary condition at the soil surface
to account for volatile chemicals. The imposed upper ‘Stagnant BC for Volatile Solutes’
allows the diffusion of gas (oxygen in this case) from the atmosphere into the soil
profile. The upper ‘Temperature BC’ allows the conduction of heat from theatmosphere to the soil profile that is initially at 20 oC. Diurnal temperature variations(daily cycles) are simulated using a sinus wave with an amplitude of 5 oC. The imposed
initial and boundary conditions that result in no-flow conditions restrict heat and solutetransport processes to be solely conductive and diffusive, respectively. Oxygen that
diffuses into the soil profile is subsequently consumed via a first-order reaction process
at a rate of 0.1 d-1 (the soil then consumes approximately 10% of the diffused gas perday).
Note that the standard HYDRUS module is intended mainly to simulate species forwhich the primary variable is the concentration in the liquid phase. Since molecular
diffusion in the gaseous phase can be very dominant for volatile species, for which the
majority of mass may be in the gaseous phase (e.g., oxygen, carbon dioxide), some
restrictions on time steps or the numerical scheme (an implicit solution in this case) are
required to achieve a stable numerical solution. Other HYDRUS modules, such as
UnsatChem or HP1, in which the primary variable can be the concentration in the gas
phase, have much more stable solutions.

Pre-Processing
Project Manager- Select Projects Tab- Button “New”- Name: Example 7- Description: Combined heat and solute transport- Button “OK”
- Select Example 8- Button “Open”
Main Processes- Heading: Combined heat and solute transport- Uncheck "Water Flow" (when water flow is unchecked, the initial condition andcorresponding water content and fluxes will be constant during the simulation).
- Check “Solute Transport”- Check “Heat Transport"- Button “Next”
Geometry Information- Leave as Default- Button “Next”
139Time Information
Time Discretization:- Final Time (day): 4- Maximum Time Step (day): 0.05 [refer to Appendix I.1 for explanation]- Check Time-Variable Boundary Conditions- Number of Time-Variable Records 1- Button “Next”
Print Information- Print Options: Leave as Default- Number of Print Times: 3- Button “Select Print Times”- Print Times (day): 1, 1.5, 4- Button “Next”
Iteration Criteria- Leave as Default- Button “Next”
Water Flow – Soil Hydraulic Model- Leave as Default- Button “Next”
Water Flow – Soil Hydraulic Parameters- Soil Catalogue: Loam
- Button “Next”
Water Flow – Boundary Conditions- Upper Boundary Condition: Atmospheric BC with Surface Runoff- Lower Boundary Condition: Constant Pressure Head
- Button “Next”
Solute Transport
Time Weighting Scheme- Check “Implicit Scheme”- Button “Next”
Solute Transport Parameters- Disp. = 5- Diffus. W. = 1.7- Diffus. G. = 17000- Button “Next”
Solute Transport and Reaction Parameters – Solute 1- Kd = 0- Nu = 0- Beta = 1.0- Henry = 32
- SinkGas1 = 0.1
140

- Button “Next”
Solute Transport Boundary Conditions- Upper Boundary Condition: Stagnant BC for Volatile Boundary- Lower Boundary Condition: Zero Concentration Gradient- Stagnant Boundary Layer: 1 cm
- Concentration in the Atmosphere: 0.3 mg/cm3 [0.2 (an oxygen fraction in air) ×
1.43 mg/cm3 (the density of oxygen)]- Initial conditions; Leave as Default- Button “Next”

Figure 8.1. The Solute Transport Boundary Conditions dialog window.

Heat Transport Parameters
You must now select heat transport parameters using the dialogue window shown in
Figure 8.2. Temperature fluctuations may either be explicitly provided as time-variableboundary conditions or may be simulated using a sinus wave having a specified
amplitude and a time period. The latter is usually 1 day (or 1 year) thus simulating
diurnal (or annual) temperature variations. Note that when the “Time-Variable
Boundary Conditions” option is selected, the temperature amplitude is specified in the“Time-Variable Boundary Conditions” dialog window where different values can be
specified on each record. When the soil textural class is selected (e.g., clay, loam, orsand), default heat transport parameters (b1, b2, and b3) are provided (in general in theunits of [WL-1K-1; MLT-3K-1], but in this particular example in [kg cm d-3 K-1]). Users
should specify the volumetric fractions of solid (Solid) and organic matter (Org.M.),
along with the longitudinal thermal dispersivity (Disp.). The units of the default values
for the “Volumetric heat capacity” are in general in [JL-3K-1; ML-1T-2K-1], except in thisparticular example where we use [kg/(cm.day2K)] and hence larger numbers. Note thatthe mass units used for solute concentrations, the bulk density (and the distribution
coefficient), or thermal parameters can be different.
141

Figure 8.2. The Heat Transport Parameters dialog window.
- Temperature Amplitude (C): 0- Interval for one temp. cycle (day): 1- Thermal Conductivity: Select “Loam” (this will provide values for b1, b2, andb3)- Volumetric Heat Capacities: Button Default (this will provide values for Co, Cn,
and Cw)- Solid: 0.59
- Org.M.: 0.15
- Disp.: 5- Button “Next”
Heat Transport Boundary Conditions- Upper Boundary Condition: Temperature BC- Lower Boundary Condition: Zero Gradient- Button “Next”
Time-Variable Boundary Conditions- Time (day): 4- Precip. (cm/day): 0- Evap. (cm/day): 0- hCritA: 50000- tTop (oC): 30- tBot (oC): 20- Ampl (oC): 5- cTop: 0 (not considered)- cBot: 0 (not considered)
Soil Profile – Graphical Editor or Toolbar- Menu: Conditions->Profile Discretization
142

- Button (on left-side Edit Bar) “Insert Fixed”: Select node at the 25-cm depth- Button (on left-side Edit Bar) “Density”: Select surface node
- Nodal Density (Lower Density) = 0.1- Menu: Conditions->Observation Points or Toolbar- Button (on left-side Edit Bar) “Insert”: Select nodes at depths 0 cm, 0.88 cm,
and 25 cm- Menu: Conditions->Initial Conditions->Pressure Head or Toolbar- Button (Edit Bar) “Edit Condition”: select with the Mouse the entire domain, left
click
- Pressure Head Distribution: Top value = -100, Bottom value = 0, uncheck "Usetop value for both"- Menu: Conditions->Initial Conditions->Temperature (Leave as default 20oC)- Menu: Conditions->Initial Conditions->Concentration (Leave as default 0)- Button “OK”
- Menu: File->Save Data or Toolbar- Menu: File->Exit
Soil Profile – Summary- Button “Next”
Execute HYDRUS- Menu: Calculation->Execute HYDRUS or a Toolbar button.

Post-Processing
Observation Points
Figures 8.3 shows liquid concentrations of the diffused gas at depths of 0, 0.88 and 25
cm. Note that the concentration at the soil surface (0 depth) attains almost
instantaneously the maximum (saturation) concentration of 0.009 mg/cm3 (inequilibrium with the atmospheric gas concentration). This apparent artifact is a resultof very fast molecular diffusion across a relatively narrow (0.5 cm) stagnant surface
layer. Had a smaller initial time step been used, the surface concentration would haveincreased more gradually before reaching an equilibrium value. Note how theconcentration at the very small depth of 0.88 cm starts from a smaller value (red line in
Figure 8.3) and then gradually attains the saturation concentration. The transport ofoxygen is slower in the soil due to tortuosity (restricted pore continuity) and retardation(equilibration between gaseous and liquid phases). After 4 days, the concentration at a
depth of 25 cm does not yet reach the saturation limit.

143

Figure 8.3. Solute concentrations at the soil surface and at depths of 0.88 and 25 cm.

Figure 8.4 shows the applied temperature boundary condition, which represents diurnal,
day-night fluctuations (with temperatures of 30 ± 5 oC and maximum temperatures at 1
p.m.). The temperature at 25 cm starts increasing in response to the temperature
gradient. Daily fluctuations are apparent also at a depth of 25 cm, albeit at a much lower
amplitude. If the simulation were extended to much larger times, one could see that at
depths larger than 60 cm the fluctuations gradually disappear.

Figure 8.4. Temperatures at the soil surface and at a depth of 25 cm.

144

Profile Information
Figure 8.5 shows simulated liquid phase oxygen concentration profiles at different
times. The downward advance of oxygen concentrations was mainly due to molecular
diffusion in the air phase (note that the gaseous diffusion coefficient “Diff. G.” is about
10,000 times larger than the liquid phase diffusion coefficient "Diff. W."). This implies
that the concentration advance always will be much faster in drier soils. Nevertheless,the concentration advance is the same in both gaseous and liquid phases due to theassumption of equilibrium partitioning between the two phases (the Henry's law).

Figure 8.5. Liquid phase oxygen concentration profiles at different times.

Figure 8.6 shows the temperature profiles at various times. Note how the imposedboundary temperature (at the soil surface) fluctuates between 25 and 35 oC (30 ± 5 oC).
At a depth of 100 cm, the temperature slowly increases above the initial 20oC; if you
were to extend the simulation to much larger times, a steady state situation equal to
30oC would be reached (i.e., diurnal fluctuations do not affect the soil temperature at
this depth).
Solute Fluxes
Figure 8.7 shows the cumulative solute fluxes, including the mass of oxygen diffusedthrough the soil surface as well as the mass of oxygen consumed by the soil (the firstorder reaction).

145

Figure 8.6. Temperature profiles at different times.

Figure 8.7. Cumulative solute fluxes.

146

Mass Balance Information
The mass balance information shown in Figure 8.8 refers to the following:
Length:Depth of the soil profile
W-volume: Total volume of soil water in the soil profile or selected subregions
In-flow:
Sum of fluxes flowing into and out of the soil profile (or subregions)hMean:
Average pressure head in the soil profile (or subregions)
HeatVol:
Total heat energy stored in the soil profiletMean:
Average temperature in the soil profile (or subregions)
ConcVol: Total solute mass in the soil profile
cMean:
Mean concentration
Top Flux: Flux into (at the top) the soil profile
Bot Flux: Flux out of (at the bottom) the soil profile

Figure 8.8. The mass balance information at Time=0.

Sample calculation for ‘HeatVol’ at the beginning of the simulation
Initial temp. = 20 oC = 293.15 K
Domain volume = 100×1×1 cm3 (assuming a unit width and length in the horizontal
directions)
Volumetric heat capacity of the solid phase, Cn = 1.433 × 1014 kg/(cm.day2K)
Volumetric heat capacity of the organic matter, Co = 1.874 × 1014 kg/(cm.day2K)
Volumetric heat capacity of the liquid phase, Cw = 3.12 × 1014 kg/(cm.day2K)
Solid fraction = 0.59
Organic matter fraction = 0.15
Water content = 0.243 (initial water content)
147

Total Volumetric Heat capacity, Csoil = (1.433×0.59 + 1.874×0.15 + 3.12×0.243) ×
1014 = 1.88 × 1014 kg/(cm.day2K)
Heat volume stored in the soil profile at initial time = Volumetric heat capacity × Soilvolume × Temperature (=Csoil ×Vsoil × T) = 1.88 × 1014 × 100 × 293.15 = 0.552 × 1019
kg.cm2/day2 = 0.552 × 1019 × 10-4 / (60×60×24)2 = 73,945 kg.m2/sec2 (Joules)
Note that HYDRUS-1D reports this value for a one-dimensional soil profile, i.e., in
kg/day2 [MT-2]. This value needs to be multiplied by the land surface area to obtain
the heat storage in three dimensions.

148

Appendix I: Notes on Spatial and Temporal Discretization
Temporal Discretization
Four different time discretizations are used in HYDRUS-1D: (1) time discretizations
associated with the numerical solution, (2) time discretizations associated with theimplementation of boundary conditions, (3) time discretizations which provide printed
output of the simulation results (e.g., nodal values of dependent variables, water, solute
mass balance components, and other information about the flow regime), and (4) time
discretizations associated with data defining the objective function in the inverse modeof HYDRUS-1D (e.g., measured water contents, pressure heads, concentrations, and/or
water or solute fluxes). Discretizations 2, 3, and 4 are mutually independent; theygenerally involve variable time steps as described in the input data file. Discretization
1 starts with a prescribed initial time increment, Δtinit. The time increment, Δt, is
automatically adjusted at each time level according to the following rules:
•
•
•

•

Discretization 1 must coincide with time values resulting from time
discretizations 2, 3, and 4.
Time increments cannot become less than a preselected minimum time step,Δtmin, nor exceed a maximum time step, Δtmax (i.e., Δtmin ≤ Δt ≤ Δtmax).
If, during a particular time step, the number of iterations necessary to reachconvergence is ≤3 (Itmin), the time increment for the next time step is increasedby multiplying Δt by a predetermined constant >1 (k1, usually between 1.1 and
1.5). If the number of iterations is ≥7 (Itmax), Δt for the next time level ismultiplied by a constant <1 (k2, usually between 0.3 and 0.9).
If, during a particular time step, the number of iterations at any time level
becomes greater than a prescribed maximum (Itcrit, usually between 10 and 50),
the iterative process for that time level is terminated. The time step issubsequently reset to Δt/3, and the iterative process restarted.

Parameters Δtinit, Δtmin, Δtmax, Itcrit, Itmin, Itmax, k1, and k2 are specified by a user at input.
The recommended values for these parameters are as follows:149

Parameter
Δtinit

Δtmin

Δtmax

Itcrit

Itmin
Itmax
k1

k2

Recommended
Commentvalue
1 s (15
The recommended value for the initial time step depends onminutes)
the type of simulation and boundary conditions used. Whensimulating a process that starts with a large initial pressure
head or concentration gradient at a particular boundary (e.g.,ponded infiltration or a sudden change in the boundary
concentration), use a small value of the initial time step(e.g., 1 s). When simulating a long-term process withvariable boundary conditions (e.g., seasonal or multiyear
simulation), start with a larger time step (e.g., 15 min). Thisbecause this initial time step is used whenever time-variableboundary conditions significantly change. If needed (if there
is no convergence for Δtinit), the program will still use asmaller time step than Δtinit, but starting with a larger Δtinit
leads to more efficient calculations. In general, smaller
initial time steps must be used for soils with more nonlinear
soil hydraulic properties (e.g., coarse textured soils) andlarger initial time steps can be used for soils with lessnonlinear soil hydraulic properties (e.g., loam)
1s
Always specify a small minimum allowed time step, on theorder of 1 s. This value may never be used, but it provides
the code with the flexibility when it may be needed, e.g.,when there is a sudden change in boundary fluxes and
HYDRUS-1D may not converge with larger time steps.
largeThis is a relatively unimportant parameter, and a large value
may be specified. Since HYDRUS-1D automatically selects
its optimal time step, there is usually no need to constrain
that. The only time when there is a need to constrain thetime step is likely for cases when HYDRUS-1D is asked togenerate internally intra-daily variations in temperature, orin evaporation and transpiration fluxes. Then there is a needto have time steps smaller (e.g., 1 h) so that these daily
variations can be properly modeled.
10
It is usually not helpful to use a larger value than 10. If
HYDRUS-1D does not converge in 10 iterations, it is not
likely that it will do so with more iterations. Even if it does,
it is much more efficient to reduce the time step and attempt
to find the solution with smaller time steps, which is done
automatically by the program when Itcrit is reached.
3
The optimal value in most cases.
7
The optimal value in most cases.
1.3
The optimal value in most cases. Only when there is asaturated zone in the profile, e.g., a perched water layer, thenumerical solution may be more stable with smaller k1 (e.g.,
1.1).
0.7
The optimal value in most cases.

150

Spatial Discretization
The finite element mesh is constructed by dividing the soil profile into linear elements
whose sizes are defined by the z-coordinates of the nodes that form the element sides.
Finite element dimensions must be adjusted to a particular problem:
1. They should be made relatively small at locations where large hydraulicgradients are expected. This often happens close to the soil surface where highlyvariable meteorological factors can cause rapid changes in the soil water content
and corresponding pressure heads. Similarly, regions with sharp gradients canbe located in the vicinity of internal sources or at the interface of two contrastingsoil types. Hence, we recommend using relatively small elements at and near
the soil surface. The size of elements can gradually increase with depth to reflect
the generally much slower changes in pressure heads at deeper depths. To
decrease numerical errors, we also recommend using very similar sizes forneighboring elements. The ratio of the sizes of two neighboring elements is notrecommended to exceed about 1.5.
2. The required size of finite elements close to the soil surface also depends verymuch on how boundary conditions are specified. When boundary conditions are
specified for daily or shorter time intervals, they often involve short-duration
fluxes of a large magnitude. The local spatial discretization needs then to be
finer (on the order of cm or even less) than when boundary conditions are
specified for longer time intervals (e.g., weekly or monthly).
3. The element dimensions should also depend upon the soil hydraulic properties.
For example, coarse-textured soils having relatively high n and α values
generally require a finer discretization than medium-textured soils. That is
because the soil hydraulic functions of coarse-textured soils are more nonlinear
and hence may cause less stable numerical solutions. This is reflected by sharperpressure head fronts for sandy soil profiles compared to medium- and manyfine-textured soil profiles. Hence, we require more FE nodes to be able todescribe the sharp moisture fronts in coarse-textured soils using elements on theorder of 1 cm or less.
One generally cannot have finite elements for variably-saturated flow problems thatare as large as for saturated problems (e.g., when using MODFLOW). While one canhave large elements in the saturated zone (where the numerical solution is often basedon the linear Boussinesq equation), one cannot have such elements in the unsaturated
zone when the solution is based on the nonlinear Richards equation, especially close
to the upper boundary where often sharp infiltration fronts are encountered (such as
for infiltration into dry soils). The HYDRUS-1D model may be used to quickly figureout what discretization one can use for specific soil hydraulic properties and boundary
conditions, and then use this discretization for two- or three-dimensional transportdomains.

151

152

Appendix II. Boundary Conditions
Specifying appropriate boundary conditions (BCs) is one of the most critical tasks whenconstructing a numerical model. A water flow BC is usually represented by a known
value of the water flux or the water pressure head along the outer boundary of the FE
mesh. The external boundary of the selected flow domain provided the interface
between the soil and the outside environment. A solute transport BC is usually
represented by the solute concentration or the solute flux. Solving the governing
equations for saturated/unsaturated flow/transport, which means finding new pressure
heads and/or concentrations at each node of the FE mesh in a time-marching scheme,
requires knowledge of the BCs. Otherwise the problem becomes mathematically not
fully defined. The BCs in a model must simulate real-life conditions and hence must be
selected with extreme care.
In HYDRUS-1D, boundary conditions are categorized as follows:• System-dependent BCs: or dynamic BCs, meaning that they can change during thesimulation (i.e., they depend on the solution at the end of each time step). They may
depend on saturation conditions (as at a seepage face or a tile drain), and on soil
hydraulic properties and/or climate conditions (as at the soil/atmosphere interface).
• System-independent BCs: This type of BCs is entirely known a priori, isimplemented by the user, and is independent of the simulation results. Systemindependent water flow BCs include:• A known pressure head (as in Constant Pressure Head and Time-Variable
Pressure Head BCs)
• A known flux (as in No Flux, Constant Flux, Time-Variable Flux, and Deep
Drainage BCs)
• A known gradient (as in a Free Drainage BC)

153

Water Flow Boundary Conditions
In HYDRUS-1D, various water flow boundary conditions are selected as indicated in
Figure II.1.

Figure II.1. The Water Flow Boundary Conditions dialogue window.
Constant Pressure Head BC
This BC refers to a constant pressure head during the entire simulation. Values of thepressure head (positive or negative) are entered as the initial pressure head for theuppermost node (for the upper BC) and/or the initial pressure head for the lowermost
node (for the lower BC). This BC hence may be applied at the upper and lower
boundaries of the model. Figure II.2 shows a variety of cases where a constant pressure
head boundary condition is applicable:•••

BC-1: A constant pressure head from a tension disc infiltrometer where thepressure head = h1–h2 (negative).
BC-2, BC-3 and BC-4: Constant pressure heads from an overflowing tank, a
Mariotte Bottle assembly, and a stream where the pressure head = h (positive).
BC-5: A constant pressure head from a hanging water column where thepressure head = h (negative).

154

Figure II.2. Various applications for a constant head boundary condition.

155

Constant Flux BCs
This BC refers to a constant flux during the entire simulation. The flux (q) has the units ofLength/time (L/T), similarly as the rainfall rate (cm/day). Fluxes are positive in the direction of
the spatial coordinate (i.e., positive fluxes are upward for vertical flow or from left to right forhorizontal flow, and vice versa, negative downward or from right to left). In the case shown in
Figure II.3, the upper boundary flux is -1 cm/day (a downward influx into the domain) and thelower boundary flux is 1 cm/day (an upward influx into the domain).

Figure II.3. The Constant Boundary Fluxes dialogue window.
Variable Pressure Head BCs
This BC refers to a time-variable pressure head that is entered as a time series. The time series has
as many records as the number of times the pressure head changes at a boundary (one may enteras many records as required). For each record in the time series, the constant pressure head iseffective during the period between the time of the previous record and the current one
To implement this BC, check the ‘Time-Variable Boundary Conditions’ check box in the 'Time
Information' dialog window (Figure II.4) and enter 3 in the edit box ‘Number of Time Variable
Boundary Records’.

156

Figure II.4. The Time Information dialogue window.

Under ‘Time-Variable Boundary Conditions’ (Figure II.5), enter the Time and hTop (an appliedpressure head):i. 0 <Time ≤10 days, h=10 cm
ii. 10 <Time ≤25 days, h=0 cm
iii. 25 <Time ≤50 days, h=-30 cm

Figure II.5. The Time Variable Boundary Conditions dialogue window.

157

Variable Pressure Head/Flux BC
This BC refers to a time-variable flux that is entered as a time series. The data is thenentered in a similar manner as the ‘Variable Head’ boundary condition, but instead of thepressure head (hTop) enter the boundary flux (FluxTop) in the time series (Figure II.6).
An example is irrigation with a time-variable intensity.
To implement this BC, press the ‘Time-Variable Boundary Conditions’ check box in the'Time Information' dialog window and enter 3 in the edit box ‘Number of Time Variable
Boundary Records’.
The variable ‘KodTop’ indicates whether a pressure head BC or a flux BC is used.
KodTop =+1 indicates a pressure head BC and KodTop =-1 indicates a flux BC. In thecase shown below, the BC at the soil surface is:
i.
ii.
iii.

0 <Time ≤ 10 days: the Pressure Head BC with h=10 cm
10 < Time ≤ 25 days: the Flux BC with a flux=-2 cm/day (influx at the surface)
25 < Time ≤ 50 days: The Pressure Head BC with h=-5 cm

Figure II.6. The Time Variable Boundary Conditions dialogue window.

158

Atmospheric Boundary Condition
This boundary condition allows you to incorporate climatic conditions such as rainfall(precipitation), evaporation, and transpiration (root uptake) by plants. The latter is only activatedif root water uptake is checked in the ‘Main Processes’ dialog window as in Figure II.7.

Figure II.7. The Main Processes dialogue window.
With this BC, you can control whether excess precipitation (i.e., when the precipitation rates
exceed the saturated hydraulic conductivity of the soil surface) runs off instantly or ponds at thesoil surface up to a certain height before it starts running off. When ‘Atmospheric BC with Surface
Layer’ is checked, then you must specify the maximum pressure head allowed at the soil surface(Figure II.8). Once this pressure head (a thickness of the surface water layer) is reached water starts
running off.

159

Figure II.8. The Water Flow Boundary Conditions dialogue window.
When an atmospheric flux boundary condition is used and the specified flux into the soil is higher
than the infiltration capacity of the soil, HYDRUS-1D switches the BC from a flux BC to apressure head BC. The pressure head at the boundary becomes zero when ‘Atmospheric BC with
Surface Run off’ is checked and becomes equal to the thickness of the soil surface water layers(up to the ‘Max h at Soil Surface’ thickness) when ‘Atmospheric BC with Surface Layer’ ischecked. The amount of water infiltrating across the soil surface is subsequently calculated
according to the hydraulic conductivity and pressure head gradient at the soil surface (using DarcyBuckingham’s law), while excess water is instantly removed as surface runoff.
The window below appears only when ‘Time-Variable Boundary Conditions’ is checked in the‘Time Information’ window. The data are then entered as a time series. The number of records inthe series depends on the variability of the data, not the length of the simulation. For example, you
may have a 5-day simulation where the precipitation is recorded at an average rate of 1 cm/day
throughout the 5-day period, in which case you will need only one record as shown here (FigureII.9).

160

Figure II.9. The Water Flow Boundary Conditions dialogue window.
However, you may also have rainfall data obtained from a 0.5 mm tipping bucket for a 12-hour
duration. The raw data are shown in columns 1 and 2 of Table II.2. Precipitation in HYDRUS-1Dis entered as a rate (having the units chosen for the simulation, cm/day in this example). Using thetipping-bucket data, we must calculate a rate during each tipping period as follows (refer to Table
II.2):•••••

Column 3: calculate the duration between consecutive tips. Convert ‘min’ to ‘day’ in
Column 4.
Column 5: convert rainfall to cm.
Column 6: calculate cumulative time (simulation time).
Column 7: calculate the rainfall rate (Column 5 / Column 4)
Columns 6 and 7 are ready to be entered as ‘Time Variable Boundary Conditions’ under‘Time’ and ‘Precip.’.

In many cases (especially for very long simulations) the weather data are best entered from a
spreadsheet. The easiest way to enter the data is as follows:••Enter the exact number of Time-Variable Boundary Records (8 in this case).
Open the ‘Time-Variable Boundary Conditions’ dialog window (Figure II.10) and copy
the 8 cells from Columns 6 and 7 (shaded cells in Table II.2). The table should contain
8 rows. Select the first cell in Row 1 under column ‘Time’, and Press “Ctrl_v” on the
keyboard; the data should appear as in the Window shown here.

161

Table II.2. Rainfall data.
Time
5:36:00
5:42:00
5:48:00
6:00:00
6:12:00
6:36:00
8:18:00
14:00:00
17:36:00

Rainfall(mm)

Duration(min)

Duration(day)

Rainfall(cm)

Time (day)Cumulative

Rate(cm/day)

1
0.5
0.5
0.5
0.5
0.5
0.5
0.5

6
6
12
12
24
102
342
216

0.004170.004170.00833
0.00833
0.016670.07083
0.237500.15000

0.1
0.05
0.05
0.05
0.05
0.05
0.05
0.05

0.00420.00830.01670.02500.04170.11250.35000.500024126630.70590.21050.3333Figure II.10. The Time Variable Boundary Conditions dialogue window.

Triggered Irrigation
When this BC is selected, you can model ‘Triggered Irrigation’ whereby a constant flux (irrigation)is introduced to the soil surface when a certain pressure head is reached at a certain node (anObservation Node). To activate this option, check the ‘Triggered Irrigation’ box in the ‘Water
Flow Boundary Conditions’ dialog window as shown below (Figure II.11).

162

Figure II.11. The Water Flow Boundary Conditions dialogue window.
The required parameters are (Figure II.12):•••

The number of the ‘Observation point’ triggering irrigation (the observation point shouldbe specified in the Graphical Editor).
The pressure head threshold at the observation point when irrigation is triggered.
Irrigation rate, irrigation duration, and the lag time (the time interval between times when
the pressure head triggering irrigation is reached and when irrigation starts).

Figure II.12. The Triggered Irrigation dialogue window.
Free Drainage
This boundary condition specifies a unit total gradient at the lower boundary (outflow, drainage)of the soil profile. It is applicable for cases when the water table is located far below the domain
of interest and thus does not affect flow in the transport domain. This BC assumes a unit total
163

vertical hydraulic gradient, that is, gravity flow with the pressure head gradient being zero. Thisboundary condition should be used only at the bottom of the domain.

Deep Drainage
This is a time-variable flux BC based on an observed experimental relationship between thepressure head at the bottom of the soil profile (usually positive) and discharge. Hence, it can only
be used when such information is available. This BC relates flow at the bottom of the soil profile
to the position of the groundwater table, which must be above the bottom of the soil profile (i.e.,
within the transport domain). The BC is based on an expression proposed by Ernst and Feddes(1979) (as reported by Hopmans and Stricker, 1989):

q = aeb h − ho(II.1)

where q is the discharge rate (or recharge rates), a and b are fitting parameters (obtained fromexperimental observations), h is the position of the groundwater level above the bottom of the soil
profile, and h0 is some equilibrium level of the groundwater table. The parameters a and b accountfor the effects of soil hydraulic properties (of the deeper layers), regional flow, and the position of
the groundwater table on discharge.
The nodal flux at each time during the simulation is calculated using Equation II.1. This flux
depends upon the pressure head, which varies during the simulation. Three values are required(Figure II.13):

Figure II.13. The Deep Drainage Boundary Conditions dialogue window.
••

The parameters ‘a’ and ‘b’ (of Equation II.1).
The reference position of the groundwater table, which is like a datum value.

164

The calculated flux is assigned to the bottom node at each time step. Applications of this BC arelimited because it requires field data to evaluate the fitting parameters. A recent example is given
by Neto et al., 2016).
Seepage Face
This is a dynamic outflow (or drainage) BC that changes according to the flow conditions during
a simulation. This BC is used at the bottom of lysimeters, laboratory soil columns, or tile drains
where water freely drains through a soil surface exposed to the outer atmosphere. It only allows
flow when saturation is attained, the pressure head at which saturation is achieved can be controlled(the default value is zero, Figure II.14).

Figure II.14. The Water Flow Boundary Conditions dialogue window.

Solute Transport Boundary Conditions
Available solute transport boundary conditions include (Figure II.15):• Concentration boundary conditions whereby users specify the liquid phase concentration
at the boundary (for the upper and lower boundaries)• Concentration flux boundary conditions whereby users specify the concentration of theinfiltrating water (for upper and lower boundaries)• A stagnant boundary layer for volatile solutes (applicable to the upper boundary)• Zero concentration gradient (free drainage) (for the lower boundary)165

Third-type (Cauchy or solute flux) boundary conditions prescribe solute fluxes across (notconcentrations at) a boundary. Due to mixing with water (and solute) initially present in the
profile, one will not get immediately the input (C0) concentration on the boundary. Since one
prescribes a solute flux, one has complete control over the mass balance and how much soluteenters into the transport domain.
When a first-type (Dirichlet or concentration) boundary condition is used, one prescribes theconcentration at the boundary (not the flux into the domain). Since solute flux consists of two
components, i.e., advective and dispersive fluxes, both of these are active on the boundary at thesame time. Initially, there is a large concentration gradient, and thus the dispersive solute flux
can be very large. In this case, one has a much larger flux into the domain than if one uses thethird-type solute flux boundary condition.
We recommend using always a third-type boundary condition since that condition is physicallythe most realistic. Dirichlet conditions in most cases do not provide a conserve mass and hence
do not provide a good alternative. You should use a Dirichlet condition only if, for example, alarge reservoir of contaminant is in contact with the transport domain, in which case one canassume that the boundary concentration is fixed (e.g., for upward flow from a water table withgroundwater having a certain concentration). Several have reviewed the implications of using
different boundary conditions for various experimental conditions (e.g., van Genuchten andParker, 1984; Skaggs and Leij, 2002)

Figure II.15. The Water Flow Boundary Conditions dialogue window.

166

Appendix III: Unsaturated Soil Hydraulic Properties
Water Retention Properties
The fundamental relationship between the soil water content (usually calculated on a volumetricbasis) and the pressure head (or the matric suction) is referred to as the soil-water retention curve(SWRC) or the soil water characteristic curve. The matric suction is the difference between thepore air pressure head and the pore water pressure head. If a change in the water content occurs at
any point in a mass of soil initially in an equilibrium state, the soil matric suction (or the negative
pore water pressure) at that point changes (assuming atmospheric conditions, the pore air pressure
has a value of 0). Following this change, there will be a difference in the energy level between the
point where the change occurred and the surrounding soil, which will cause water to flow into thesoil until new equilibrium conditions are established (Croney, 1952). Knowledge of the water
retention curve is essential when investigating how soil moisture reacts in response to changes in
the pressure head, resulting from events such as infiltration, evaporation, and drainage. Due to thehysteretic behavior of soils with respect to wetting and drying processes, different SWRC curves
may be obtained depending upon whether the soil is undergoing a drying cycle, a wetting cycle,
or a re-drying cycle (Collis-George, 1955).
Shape of the Soil-Water Retention Curve
The general shape of the SWRC under wetting and drying conditions is shown in Figure III.1 ona semi-log scale. The following parameters are identified:••••

The saturated volumetric water content, θs, which theoretically corresponds to the soil’s
porosity, but practically is often about 10-25% less because of dissolved and/or entrapped air.
The air-entry value, AEV, is a critical value of suction, at which the largest pores in the soil
matrix start losing water.
The residual water content, θr, is the water content beyond which a further increase in the soil’s
suction results in only marginal changes in the water content.
At very high suctions (negative pressure heads), the water content eventually becomes zero.
This value is about 106 kPa (107 cm) and has been supported experimentally (e.g., Croney andColeman, 1961; Koorevaar et al., 1983), as well as by thermodynamic considerations(Richards, 1965).

167

0.5

Volumetric water content

θs

AEV

Drying curve

0.4

Wetting curve

0.3

0.2

0.1

θr

0
0.1

1

10

100

1000

10000

1000001000000

Matric suction (kPa)

Figure III.1. The typical shape of a soil water retention curve.
The residual water content (θr) and the corresponding suction, sometimes termed the residualsuction, have been defined in a number of ways in the literature. However, the accepted use of
these terms is not always agreed upon. Tinjum et al. (1997) defined it as the water content
corresponding to the asymptote of the SWRC at low degrees of saturation. van Genuchten (1980)
defined θr as the water content at which its rate of change with respect to suction becomes zero,or when the liquid flow of water ceases (excluding the region near saturation).
Classically, moisture in the unsaturated zone is separated into three components: gravitational,
pellicular, and hygroscopic water (e.g., Everett, 1993). Gravitational water is the moisture in a soilthat can be drained by gravitational forces. Pellicular water is moisture in a soil that cannot be
drained by gravity forces but can be lost by evaporation. Finally, the hygroscopic water is moisturethat will never be lost through the above natural forces. The residual water content may be viewed
as the water content at the limit of liquid water extraction. De-saturation beyond residual
conditions occurs primarily as a result of vapor flow up to the point where the soil water content
is in equilibrium with the vapor pressure of its surrounding. Brooks and Corey (1966) pointed outthat it is difficult to determine the residual saturation of clayey materials; they related the residual
saturation to the clay content of the soil. The well-defined, residual state for sands, silts, and their
mixtures can reasonably well be predicted from the SWRC as shown in Figure III.1.

Models Describing the SWRC
A mathematical representation of the SWRC is essential when using numerical models for
simulating fluid flow and mass transfer in the unsaturated zone. HYDRUS-1D allows users to
select from five analytical models to describe the soil hydraulic properties: van Genuchten (1980),
modified van Genuchten (Vogel and Cislerova, 1988), Brooks and Corey (1964), Kosugi (1996),
and Durner (1994). The first three of these models are described briefly in the following section.
168van Genuchten (VG) Modelvan Genuchten (1980) proposed a mathematical representation of the SWRC, which depicts an Sshaped curve when plotted on a semi-log scale. His equation is given by:

S=
e






n −m

1 + α h 

(III.1)

where α, m, and n are mostly semi-empirical fitting parameters (usually m=1-1/n), h is thepressure head, and Se is the normalized volumetric water content (often referred to as effective
saturation):S e ( h) =θ ( h) − θ r
θs − θr(III.2)in which θ is the volumetric water content at pressure head h, and θs and θr are the saturated andresidual water contents, respectively.
Each of the three fitting parameters (α, m, and n) control part of the S-shape curve. Figure III.2
shows the effect of α (Alpha), which is closely related to the reciprocal of the air-entry value(AEV) of a soil. Coarse-grained soils have a low AEV and a high α value, whereas fine-textured
soils have a lower α.
0.35Volumetric water content

0.3

0.25

0.2

0.15

Alpha = 0.015
Alpha = 0.05

0.1

Alpha = 0.1

0.05
1

10

100

1000

Suction (cm)Figure III.2. The effect of α in the van Genuchten model on the SWRC.
Figure III.3 shows the effect of the n parameter, which controls the slope of the curve and reflects
the width of the particle size distribution. The plot indicates that a unit increase in suction at ornear the steepest part of the curve causes more water to be extracted from coarse-textured soils(high n) than from fine-textured soils (low n), but not at other places of the curve.

169

0.35

Volumetric water content

0.3
n=1.5
0.25

n=3
n=5

0.2

0.15

0.1

0.05
10

100

1000

Suction (cm)Figure III.3. The effect of n in the van Genuchten model on the SWRC.
When the parameters n and m are correlated (as when m=1-1/n), a relatively simple equation canbe derived for the hydraulic conductivity function, K(h), when the retention function is combinedwith the theoretical pore-size distribution model of Mualem (1976) for K (see below). Figure III.4
shows that when the two parameters are correlated, the curvatures (A and B) of the S-shape curve
become correlated. When individual values are assigned to n and m the curvatures in the wet anddry regions of the retention curve become independent (van Genuchten and Nielsen, 1985).
0.35
m = 1-1/n = 0.66

A

Volumetric water content

0.3

m = 0.3

Curvatures A &B are dissimilar

Alpha = 0.015
n=3

0.25

0.2

0.15

Curvatures A & B are similar(correlated; m=1-1/n)

B

0.1

0.05
10

100

1000
Matric suction (cm)Figure III.4. The relation of m and n in the van Genuchten model.

170

10000

van Genuchten Model with AEV of -2 cm
The unsaturated soil hydraulic properties are often described using the Mualem–van Genuchten(MVG) type analytical functions (van Genuchten, 1980). Recent studies suggest several
shortcomings of these functions near saturation, notably the lack of the second-order continuity ofthe soil water retention function at saturation and the inability of the hydraulic conductivity
function to account for macroporosity. Vogel et al. (1985) and Schaap and van Genuchten (2006)
showed that a modified MVG formulation with a small but constant air-entry pressure in the water
retention curve substantially improves the description of the hydraulic conductivity nearsaturation. It is thus recommended, especially for fine-textured soils with small values of the nparameter (n < 1.1-1.2), to use a modified van Genuchten-Mualem model with an air-entry valueof -2 cm. In addition to improving the description of the hydraulic conductivity function, thismodified MVG model also dramatically improves the convergence of the numerical solution.
The van Genuchten (VG) model that implements a small non-zero air-entry value (AEV) of -2 cm
shows almost no effect on the macroscopic description of the water retention function of finetextured soils compared to the original VG formulation. Selecting this option will introduce a small
correction (like an air-entry value) in the water retention function to force the slope of the retentionfunction (dθ/dh) to become zero when approaching saturation at h=0 (see Figure III.5; note thatthe water content remains at saturation for suctions below 2 cm).
0.45
Modified VG; AEV=-2 cm
VG

Volumetric water content

Volumetric water content

0.35

0.25

0.15

0.43

0.42
Modified VG
VG

0.41
1

10

Suction (cm)

0.05
1

10

100

1000

Suction (cm)Figure III.5. The van Genuchten model with AEV of –2 cm.

We will show later that implementing this low non-zero AEV has a significant effect on thehydraulic conductivity function. This option is recommended for fine-textured soils (e.g., silty
clays and clays), but should not be used for medium-, and especially coarse-textured, soils whenthe value of the parameter α becomes relatively large (say >0.02 or 0.05 cm). This rule applies
171

unless n is again less than about 1.1 or 1.2, which generally is not the case for medium- and coarsetextured soils. The -2 cm AEV should then be avoided since it introduces a Brooks and Corey(1964) type air-entry correction into the retention function to make the function less accurate formany soils (Schaap and van Genuchten, 2006).

Modified van Genuchten Model (Vogel and Cislerova)
Vogel and Cislerova (1988) modified the van Genuchten model by incorporating a non-zero AEV
into the model. The modification was implemented by introducing a fictitious water content θmthat is higher than θs and replaces θs in the van Genuchten model. This fictitious water content isused only when h<AEV. Above AEV the water content is equal to θs. Their retention model isgiven byθm − θ aθ a +n m
1+ α h
θ =

θs


(

)

h < AEV(III.3)h ≥ AEV

Figure III.6 shows how the model is implemented. The data points at suctions <17 cm are only ahypothetical continuation to the S-shape curve. That is, when the pressure head is below thisvalue, the water content is fixed at θs.
0.5
Theta_m = 0.48
Theta_s = 0.43

Volumetric water content

0.4

0.3

Vogel & Cislerova

0.2

van Genuchten

0.1

AEV (hs) =17 cm

0
1

10

100
Suction (cm)Figure III.6. The modified van Genuchten model.

172

1000

Brooks and Corey (BC) Model
Brooks and Corey (1964, 1966) concluded from comparisons of a large number of experimental
data that the SWRC could be described by the following formula:
given by(

Se = AEVh)n(III.4)

where AEV is the air entry value (corresponds to 1/α in the HYDRUS-1D GUI), and n is a soil
characteristic parameter, called the pore-size-distribution index. Experimental data of suction andwater content are often plotted as logh versus log θ. Linear regression can then be performed on
the straight-line relationship. The slope and intercept of the best-fit line then correspond to 1/n and
AEV, respectively. Figure III.7 shows a comparison of the BC model with the VG model using theparameters listed in Table III.1.
When using the Brooks and Corey functions in HYDRUS-1D, we recommended to either disable
the internal interpolation tables or to set the lower limit above the air-entry value.

Table III.1. Parameters for two SWRCs in Figure III.7.
van Genuchten
Brooks & Corey

θr

θs

α

0.078
0.027

0.43
0.434

0.036
0.0897

n
1.56
0.293

0.45

Volumetric water content

0.35

0.25

Brooks and Corey

0.15

van Genuchten

0.05
1

10

100

1000

Suction (cm)Figure III.7. Comparison of the Brooks and Corey and van Genuchten models.

173

Hydraulic Conductivity Functions
van Genuchten-Mualem Model
Mualem (1976) developed the following equation for predicting the relative hydraulic conductivity
Kr from the knowledge of the SWRC:
K r = Se lSe

1

∫0 h (τ )dτ ∫0 h (τ )dτ−1−1(III.5)where h is the pressure head, given as a function of the dimensionless water content Se defined by
Equation (III.2).
For the special case where m=1-1/n (the fitting parameters of Equation (III.1)), a closed-form
solution was obtained by van Genuchten (1980) as follows:mK=(h) K s S 1 − (1 − Se1/ m ) 


l
e

2

(III.6)

where l is a pore connectivity parameter established to be 0.5 by Mualem (1976) and -1 by Schaapand Leij (2000). Equation (III.6) in combination with the VG retention model has been usedsuccessfully in an extremely large number of studies. They are the default functions in HYDRUS1D.
van Genuchten Model with AEV of -2 cm
Choosing the approach by Vogel and Cislerova (1988) for the retention function substantiallychanges the shape of the predicted hydraulic conduction function, K(h), based on Equation (III.5),
as compared to the curve predicted using the standard VGM statistical pore-size distribution model(Vogel et al., 2001). The predicted K(h) curve is then far less nonlinear and much more accurate
for fine-textured soils (typically clays and silty clays). Using this option may also result in fewer
numerical problems when modelling ponded infiltration or other variably-saturated flow problems.
The -2 cm correction was first suggested by Vogel and Cislerová (1988) and worked out further
by Vogel et al. (2001). Schaap et al. (2006) later showed that a -4 cm correction would be slightly
better than -2 cm, but the difference is minimal.
Figure III.8 shows the impact of using this option. Note that the hydraulic conductivity functions
are usually plotted on the log-log scale, but in this case we chose to plot it on a semi-log scale toemphasize the dramatic differences near saturation.

174

25

Hydraulic conductivity (cm/day)van Genuchten
van Genuchten with AEV=-2 cm

20

Vogel & Cislerova

15

10

5
LOAM
0
1

10

100

1000

Suction (cm)Figure III.8. The effect of the non-zero AEV on the hydraulic conductivity function.

Modified van Genuchten Model (Vogel and Cislerova)
When using this model, the hydraulic conductivity function is further discretized into three zones
as shown in Figure III.9:h >= hs (AEV)hk < h < hs
h <= hk

K(h) = Ks
K(h) is linear value between Ks and Kk
K(h) is non-linear (Mualem model)

where hk is the head corresponding to a user-defined water content θk (Qk in the HYDRUS GUI),θk is the water content below saturation where a measurement of the hydraulic conductivity Kk isavailable. The mathematical formulation for this model is found in the online Help and will not be
repeated here. A comparison between the results of this model and the original van Genuchten
model is shown in Figure III.10. Figure III.8 shows the significant differences near saturation. The
modification was made to enable relatively large increases in K near saturation that are often
observed in structured (macroporous) soils.

175

0.5

25

Conductivity versus head

20

0.4

kk = 15
15

Qk = 0.38

10
0.3

Volumetric water content

Hydraulic conductivity (cm/day

Water content versus head5hs ; AEV

hk

0
0

25

50

75

100

125

0.2
150

Head (-cm)

Figure III.9. The Vogel and Cislerova hydraulic conductivity function.

Hydraulic conductivity (cm/day

1.0E+01

1.0E-01

1.0E-03

van Genuchten
Vogel & Cislerova

1.0E-05

1.0E-07

1.0E-09
1

10

100

1000

10000

Suction (cm)Figure III.10. The van Genuchten and Vogel and Cislerova hydraulic conductivity functions.

Brooks and Corey Model
For heads <AEV, the hydraulic conductivity is defined as follows:K ( h) = K s

( AEV
h )
176

2+3n

(III.7)

A comparison of the results obtained with this model and those obtained with the van Genuchten
model is shown in Figure III.11. The hydraulic parameters listed in Table III.1 were used for thispurpose.

1.0E+03
van Genuchten
Hydraulic conductivity (cm/day)

1.0E+01

Brooks & Corey

1.0E-01

1.0E-03

1.0E-05

1.0E-07

1.0E-09
1

10

100

1000

10000

Suction (cm)Figure III.11. The van Genuchten and Brooks and Corey hydraulic conductivity functions.

177

178

References
Brooks, R. H., and A. T. Corey, Hydraulic properties of porous media, Hydrol. Paper No. 3,
Colorado State Univ., Fort Collins, CO, 1964.
Brooks, R.H. and Corey, A.T. Properties of porous media affecting fluid flow. J. Irrig. Drain. Div.,ASCE, 92(IR2): 61-88, 1966.
Carsel, R.F., and R. S., Parrish, Developing joint probability distributions of soil water retention
characteristics, Water Resour. Res, 24, 755-769, 1988.
Collis-George, N. Hysteresis in moisture content-suction relationships in soils. Proc. Nat.
Academy of Science (India), 24a: 80-85, 1955.
Constantz, J., Temperature dependence of unsaturated hydraulic conductivity of two soils. Soil
Sci. Soc. Am. J., 46(3), 466-470, 1982.
Croney, D. The movement and distribution of water in soils. Geotechnique, 3(1): 1-16, 1952.
Croney, D., and Coleman, J.D. Pore pressure and suction in soils. Proc. Conf. on Pore Pressureand Suction in Soils. Butterworths, London: 31-37, 1961.
Dash, Ch.J., A. Sarangi, D.K. Singh, A.K. Singh, and P.P. Adhikari, Prediction of root zone waterand nitrogen balance in an irrigated rice field using a simulation model, Paddy Water Environ.,13, 281–290, doi:10.1007/s10333-014-0439-x, 2015.
Durner, W., Hydraulic conductivity estimation for soils with heterogeneous pore structure, Water
Resour. Res., 32(9), 211-223, 1994.
Ebrahimian, H., A. Liaghat, M. Parsinejad, E. Playan, F. Abbasi, and M. Navabian, Simulation of
1D surface and 2D subsurface water flow and nitrate transport in alternate and conventional
furrow fertigation, Irrig. Sci., 31, 301–316, doi: 10.1007/s00271-011-0303-3, 2013a.
Ebrahimian, H., A. Liaghat, M. Parsinejad, E. Playan, F. Abbasi, M. Navabian, and B. Lattore.
Optimum design of alternate and conventional furrow fertigation to minimize nitrate loss, J.
Irrig. Drain. Eng., 139, 911–921, doi: 10.1061/(ASCE)IR.1943-4774.0000635, 2013b.
Eching, S. O., and J. W. Hopmans, Optimization of hydraulic functions from transient outflow andsoil water pressure data, Soil Sci. Soc. Am. J. 57, 1167-1175, 1993.
EFSA, European Food Safety Authority. Scientific Opinion of the Panel on Plant Protection
Products and their Residues on a request from EFSA related to the default Q10 value used todescribe the temperature effect on transformation rates of pesticides in soil. The EFSA Journal(2007) 622, 1-32, 2008.
Everett, EG. Vadose zone monitoring. In Geotechnical Practice for Waste Management, Ed. D.E.
Daniel, Chapman and Hall, London: 651-675, 1993.
Feddes, R. A., P. J. Kowalik, and H. Zaradny, Simulation of Field Water Use and Crop Yield, John
Wiley & Sons, New York, NY, 1978.
Flury, M., and T. F. Gimmi, Solute Diffusion. Methods of Soil Analysis: Part 4 Physical Methods.
SSSA Book Series, 5.4:1323-1351. doi:10.2136/sssabookser5.4.c55, 2002.
Freeze, R.A., and J.A., Cherry, Groundwater: Englewood Cliffs, NJ, Prentice-Hall, 604 p., 1979.
Green, T. W., Z. Paydar, H. P. Cresswell, and R. J. Drinkwater, Laboratory Outflow Techniquefor Measurement of Soil Water Diffusivity and Hydraulic Conductivity, Technical Report No.
12/98, CSIRO, 30 pp., 1998.
Hopmans, J. W., and J. N. M. Stricker, Stochastic analysis of soil water regime in a watershed, J.
Hydrol., 105, 57-84, 1989.

179

Hopmans, J. W., J. Šimůnek, N. Romano, and W. Durner, Inverse Modeling of Transient Water
Flow, In: Methods of Soil Analysis, Part 1, Physical Methods, Chapter 3.6.2, Eds. J. H. Daneand G. C. Topp, Third edition, SSSA, Madison, WI, 963-1008, 2002.
Jacques, D., J. Šimůnek, User Manual of the Multicomponent Variably-Saturated Flow andTransport Model HP1, SCK•CEN-BLG-998, 2005.
Jacques, D., J. Šimůnek, D. Mallants, M. Th. van Genuchten, Operator-splitting errors in coupled
reactive transport codes for transient variably saturated flow and contaminant transport in layer
soil profiles, J. Contam. Hydrol., 88, 197–218, 2006.
Jacques, D., J. Šimůnek, D. Mallants and M. Th. van Genuchten, The HPx code for coupled
reactive transport under variable saturated flow conditions: Recent developments andapplications, Journal of Hydrology and Hydromechanics, 66(2), 211-226, doi: 10.1515/johh2017-0049, 2018.
Jensen, D. T., G. H. Hargreaves, B. Temesgen, and R. G. Allen, Computation of Eto under nonideal
conditions, J. Irrig. Drainage, 123(5), 394-400, 1997.
Koorevaar, P., Menelik, G., and Dirksen, C. Elements of soil physics, Elsevier Science
Publications, Amsterdam, The Netherlands, 1983.
Kosugi, K., Lognormal distribution model for unsaturated soil hydraulic properties, Water Resour.
Res., 32(9), 2697-2703, 1996.
Kumar, R., V. Shankar, and M. K. Jat, Evaluation of root water uptake models – a review, ISH
Journal of Hydraulic Engineering, 1-10, doi: 10.1080/09715010.2014.981955, 2014
Lenhard, R. J., and J. C. Parker, Modeling multiphase fluid hysteresis and comparing results tolaboratory investigations, In, Proc. Intl. Workshop on Indirect Methods for Estimating theHydraulic Properties of Unsaturated Soils, edited by M. Th van Genuchten, F. J. Leij, and L.
J. Lund, University of California, Riverside, pp. 233-248, 1992.
Li Y.-H. and S. Gregory, Diffusion of ions in sea water and in deep-sea sediments, Geochemica et
Cosmochimica Acta, 38, 703-714, 1974.
Luckner, L., M. Th. van Genuchten, and D. R. Nielsen. A consistent set of parametric models forthe two-phase flow of immiscible fluids in the subsurface. Water Resour. Res., 25(10):2187-2193, 1989.
Maas, E. V., Crop salt tolerance, In "Agricultural Salinity Assessment and Management", K. K.
Tanji (ed.), ASCE Manuals and Reports on Engineering practice, No 71, New York, 1990.
Mallants, D., J. Šimůnek, M. Th. van Genuchten, and D. Jacques, Using Hydrus and its modules
to simulate the fate and transport of coal seam gas chemicals in variably-saturated soils, Water,
9(6), 385, doi: 10.3390/w9060385, 2017.
Mallants, D., M. Vanclooster, J. Feyen, A transect study on solute transport in a macroporous soil.
Hydrol Process 10, 55–70, 1996a.
Mallants, D., D. Jacques, M. Vanclooster, J. Diels, and J. Feyen, A stochastic approach to simulate
water flow in a macroporous soil, Geoderma, 70, 299–324, 1996b.
Mallants, D., M. Th. van Genuchten, J. Šimůnek, D. Jacques, and S. Seetharam, Leaching ofcontaminants to groundwater – Chapter 18, In (F.Swartjens, Editor) Dealing with
Contaminated Sites. From Theory towards Practical Application, ISBN 978-90-481-9756-9,
Springer Verlag, Dordrecht, The Netherlands, 2011, pp. 787-850. doi: 10.1007/978-90-4819757-6_1. 2011.
Marrero, T.R., and E.A. Mason, Gaseous diffusion coefficients, J. Phys. Chem. Ref. Data, 1(1), 3118, 1972.

180

Millington, R. J., and J. M. Quirk, Permeability of porous solids, Trans. Faraday Soc., 57, 12001207, 1961.
Mualem, Y. A new model for predicting the hydraulic conductivity of unsaturated porous media.
Water Resources Engineering, 12: 513-522, 1976.
Perrochet, P., and D. Berod, Stability of the standard Crank-Nicolson-Galerkin scheme applied tothe diffusion-convection equation: some new insights, Water Resour. Res., 29(9), 3291-3297,
1993.
Moldrup, P., T. Olesen, J. Gamst, P. Schjønning, T. Yamaguchi, and D. E. Rolston, Predicting thegas diffusion coefficient in repacked soil: water-induced linear reduction model, Soil Sci. Soc.
Am. J., 64, 1588-1594, 2000.
Monteith, J. L., Evaporation and surface temperature, Quarterly J. Royal Meteo. Soc., 107, 1-27,
1981.
Monteith, J. L. and M. H. Unsworth, Principles of Environmental Physics, Edward Arnold, London,
1990.
Morales I., J. Cooper, J. A. Amador, and T. B. Boving, Modeling nitrogen losses in conventional
and advanced soil-based onsite wastewater treatment systems under current and changing
climate conditions, PLoS ONE, 11(6), e0158292, doi: 10.1371/journal.pone.0158292, 2016.
Mualem, Y., A new model for predicting the hydraulic conductivity of unsaturated porous media,
Water Resour. Res., 12(3), 513-522, 1976.
Neto, D. C., H. K. Chang, and M. Th. van Genuchten, A mathematical view of water tablefluctuations in a shallow aquifer in Brazil, Ground Water, doi: 10.1111/gwat.12329, 2016.
Passioura, J. B., Determining soil water diffusivities from one-step outflow experiments, Aust. J.
Soil Res., 15, 1-8, 1976.
Philip, J. R., and D. A. de Vries, Moisture movement in porous media under temperature gradients,
Eos Trans. AGU, 38(2), 222-232, 1957.
Radcliffe, D., and J. Šimůnek, Soil physics with HYDRUS: Modeling and Applications, CRC Press,
Taylor & Francis Group, Boca Raton, FL, 2010.
Schaap, M. G., & Leij, F. J. Improved prediction of unsaturated hydraulic conductivity with theMualem-van Genuchten model. Soil Science Society of America Journal, 64(3), 843-851,
2000.
Schaap, M. G., Leij, F. J., and van Genuchten, M. Th., Rosetta: a computer program for estimating
soil hydraulic parameters with hierarchical pedotransfer functions, Journal of Hydrology, 251,
163-176, 2001.
Schaap, M. G., and M. Th. van Genuchten. A modified Mualem-van Genuchten formulation forimproved description of the hydraulic conductivity near saturation, Vadose Zone J., 5, 27-34,
2006.
Sierra, C.A., M. Müller, and S. E.Trumbore, Models of soil organic matter decomposition: The
SoilR package, version 1.0, Geosci. Model Dev., 5, 1045–1060, 2012.
Šimůnek, J., and J. W. Hopmans, Parameter optimization and nonlinear fitting, In: Methods of Soil
Analysis, Part 1, Physical Methods, Chapter 1.7, Eds. J. H. Dane and G. C. Topp, Third
edition, SSSA, Madison, WI, 139-157, 2002.
Šimůnek, J., D. Jacques, J. W. Hopmans, M. Inoue, M. Flury, and M. Th. van Genuchten, Solute
Transport During Variably-Saturated Flow - Inverse Methods, In: Methods of Soil Analysis,
Part 1, Physical Methods, Chapter 6.6, Eds. J. H. Dane and G. C. Topp, Third edition, SSSA,
Madison, WI, 1435-1449, 2002.

181

Šimůnek, J., M. Sejna, M. Sakai, H. Saito, and M. Th. van Genuchten, The HYDRUS-1D software
package for simulating the one-dimensional movement of water, heat, and multiple solutes invariably-saturated media, Version 4.0x, HYDRUS-1D Series 3, Department of Environmental
Sciences, University of California Riverside, Riverside, CA, USA, 2008.
Šimůnek, J., M. Th. van Genuchten, and M. Šejna, Development and applications of the
HYDRUS-1D and STANMOD software packages, and related codes, Vadose Zone J., doi:
10.2136/VZJ2007.0077, Special Issue “Vadose Zone Modeling”, 7(2), 587-600, 2008.
Šimůnek, J. and J. W. Hopmans, Modeling compensated root water and nutrient uptake, EcologicalModeling, doi: 10.1016/j.ecolmodel.2008.11.004, 2008.
Šimůnek, J., and M. Th. van Genuchten, Modeling nonequilibrium flow and transport processes using
HYDRUS-1D, Vadose Zone J., 7, 782–797, 2008.
Šimůnek, J., M. Th. van Genuchten, and M. Šejna, Recent developments and applications of the
HYDRUS-1D computer software packages, Vadose Zone Journal, 15(7), pp. 25, doi:
10.2136/vzj2016.04.0033, 2016.
Skaggs, T. H., and F. J. Leij, Solute transport: Theoretical background, In: Methods of Soil
Analysis, Part 1, Physical Methods, Chapter 1.7, Eds. J. H. Dane and G. C. Topp, Third
edition, SSSA, Madison, WI, 1353-1380, 2002.
Søgaard-Hansen J. and A. Damkjaer, Determining 222Rn diffusion lengths in soils and sediments,
Health Physics, 53(5), 455-459, 1987.
Tinjum, J.M., Benson, C.H., and Blotz, L.R. “Soil-water characteristic curves for compacted
clays.” Journal of Geotechnical and Geoenvironmental Engineering, ASCE, 123(11): 10601069, 1997.
Twarakavi, N. K. C, M. Sakai, and J. Šimůnek, An objective analysis of the dynamic nature offield capacity, Water Resources Research, 45, W10410, doi:10.1029/2009WR007944, 9 pp.,2009.
van Dam, J. C., J. N. M. Stricker, and P. Droogers, Inverse method for determining soil hydraulicfunctions from one-step outflow experiment, Soil Sci. Soc. Am. J., 56, 1042-1050, 1992.
van Dam, J. C., J. N. M. Stricker, and P. Droogers, Inverse method to determine soil hydraulicfunctions from multistep outflow experiment, Soil Sci. Soc. Am. J., 58, 647-652, 1994.
Vanderborght, J., and H. Vereecken, Review of dispersivities for transport modeling in soils,
Vadose Zone J., 6, 29–52, 2007.
van Genuchten, M. Th., and P. J. Wierenga, Mass transfer studies in sorbing porous media: I.
Analytical solutions. Soil Sci. Soc. Am. J., 40(4): 473-480, 1976.
van Genuchten M. Th., and D. R. Nielsen, On describing and predicting the hydraulic propertiesof unsaturated soils, Annales Geophysicae, 3(5): 615-628, 1985.
van Genuchten, M. Th., and P. J. Wierenga, Solute dispersion coefficients and retardation factors,
In: A. Klute (ed.), Methods of Soil Analysis, part 1, Physical and Mineralogical Methods,
Agronomy, 9(1), 1025-1054. 2nd ed., Am. Soc. Agron., Madison, WI, 1986.
van Genuchten, M. Th., A closed-form equation for predicting the hydraulic conductivity ofunsaturated soils, Soil Sci. Soc. Am. J., 44, 892-898, 1980.
van Genuchten, M. Th., A numerical model for water and solute movement in and below the rootzone. Unpublished Report, U.S. Salinity Laboratory, USDA, ARS, Riverside, California,
1987.
van Genuchten, M. Th., and J. C. Parker, Boundary conditions for displacement experiments
through short laboratory soil columns, Soil Sci. Soc. Am. J., 48(4): 703-708, 1984.

182

van Genuchten, M. Th., F. J. Leij, and S. R. Yates. The RETC code for quantifying the hydraulicfunctions of unsaturated soils. Report No. EPA/600/2-91/065. R. S. Kerr Environmental
Research Laboratory, U. S. Environmental Protection Agency, Ada, OK. 85 p., 1991
Vogel, T., M. Cislerova, and M. Šir, Reliability of indirect estimation of soil hydraulicconductivity, Vodohospodarsky Casopis, 33(2), 204–224, 1985.
Vogel, T., and M. Císlerová, On the reliability of unsaturated hydraulic conductivity calculated
from the moisture retention curve, Transport in Porous Media, 3, 1-15, 1988.
Vogel, T., K. Huang, R. Zhang, and M. Th. van Genuchten, The HYDRUS code for simulating
one-dimensional water flow, solute transport, and heat movement in variably-saturated media,
Version 5.0, Research Report No 140, U.S. Salinity Laboratory, USDA, ARS, Riverside, CA,
1996.
Vogel, T., van Genuchten, M. Th., and Cislerova. Effect of the shape of the hydraulic functions
near saturation on variably-saturated flow predictions. Advances in Water Resources, (24),133-144, 2001.
Walker, A., A simulation model for prediction of herbicide persistence, J. Environ. Quality, 3(4),396-401, 1974.
Wang, L., E. Martens, D. Jacques, P. De Cannière, J.A. Berry, and D. Mallants, Review of sorption
values for the cementitious near field of a near-surface radioactive waste disposal facility,
NIROND-TR 2008-23 E, NIRAS-ONDRAF, Belgium, 2009.

183

CONTACT US

FOR FURTHER INFORMATION

t 1300 363 400+61 3 9545 2176
e csiroenquiries@csiro.au
w www.csiro.au

David Rassam
CSIRO Land and Water
Dutton Park, Brisbanet +61 7 3833 5586
e david.rassam@csiro.au

AT CSIRO, WE DO THE
EXTRAORDINARY EVERY DAY
We innovate for tomorrow and help
improve today – for our customers, allAustralians and the world.
Our innovations contribute billions ofdollars to the Australian economy
every year. As the largest patent holder
in the nation, our vast wealth ofintellectual property has led to morethan 150 spin-off companies.
With more than 5,000 experts and aburning desire to get things done, we are
Australia’s catalyst for innovation.
CSIRO. WE IMAGINE. WE COLLABORATE.
WE INNOVATE.

4 | Tutorial

Dirk Mallants
CSIRO Land and Water
Waite Campus, Adelaidet +61 8 8303 8595
e dirk.mallants@csiro.au