Czech University of Life Sciences, Prague
Faculty of Agrobiology, Food and Natural Resourcesand

PC-Progress, Ltd.
The Second HYDRUS Workshop
March 28th 2008
Prague, Czech Republic

Proceedingsedited byJirka Šimůnek and Radka Kodešová

The texts of various papers in this volume were set individually by the authors or under theirsupervision and were not subject to additional editorial work.

ISBN: 978-80-213-1783-3

The Second HYDRUS Workshop
Content
1. Šimůnek, J. and R. Kodešová, Preface

1

2. Šimůnek, J., M. Šejna, and M. Th. van Genuchten, New features and developments
in HYDRUS software packages.

3

3. Ramos, T. B., M. C. Gonçalves, A. Prazeres, and J. C. Martins, N. L. Castanheira,Multicomponent solute transport in two multifactorial experiments.

13

4. Pot, V., P. Benoit, M. Le Menn, O. M. Eklo, T. Sveistrup, J. Kværner, Is thedispersion coefficient of sorbing tracers modified by retardation? Inverse modelingwith HYDRUS-1D.

21

5. Cheviron B. and Y. Coquet, and, Sensitivity analysis of HYDRUS-1D to transientMIM parameters: a case study related to pesticide fate in soils.

27

6. De Wilde, T., J. Mertens, J. Šimůnek, K. Sniegowksi, P. Spanoghe, J. Ryckeboer, P.
Jaeken, and D. Springael, Evaluating sorption and degradation characteristics ofpesticides using column displacement experiments.
35
7. Cheyns, K., J. Mertens, J. Šimůnek, E. Smolders, and D. Springael, Describing
pesticide degradation using Monod kinetics: information content of batch andcolumn data.

41

8. Twarakavi, N. K. C., J. Šimůnek, and S. Seo, A HYDRUS based approach forcoupled modeling of vadose zone and ground water flow at different scales.

47

9. Zhao, Y., S. Peth, and R. Horn, Modeling of coupled water and heat fluxes in both
unfrozen and frozen soils.

55

10. Jiménez-Martínez, J., T.H. Skaggs, M.Th. van Genuchten, L. Candela, HYDRUS1D Modeling of an Irrigated Agricultural Plot with Application to Aquifer Recharge
Estimation.
61
11. Matisoff, G., L. Vitko, P. Whiting, M. Ketterer, J. W. Mietelski, E. Lokas, K. Rosen,
H. Persson, Downward migration of Chernobyl-derived radionuclides in soils in
Poland and Sweden.
67
12. Langergraber, G., Overview of applications of the HYDRUS constructed wetland
module.

71

13. Jakimavičiūtė-Maselienė, V., J. Mažeika, Ž. Skuratovič, Water flow in the different
types of unsaturated soils in Lithuania.
79
14. Jacques, D., J. Šimůnek, D. Mallants, and M. Th. van Genuchten, Simulation of theeffect of time-variable water flow on the speciation and mobility of reactive
elements in soil.
85
15. Sakai, M., J. Šimůnek, and H. Saito, Evaluating Evapotranspiration from
Meteorological Data Using HYDRUS-1D.

93

16. Kodešová, R., M. Kočárek, J. Kozák, V. Kodeš, J. Šimůnek, N. Vignozzi, and M.
Pagliai, Impact of soil micro-morphological features on character of water flow andsolute transport.
99
17. Szegedi, K., D. Vetterlein, R. Jahn, H.-U. Neue, What comes before upscaling? Modeling transport and speciation in the rhizosphere with the new tool RhizoMath. 101
18. Slama, F., R. Bouhlila, and J. Tarhouni, Field data and modeling of salt transport ina coastal irrigated plain.
103
19. Barão, A. L. B., Pedro Chambel Leitão, Ramiro J.Neves, Maria Conceição
Gonçalves, Tiago B. Ramos, Nádia Castanheira, Simulation of water dynamics in
two irrigated soils.
105Preface
The first HYDRUS workshop was held on October 19, 2005 at the Department of Earth
Sciences of the Utrecht University in Utrecht, the Netherlands. 22 contributions were
presented at the workshop and assembled in the workshop proceedings (Torkzaban and
Hassanizadeh, 2005).
Since then the community of HYDRUS users has been continuously growing.
HYDRUS codes have been downloaded over two thousand times and HYDRUS web
pages are visited on average by about seven hundred individual visitors daily. Hundredsof research papers, in which HYDRUS codes have been used, have appeared in thepeer-reviewed literature. Also two major releases of new versions of HYDRUS
software packages have occurred. While the HYDRUS (2D/3D) software package wasreleased in 2006 as a complete rewrite of HYDRUS-2D and its extensions for two- andthree-dimensional geometries, version 4.0 of HYDRUS-1D was released in 2008.
To give HYDRUS community an opportunity to meet and share their experience with
HYDRUS codes and to learn about new features and developments in HYDRUS codes,
The Second HYDRUS Workshop was organized by the Department of Soil Scienceand Geology at Czech University of Life Sciences, Prague, Czech Republic on March
28, 2008, following a HYDRUS short course that was held during two previous days.
The purpose of the workshop was to bring together the users and developers of the
HYDRUS software packages, to present the latest innovations in the model applications,and to discuss the capabilities and limitations of HYDRUS. This was a unique
opportunity for all HYDRUS users to give a brief overview of their work and discuss
future needs and directions for obtaining more reliable modeling results. Over 40
scientists from four continents participated at this workshop.
These proceedings contain the collection of papers presented at the workshop. This
collective work includes contributions by many users of the HYDRUS software
packages ranging from the very fundamental to the most compelling and important
applications.
The editors:
Jirka Šimůnek
Radka Kodešová
Please refer to papers from these proceedings as follows:
Authors, Paper title, In: J. Šimůnek and R. Kodešová (eds.), Proc. of The Second
HYDRUS Workshop, March 28, 2008, Dept. of Soil Science and Geology, Czech
University of Life Sciences, Prague, Czech Republic, ISBN 978-80-213-1783-3, pp. ????, 2008.

1

2

New Features and Developments in HYDRUS Software Packages
J. Šimůnek
Department of Environmental Sciences, University of California Riverside, Riverside, CA 92521, USA

M. Šejna
PC-Progress, s.r.o., Anglicka 28, Prague 120 00, Czech Republic

M. Th. van Genuchten
U.S. Salinity Laboratory, USDA, ARS, Riverside, CA 92507, USA

ABSTRACT: The first HYDRUS workshop was held on October 19, 2005 at the Department ofEarth Sciences of the Utrecht University in Utrecht, the Netherlands. Two major developments
related to HYDRUS software packages have occurred since then. The main development
undoubtedly was the replacement of HYDRUS-2D with HYDRUS (2D/3D) in 2006. The
HYDRUS (2D/3D) software package is an extension and replacement of HYDRUS-2D and
SWMS_3D. This software package is a complete rewrite of HYDRUS-2D and its extensions fortwo- and three-dimensional geometries. The second major development was the release of the
new version (4.0) of HYDRUS-1D in 2008. In addition to many new features, GUIs of
HYDRUS-1D and HYDRUS (2D/3D) support also the biogeochemical flow and transportmodel HP1 and the constructed wetland module CW2D, respectively. Both software packages
represent major upgrades of previous versions, with many new processes considered and withsignificantly improved graphical user interfaces and more detailed online helps. Additionally,
HYDRUS-1D was significantly simplified and incorporated as the HYDRUS package into thegroundwater flow model MODFLOW.
1 INTRODUCTION
The first HYDRUS workshop was held on October 19, 2005 at the Department of Earth Sciences of the Utrecht University in Utrecht, the Netherlands. 22 contributions were presented at
the workshop and assembled in the workshop proceedings (Torkzaban and Hassanizadeh, 2005).
Since then two major releases of new versions of HYDRUS software packages has occurred.
While the HYDRUS (2D/3D) (Šimůnek et al., 2006; Šejna and Šimůnek, 2007) software package was released in 2006 as a complete rewrite of HYDRUS-2D (Šimůnek et al., 1999) and its
extensions for two- and three-dimensional geometries, version 4.0 of HYDRUS-1D (Šimůnek etal., 2008a) was released in 2008. In the text below we summarize new features that were implemented in both software packages. We also briefly discuss other modeling developments related
to HYDRUS family of models.
2 HYDRUS (2D/3D)
The HYDRUS (2D/3D) software package (Šimůnek et al., 2006; Šejna and Šimůnek, 2007)(Figure 1) is an extension and replacement of HYDRUS-2D (version 2.0) (Šimůnek et al., 1999)
and SWMS_3D (Šimůnek et al., 1995). This software package is a complete rewrite of HYDRUS-2D and its extensions for two- and three-dimensional geometries.

32.1 New features and processes in computational modules
In addition to features and processes available in HYDRUS-2D and SWMS_3D, the new computational modules of HYDRUS (2D/3D) consider (a) water flow and solute transport in a dualporosity system, thus allowing for preferential flow in fractures or macropores while storing water in the matrix (Šimůnek et al., 2003), (b) root water uptake with compensation, (c) the spatialroot distribution functions of Vrugt et al. (2001), (d) the soil hydraulic property models of Kosugi (1996) and Durner (1994), (e) the transport of viruses, colloids, and/or bacteria using an attachment/detachment model, filtration theory, and blocking functions (e.g., Bradford et al.,2002), (f) a constructed wetland module (only in 2D) (Langergraber and Šimůnek, 2005, 2006),(g) the hysteresis model of Lenhard et al. (1991) to eliminate pumping by keeping track of historical reversal points, (h) new print management options, (i) dynamic, system-dependent
boundary conditions, (j) flowing particles in two-dimensional applications, and (k) calculationsof actual and cumulative fluxes across internal meshlines.
2.2 Wetland module
A multi-component reactive transport model CW2D (Constructed Wetlands 2D) (Langergraberand Šimůnek 2005, 2006) was developed to model the biochemical transformation and degradation processes in subsurface-flow constructed wetlands. The model was incorporated into the
HYDRUS (2D/3D) variably-saturated water flow and solute transport software package. Constructed wetlands have become increasingly popular for removing organic matter, nutrients,
trace elements, pathogens, or other pollutants from wastewater and/or runoff water. Such wetlands involve a complex mixture of water, substrate, plants, litter, and a variety of microorganisms to provide optimal conditions for improving water quality. The water flow regime in subsurface-flow constructed wetlands can be highly dynamic and requires the use of transientvariably-saturated flow model. The biochemical components defined in CW2D include dissolved oxygen, three fractions of organic matter (readily- and slowly-biodegradable, and inert),
four nitrogen compounds (ammonium, nitrite, nitrate, and dinitrogen), inorganic phosphorus,
and heterotrophic and autotrophic micro-organisms. Organic nitrogen and organic phosphorus
were modeled as part of the organic matter. The biochemical degradation and transformation
processes were based on Monod-type rate expressions, such as for NO3-based growth of heterotrophs on readily biodegradable COD (denitrification):r = μDN

K DN,O2

cNO3

K DN,NO2

K DN,O2 + cO2 K DN,NO3 + cNO3 K DN,NO2 + cNO2

(1)

cCR
f N,DN cXH
K DN,CR + cCR

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

2.3 New features in the graphical user interface
New features of the Graphical User Interface of HYDRUS (2D/3D) include, among other things,(a) a completely new GUI based on Hi-End 3D graphics libraries, (b) the MDI (multi document
interface) architecture with multiple projects and multiple views, (c) a new organization of geometric objects, (d) a navigator window with an object explorer, (e) many new functions improving the user-friendliness, such as drag-and-drop and context sensitive pop-up menus, (f) im4

proved interactive tools for graphical input, (g) options to save cross-sections and mesh-lines forcharts within a given project, (h) a new display options dialog where all colors, line styles, fonts
and other parameters of graphical objects can be customized, (i) extended print options, (j) extended information in the Project Manager (including project previews), and (k) an option to export input data for the parallelized PARSWMS code (Hardelauf et al., 2007).

Figure 1. Main window of the HYDRUS (2D/3D) software package. Input and output data are accessibleusing the data tree at the Navigator Bar on the left. The computational domain with its finite element discretization, various domain properties, initial and boundary conditions, and results are displayed in one or
multiple View Windows in the middle. Various tools for manipulating data in the View Window areavailable on the Edit Bar on the right. The tabs in the View Window allow for fast access to different
types of data (Šimůnek et al., 2008b).

2.4 Website and documentation
The HYDRUS web site hosts a discussion forum for HYDRUS (2D/3D) (as well as other related programs) where users, after registering, can submit questions about the different software
packages and how to use them for their particular applications. Users there can also discussvarious topics related to modeling, or respond to questions posted by other users. The large
number of users of these discussion forums has made the forums nearly self-supporting in termsof software support and feedback.
The HYDRUS website also provides tutorials (Figure 2), including brief downloadable videos in which these tutorials are carried out step by step, thus allowing software users to teachthemselves interactively about the basic components of the software, including the process ofdata entry and display of calculated results.
We have also dramatically extended the documentation for HYDRUS (2D/3D). The installation of the latest HYDRUS (2D/3D) is accompanied with 240 pages of information in the technical manual, 200 pages of user manual, and over thousand pages of online context-sensitive
help. The software package furthermore comes with a suite of test problems, some of which aredescribed in detail in the technical manual.
Since previously published studies in which the program has been used can be a major sourceof information for new users, we are continuously updating the list of such publications at5http://www.pc-progress.cz/Pg_Hydrus_References.htm for HYDRUS-2D (or 2D/3D) and its
predecessors. Similar information is collected for HYDRUS-1D and related software packages
at http://www.pc-progress.cz/Pg_Hydrus1D_References.htm.

Figure 2. HYDRUS web page with HYDRUS tutorials and brief downloadable videos, demonstrating
step by step the use of the software package.

3 HYDRUS-1D
New features in version 4.0 of HYDRUS-1D (Šimůnek et al., 2008a) as compared to version 3.0(Šimůnek et al., 2005) include a) coupled water, vapor, and energy transport, b) dualpermeability type water flow and solute transport, c) dual-porosity water flow and solute transport, with solute transport subjected to two-site sorption in the mobile zone, d) option to calculate potential evapotranspiration the Penman-Monteith combination equation or with Hargreaves
equation, e) daily variations in the evaporation, transpiration, and precipitation rates, and f) support for the HP1 code, which was obtained by coupling HYDRUS with the PHREEQC biogechemical code (Parkhurst and Appelo, 1999). Selected new features are briefly discussed below.
3.1 Coupled water, vapor, and energy transport
Version 3.0 of HYDRUS-1D numerically solved the Richards equation that considered only water flow in the liquid phase and ignored the effects of the vapor phase on the overall water massbalance. While this assumption is justified for the majority of applications, a number of problems exist in which the effect of vapor flow can not be neglected. Vapor movement is often an
important part of the total water flux when the soil moisture becomes relatively low. Version 4.0
of HYDRUS-1D offers an option to simulate nonisothermal liquid and vapor flow, closely coupled with the heat transport (Saito et al., 2006):6∂θ ( h ) ∂ ⎡
∂h
∂T ⎤
=
( K + K vh ) ⎛⎜ + 1⎞⎟ + ( K LT + K vT ) ⎥ - S (h )
⎢
∂t
∂x ⎣
∂x ⎦
⎝ ∂x
⎠

(2)

∂T
∂θ
∂ ⎛ ∂T ⎞∂T∂q T∂q+ L0 v = ⎜ λ− Cv v − L0 v
⎟ - Cw q∂t∂t∂x ⎝ ∂x ⎠
∂x
∂x
∂x

(3)

C (θ )

In eq. (2), θ = total volumetric water content, being the sum (θ=θl +θv) of the volumetric liquidwater content, θl, and the volumetric water vapor content, θv (both expressed as an equivalentwater content), h = pressure head [L], T = temperature [K], K = isothermal hydraulic conductivity for the liquid phase [LT-1], KLT = thermal hydraulic conductivity for the liquid phase [L2K-1T1], Kvh = isothermal vapor hydraulic conductivity [LT-1], KvT = thermal vapor hydraulic conductivity [L2K-1T-1], and S = sink term representing root water uptake [T-1]. Overall water flow in(2) is given as the sum of isothermal liquid flow, isothermal vapor flow, gravitational liquid
flow, thermal liquid flow, and thermal vapor flow. Since several terms of (2) are a function oftemperature, this equation should be solved simultaneously with the heat transport equation (3)to properly account for temporal and spatial changes in soil temperature.
In eq. (3), λ = apparent thermal conductivity of the soil [MLT-3K-1] (e.g. Wm-1K-1); C(θ) and Cw =volumetric heat capacities [ML-1T-2K-1] (e.g. Jm-3K-1) of the porous medium and the liquid phase, respectively, q = fluid flux density [LT-1], L0 = volumetric latent heat of vaporization of liquid water[ML-1T-2] (e.g., Jm-3), and qv = vapor flux density [LT-1]. In equation (3), the total heat flux density is defined as the sum of the conduction of sensible heat as described by Fourier’s law (thefirst term on the right side), sensible heat by convection of liquid water (the second term) andwater vapor (the third term), and of latent heat by vapor flow (the forth term).
3.2 Physical and chemical nonequilibrium models
3.2.1 Physical nonequilibrium models
Version 4.0 of HYDRUS-1D implements several physical nonequilibrium water flow and solutetransport models. While mathematical description of these models is given in detail in Šimůnek
et al. (2008c) or the HYDRUS-1D manual, here we will present only the conceptual description.
A hierarchical set of physical nonequilibrium flow and transport models can be derived from the
Uniform Flow Model (Figures 3a). The equilibrium flow and transport model can be modified
by assuming that the soil particles or aggregates have their own micro-porosity and that water
present in these micropores is immobile (the Mobile-Immobile Water Model in Figure 3b).
While the water content in the micropore domain is constant in time, dissolved solutes can moveinto and out of this immobile domain by molecular diffusion. This simple modification leads tophysical nonequilibrium solute transport while still maintaining uniform water flow.
The mobile-immobile water model can be further expanded by assuming that both water andsolute can move into and out of the immobile domain (Šimůnek et al., 2003), leading to theDual-Porosity Model in Figure 3c. While the water content inside of the soil particles or aggregates is assumed to be constant in the Mobile-Immobile Water Model, it can vary in the DualPorosity Model since the immobile domain is now allowed to dry out or rewet during drying
and wetting processes. Water flow into and out of the immobile zone is usually described using
a first-order rate process. Solute can move into the immobile domain of the Dual-Porosity
Model by both molecular diffusion and advection with flowing (exchanging) water. Since water
can move from the main pore system into the soil aggregates and vice-versa, but not directly between the aggregates themselves, water in the aggregates can be considered immobile from alarger scale point of view.
The limitation of water not being allowed to move directly between aggregates is overcome
in Dual-Permeability Models (e.g., Gerke and van Genuchten, 1993). Water and solutes in suchmodels move also directly between soil aggregates as shown in Figure 3d. Dual-Permeability
Models assume that the porous medium consists of two overlapping pore domains, with waterflowing relatively fast in one domain (often called the macropore, fracture, or inter-porosity domain) when close to full saturation, and slow in the other domain (often referred to as the mi-7cropore, matrix, or intra-porosity domain). Like the Dual-Porosity Model, the Dual-Permeability
Model allows the transfer of both water and solutes between the two pore regions.
Finally, the Dual-Permeability Model can be further refined by assuming that inside of thematrix domain an additional immobile region exists into which solute can move by molecular
diffusion (the Dual-Permeability Model with MIM in Figure 3e) (Pot at al., 2005; Šimůnek etal., 2008c).
Equilibrium Model -----------------------------------------Non-Equilibrium Models ----------------------------------------a. Uniform Flowb. Mobile-Immobile Water

Water

Water

c. Dual-Porosity

Water
Immob.

SoluteSoluteImmob.

θS = θcd. Dual-Permeability e. Dual-Permeability with MIM

Water

MobileSlow

SoluteMobileImmob.

Water
Fast

Slow

SoluteMobileSlow

Fast

SoluteIm.

Fastθ = θ im + θ mo

θ = θ im + θ mo

θ = θM + θFS = θ im cim + θ mo cmo

S = θ im cim + θ mo cmo

S = θ M cm + θ F c f

Slow

Fast

θ = θ M ,im + θ M ,mo + θ F
S = θ M ,im cM ,im + θ M ,mo cM ,mo + θ F c f

Figure 3. Conceptual physical nonequilibrium models for water flow and solute transport. In the plots, θ
is the water content, θmo and θim in (b) and (c) are water contents in the mobile and immobile flow regions, respectively; θM and θF in (d) are water contents in the matrix and macropore (fracture) regions, respectively, and θM,mo, θM,im, and θF in (e) are water contents in the mobile and immobile flow regions of
the matrix domain, and in the macropore (fracture) domain, respectively; c are concentrations in corresponding regions, with subscripts having the same meaning as for water contents, while S is the total solute content of the liquid phase (Šimůnek et al., 2008c).

3.2.2 Chemical nonequilibrium models
Chemical nonequilibrium models implemented into HYDRUS-1D are schematically shown in
Figure 4. The simplest chemical nonequilibrium model assumes that sorption is a kinetic process(the One Kinetic Site Model in Figure 4a), usually described by means of a first-order rate equation. The one-site kinetic model can be expanded into a Two-Site Sorption model by assumingthat the sorption sites can be divided into two fractions. The simplest two-site sorption model
arises when sorption on one fraction of the sorption sites is assumed to be instantaneous, while
kinetic sorption occurs on the second fraction (Two-Site Model in Figure 4b). This model can befurther expanded by assuming that sorption on both fractions is kinetic and proceeds at different
rates (the Two Kinetic Sites Model in Figure 4c). The Two Kinetic Sites Model reduces to the
Two-Site Model when one rate is so high that it can be considered instantaneous, to the One Kinetic Site Model when both rates are the same, or to the chemical equilibrium model when bothrates are so high that they can be considered instantaneous.
a. One Kinetic Site Modelθb. Two-Site Modelseθc. Two Kinetic Sites Model

s1kd. Dual-Porosity Model withOne Kinetic Sitee. Dual-permeability Modelwith Two-Site Sorption

Immob. Mob.

Slow Fast

θ

θim

θmo

cim

cmo

smoe

sme

smok

smk

θM

θF

cm

cf

sfe

sime

sk
c

sk

c

s2k

c

sfk

Figure 4. Conceptual chemical nonequilibrium models for reactive solute transport. In the plots, θ is thewater content, θmo and θim in (d) are water contents of the mobile and immobile flow regions, respectively; θM and θF in (e) are water contents of the matrix and macropore (fracture) regions, respectively; c
are concentrations of the corresponding regions, se are sorbed concentrations in equilibrium with the liquid concentrations of the corresponding regions, and sk are kinetically sorbed solute concentrations of thecorresponding regions (Šimůnek et al., 2008c).

8

3.2.3 Physical and chemical nonequilibrium models
The combined physical and chemical nonequilibrium approach may be simulated with HYDRUS-1D using the Dual-Porosity Model with One Kinetic Site (Figure 4d). This model considers water flow and solute transport in a dual-porosity system (or a medium with mobileimmobile water), while assuming that sorption in the immobile zone is instantaneous. However,following the two-site kinetic sorption concept, the sorption sites in contact with the mobile
zone are now divided into two fractions, subject to either instantaneous or kinetic sorption.
Since the residence time of solutes in the immobile domain is relatively large, equilibrium likely
exists between the solution and the sorption complex here, in which case there is no need toconsider kinetic sorption in the immobile domain. The model, on the other hand, assumes thepresence of kinetic sorption sites in contact with the mobile zone since water can move relatively fast in the macropore domain and thus prevent chemical equilibrium (Šimůnek et al.,2008c).
Finally, chemical nonequilibrium can also be combined with the Dual-Permeability Model.
This last nonequilibrium option implemented into HYDRUS-1D (the Dual-Permeability Modelwith Two-Site Sorption in Figure 4e) assumes that equilibrium and kinetic sites exist in both themacropore (fracture) and micropore (matrix) domains. Applications of this transport model thatconsiders simultaneously both physical and chemical nonequilibrium has recently been presented by Pot et al., (2005), Köhne et al. (2006), and Kodešová et al. (2008).
3.3 Calculation of potential evapotranspiration
Potential evapotranspiration may be calculated in HYDRUS-1D using either the FAO recommended Penman-Monteith combination equation for evapotranspiration (ET0) (FAO, 1990) or
the Hargreaves equation (Jensen et al., 1997). With the Penman-Monteith approach, ET0 is determined using a combination equation that combines the radiation and aerodynamic terms asfollows [FAO, 1990]:
ET0 = ETrad + ETaero =ρ c p (ea - ed ) / ra ⎤1 ⎡ Δ ( Rn - G )+⎢⎥λ ⎣ Δ + γ (1 + rc / ra ) Δ + γ (1 + rc / ra ) ⎦(4)

where ET0 is the evapotranspiration rate [mm d-1], ETrad is the radiation term [mm d-1], ETaero is the
aerodynamic term [mm d-1], λ is the latent heat of vaporization [MJ kg-1], Rn is net radiation at
surface [MJ m-2d-1], G is the soil heat flux [MJ m-2d-1], ρ is the atmospheric density [kg m-3], cp isthe specific heat of moist air [i.e., 1.013 kJ kg-1 oC-1], (ea-ed) is the vapor pressure deficit [kPa], ea isthe saturation vapor pressure at temperature T [kPa], ed is the actual vapor pressure [kPa], rc is thecrop canopy resistance [s m-1], and ra is the aerodynamic resistance [s m-1].
The potential evapotranspiration can also be evaluated using the much simpler Hargreaves
formula (e.g., Jensen et al., 1997):(5)
ETp = 0.0023Ra (Tm + 17.8 ) TR-1
-2 -1
where Ra is extraterrestrial radiation in the same units as ETp [e.g., mm d or J m s ], Tm is thedaily mean air temperature, computed as an average of the maximum and minimum air temperatures [oC], TR is the temperature range between the mean daily maximum and minimum air temperatures [oC].
3.4 Daily variations in the evaporation, transpiration, and precipitation rates
Variations in potential evaporation and transpiration during the day can be generated with HYDRUS-1D using the assumptions that hourly values between 0-6 a.m. and 18-24 p.m. represent
1% of the total daily value and that a sinusoidal shape is followed during the rest of the day(Fayer, 2000), i.e.,

Tp (t ) = 0.24Tp

t < 0.264d, t > 0.736d⎛ 2π t π ⎞Tp (t ) = 2.75Tp sin ⎜− ⎟ t ∈ (0.264d, 0.736d)⎝ 1day 2 ⎠9(6)where Tp is the daily value of potential transpiration (or evaporation). Similarly, variation of
precipitation can be approximated using a cosine function as follows:⎡⎛ 2π t⎞⎤P(t ) = P ⎢1 + cos ⎜− π ⎟⎥
⎝ Δt
⎠⎦
⎣

(7)

where P is the average precipitation rate of duration Δt.
3.5 GUI support for the HP1 code
Graphical User Interface of HYDRUS-1D provides a support for the biogeochemical transport
code HP1. This is a complex modeling tool that was recently developed by coupling HYDRUS1D with the PHREEQC geochemical code (Parkhurst and Appelo, 1999). This coupling resulted
in a new comprehensive simulation tool, HP1 (acronym for HYDRUS1D-PHREEQC) (Jacquesand Šimůnek, 2005; Jacques et al., 2006, 2008ab). The combined code contains modules simulating (1) transient water flow in variably-saturated media, (2) the transport of multiple components, (3) mixed equilibrium/kinetic biogeochemical reactions, and (4) heat transport. HP1 is asignificant expansion of the individual HYDRUS-1D and PHREEQC programs by combining
and preserving most of their original features and capabilities into a single numerical model. The
code still uses the Richards equation for variably-saturated flow and advection-dispersion typeequations for heat and solute transport. However, the program can now simulate also a broad
range of low-temperature biogeochemical reactions in water, the vadose zone and in ground water systems, including interactions with minerals, gases, exchangers, and sorption surfaces,
based on thermodynamic equilibrium, kinetics, or mixed equilibrium-kinetic reactions.
4 HYDRUS PACKAGE FOR MODFLOW
Although computer power has increased tremendously during the last few decades, large scale
three-dimensional applications evaluating water flow in the vadose zone are often still prohibitively expensive in terms of computational resources. To overcome this problem, Seo et al.
(2007) developed a computationally efficient one-dimensional unsaturated flow HYDRUS
package and linked it to the three-dimensional modular finite-difference ground water model
MODFLOW-2000 (Harbaugh et al. 2000). The HYDRUS unsaturated flow package used HYDRUS-1D to simulate one-dimensional vertical variably-saturated flow. MODLOW zone arrays
were used to define the cells to which the HYDRUS package was applied. MODFLOW used thetime-averaged flux from the bottom of the unsaturated zone as recharge, and calculated a watertable depth which was then used as a pressure head bottom boundary for HYDRUS. Twarakavi
et al. (2008) provided a comparison of the HYDRUS package to other MODFLOW packages
that evaluate processes in the vadose zone and presented a field application demonstrating thefunctionality of the package.

5 CONCLUSSIONSOver the last 15 years the close collaboration between the University of California Riverside,and the U.S. Salinity Laboratory, and more recently with PC-Progress in Prague, Czech Republic, and SCK•CEN, Mol, Belgium, resulted in the development of a large number of computer
tools that are currently being used worldwide for a variety of applications involving the vadosezone. The need for codes such as HYDRUS is reflected by the frequency of downloading fromthe HYDRUS web site. For example, HYDRUS-1D was downloaded more than 200 times inMarch of 2007 by users from 30 different countries, and over one thousand times in 2006. The
HYDRUS web site receives on average some 700 individual visitors each day. We hope that the
HYDRUS family of models will remain as popular in the future as it is now.
10

REFERENCES
Bradford, S. A., S. R. Yates, M. Bettehar, and J. Šimůnek, Physical factors affecting the transport and fate ofcolloids in saturated porous media, Water Resour. Res., 38(12), 1327, doi:10.1029/2002WR001340,
63.1-63.12, 2002.
Durner, W., Hydraulic conductivity estimation for soils with heterogeneous pore structure, Water Resour.
Res., 32(9), 211-223, 1994.
Fayer, M. J. UNSAT-H version 3.0: Unsaturated soil water and heat flow model, theory, user manual, andexamples, Rep. 13249, Pacific Northwest National Laboratory, Richland, Washington, 2000.
Food and Agriculture Organization of the United Nations, Expert consultation on revision of FAO methodologies for crop water requirements, ANNEX V, FAO Penman-Monteith Formula, Rome Italy, 1990.
Gerke, H. H., & van Genuchten, M. Th., A dual-porosity model for simulating the preferential movementof water and solutes in structured porous media, Water Resour. Res. 29, 305-319, 1993.
Harbaugh, A.W., E. R., Banta, M. C. Hill, and M. G. McDonald, MODFLOW-2000, the U.S. Geological
Survey modular ground-water model user guide to modularization concepts and the ground-water flow
process, Denver, CO, Reston, VA, U.S. Geological Survey, 2000.
Hardelauf, H., M. Javaux, M. Herbst, S. Gottschalk, R. Kasteel, J. Vanderborght, and H. Vereecken,
PARSWMS: a parallelized model for simulating 3-D water flow and solute transport in variably saturated soils, Vadose Zone Journal, 6(2), 255-259, 2007.
Jacques, D., and J. Šimůnek, User Manual of the Multicomponent Variably-Saturated Flow and Transport
Model HP1, Description, Verification and Examples, Version 1.0, SCK•CEN-BLG-998, Waste andDisposal, SCK•CEN, Mol, Belgium, 79 pp., 2005.
Jacques, D., J. Šimůnek, D. Mallants, and M. Th. van Genuchten, Operator-splitting errors in coupled reactive transport codes for transient variably saturated flow and contaminant transport in layered soil
profiles, J. Contam. Hydrology, 88, 197-218, 2006.
Jacques, D., J. Šimůnek, D. Mallants, and M. Th. van Genuchten, Modeling coupled hydrological andchemical processes in the vadose zone: A case study on long term uranium migration following mineral phosphorus fertilization, Vadose Zone Journal, Special Issue “Vadose Zone Modeling”, 2008a (inpress).
Jacques, D., J. Šimůnek, D. Mallants and M. Th. van Genuchten, Modelling coupled water flow, solutetransport and geochemical reactions affection heavy metal migration in a Podzol soil, Geoderma,
2008b (in press).
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.
Kodešová, R., M. Kočárek, V. Kodeš, J. Šimůnek, and J. Kozák, Impact of soil micromorphological features on water flow and herbicide transport in soils, Vadose Zone Journal, Special Issue “Vadose
Zone Modeling”, 2008, (in press).
Köhne, J. M., S. Köhne, and J. Šimůnek, Multi-process herbicide transport in structured soil columns:
Experiment and model analysis, J. Contam. Hydrology, 85, 1-32, 2006.
Kosugi K., Lognormal distribution model for unsaturated soil hydraulic properties, Water Resour. Res.,32(9), 2697-2703, 1996.
Langergraber, G. and J. Šimůnek, Modeling variably-saturated water flow and multi-component reactive
transport in constructed wetlands, Vadose Zone Journal, 4, 924-938, 2005.
Langergraber, G., and J. Šimůnek, The Multi-component Reactive Transport Module CW2D for Constructed Wetlands for the HYDRUS Software Package, Manual – Version 1.0, HYDRUS Software Series 2, Department of Environmental Sciences, University of California Riverside, Riverside, CA, 72 pp.,
2006.
Lenhard, R. J., J. C. Parker, and J. J. Kaluarachchi, Comparing simulated and experimental hysteretic
two-phase transient fluid flow phenomena, Water Resour. Res., 27(8), 2113-2124, 1991.
Parkhurst, D. L., and C. A. J. Appelo, User’s guide to PHREEQC (version 2) – A computer program forspeciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations, Water Resources Investigation, Report 99-4259, Denver, Co, USA, 312 pp, 1999.
Pot, V., J. Šimůnek, P. Benoit, Y. Coquet, A. Yra, and M.-J. Martínez-Cordón, Impact of rainfall intensity on the transport of two herbicides in undisturbed grassed filter strip soil cores, J. of ContaminantHydrology, 81, 63-88, 2005.
Saito, H., J., J. Šimůnek, and B. Mohanty, Numerical analyses of coupled water, vapor and heat transport in
the vadose zone, Vadose Zone Journal, 5, 784–800, 2006.
Šejna, M., and J. Šimůnek, HYDRUS (2D/3D): Graphical User Interface for the HYDRUS Software
Package Simulating Two- and Three-Dimensional Movement of Water, Heat, and Multiple Solutes in
Variably-Saturated Media, published online at www.-pc-progress.cz, PC-Progress, Prague, Czech Republic, 2007.
Seo, H. S., J. Šimůnek, and E. P. Poeter, Documentation of the HYDRUS Package for MODFLOW11

2000, the U.S. Geological Survey Modular Ground-Water Model, GWMI 2007-01, International
Ground Water Modeling Center, Colorado School of Mines, Golden, Colorado, 96 pp., 2007.
Šimůnek, J., K. Huang, and M. Th. van Genuchten, The SWMS_3D code for simulating water flow and solutetransport in three-dimensional variably saturated media. Version 1.0, Research Report No. 139, U.S. Salinity Laboratory, USDA, ARS, Riverside, California, 155 pp., 1995.
Šimůnek, J., M. Šejna, and M. Th. van Genuchten, The HYDRUS-2D software package for simulating twodimensional movement of water, heat, and multiple solutes in variably saturated media. Version 2.0,
IGWMC - TPS - 53, International Ground Water Modeling Center, Colorado School of Mines, Golden,
Colorado, 251 pp., 1999.
Šimůnek, J., N. J. Jarvis, M. Th. van Genuchten, and A. Gärdenäs, Review and comparison of models fordescribing non-equilibrium and preferential flow and transport in the vadose zone, Journal of Hydrology, 272, 14-35, 2003.
Šimůnek, J., M. Th. van Genuchten, and M. Šejna, The HYDRUS-1D software package for simulating the
one-dimensional movement of water, heat, and multiple solutes in variably-saturated media. Version 3.0,
HYDRUS Software Series 1, Department of Environmental Sciences, University of California Riverside,
Riverside, CA, 270 pp., 2005.
Šimůnek, J., M. Th. van Genuchten, and M. Šejna, The HYDRUS Software Package for Simulating Twoand Three-Dimensional Movement of Water, Heat, and Multiple Solutes in Variably-Saturated Media,
Technical Manual, Version 1.0, PC Progress, Prague, Czech Republic, pp. 241, 2006.
Šimůnek, J., M. Šejna, and M. Th. van Genuchten, The HYDRUS Software Package for Simulating Twoand Three-Dimensional Movement of Water, Heat, and Multiple Solutes in Variably-Saturated Media,
User Manual, Version 1.0, PC Progress, Prague, Czech Republic, pp. 161, 2006.
Šimůnek, J., M. Šejna, H. Saito, M. Sakai, and M. Th. van Genuchten, The HYDRUS-1D Software Package for Simulating the Movement of Water, Heat, and Multiple Solutes in Variably Saturated Media,
Version 4.0, HYDRUS Software Series 3, Department of Environmental Sciences, University of California Riverside, Riverside, California, USA, pp. 315, 2008a.
Šimůnek, J., M. Th. van Genuchten, and M. Šejna, Development and applications of the HYDRUS andSTANMOD software packages, and related codes, Vadose Zone Journal, 5(2), Special Issue “Vadose
Zone Modeling”, 2008b (in press).
Šimůnek, J. and M. Th. van Genuchten, Modeling nonequilibrium flow and transport with HYDRUS,
Vadose Zone Journal, 5(2), Special Issue “Vadose Zone Modeling”, 2008c (in press).
Torkzaban, S. and S. M. Hassanizadeh (eds.), Proc. of Workshop on HYDRUS Applications, October 19,
2005, Department of Earth Sciences, Utrecht University, The Netherlands, ISBN 90-39341125, pp.
94, 2005.
Twarakavi, N. K. C., J. Šimůnek, and H. S. Seo, Evaluating interactions between groundwater and vadosezone using HYDRUS-based flow package for MODFLOW, Vadose Zone J., Special Issue “Vadose
Zone Modeling”, 2008 (in press).
Vrugt, J.A., M. T. van Wijk, J. W. Hopmans, and J. Šimůnek, One-, two-, and three-dimensional root water uptake functions for transient modeling, Water Resour. Res., 37(10), 2457-2470, 2001.

12

Multicomponent Solute Transport in Two Multifactorial
Experiments
T. B. Ramos, M. C. Gonçalves, A. Prazeres & J. C. Martins
Department of Soil Science, Estação Agronómica Nacional, Oeiras, Portugal

N. L. Castanheira
Department of Rural Engineering, University of Évora, Portugal

ABSTRACT: In this study we use HYDRUS-1D to analyze multicomponent solute transport intwo soils with medium and coarse textures, under field conditions, in a three year experiment, atthree depths (20, 40 and 60 cm). A triple emitter source irrigation system was used for water,
salt and nitrogen application. With this system it is intended to obtain two gradients, one for thefertilizer and the other for salinity, and to observe its effects on soil and culture. HYDRUS-1D
successfully predicted field measurements of the water content (R2>0.96; b>0.86), overall salinity (R2>0.71; b<1.17), and the concentration of soluble Na+ (R2>0.71; b>0.91). The correspondence between measured and simulated values for soluble Ca2+ and Mg2+ was lower (R2>0.571;
b<1.61), affecting the predictions of the sodium adsorption ratio (R2>0.72; b>0.61). The exchangeable sodium percentage presented the poorest relation between measured and simulated
values (R2>0.34; b<1.50). The simulations resulted in a closer correspondence in the medium
texture soil.
1 INTRODUCTION
The past several decades have seen considerable progress in the conceptual understanding andmathematical description of water flow and solute transport processes in the unsaturated zone
(van Genuchten & Šimůnek 2004). A variety of analytical and numerical models are now available to predict water and/or solute transfer processes between the soil surface and the groundwater table.
The HYDRUS-1D software package (Šimůnek et al. 2005) is one of the most popular models
and is based on the numerical solution of the Richards equation for variably saturated flow, andanalytical or numerical solutions of the Fickian-based convection-dispersion equation for solutetransport. The model also incorporates the major ion chemistry module developed for the UNSATCHEM model (Šimůnek et al. 1996, Šimůnek & Suarez 1994) which consider kinetic
chemical reactions and can be used to evaluate chemical reactions for solutions having very highionic strength.
The effectiveness of the major ion chemistry module of HYDRUS-1D has been recently
evaluated in a multi-year experiment conducted in soil lysimeters under field conditions (Gonçalves et al. 2006). The authors compared model simulated and experimental results on soil water content, concentrations of Na+, Ca2+, and Mg2+, electrical conductivity (EC) of the soil solution, the sodium adsorption ration (SAR) and the exchangeable sodium percentage (ESP), with
R2 ranging from 0.60 for the water content to 0.87 for the sodium adsorption ratio.
The objective of this study was to test the effectiveness of the major ion chemistry module ofthe HYDRUS-1D software package in multifactorial experiments under field conditions, in soils
with different textures, and using maize as a crop.

13

2 MATERIAL AND METHODS
2.1 Field experiment design
The field plot experiments were conducted from 2004 to 2006 using maize as crop, on a medium
textured Eutric Fluvisol at the Alvalade Experimental Station, and on a coarse texture Hortic
Antrosol in the Campus of Mitra of the University of Évora, both located in the Alentejo regionof southern Portugal.
A triple emitter source irrigation system was used for water, salt (Na+), and fertilizer (N) application. This system, adapted from Beltrão et al. (2002), consists of three trickle laterals connected together in order to form a triple joint lateral. The first of the laterals is connected to thesalt stock solution while the second one is connected to the nitrogen reservoir. The third lateralcontains fresh water and is used to obtain a constant water application rate for each dripping
point along the triple joint lateral. Gradients of applied salt (Na+) and nitrogen (N) are then possible by having different emitters along the laterals and varying their discharges to obtain various mixings between the three lines while maintaining constant application rates in each dripping point.

Figure 1. Layout of the triple emitter source design.

Each experimental field (Figure 1) was divided into four groups (I to IV) with three triple
joint laterals each, and the N gradient decreasing from group I to IV. Each group was then divided into 3 treatments (A to C), with 6.75 m2 (2.25 m wide x 3 m long) each, and the Na+ gradient decreasing from treatment A to C. The discharges of the emitters along the triple joint laterals are described in Table 1. Three sources of water supply together a constant amount of 18
L/h in each dripping point and per linear meter of crop. The dripping points were distanced one
meter from the next, totalizing 9 emitters in each of the 12 treatments. Two laterals of fresh water bordered the different groups while each treatment was bordered with earthen ridges to prevent surface runoff during rainfall and irrigation.
Table 1. Discharge rates of the laterals applying salt (Na+), nitrogen (N) and fresh water (W) in each experimental plot.
Application Rates (L/h)
TreatmentGroup IGroup IIGroup III
Group IV
N
W
Na+
N
W
Na+
N
W
Na+
N
W
Na+
A
12
6
0
12
4
2
12
2
4
12
0
6
B
6
6
6
6
4
8
6
2
10
6
0
12
C
0
6
12
0
4
14
0
2
16
0
0
18

The N dynamics will not be referred any further in this manuscript as it was not handled with
HYDRUS-1D. The total amount of water applied during the three years of the experiment is described in Table 2. The concentration of saline waters applied in all emitters of the salt laterals
14was constant in each year and is described in Table 3. However, different discharges of emittersalong the laterals, made possible to apply different amounts of Na+ in the three treatments (A, B,
C). The total amount of Na+ applied in the saline water is described in Table 4.
Table 2. Total amount of water applied (fresh water + saline water + water with fertilizer) in Alvalade andMitra during the three irrigation cycles of the experiment.
Water applied (mm)Experimental field
2004
2005
2006
Alvalade
997
1012
1028
Mitra
1067
725
729
Table 3. Ion concentrations of irrigation waters.

AlvaladeFresh water
Saline water (2004-2005)
Saline water (2006)
Mitra
Fresh water
Saline water (2004-2006)

Ca2+

Mg2+

Na+
K+
-1
(meq L )

Cl-

1.64
2.05
2.1

2.25
2.28
1.84

4.16
56.13
100.0

0.06
0.24
0.2

8.11
60.70
104.14

1.90
2.80

1.36
2.68

1.00
50.00

0.04
0.18

4.30
55.66

Table 4. Total amount of salts (Na+) applied in Alvalade and Mitra in each Treatment of the four Groups
during the three irrigation cycles of the experiment.
Salt (g m-2)
Treatment
2004
2005
2006
Alvalade Mitra Alvalade Mitra Alvalade Mitra
A
1365
1352
2055
962
2792
1229
B
683
676
1027
481
1396
614
C
0
0
0
0
0
0

Soil solutions were collected from ceramic cups from Group IA (highest N and Na+ application) and Group IVC (lowest N and Na+ application), at three depths (20, 40 and 60 cm), twice aweek during the irrigation periods and once a week during the remaining months, to measure thesoluble ions (Na+, Ca2+, Mg2+, K+) and the electrical conductivity (EC). Soil water contents weremonitored using a TDR system at the same locations, depths and dates. At the beginning of theexperiment, at the end of each irrigation period, and after the following winters, soil samples
were collected at 3 depths (0-20, 20-40, 40-60 cm) to measure the exchangeable cations.
2.2 Input data
The following parameters were determined during the experiment in order to satisfy the considerable input data needed to run the major ion chemistry module of the HYDRUS-1D software
package (Šimůnek et al. 2005):− Crop evapotranspiration (ETc) computed, between April 2004 and February 2007, from theproduct of the reference evapotranspiration (using weather data and according to PenmanMonteith method) and a crop coefficient. The Leaf Area Index and the corresponding Soil
Cover Factor were used to account for different stages of the soil cover and to divide the ETcdaily values into crop transpiration and soil evaporation rates as required by HYDRUS-1D;− Soil hydraulic properties (Table 5) measured in the laboratory in the beginning of the experiment, and described with the Mualem-van Genuchten equation using RETC code (vanGenuchten et al. 1991);− Solute transport parameters (dispersivities) obtained by fitting analytical solutions of the
CDE to observed chloride breakthrough data using the nonlinear parameter estimation codeCXTFIT 2.1 (Toride et al. 1995) (Table 6);
Table 5. Mualem-van Genuchten parameters for soil hydraulic functions.
15

Depth (cm)θrθsα (cm-1)
ηA

Ks (cmd-1)

0-30
0
0.372
0.245
1.146
-7.333
41.4

Mitra
30-50
0
0.378
0.141
1.135
-3.14342.0

Alvalade
0-30
30-75
0.05
0.108
0.380
0.380
0.027
0.11481.205
1.186-4.411-5.366
16.6
84.4

50-90
0
0.332
0.079
1.184
0
101.6

Table 6. Dispersivity values for the two studied soils.
Mitra
Depth (cm)
0-30
30-50
Bulk density (g/cm3)
1.51
1.70
Dispersivity (cm)
273.2
54.2

50-90
1.69
209.2

Alvalade
0-30
30-75
1.49
1.51
5.36
0.68

− Physical and chemical soil characteristics to describe the initial soil conditions, namely thesoluble and exchangeable cations (Table 7);
Table 7. Physical and chemical soil characteristics (initial conditions) in two experimental fields.
Mitra
Alvalade
Treatment A
Treatment C
Treatment A
Treatment C
Depth (cm)0-30 30-50 50-90 0-30 30-50 50-90 0-30 30-75 0-30 30-75
Soluble
1.06
0.71
1.40
0.91
0.63
2.96
2.94
2.75
0.94
Ca2+ 1.23
cations
Mg2+ 0.98
0.83
0.59
1.10
0.68
0.53
1.28
1.18
1.33
0.56
(meq/L)
Na+
3.05
2.44
2.98
2.54
2.20
2.94
0.60
0.80
0.66
0.62
+
0.44
0.60
0.09
0.58
0.40
0.08
0.71
0.57
0.55
0.30
K
5.69
4.92
4.36
5.62
4.19
4.17
5.55
5.49
5.29
2.42
ClExchang.
Ca2+ 60.01 61.64 62.77 62.51 64.18 66.30 98.10 78.10 72.25 50.24
cations
Mg2+ 17.50 18.51 20.35 17.37 18.01 21.48 13.16 12.14 11.34 8.79
(meq/kg)
Na+
2.53
2.28
3.49
2.17
2.41
3.71
0.99
1.03
0.90
1.03
+
K
5.71
7.10
2.95
6.35
5.86
2.95
6.16
4.68
4.06
3.54
CEC (meq/kg)
85.75 89.53 89.56 88.40 90.46 94.44 118.4 95.95 88.55 63.60− Gapon selectivity coefficients (White & Zelazny, 1986) for describing the partition between
the solid exhange phase and the solution phase (cation exhange reaction) (Table 8);
Table 8 - Gapon selectivity coefficients for Ca2+, Na2+, Mg2+ and K+ relations.
Alvalade
Mitra
K
Treatment
0-30
30-75
75-100
0-30
30-50
A
0.326
0.339
0.354
0.204
0.246
KMg/Ca
C
0.314
0.325
0.353
0.226
0.227
A
2.923
2.872
2.855
1.540
1.582
KCa/Na
C
2.767
2.747
2.972
1.425
1.395
A
0.112
0.108
0.096
0.090
0.096
KK/Ca
C
0.111
0.098
0.088
0.089
0.088

50-90
0.232
0.213
1.382
1.534
0.092
0.079

− Concentration of the irrigation water (fresh and saline waters) applied during the irrigation
periods (Table 3);− Root development monitored using the Minirhizotron technique associated with a video
camera system and digital technique for image analyses. Maize reached about 60 cm depth.
2.3 Statistical analysis
In addition to a visual check, field-measured values were compared with results of the HYDRUS-1D predictions using simple regression analysis (y=bx, with zero intercept). Values close
to 1.0 for the determination (R2) coefficient, and for the regression coefficient (b) represent a
good agreement. Regression analyses between measured and simulated data were carried out for
16

water contents, electrical conductivities of the soil solution (representing the overall salinity),
for the individual ion concentrations (Na+, Ca2+, Mg2+), and for the SAR and ESP relationships.
3 RESULTS AND DISCUSSION
The experiments started on 22 and 27 April 2004 for Alvalade and Mitra, respectively, andlasted for a period of 3 years. As initial conditions for the HYDRUS-1D simulations we used, inboth locations, a water content value of 0.25 cm3 cm-3 for the entire soil profile. After a small
period of adjustment, simulated water content values in HYDRUS-1D closely mirrored measured values at the three monitored depths. In Alvalade, an overall determination coefficient (R2)of 0.986 and a regression coefficient (b) of 0.997 was obtained for n=282 observations (Table
9). In Mitra, R2 was 0.962 and b was 0.862 for n=379 observations. Figures 2 and 3 shows anexample of the soil observed and simulated water content data in Alvalade and Mitra, respectively, at 20 cm depth. When comparing measured and simulated values, one must realize thatwater contents were measured only sporadically and never during or immediately after rainfall
events, thus avoiding the highest water content values. On the other hand, HYDRUS-1D calculates the entire process and predicts continuous water contents during all events, including possible extremes. These differences are more visible in the coarse textured soil of Mitra, mainly
because its higher water content dynamics resulted in greater variations of the water content in
the first hours after irrigation.
HYDRUS-1D simulations of the total soil salinity resulted in a generally good agreementwith the observed patterns in both locations. Values of 0.819 for R2 and 1.170 for b (n = 247 observations) were obtained for the relation between measured and simulated EC values in Alvalade. In Mitra, R2 was 0.710 and b was 1.160 for n=363 observations (Table 9). Figures 2 and3 show examples of observed and simulated values of the electrical conductivity of the soil solution in Alvalade and Mitra, respectively, at a 40-cm depth.
Measured and simulated concentrations of soluble Na+ resulted also in a generally good
agreement in both locations. In Alvalade, values of 0.809 for R2 and 1.046 for b were obtained
between measured and simulated sodium concentrations for the 242 observation data points. In
Mitra, the agreement between simulated and measured sodium concentrations presented valuesof 0.710 for R2 and 0.907 for b for the 356 observation data points.
Differences between measured and simulated concentrations for the soluble Ca2+ and Mg2+
resulted in lower determination coefficients than for Na+ and EC. Also, b was further from 1.0
than before. The R2 for the soluble Ca2+, in both fields, was not higher than 0.644 and b gave
values higher than 1.212, while for the soluble Mg2+ the best value found was 0.659 for R2 and1.375 for b. These lower R2, and b values different from 1.0, were certainly due to its occurrencein small quantities (they were not added to the irrigation waters) and due to erroneous readingsof these cations in the atomic absorption spectrophotometer equipment in 2006. Consequently,
the SAR determination was also affected in 2006 (Table 9, and Figures 2, and 3).
ESP calculated from HYDRUS-1D simulations of the exchangeable cations resulted in poor
relations with the ESP determined from the measured exchangeable cations for Alvalade and especially for Mitra.
Table 9. Results of the regression analysis between measured and simulated soil water contents (θ), electrical conductivity of the soil solution, soluble Na+, Ca2+, and Mg2+ concentrations, sodium adsorption ratios, and exchangeable sodium percentage.
Alvalade
Mitra
Parameters
n
b
R2
n
b
R2
θ
0.997
0.986
282
0.862
0.962
379
EC
1.170
0.819
247
1.160
0.710
363
1.046
0.809
242
0.907
0.710
356
Na+
1.212
0.644
227
1.396
0.638
345
Ca2+
1.375
0.659
233
1.610
0.571
343
Mg2+
SAR
0.726
0.755
219
0.605
0.717
327
ESP
1.498
0.773
36
1.439
0.338
36

17

Alvalade
20 cm

0.4

3

-3

Water content (cm cm )

0.5

0.3
0.2
0.1
I

R

I

18-12-04

06-07-05

R

I

R

0.0
01-06-04

22-01-06

10-08-06

26-02-07

Time

Alvalade
40 cm

20

EC (dS/m)

15
I

R

I

R

I

R

10
5
0
01-06-04

18-12-04

06-07-05

22-01-06
Time

10-08-06

26-02-07

Alvalade
20 cm

140

100
80

I

R

01-06-04

18-12-04

60

I

R

I

R

+

[Na ] (mmol(c)/L)

120

40
20
0

06-07-05

22-01-06

10-08-06

26-02-07

Time

Alvalade
60 cm

SAR (meq/L)

0.5

40
30
I

R

I

R

I

20

R

GI A - HYDRUS
GI A - Measured
GIV C - HYDRUS
GIV C - Measured

10001-06-04

18-12-04

06-07-05

22-01-06
Time

10-08-06

26-02-07

Figure 2. Example of the results obtained in Alvalade for the water content, electrical conductivity (EC)of the soil solution, soluble sodium concentration and sodium adsorption ration. I and R correspond to theirrigation and rainfall periods, respectively.
18

Mitra
20 cm

0.4

3

-3

Water content (cm cm )

0.5

0.3
0.2
0.1
I
0.0
01-06-04

IR
18-12-04

R

06-07-05
22-01-06
Time

I

R

10-08-06

26-02-07

Mitra
40 cm

15

EC (dS/m)

12
I

9

R

I

R

I

R

22-01-06

10-08-06

6
3
0
01-06-04

18-12-04

06-07-05

26-02-07

Time

Mitra
60 cm

60
[Na+] (mmol(c)/L)

50
40
I

30

R

I

R

I

R

20
10
0
01-06-04

18-12-04

06-07-05
22-01-06
Time

0.5

26-02-07

Mitra
40 cm

20

SAR (meq/L)

10-08-06

GI A - HYDRUS
GI A - Measured
GIV C - HYDRUS
GIV C - Measured

15

10
I

R

I

18-12-04

06-07-05

RI

22-01-06

10-08-06

R

5

0
01-06-04

26-02-07

Time

Figure 3. Example of the results obtained in Mitra for the water content, electrical conductivity (EC) of
the soil solution, soluble sodium concentration and sodium adsorption ration. I and R correspond to the irrigation and rainfall periods, respectively.
19

4 CONCLUSIONHYDRUS-1D successfully simulated the water regime, as well as the effects of different irrigation waters on the geochemistry of the studied Fluvisol and Antrossol. The correspondence between observed and simulated variables was good, considering that simulations were carried out
to predict field measurements over a considerable time period (3 years) without any calibration,
and with all input variables (i.e., the soil hydraulic properties, solute transport parameters, atmospheric demand, Gapon constants, physical and chemical characteristics of the soil, LAI androot depth) measured independently.
In Alvalade, the agreement between measured and simulated values was best for the watercontent (R2=0.97; b=1.00), EC (R2=0.82; b=1.17) and soluble sodium concentrations (R2=0.81;
b=1.05). In Mitra, the best agreement was obtained also for the water content (R2=0.96; b=0.86),
EC (R2=0.71; b=1.16) and soluble sodium concentrations (R2=0.71; b=0.91), but with slightly
lower coefficients.
Soluble calcium (R2=0.64; b=1.21) and magnesium concentrations (R2=0.66; b=1.38), SAR
2(R =0.76; b=0.73), and ESP (R2=0.77; b=1.50) were predicted less well in Alvalade. In Mitra,
the same parameters provided again slightly worse coefficients.
In spite of the considerable input data demands, HYDRUS-1D proved to be an effective toolthat may become very useful for irrigation management, and for predicting the effects of irrigation water quality on soil and groundwater quality.
ACKNOWLEDGMENTS
This work was financed by project PTDC/AGR-AAM/66004/2006 of the Fundação para a Ciência e a Tecnologia.
REFERENCES
Beltrão, J., Jesus, S.B., Silva, V., Sousa, P.B., Carvalho, I., Trindade, D., Rodrigues, M.H. & Machado,
A. 2002. Efficiency of triple emitter source (TES) for irrigation experiments of horticultural crops.
Acta Hort. (ISHS) 573:183-188.
Gonçalves, M.C., Šimůnek, J., Ramos, T.B., Martins, J.C., Neves, M.J. & Pires, F.P. 2006. Multicomponent solute transport in soil lysimeters irrigated with waters of different quality. Water Resour. Res.,
42 (W08401, doi:10.1029/2005WR004802).
Šimůnek, J. & D.L. Suarez, D.L. 1994. Two-dimensional transport model for variably saturated porous
media with major ion chemistry. Water Resources Research, 30, 1115-1133.
Šimůnek, J., Suarez, D.L. & Šejna, M. 1996. The UNSATCHEM software package for simulating onedimensional variably saturated water flow, heat transport, carbon dioxide production and transport, andmulticomponent solute transport with major ion equilibrium and kinetic chemistry, Version 2.0, Research
Report No. 141, U.S. Salinity Laboratory, USDA, ARS, Riverside, California, 186pp.
Šimůnek, J., van Genuchten, M. Th. & Šejna M. 2005. The HYDRUS-1D software package for simulating the
one-dimensional movement of water, heat, and multiple solutes in variably-saturated media. Version 3.0,
HYDRUS Software Series 1. Department of Environmental Sciences, University of California Riverside,
Riverside, CA, 270pp.
Toride, N., Leij, F.J. & van Genuchten M. Th. 1995. The CXTFIT code for estimating transport parameters from laboratory or field tracer experiments, version 2.0. Res. Rep. 137, U.S. Salinity Lab., Riverside, Calif.
Van Genuchen, M. Th. & Šimůnek, J. 2004 Integrated modeling of vadose zone flow and transport processes. In R.A. Feddes, G.H. de Rooij & J.C. van Dam (eds), Unsaturated Zone Modelling: Progress,
Challenges and Applications. Proc. Wageningen, The Netherlands, October 3-5, pp. 37-69.
van Genuchten, M. Th., Leij, F.J. & Yates, S.R. 1991, The RETC code for quantifying the hydraulicfunctions of unsaturated soils, Rep. EPA/ 600/2-91-065, U.S. Environ. Prot. Agency, Washington,
D.C.
White, N.L. & Zelazny L.M. 1986. Charge properties in soil colloids. In D.L. Sparks (ed) Soil Physical
Chemistry, pp. 39– 81, CRC Press, Boca Raton, Fla.

20

Is the dispersion coefficient of sorbing tracers modified bysorption? Inverse modeling with HYDRUS-1D
V. Pot, P. Benoit, M. Le Menn,
UMR INRA AgroParisTech, Environment and Arable Crops, F-7885,Thiverval-Grignon, France

O.M. Eklo,
Bioforsk – Norwegian Institute for Agricultural and Environment Research, Plant health and Plant
Protection Division, Høgskoleveien 7, N-1432 Ås, Norway

T. Sveistrup, J. Kværner
Bioforsk – Norwegian Institute for Agricultural and Environment Research, Soil and Environment
Division, Frederic A. Dahls vei 20, N-1432 Ås, Norway

ABSTRACT: Dispersion is a key process describing spreading of solutes. Aris (1959) demonstrated that retardation (equilibrium sorption) does affect the value of the dispersion coefficientand increases it as compared to non-sorbing solute in the case of cylindrical pores. This phenomenon is named Fickian sorption-enhanced dispersion. This paper presents an attempt toevaluate the presence of the latter process in the complicated geometry of soil pore space. To doso, Fickian sorption-enhanced dispersion and non-equilibrium two-site sorption were evaluatedto describe metribuzin transport in undisturbed soil cores sampled in the A, B and C horizons ofa sandy loam soil. Inverse modelling of the water tracer (bromide) and metribuzin breakthroughcurves (BTCs) was performed in the framework of the CDE, with HYDRUS-1D.
1 INTRODUCTION
Transport of non-sorbing tracer is usually correctly described by the convection dispersion
equation (CDE), in the absence of physical non-equilibrium transport processes.
In the case of viscous Poiseuille flow in cylindrical tubes and parallel plates, Taylor and Aris
explained the hydrodynamic dispersion in the longitudinal direction as a combined effect of molecular diffusion in the liquid phase and velocity variations with the radius (thickness) of the
tube (parallel plate). The following expression for dispersion coefficient upscales the two effects
for asymptotic times and is further referred to as Taylor dispersion (Taylor, 1953):

D Lns = Dm + κ

U 2a2
Dm

(1)

where Dm is the molecular diffusion coefficient (cm2 h-1), U is the average velocity in the tube
(cm h-1), a is the radius (thickness) of the tube (parallel plate) (cm) and κ is a constant dependingon the velocity profile in the medium (-). The constant κ is 1/48 and 1/210 for the tube and parallel plate, respectively. Using Taylor dispersion (1) in the CDE exactly describes the concentration distribution of a non-sorbing tracer in these geometries, that is, Taylor dispersion is truly
Fickian.
Transport of sorbing tracers is further described by an additional term to the CDE. One usually assumes that the dispersion coefficient of the sorbing tracer is identical to that of the nonsorbing tracer, provided for a correction of the molecular diffusion coefficient (e.g. Prata et al.,2003). However, Aris (1959) established a general expression of the variation in time of the second moment, σ2, of the spatial distribution of the sorbing solute in a fluid flowing through a tube
applied to the case in which the solute can also pass into another fluid phase flowing in an annular region around the first. Matrix diffusion, instantaneous and kinetic sorptions were taken into
account. Reducing the Aris’s general expression to the case of one sorbing wall, one fluid, laminar flow, isotropic and uniform diffusion coefficient gives (Harvey and Marsily, 1995):

21

(R − 1)2 U 2 a
1 dσ 2 1 ⎛U 2a2 ⎞ R −1⎟⎟ += ⎜⎜ Dm + κ 1
Ds +
2 dt
R⎝
Dm ⎠
R
R 3 2k d α k(2)where Ds is the molecular diffusion in solid phase (cm2 h-1), kd is the partition coefficient between solid and liquid phases (L kg-1), αk is the mass transfer rate coefficient (h-1) and R is theretardation factor.
The first term of the right-hand side of equation (2) represents the Taylor-Aris dispersion
mechanism due to the combined effect of molecular diffusion in the liquid phase and velocity
profile in the liquid phase. It is similar to equation (1) except for the constant κ1 that differs fromthe constant κ due to the retardation factor:⎛

κ 1 = κ ⎜ A1 +⎝A2 A3 ⎞+⎟
R R2 ⎠(3)

where A1, A2 and A3 are constants depending on the geometry. For a tubular geometry (A1 = 11,
A2 = -16, A3 = 6), for a parallel plate (A1 = 78, A2 = -147, A3 = 70).
The second term of the right-hand side of equation (2) represents the effect of solute diffusion into the solid phase. This term is generally neglected when considering sorption process.
The third term of the right-hand side of equation (2) represents the effect of rate-limited masstransfer. This latter term results in tailing of the concentration distribution due to the rate-limited
remobilization of solute from either an immobile phase or kinetic sorption sites. This process is
described in the CDE by a source/sink term and is thus not of a Fickian form. From equation (6),Aris demonstrated that retardation does affect the value of the dispersion coefficient and increases it as compared to non-sorbing solute. The dispersion coefficient of a sorbing tracer
writes thus:

DLs = Dm + κ 1

U 2a2
Dm

(4)

Introducing (4) in the CDE of a sorbing tracer correctly described the concentration distribution of a sorbing tracer in parallel plate as numerically shown by Pot and Genty, 2005. The enhanced dispersion was also correctly predicted by experimental investigation of solute exchange
between phases in capillary tube (Grosser et al., 1991). This sorption-enhanced dispersion isthus a truly Fickian mechanism and is hereafter named as Fickian sorption-enhanced dispersion.
From equations (2) and (3), retardation is expected to be the predominant mechanism in increasing the dispersion (by a factor up to A1) for highly sorbing solutes, while rate-limited masstransfer becomes negligible.
In this work, we attempted to evaluate the presence of Fickian sorption-enhanced dispersion
in the complicated geometry of soil pore space.
2 MATERIALS AND METHODS
We analyzed sorption and transport of metribuzin [4-amino-6-tert-butyl-4,5-dihydro-3methylthio-1,2,4-triazin-5(4H)-one], a fairly mobile herbicide, by the means of transport column
experiments. Displacement experiments of Br- and metribuzin were realized in undisturbed soilcores (12 cm inner diameter) sampled in 2002 at Grue in Southern Norway in the A, B and C
horizons of a sandy loam soil from the depths 0-25 cm, 28-43 cm and 50-75 cm respectively.
The experimental transport conditions were established using fixed and steady-state water potential conditions (-10 cm), in order to achieve transport in the same range of pore size for all
soil columns. A pulse of CaCl2 water solution (5 x 10-3 mol L-1) containing the water tracer, Br-,and the herbicide metribuzin was injected.
Inverse modelling of Br- and metribuzin breakthrough curves (BTCs) was performed in theframework of the CDE with HYDRUS-1D (Šimůnek et al., 1998, 2003). The metribuzin data
were first analyzed using the chemical equilibrium transport model written for a sorbing and decaying solute. The metribuzin dispersion was either fixed to the value obtained from Br- BTCs
(hypothesis 1, model CDE I) or estimated (hypothesis 2, model CDE II). Then the same analysis
was repeated using the chemical non-equilibrium transport model based on the two-site sorption

22

model (Cameron and Klute, 1977) using hypothesis 1 (two-site I model) and hypothesis 2 (Twosite II model).
The soil profile was assumed to be initially solute-free. A flux concentration (third-type) anda zero concentration gradient were used as the upper and lower boundary conditions, respectively. The transport models were calibrated against measured BTCs using flux concentrations at
the lower boundary.
3 RESULTS
All measured bromide BTCs were well described by the classical physical equilibrium transportmodel based on the CDE (Figure 1). The fitted dispersivity values were found smaller than 1 cm
for all soil columns indicating thus that spatial heterogeneity was small.
The non-equilibrium chemical models generally performed better than the equilibrium
chemical models in the order: Two-site II > Two-site I > CDE II >> CDE I (Figure 1). The
range of the values of the distribution coefficient, kd was lower than the kd values measured inbatch conditions (Benoit et al., 2007; Stenrød et al., 2007). For the surface horizon, the half-life,of about 4 days, was shorter than corresponding static laboratory incubations of 26 days at 28°C
(Stenrød et al., 2007). This high degradation rate could be attributed to irreversible sorption (defined as fast attachment and much slower detachment) at the time scale of transport (Beigel and
Di Pietro, 1999; Pot et al., 2005). The metribuzin dispersion coefficients in the CDE II and Twosite II models were larger than the bromide dispersion coefficients but generally to a lesser extent when using the Two-site II model than when using the CDE II model.
For the Grue soil C horizon, only the CDE II model yielded sound transport parameters. This
result would support the hypothesis that retardation caused a Fickian increase of the dispersion
coefficient during transport of metribuzin in this horizon. The mean ratio of metribuzin to Brdispersion coefficients was of 2.0 ± 0.5 for a mean retardation factor value of 1.4 ± 0.1. This ratio is in good agreement with the one found in numerical investigation of Fickian sorptionenhanced dispersion in a natural rough fracture (ratio of about 2 for a retardation factor value of
1.6) (Pot and Genty, 2007). However, this increased dispersion coefficient may also reflect additional non-equilibrium processes not taken into account by non-equilibrium two-site sorption.
For instance, slow diffusion of metribuzin to the sorption places or intraparticle diffusion could
result in an effective dispersion coefficient higher than the Br- dispersion coefficient (Brusseau,
1993).
For the A and B horizons, the Two-site II model yielded larger fraction of equilibrium sites
together with smaller αk values compared to the ones obtained by the Two-site I model, so thatthe overall rate of reaction predicted by both models is more or less the same. It is thus to notethat both non-equilibrium models (Two-site I and Two-site II models) yielded similar nonequilibrium degree, indicating that the increase of dispersion was not compensated by lower
non-equilibrium. This is indicative of robustness of the models. The Two-site II model thus described metribuzin transport in the A and B horizons of Grue by assuming high dispersion coefficients and kinetic sorption. The mean ratio of metribuzin to bromide dispersion coefficients of
2.1 ± 0.6 (R = 1.8 ± 0.2) was similar to the one found for the C horizon of Grue soil and thenatural rough fracture (Pot and Genty, 2007). However, as discussed for the CDE II model, additional non-equilibrium (like intraparticle diffusion or film diffusion) may also result in a similar effective increase of the metribuzin dispersion coefficient.
4 CONCLUSIONSFor the sandy loam soil investigated in this work, the classical physical equilibrium transportmodel based on the CDE satisfactorily described bromide transport. The non-equilibrium twosite sorption combined with Fickian sorption-enhanced dispersion (Two-site II model) better described the metribuzin fate during transport in the A and B horizons of Grue soil. Fickian sorption-enhanced dispersion (CDE II model) better described the metribuzin fate during transport in
the C horizon of Grue soil. In our opinion, we can only conclude from this study that either retardation may influence dispersion and increase it in porous medium or that other nonequilibrium processes not included in the models we used occurred and resulted in higher effective dispersion. We thus emphasize that more theoretical and modelling studies are needed to23

elucidate the potential dependency relation of a Fickian form between dispersion and retardationin porous medium and macrodispersion regime.
REFERENCES
Aris R., 1959. On the dispersion of a solute by diffusion, convection, and exchange between phases.
Proc. Roy. Soc. A, 252, 538-550.
Beigel C., Di Pietro L., 1999. Transport of triticonazole in homogeneous soil columns : influence of nonequilibrium sorption. Soil Sci. Soc. Am. J., 63, 1077-1086.
Benoit P., Perceval J., Stenrød M., Moni C., Eklo O-M., Barriuso E., Sveistrup T., Kværner J., 2007.
Availability and biodegradation of metribuzin in alluvial soils as affected by temperature and soil
properties. Weed Research, 47, 517-526.
Brusseau M.L., 1993. The influence of solute size, pore-water velocity, and intraparticle porosity on solute dispersion and transport in soil. Water Resources research, 29, 1071-1080.
Cameron D.R., Klute A., 1977. Convective-dispersive solute transport with a combined equilibrium andkinetic adsorption model. Water Resour. Res., 13, 183-188.
Grosser K.A., Erickson K.L., Carbonell R.G., 1991. Enhanced dispersion resulting from solute exchange
between phases. AIchE J., 37, 512-526.
Harvey C.F., de Marsily G. 1995. Are the local dispersion coefficients of a sorbing and a non-sorbing
tracer identical? An analysis of Taylor-Aris dispersion and experimental evidence. Draft for publication, pers. comm.
Pot V., Šimůnek J., Benoit P., Coquet Y., Yra A., Martínez-Cordón M.J., 2005. Impact of rainfall intensity on the transport of two herbicides in undisturbed grassed filter strip soil cores. J. ContaminantHydrol., 81, 63-88.
Pot V., Genty A., 2005. Sorbing and non-sorbing solute migration in rough fractures with a multi-species
LGA model: dispersion dependence on retardation and roughness. Transp. Porous Media, 59, 175196.
Pot V., Genty A., 2007. Dispersion dependence on retardation in a real fracture geometry using latticegas cellular automaton. Adv. Water Resour., 30, 273-283.
Prata F., Lavorenti A., Vanderborght J., Burauel P., Vereecken H., 2003. Miscible displacement, sorption
and desorption of atrazine in a Brazilian Oxisol. Vadose Zone Journal, 2, 728-738.
Šimůnek J., Šejna M., van Genuchten M.Th., 1998. The HYDRUS-1D software package for simulating
the one-dimensional movement of water, heat, and multiple solutes in variably-saturated media. Version 2.0, ICWMC-TPS-70. International Ground Water Modeling Center, Colorado School of Mines,
Golden, Colorado. 202 pp.
Šimůnek J., Jarvis N.J., van Genuchten M.Th., Gärdenäs A., 2003. Review and comparison of models fordescribing non-equilibrium and preferential flow and transport in the vadose zone. J. Hydrol., 272,
14-35.
Stenrød M., Perceval J., Benoit P., Almvik M., Bolli R. I., Eklo O.M., Sveistrup T. E., Kværner J., 2007.
Cold climatic conditions; effects on bioavailability and leaching of the mobile pesticide metribuzin in
a
silt
loam
soil
in
Norway.
Cold
Regions
ScienceandTechnology,
doi:10.1016/j.coldregions.2007.06.007.
Taylor G., 1953. Dispersion of soluble matter in solvent flowing slowly through a tube. Proc. Roy. Soc.
A, 219, 186-203.

24

Figure 1. Experimental bromide and metribuzin- BTCs, presented as the relative concentration
against time, and best fit obtained with the CDE model for bromide and the CDE I, CDE II,
Two-site I and Two-site II models for metribuzin.
25

26

Sensitivity Analysis of HYDRUS-1D to transient-MIM
parameters: a case study related to pesticide fate in soil
B. Cheviron§, Y. Coquet
UMR 1091 INRA / AgroParisTech Environment and Arable Crops, BP 01, 78850 Thiverval-Grignon,
France§currently at UMR 1221 INRA/IRD/ Montpellier SupAgro LISAH, 2 place Viala, 34060 Montpellier
cedex 1, France

ABSTRACT: We evaluated version 3.0 of HYDRUS-1D for its capacity to simulate pesticidetransport in soil under transient MIM flow. We present a one-parameter-at-a-time sensitivity
analysis of HYDRUS-1D based on 3 base case scenarios corresponding to three soil profiles, forwhich a comprehensive set of pesticide fate parameters were available, including MIM parameters. The sensitivity analyses were performed using the SENSAN software and were based on
the maximum instantaneous and annually-averaged concentration in soil solution at 1 m depth.
The most sensitive parameters for the three soil profiles were MIM parameters θm,s and θim,s, thevan Genuchten shape parameter n and, more classically, pesticide sorption coefficient and degradation rate. The high sensitivity of pesticide concentration to MIM parameters pleads for the
inclusion of transient MIM processes in pesticide fate models.

1 INTRODUCTION
The second half of the twentieth century has seen an important degradation of groundwater ressources in Europe as a result of the development of intensive agriculture. A large part of thisdegradation is due to contamination by pesticides (Leistra & Boesten, 1989). The situation in
France is no exception. The quality of the Beauce groundwater, a unique national groundwater
resource with a volume estimated to 20 billions of m3 at full capacity, is also degraded by physicochemical contaminations, especially by agricultural pesticides (IFEN, 2004; FREDON,
2007). The recharge zone of the Beauce groundwater aquifer covers 9000 km2 over 6 French
departments (Eure-et-Loir, Loir-et-Cher, Loiret, Yvelines, Essonne, Seine-et-Marne). Annual
uses of this groundwater are estimated to 100-150 millions of m3 for drinking water supply(DWS) and industrial uses, and to 300-500 millions of m3 of irrigation in agriculture.
The research project ESHEL, supported by the French Ministry for Ecology, Sustainable Development and Land Planning, had for main objective the understanding of the processes at theorigin of Beauce groundwater contamination by pesticides. It extended an action plan initiatedby the GREPPPES (Regional Group for the study of pollution par phytopharmaceutical products
in waters and soils for the Centre region) and the Chambre d’Agriculture d’Eure-et-Loir, which
aimed at limitating the contamination of the DWS of the Ouarville commune by agricultural
pesticides. This action plan included a monitoring of pesticide concentrations in many wells andpiezometers located in the water basin of Ouarville (20 km2) from 1999 to 2003, and showed apermanent contamination by a variety of molecules (atrazine, isoproturon, chlortoluron, bentazone, fluroxypyr).
A preliminary work on solute transport in the loamy soils of Beauce had shown that preferential flow through macropores (earthworm or root channels) was ineffective due to the high hydraulic conductivity of the soil matrix. However, flow rates of conservative inert tracers (bromide anion) were not consistent with a movement of all the water contained in the soil. The
research results from the ESHEL project have demonstrated the existence of a « mobile/immobile water » (MIM)-type of transport in Beauce soils (Alletto et al., 2006a). Thesesoils are characterized by a large immobile water fraction, of the order of 50 % of the soil watercontent when close to saturation. This leads to larger solute flow rate than if all the soil water27

was participating to the flow. In consequence, the MIM flow process has to be accounted for inmodelling the risk of pesticide leaching to groundwater in the Beauce region. For that purpose,
we intend to use version 3.0 of the HYDRUS-1D model (Šimůnek et al., 2005), which allows tosimulate MIM solute transport under transient conditions. However, no information is currently
available about the sensitivity of HYDRUS-1D to transient MIM parameters. This paper presents the results of a one-parameter-at-a-time sensitivity analysis of HYDRUS-1D to soil parameters determining pesticide fate. This sensitivity analysis is based on 3 base case scenarios
corresponding to 3 Beauce soils studied in the ESHEL project, and for which a comprehensive
set of pesticide fate parameters, including MIM parameters, have been measured.
2 METHODS AND TOOLS
2.1 HYDRUS-1D
Two ways of modeling transient MIM flow in soils are proposed in HYDRUS-1D (Šimůnek etal., 2003). While the darcian water flow is limited to the mobile region, the exchange of water
between the mobile and the immobile regions, which is compulsory to account for the dynamics
of the partitioning between the two zones under transient conditions, may be described as afunction of local non equilibrium expressed either in terms of effective saturation or in terms ofwater potential.
The effective saturation-based MIM model is written:⎧ ∂θ m∂ ⎡⎛ ∂h ⎞⎤= ⎢ K (h )⎜ + 1⎟⎥ − ω (Sem − Seim )
⎪
⎪ ∂t
∂z ⎣
⎝ ∂z ⎠⎦(1)⎨⎪ ∂θ im = ω (Se − Se )
m
im
⎪⎩ ∂t
where θm and θim [-] are the volumetric water contents in the mobile and immobile water regions, K [L T-1] is the hydraulic conductivity, h [L] is the soil water pressure head, ω [T-1] is thewater transfer rate between the mobile and the immobile regions, Sem and Seim [-] are the effective water saturations in the mobile and immobile regions, z [L] denotes the position on the vertical axis taken positive upwards, and t [T] is the time. The effective water saturations in themobile and immobile regions are defined as:θ m − θ m ,r
⎧
Se
=
m
⎪
θ m , s − θ m ,r⎪⎨⎪Se = θ im − θ im,r⎪ im θ im ,s − θ im ,r⎩(2)

where θm,r and θm,s [-] are the residual and saturated volumetric water contents in the mobile region while θim,r and θim,s [-] are the residual and saturated volumetric water contents in immobile
region. Effective saturation of the mobile region and soil hydraulic conductivity are controlled
by soil water potential according to the Mualem-van Genuchten formulation (van Genuchten,1980):[]⎧Se = 1 + (α h )n − (1−1 / n )⎪ m⎨
2⎪⎩ K (h ) = K s Sem l 1 − (1 − Sem )1 / (1−1 / n ) (1−1 / n )

[(

)

]

(3)

where α [L-1] and n [-] are van Genuchten’s shape factors, Ks [L T-1] is the soil saturated hydraulic conductivity and l [-] is the Mualem exponent related to pore tortuosity and connectivity.
The water potential-based MIM model is written:28⎧ ∂θ m∂ ⎡⎛ ∂h⎞⎤
= ⎢ K (hm )⎜ m + 1⎟⎥ − ω * (hm − him )
⎪
⎪ ∂t
∂z ⎣
⎝ ∂z
⎠⎦
⎨
⎪ ∂θ im = ω * (h − h )
m
im
⎪⎩ ∂t

(4)

where hm and him [L] are the water potential head in the mobile and immobile water regions, andω* [T-1] is the water transfer rate between the mobile and the immobile regions. In this formulation, the water exchange between the immobile and the mobile regions is driven by difference ofwater potential between the two regions. This implies that the two regions have different water
retention characteristics:[]⎧Se = 1 + (α h )nm − (1−1 / nm )
m m
⎪ m
2
⎪
l
1 / (1−1 / nm ) (1−1 / nm )
(
)
(
)
K
h
=
K
Se
−
−
Se
1
1
⎨
m
s
m
m
⎪
nim −(1−1 / nim )⎪Seim = 1 + (α im him )
⎩

[

[(

)

]

]

(5)

Depending on the choice of the formulation, the modelling of MIM hydraulic flow in soil will
necessitate:- 9 parameters for the Se-based formulation: θm,s, θm,r, θim,s, θim,r, ω, α, n, Ks, l- 11 parameters for the h-based formulation: θm,s, θm,r, θim,s, θim,r, ω*, αm, nm, αim, nim,
Ks, lLike Köhne et al. (2004), we decided to use the formulation that required less parameters: theSe-based one.
For modelling the fate of pesticides in soil, we had to consider solute sorption and decay inaddition to transport:⎧ ∂ (θ m cm ) ∂ (ρfsm ) ∂ ⎛∂c ⎞ ∂q c+= ⎜θ m Dm m ⎟ − m m − μθ m cm − Γs
⎪
∂t
∂z ⎝
∂z ⎠
∂z
⎪ ∂t
⎪⎪ ∂ (θ im cim ) ∂ (ρ (1 − f )sim )+= − μθim cim − Γs
⎨
∂t
⎪ ∂t
⎪
⎧Γwcm si Γw > 0⎪Γs = ωs (cm − cim ) + ⎨⎪⎩⎩Γwcim si Γw < 0(6)where cm and cim [ML-3] are the solute concentration in the mobile and immobile water regions,ρ [ML-3] is the bulk density, sm and sim [-] are the solute concentration adsorbed to the solid
phase in the mobile and immobile water regions, f [-] is the fraction of the sorption sites locatedin the mobile region, qm [L T-1] is the flux density in the mobile region, Γs [ML-3T-1] is the solutetransfer rate between the mobile and the immobile regions, μ [T-1] is the first-order solute degradation rate, ωs [T-1] is the first-order solute transfer rate coefficient between the mobile and immobile region due to diffusion, and Dm [L2T-1] is the dispersion coefficient in the mobile region
defined as:
7

λ q
θ 3
Dm = m m + D0 m 2
θm
θ m ,s(7)

where λm [L] is the dispersivity of the mobile region, and D0 [L2T-1] is the molecular diffusion
coefficient.
Due to numerical convergence problems, we had to consider a linear sorption isotherm instead of a Freundlich isotherm:s = kd ⋅ c(8)29where kd [L3M-1] is the linear sorption coefficient. Sorption coefficient was assumed to be similar in both mobile and immobile regions. Pesticide degradation was described by a first orderkinetics. Pesticide degradation rates, μ, were considered to be depending on soil temperature andwater content. While temperature effects were modelled by the Arrhenius equation, water content effects were implemented through Walker (1974)’s equation, which assumes that μ is proportional to water content:⎛ θ ⎞⎟μ = μ ref ⎜⎜⎟θ⎝ ref ⎠b(9)

where μref [L3M-1] is the degradation rate measured at some reference water content θref., and b isan empirical parameter.
2.2 Base case scenarios
All the horizons of three soil profiles from the Ouarville catchment (Table 1) were fully characterized for their hydraulic, MIM and pesticide fate parameters through independent measurements.
Table 1. Description of the three Beauce soil profiles.

_____________________________________________________________________________

Site 1Site 2_________________________________

Site 3_________________

HorizonDepth Horizon Depth Horizon Depth
(cm)
(cm)
(cm)
_____________________________________________________________________________
Seedbed
LA1-1
0-8
LA1-2
0-8
LA1-3
0-8
Plough
layer
LA2-1
8-28
LA2-2
8-28
LA2-3
8-28
_____________________________________________________________________________
Clay-enriched
BT1-1
28-48 BT-2
28-68 BT-3
28-68
layers
BT2-1
48-78
B-3
68-90
Weatheredlimestone
C-1
78C-2
68C-3
90_____________________________________________________________________________

Hydraulic parameters θs, θr, α, n, Ks, l were estimated using Wind evaporation experiments,
while setting l values to 0.5. θm,s, θm,r, θim,s, θim,r were calculated from θs, θr values assuming θim,rto be equal to zero and using the mobile water fraction θ/θm measured at h = -1 cm with tension
infiltrometers and anionic tracers (Alletto et al., 2006a) to calculate θm,s and θim,s. Water exchange parameter ω was taken from Köhne et al. (2004). Bulk density was measured by the cylinder method. Solute exchange coefficient ωs was measured using the technique of Jaynes et al.
(1995) (Alletto et al., 2006a). Linear sorption coefficients of isoproturon, a substituted-urea herbicide used on cereal crops, were obtained from batch experiments and degradation rates fromlaboratory incubations (Alletto et al., 2006b). Degradation experiments made at different temperatures and water contents were used to derive the parameters of the Arrhenius and Walker
equations (Alletto et al., 2006b). Dispersivity λm was calculated from a field bromide tracer experiment and was assumed to be identical for all the horizons of each soil profile. All parametervalues are listed in Table 2.

30

Table 2. Reference values of HYDRUS-1D parameters for the horizons of the three soil profiles.
parametercode

theta m,rtheta m,stheta im,rMIM
theta im,somega m->im
mobile fractionalpha vG
VAN
n vG
GENUCHTEN
L vG
Ks
diffusivity in water
kd
beta
mu
SOLUTEalpha
walker bea (kd)ea (beta)
ea (mu)bulk density
dispersivity
PHYSICS
solid fraction
organic fraction
thermal dispersivity

m1
m2
m3
m4
m5
m6
v1
v2
v3
v4
s1
s2
s3
s4
s5
s6
s7
s8
s9
p1
p2
p3
p4
p5

site 10,00000,47840,00000,1716
1,30
0,74
0,42
1,23
0,5
156,0
1,797
5,00
1
0,06080,0192
0,68
45000
45000
45000
1,30
2,30
0,35
0,01
5,0

layer 1
site 20,00000,23580,00000,3832
1,30
0,38
0,60
1,21
0,5
213,6
1,797
5,00
1
0,06080,0192
0,68
45000
45000
45000
1,10
1,5
0,38
0,02
5,0

site 30,00000,24840,00000,4416
1,30
0,36
0,37
1,34
0,5
528,0
1,797
5,00
1
0,06080,0192
0,68
45000
45000
45000
1,20
2,10
0,30
0,01
5,0

layer 2
site 20,00000,13200,00000,2680
1,30
0,33
4,80
1,10
0,5
633,6
1,797
4,00
1
0,05780,1320
0,68
45000
45000
45000
1,20
1,5
0,43
0,01
5,0

site 10,00000,32900,00000,1410
1,30
0,65
2,00
1,25
0,5
160,0
1,797
4,00
1
0,05780,1320
0,68
45000
45000
45000
1,35
2,30
0,50
0,01
5,0

site 30,00000,14850,00000,3015
1,30
0,33
0,29
1,42
0,5
528,0
1,797
4,00
1
0,05780,1320
0,68
45000
45000
45000
1,45
2,10
0,55
0,01
5,0

site 10,00000,27050,00000,2295
4,50
0,55
0,21
1,32
0,5
158,4
1,797
0,90
1
0,00620,0360
0,53
45000
45000
45000
1,41
2,30
0,53
0,01
5,0

layer 3site 20,14000,25900,06000,1110
4,50
0,70
2,49
1,25
0,5
604,8
1,797
1,00
1
0,00420,0288
0,53
45000
45000
45000
1,35
1,5
0,47
0,00
5,0

site 30,17520,33100,09480,1790
4,50
0,65
1,62
1,24
0,5
233,0
1,797
1,00
1
0,00620,0240
0,53
45000
45000
45000
1,45
2,10
0,49
0,01
5,0

site 1
0,07600,18800,11400,2820
4,50
0,40
2,00
1,16
0,5
600,0
1,797
0,70
1
0,0042
0,0552
0,53
45000
45000
45000
1,41
2,30
0,53
0,01
5,0

layer 4
site 20,07600,20120,0094
0,2489
4,50
0,45
0,48
1,30
0,5
1106,41,797
0,30
10,0034
0,6528
0,38
45000
45000
45000
1,48
1,5
0,55
0,00
5,0

site 30,2506
0,36890,10940,1611
4,50
0,70
2,78
1,16
0,5
205,0
1,797
0,60
1
0,0047
0,0288
0,53
45000
45000
45000
1,20
2,10
0,47
0,00
5,0

layer 5site 1
site 30,0300 0,0338
0,0990 0,0783
0,0100 0,1562
0,3280 0,3617
4,50
4,50
0,23
0,28
3,03
5,45
1,22
1,20
0,5
0,5
1555,2 383,0
1,797
1,797
0,50
0,30
1
1
0,0065 0,0007
0,5808 0,1440
0,38
0,38
45000 45000
45000 45000
45000 45000
1,49
1,40
2,30
2,10
0,53
0,56
0,01
0,00
5,0
5,0

Units : ω (d-1), α (cm-1), Ks (cm d-1), D0 (cm2 day-1), kd (g cm-3), μ (d-1), ωs (solute alpha) (d-1), Ea (activation energy in Arrhenius equation) (J), bulk density (g cm-3), λm (cm), thermal dispersivity (cm).

2.3 Sensitivity analyses
We used the SENSAN software (Doherty, 2004) to predefine parameter sets for HYDRUS andautomate HYDRUS runs until all parameter sets had been tested. Batch files were created to activate FORTRAN routines to count, record and reformat HYDRUS runs and outputs to ad hoc
formats and files, and extracting concentration values simulated at 1 m depth, which was the
only output variable considered in the sensitivity analyses. Some constrains were put on parameter combinations for physical reasons:⎧θ m ,r < θ m,s⎨⎩θ im ,r < θ im,s

(10)

Range of variations were determined according to measured local variability of the parameters
or based on a priori estimation of this variability or based on measurement accuracy (Table 3).
Table 3. Range of variation of HYDRUS-1D parameters.
parametercode

theta m,rtheta m,stheta im,rtheta im,somega m->im
mobile fractionalpha vG
n vG
LvG
Ks
diffusivity in water
kd
beta
mu
alpha
walker bea (kd)ea (beta)
ea (mu)bulk density
dispersivity
solid fraction
organic fraction
thermal dispersivity

m1
m2
m3
m4
m5
m6
v1
v2
v3
v4
s1
s2
s3
s4
s5
s6
s7
s8
s9
p1
p2
p3
p4
p5

SITE 1layers 1-2
layers 3-4-5
fixed
0.20 - 1.80
0.20 - 1.58
0.32 - 1.79
fixed
0.20 - 1.80
0.20 - 1.80
0.28 - 1.79
0.10 - 5.00
0.11 - 4.96
0.71 - 1.30
0.70 - 1.30
0.80 - 1.20
0.90 - 1.10
-4.00 - 4.00
0.25 - 4.00
0.26 - 3.97
0.65 - 1.35
0.6 - 1.40
0.21 - 1.79
fixed
0.60 - 1.40
0.09 - 1.90
0.10 - 5.00
0.11 - 4.96
0.5 - 1.5

SITE 2layers 1-2
layers 3-4-5
fixed
0.20 - 1.80
0.20 - 1.80
0.47 - 1.80
fixed
0.20 - 1.80
0.20 - 1.80
0.46 - 1.80
0.03 - 32.23
0.03 - 32.00
0.70 - 1.30
0.80 - 1.20
0.90 - 1.10
-4.00 - 4.00
0.25 - 4.00
0.65 - 1.35
0.60 - 1.40
0.20 - 1.80
fixed
0.60 - 1.40
0.13 - 1.86
0.10 - 5.00
0.11 - 1.87
0.5 - 1.5

SITE 3layers 1-2
layers 3-4-5
0.21 - 1.79
fixed
0.50 - 1.79
0.20 - 1.80
0.21 - 1.79
fixed
0.20 - 1.80
0.20 - 1.80
0.11 - 4.96
0.10 - 5.00
0.70 - 1.30
0.80 - 1.20
0.90 - 1.10
-4.00 - 4.00
0.26 - 3.97
0.25 - 4.00
0.65 - 1.35
0.21 - 1.79
0.60 - 1.40
fixed
0.09 - 1.90
0.60 - 1.40
0.11 - 4.96
0.10 - 5.00
0.5 - 1.50.80 - 1.20

0.80 - 1.20

0.80 - 1.20

0.85 - 1.15
0.20 - 5.00
0.33 - 3.00
0.90 - 1.10
0.73 - 1.27
0.67 - 1.33
0.50 - 1.50

0.85 - 1.15
0.33 - 3.00
0.90 - 1.10
0.72 - 1.28
0.50 - 1.50

0.85 - 1.15
0.33 - 3.00
0.20 - 5.00
0.90 - 1.10
0.71 - 1.29
0.50 - 1.50

Legend : 0.50 - 1.25 -> a variation between 50% and 125% of the reference value has been tested

In the sensitivity analyses, we created two groups of horizons: the topsoil group including theseedbed and the plough layer, and the subsoil group including all the horizons below the plough
layer. This choice was motivated by the well-known importance of topsoil properties in pesticide dissipation.
31

3 RESULTS
Figure 1 gives some examples of the sensitivity of the 1 m depth concentration to MIM parameters for the soil profile of site 2.
45

35
0.2 refref1.8 ref

SITE 2 - MIM Parameters
30 theta m,r (layers 3-4)
25
20

SITE 2 - MIM Parameterstheta m,s (layers 1-2)

40

0.2 refref1.8 ref35
30µg / lµg / l

25
20

15

15

10

10
5

5
days

days

0
365

730

1095

1460

1825

2190

2555

2920

0
365

730

1095

1460

1825

2190

2555

2920

100

35
0.2 refref1.8 ref

SITE 2 - MIM Parameters
30 theta im,r (layers 3-4)

SITE 2 - MIM Parameterstheta im,s (layers 1-2)

90

0.2 refref1.8 ref80

25

70µg / lµg / l

60

20

50
15

40
30

10

20
5

10

days

0
365

45

730

1095

1460

1825

2190

2555

0
365

2920
35

SITE 2 - MIM Parameters0.03 refref32.23 ref40 omega m->im (layers 1-2)
35

days

730

1095

1460

1825

2190

SITE 2 - MIM Parameters

2555

2920

0.7 refref1.3 ref

30 mobile fraction (layers 1-2)2530

µg / lµg / l

25

20

20

15

15

10

10
5

5

days

days

0
365

730

1095

1460

1825

2190

2555

2920

0
365

730

1095

1460

1825

2190

2555

2920

The most sensitive parameters are θm,s and θim,s. Water transfer rate ω is less sensitive. We
used the Ratio Of Variation (ROV) to estimate the sensitivity of our target output (maximum
concentration at 1 m depth, C) to variations of any input parameter p:

ROV =C − Crefp − pref

Cref

pref

(11)

We calculated the maximum mean annual ROV values for each site and calculated averagevalues between the three sites and the corresponding standard deviations. Table 4 shows the results for the maximum and minimum input values for all parameters, sorted according to their
sensitivities.
Table 4. Sensitivities of pesticide concentration at 1 m depth to input parameters.

32

MINIMUM ENTRIESPARAMETER LAYERS
1-2
theta im,s
3+
kd
1-2
mu
1-2
kd
3+
mu
1-2
theta m,s
1-2
bulk density
3+
theta im,s3+
theta m,s1-2alpha mim
3+bulk density
1-2
omega mim3+
dispersivity
3+
theta m,r3+
theta im,r3+
thermal disp.
1-2
solid fraction
3+
alpha vG3+solid fraction
3+
L vG
1-2
L vG
1-2
organic fraction
3+organic fraction
3+
omega mim
1-2
diffusivity in water
3+
alpha mim3+
mobile fraction
1-2
dispersivity
1-2
thermal disp.
1-2alpha vG
3+
Ks
1-2
Ks
1-2
mobile fraction
3+
walker b
1-2
walker b
1-2
ea (kd)1-2ea (mu)
3+
n vG
1-2
n vG

ROV
-2,32
-1,99
-1,93
-1,63
-1,08
-0,96
-0,65
-0,46
-0,42
-0,41
-0,21
-0,19
-0,15
-0,14
-0,09
-0,07
-0,06
-0,05
-0,05
-0,02
-0,01
0,00
0,00
0,00
0,01
0,02
0,04
0,05
0,07
0,11
0,11
0,15
0,17
0,24
0,37
0,37
1,03
1,77
2,00

MAXIMUM ENTRIESS-DEVPARAMETER LAYERS
0,97
1-2
mu
0,60
1-2
kd
0,33
1-2
theta m,s0,36
1-2
theta im,s
0,47
3+
mu
0,83
3+
theta m,s
0,17
1-2
bulk density
0,43
3+
theta im,s
0,19
3+
kd
0,29
3+
solid fraction
0,19
3+
theta m,r
1-2
0,17 diffusivity in water
0,17
3+
theta im,r
0,08
1-2
solid fraction
0,05
3+alpha vG
0,01
3+
thermal disp.
0,04
1-2
organic fraction
0,09
3+
L vG
0,05
3+
alpha mim
0,01
3+
dispersivity
0,01
1-2
alpha mim
0,00
3+
organic fraction
0,01
1-2
omega mim
0,02
3+
omega mim
0,03
1-2
L vG
0,03
3+
Ks
0,05
3+
bulk density
0,19
1-2
thermal disp.
0,04
1-2
Ks
0,12
3+
mobile fraction
0,04
1-2
alpha vG
0,02
1-2
dispersivity
0,02
1-2
mobile fraction
0,08
3+
walker b
0,20
1-2
walker b
0,16
1-2
ea (kd)0,111-2
ea (mu)
0,68
3+
n vG
1,24
1-2
n vG

ROV
-0,70
-0,63
-0,56
-0,49
-0,41
-0,35
-0,34
-0,34
-0,25
-0,11
-0,09
-0,09
-0,08
-0,05
-0,04
-0,04
-0,02
-0,02
-0,02
-0,02
-0,01
-0,01
0,00
0,00
0,01
0,03
0,04
0,05
0,06
0,07
0,11
0,11
0,18
0,19
0,35
0,41
1,26
1,37
3,95

S-DEV
0,10
0,43
0,24
0,49
0,20
0,18
0,15
0,16
0,04
0,02
0,09
0,18
0,09
0,04
0,07
0,00
0,01
0,00
0,07
0,01
0,00
0,01
0,01
0,00
0,01
0,02
0,10
0,04
0,02
0,04
0,12
0,05
0,05
0,08
0,19
0,19
0,13
0,79
0,18

The most sensitive input parameters are pesticide fate parameters kd and μ but also the n parameter of the van Genuchten function, for all horizon groups. The activation energy parameterfor the degradation rate is also very sensitive, showing the importance of temperature effects on
simulated pesticide concentrations. MIM parameters θm,s and θim,s are also very sensitive,
thereby showing the importance of MIM processes in the simulation of pesticide leaching togroundwater.
4 CONCLUSION
We presented a sensitivity analysis of HYDRUS-1D in its transient-MIM version. The framework of this analysis was the prediction of pesticide leaching to groundwater. Three base casescenarios were defined using independent measurements of the HYDRUS-1D parameters. The
results show the high sensitivity of pesticide concentrations at 1 m depth to pesticide fate parameters, but also to the n parameter of the van Genuchten water retention function and to MIM
parameters. Sensitivity to the mobile and immobile saturated water contents was particularly
important. These results show the importance of including MIM processes in the models used toevaluate pesticide leaching risks to groundwater.
REFERENCES
Alletto, L., Coquet, Y., Vachier, P., & Labat, C. 2006a. Hydraulic conductivity, immobile water content,
and exchange coefficient in three soil profiles, Soil Sci. Soc. Am. J. 70:1272-1280.
Alletto, L., Coquet, Y., Benoit, P., & Bergheaud, V., 2006b. Effects of temperature and water content on
degradation of isoproturon in three soil profiles, Chemosphere 64:1053-1061.
FREDON (Fédération régionale de défense contre les organismes nuisibles de la région Centre), 2007.
Les réseaux de surveillance des eaux souterraines et superficielles du Grepppes. Synthèse des résultats de 2000 à 2006. FREDON/GREPPPES, Orléans, France.
IFEN (Institut Français de l’Environnement), 2004. L’environnement en région Centre. IFEN, Paris,
France, 173 p.

33

Jaynes, D., Logsdon, S., & Horton, R. 1995. Field method for measuring mobile/immobile water content andsolute transfer rate coefficient. Soil Sci. Soc. Am. J. 59:352–356.
Köhne, J.M., Köhne, S., Mohanty, B.P., & Šimůnek, J. 2004. Inverse mobile-immobile modeling of transport during transient flow: effects of between-domain transfer and initial water content. Vadose Zone
J. 3:1309-1321.
Leistra, M. & Boesten, J. J. T. I. 1989. Pesticide contamination of groundwater in western Europe. Agric. Ecosyst. Envir. 26: 369-389.
Šimůnek, J., van Genuchten, M.Th., & Sejna, M., 2005. The HYDRUS-1D Software Package for Simulating the Movement of Water, Heat, and Multiple Solutes in Variably Saturated Media, Version 3.0.
HYDRUS Software Series 1, Department of Environmental Sciences, University of California Riverside, Riverside, California, USA.
Šimůnek, J., Jarvis, N.J., van Genuchten, M.Th., & Gärdenäs, A. 2003. Review and comparison of models for describing non-equilibrium and preferential flow and transport in the vadose zone. J. Hydrol.
272:14-35.
van Genuchten, M.Th. 1980. A closed form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892-898.
Walker, A. 1974. A simulation model for prediction of herbicide persistence. J. Environ. Qual. 3: 396401.

34

Evaluating sorption and degradation characteristics of pesticides
using column displacement experiments
T. De Wilde1, J. Mertens2, J. Šimůnek3, K. Sniegowski2, P. Spanoghe1, J. Ryckeboer2,
D. Springael3
1

Laboratory of Crop Protection Chemistry, Ghent University, Coupure Links 653, B-9000 Ghent,
Belgium, 2Division Soil and Water Management, K.U. Leuven, Kasteelpark Arenberg 20, B-3001
Heverlee, Belgium and 3Department of Environmental Sciences, University of California Riverside,
Riverside, CA 92521, USA

ABSTRACT: HYDRUS-1D was applied to analyze breakthrough curves (BTCs) generated incolumn displacement experiments. These column displacement experiments represent small
scale biopurificiation systems to treat pesticide contaminated water. Pesticides were retainedand degraded inside the column. For certain pesticides (i.e. metalaxyl), degradation was not a
first-order process, therefore the convection-dispersion equation was augmented with the
Monod type kinetics.

1 INTRODUCTION
40-90% of surface water contamination is attributable to point sources or direct losses (Carter,
2000; Decoin, 2003). Point source pollution results, for example, from spills at filling operations, leakages of the spray equipment, spray leftovers, spills of rinsing water from cleaning ofinternal tanks, and spills of water from external cleaning (Torstensson & Castillo, 1997; Ramwell et al., 2004). Therefore biopurification systems, introduced by Torstensson and Castillo(1997) in Sweden, have been developed to clean pesticide contaminated water which is generated on farms. This system uses natural materials (agricultural waste products or products available on-farm) as sorbents to clean the pesticide contaminated water. Retention of pesticides inthese systems does not depend solely on the structure and composition of the organic substrates,
but also on the physicochemical characteristics of the pesticide. The efficiency of these systems
is much lower in retaining and degrading mobile pesticides, such as bentazon, metsulfuronmethyl, and mecoprop, than less mobile pesticides such as linuron (De Wilde et al. 2007). At
present, biopurification systems are operated as a black box, since no research has yet been carried out to sufficiently understand processes taking place inside of these systems. In order to optimize the efficiency of these systems for a broad range of pesticides and organic substrates, thefate of pesticides, and the effects of degradation and retention processes, inside of these systems
need to be well characterized. This study focuses on the fate of four pesticides, i.e., linuron, bentazon, metalaxyl and isoproturon, in small scale biopurification systems. Transport of these pesticides was evaluated in column displacement experiments containing different mixtures of organic materials. Different matrix substrates that can be used in biopurification systems and thatwere evaluated in this study have different contents and characteristics of the organic carbon andcould thus influence pesticide sorption and degradation inside of the soil column. The onedimensional transport model HYDRUS-1D (Simunek et al., 2005) was used to identify andquantify solute transport and hydraulic parameters.

35

2 MATERIALS AND METHODS
2.1 Pesticide and susbtrate properties
Linuron was selected as a rather immobile pesticide (Koc = 410 L kg-1), metalaxyl (Koc = 47 Lkg-1) and isoproturon (Koc = 36 L kg-1) as intermediately mobile pesticides, and bentazon (Koc =13 L kg-1) as a mobile pesticide. The following organic substrates were selected for this study:
peat mix, garden waste compost, straw, sandy loam soil, cow manure and coco chips.
2.2 Displacement experiments
Laboratory columns were packed in triplicates with five different mixtures of air dried organic
substrates and sandy loam soil (Table 1). The columns were 15 cm high and had an inner diameter of 10 cm. The effluent was collected at the bottom in a fraction collector. A peristaltic pumpdelivered a constant Darcy flux of 1.74 cm d-1. Bromide in the form of KBr was used as a nonreactive tracer in the column experiments to determine physical transport parameters. The bromide solution (0.1 mM Br-) was injected into the columns as a pulse of 18.3 h. Pesticides werepumped onto the column together with the bromide solution. The pesticide solution pumped
onto the column contained 0.001 M CaCl2 and 10 mg L-1 linuron, isoproturon, metalaxyl andbentazon. Pesticides were added continuously as a step input. The effluent was collected every
2-3 days and both outflow volumes and pesticide concentrations were measured.
Table 1: Composition of the laboratory columns
Mix 1
Garden waste compost
/
Straw
50
Coco chips/
Peat mix
25
Sandy loam soil
25
Cow manure
/

Mix 2
25
50
/
/
25
/

Mix 3
/
50
/
20
25
5

Mix 4
/
25
25
25
25
/

Mix 5
/
25
25
20
25
5

2.3 Transport models
The HYDRUS-1D model that simulates one-dimensional water flow and transport of solutes insoils was used to describe the transport of pesticides (Šimůnek et al., 2005). It was assumed thatexperimental breakthrough curves (BTCs) could be described using the classical convectiondispersion transport model (Lapidus & Amundson, 1952) that neglects both physical and chemical equilibrium. Assuming that decay is active only in the liquid phase and a non-linear relationship between liquid and solid phase concentrations can be described using the Freundlich adsorption isotherm: s = Kf cn, where c [M L-3] and s [M M-1] are concentrations in the liquid andsolid phases, respectively, and Kf and n are Freundlich parameters, the transport of a reactive
solute for steady-state water flow conditions can be written as follows:

R

∂c
∂ ²c
∂c
=D
−ν
− μl c
∂t
∂z ²
∂z

(1)

where D is the dispersion coefficient [L² T-1], ν is the pore water velocity (L T-1), ν = q/θ inwhich q is the Darcian water flux [L T-1] and θ is the volumetric water content [L³ L-3], µl is thefirst-order degradation constant for the solute in the liquid phase [T-1], t [T] and z [L] are thetemporal and spatial coordinates, respectively and R is the retardation factor [-], in which ρb is
the bulk density [M L-3]:

R = 1+

ρb⋅ K f ⋅ n ⋅ c n −1θ(2)

In the model described above (referred to below as the CDE model), pesticide degradation isassumed to be a first-order process depending only on the pesticide concentration. However, thisis in contrast with a commonly observed lag phase (a start-up period, during which degradation36

is not detected) in laboratory mineralization experiments. Therefore, we hypothesize that thesimplified version of the Monod kinetics (Guimont et al., 2005) needs to be used to describe
BTCs where a lag phase was clearly present (e.g., metalaxyl BTCs) (Mertens et al., 2008). In
such model, bacteria growth can be described using the following Monod kinetics:

dX ⎛
c ⎞
= ⎜ μm⎟ X − kdecay X
dt ⎝Ks + c ⎠(3)

where X is the biomass concentration [Mb Lw-3], c is the liquid pesticide concentration [Mp Lw-3],µm is the mass growth rate [T-1], Ks is the half saturation constant [Mp Lw-3], and kdecay is the decay rate [T-1] (subscripts b, p and w refer to biomass, pesticide and water, respectively). The bacteria Monod growth function can be related with the pesticide consumption utilizing a yield coefficient, Y, where Y is a measure of the organisms formed per pesticide utilized [Mb Mp-1]. The
local change in the pesticide concentration when neglecting the effects of the transport can thenbe expressed as follows:

dc
1⎛
c ⎞
= − ⎜ μm
⎟Xdt
Y⎝
Ks + c ⎠(4)

If we can assume that pesticide concentration (c) is significantly lower than the half saturation
constant (c <<Ks), then (3) and (4) simplify as follows:dX ⎛ μm ⎞=⎜c ⎟ X − k decay X = ( μm* c ) X − k decay X
dt ⎝ K s ⎠

(5)

⎛ μ
⎞
dc
= − ⎜ m X ⎟ c = − ( μ* X ) c
dt⎝ YK s ⎠(6)

where µm/Ks is the modified mass growth rate [Lw3Mp-1T-1] and µ* has units of [Lw3Mb-1T-1]. Incorporating Eq. (6) into the one-dimensional transport equation (1), results in:

R

∂c
∂ ²c
∂c
=D
−ν
− ( μ * X )c − μl c
∂t
∂z ²
∂z

(7)

These equations (referred to below as the Monod model) have been implemented in HYDRUS-1D. The fitting parameters for this model are µ * , µm* (yield Y) and kdecay.
3 RESULTS AND DISCUSSION
The physical transport was evaluated first by analyzing the bromide BTCs. The volumetric water content, θ, and the longitudinal dispersivity, λ (D= λ v), were fitted to the observed Br BTCs.
Once the physical transport was fully characterized, chemical processes for isoproturon,
metalaxyl, linuron and bentazon were analyzed by inverting the pesticide BTCs, using thephysical transport parameters estimated for the bromide transport. In all inversions of the pesticide BTCs, the soil water content, θ, the water flux, q, the bulk density, ρ, and the dispersivity,
λ, were fixed. Estimated reaction parameters included the Freundlich parameters Kf and n, andthe liquid degradation constant µl.
3.1 Bromide
Fitted transport parameters are presented in Table 2. High water contents, likely caused by thehigh water sorbing capacity of the organic substrates such as peat mix and coco chips, were present in all columns.

37

Table 2. Estimated physical transport parameters (± 95% confidence interval).
Mix
1
2
3
4
5

θ
0.58 ± 0.02
0.54 ± 0.02
0.62 ± 0.02
0.54 ± 0.02
0.63 ± 0.03

λ

R²

cm
1.97 ± 0.28
1.41 ± 0.18
2.54 ± 0.28
4.03 ± 0.45
2.03 ± 0.39

0.90
0.92
0.89
0.92
0.77

3.2 Isoproturon
The fitted solute reaction parameters are presented in Table 3. The CDE model described wellthe BTCs of isoproturon, as is reflected by high values of the determination coefficient R² (0.98– 1.00). A significant shift of BTCs to the right compared to bromide indicated that the pesticide
is strongly adsorbed to the organic substrates, as is also reflected by the retardation coefficient R(Table 3). A comparison of the Freundlich parameters obtained using column transport experiments with those obtained using batch experiments (De Wilde et al., 2008) showed that Kf values were about 2.4-6.9 times smaller for the column experiments. This means that isoproturon ismore strongly sorbed to the organic substrates in batch experiments and that using these coefficients in transport models would overestimate sorption. Relatively small differences exist in
degradation between different columns.
Table 3. Solute reaction parameters fitted to isoproturon BTCs (± 95% confidence interval).
Mix
1
2
3
4
5

Kf
L kg-18.00 ± 0.39
3.70 ± 0.19
4.93 ± 0.03
5.11 ± 0.07
2.05 ± 0.03

µl (x 10-3)
h-1
2.23 ± 0.17
2.40 ± 0.12
1.70 ± 0.08
1.87 ± 0.15
2.01 ± 0.11

n1.02 ± 0.04
1.31 ± 0.02
1.04 ± 0.03
1.04 ± 0.07
1.32 ± 0.08

R

R²

12.07 ± 2.01
15.39 ± 3.37
7.20 ± 0.86
7.97 ± 1.72
7.36 ± 2.14

0.99
0.99
1.00
0.99
0.98

3.3 Bentazon
Experimental BTCs of bentazon were again analyzed using the CDE model (Table 4) with a determination coefficient R² from 0.72 to 0.92. In general, the retardation factors R were more orless equal to 1, indicating that hardly any retention occurred for all mixtures. The fitted n values
for most mixtures were relatively high. The n and Kf values are determined from the changing
slope of the part of the BTC curve were C/C0 increases from 0 to equilibrium conditions. Since
bentazon breakthrough occurs very fast, the number of data points in this part of the curve islimited. This results in a poor estimation of the n and Kf value. Degradation of bentazon was
about 10 times slower than for isoproturon and was almost negligible.
Table 4. Solute reaction parameters fitted to bentazon BTCs (± 95% confidence interval).
Mix
1
2
3
4
5

Kf
L kg-10.00 ± 0.00
0.02 ± 0.02
0.01 ± 0.02
0.01 ± 0.03
0.02 ± 0.02

µl (x 10-4)
h-1
3.11 ± 4.64
5.18 ± 6.43
1.82 ± 2.83
3.25 ± 1.49
0.00

n
4.78 ± 2.14
2.55 ± 0.41
3.15 ± 0.85
3.61 ± 1.06
2.02 ± 0.43

38

R

R²

1.10 ± 0.17
1.00 ± 0.00
1.03 ± 0.02
1.27 ± 0.33
1.06 ± 0.10

0.90
0.91
0.72
0.85
0.92

3.4 Metalaxyl
Although a step input was applied on the top of the columns, the BTCs for mixtures 2, 3 and 5
showed a pattern that is characteristic for a pulse input (Figure 1). We hypothesized that these
BTCs can be described using a Monod-type degradation model, that considers bacteria growthand decay. This behavior of pesticides can be described by considering the simplified Monod
kinetics in the CDE model. The simplified Monod kinetics (eq. (5) and (7)) were implemented
into HYDRUS-1D. This model has the following additional reaction parameters (µ*, µm* (Y) andkdecay) that needs to be either measured or estimated from collected breakthrough data. The classical CDE model (1) was first used to analyze Metalaxyl BTCs for all mixtures. The determination coefficient R² for the CDE model ranged from 0.64 to 0.99 (Table 6). First-order degradation constants for metalaxyl (Table 5) were of the same magnitude as for isoproturon. The
Monod kinetics parameters and the yield coefficients determined by fitting the entire experimental data curves are presented in Table 6.
Mix 3

Mix 2
1 ,0

0 ,8

0 ,8

A c a lc u la te dA e x p e rim e n ta l
B c a lc u la te d
B e x p e rim e n ta l
C c a lc u la te dC e x p e rim e n ta l

0 ,6

C/C0

C/C0

0 ,6

0 ,4

0 ,4
A c a lc u la te dA e x p e rim e n ta l
B c a lc u la te d
B e x p e rim e n ta l
C c a lc u la te dC e x p e rim e n ta l

0 ,2

0 ,2

0 ,0

0 ,0
0

5

10

15

20

0

PV

5

10

15

20

PV

Mix 4
1 ,0

0 ,8

C/C0

0 ,6

0 ,4
A c a lc u la te dA e x p e rim e n ta l
B c a lc u la te d
B e x p e rim e n ta l
C c a lc u la te dC e x p e rim e n ta l

0 ,2

0 ,0
0

5

10

15

20

PV

Figure 1: Metalaxyl BTCs measured and calculated for 3 columns with delayed degradation (in triplicate:
A, B, C) Symbols are the experimental data, full, dashed and dotted lines are the fits given by the modified convection-dispersion equation. Relative concentrations are presented against the number of porevolumes (PV).
Table 5. Solute reaction parameters fitted to metalaxyl BTCs. (± 95% confidence interval).
Mix
1
2
3
4
5

Kf
L kg-1
3.39 ± 0.37
2.26 ± 0.53
4.55 ± 0.15
1.37 ± 0.28
2.82 ± 0.44µl (x 10-3)
h-1
1.38 ± 0.12
3.23 ± 0.12
3.44 ± 0.05
1.26 ± 0.09
4.26 ± 0.23

n1.17 ± 0.05
1.02 ± 0.11
0.78 ± 0.02
1.58 ± 0.01
0.76 ± 0.09

39

R

R²

7.54 ± 1.76
6.52 ± 3.97
3.69 ± 0.60
8.51 ± 1.49
2.76 ± 0.93

0.99
0.97
0.92
0.96
0.70

Table 6. Monod kinetic parameters and the yield coefficient for the degradation of metalaxyl
Mix
2b
2c
3a
3b
3c
5a
5b
5c

µ* (x 10-7)
9.15 ± 14.02
13.43 ± 2.92
9.99 ± 0.89
10.95 ± 3.96
4.62 ± 0.95
3.38 ± 0.65
4.43 ± 0.70
1.99 ± 0.25

µm* (x 10-4)-0.18 ± 8.21
-2.11 ± 1.93
-5.27 ± 1.09
-0.98 ± 3.23
-2.42 ± 1.41
-3.95 ± 0.99
-2.69 ± 0.54
-2.51 ± 0.93

kdecay3.15·10-6 ± 7.01·10-3
1.08·10-3 ± 1.37·10-3
2.57·10-3 ± 6.22·10-4
2.30·10-5 ± 2.01·10-3
6.26·10-4 ± 1.06·10-3
5.27·10-4 ± 4.88·10-4
4.90·10-8 ± 2.14·10-4
1.12·10-7 ± 4.30·10-4

R²
0.96
0.96
0.96
0.86
0.98
0.93
0.91
0.90

Y
20.20
156.89527.3989.66
524.841168.18
607.711266.12

3.5 Linuron
Linuron is a pesticide with a higher Koc value than isoproturon, metalaxyl and bentazon, and it isthus quit immobile. No breakthrough of this pesticide was observed for 4 mixtures (out of 5)even after about 20 pore volumes. A slight breakthrough occurred only in mixture 4.
4 CONCLUSION
Pesticide displacement experiments were carried out for isoproturon, bentazon, metalaxyl andlinuron in columns filled with five different types of organic mixtures. The experimental bromide, isoproturon and bentazon BTCs were well described using the transport model based on
the convection-dispersion equation. The transport model had to be enhanced with the Monod kinetics to describe experimental BTCs of metalaxyl.
REFERENCES
Carter, A. D. How pesticides get into water - and proposed reduction measures. Pesticide Outlook 11, 149157. 2000.
De Wilde, T., Spanoghe, P., Debaer, C., Ryckeboer, J., Springael, D. & Jaeken, P. 2007. Overview of on-farmbioremediation systems to reduce the occurrence of point source contamination. Pest Management Science
63: 111-128.
De Wilde, T., Spanoghe, P., Ryckeboer, J., Jaeken, P. & Springael, D. 2008. Sorption characteristics of pesticides on organic substrates used in biopurification systems. Environmental Science and Pollution Research, In press.
Decoin, M. 2003. Où en est la Fontaine-du-Theil. Phytoma 557: 29-32.
Guimont, S., Perrin-Ganier, C., Real, B. & Schiavon, M. 2005. Effects of soil moisture and treatment volume
on bentazon mobility in soil. Agronomy for Sustainable Development 25: 323-329.
Lapidus, L. & Amundson, N. R. 1952. Mathematics of adsorption in beds: VI. The effect of longitudinal diffusion in ion exchange and chromatographic columns. Journal of Physical Chemistry 56: 984-988.
Ramwell, C. T., Johnson, P. D., Boxall, A. B. A. & Rimmer, D. A. 2004. Pesticide residues on the externalsurfaces of field-crop sprayers: environmental impact. Pest Management Science 60: 795-802.
Šimůnek, J., van Genuchten, M. Th. & Šejna, M. 2005. The HYDRUS-1D software package for simulating
the one-dimensional movement of water, heat, and multiple solutes in variably-saturated media. Version
3.0, HYDRUS Software Series 1, Department of Environmental Sciences, University of California Riverside, Riverside, CA, 270 pp.
Torstensson, L. & Castillo, M. D. 1997. Use of biobeds in Sweden to minimize environmental spillages fromagricultural spraying equipment. Pesticide Outlook June: 24-27.

40

Describing pesticide degradation using
information from batch and column dataMonodkinetics:

Cheyns1, K., J. Mertens1, J. Šimůnek2, E. Smolders1 and D. Springael1
1 Division soil- and water management, KULeuven, Arenbergpark 20, 3001 Heverlee, Belgium

2 Department of Environmental Sciences, University of California Riverside, Riverside, CA, USA

1 INTRODUCTION
Degradation and sorption characteristics of pesticides are crucial to quantify their transport
through the unsaturated zone towards the groundwater. Parameters needed for the mathematicaldescription of both processes are usually estimated from laboratory batch or column
experiments and thereafter incorporated into larger scale models. This study evaluates the use ofbatch mineralization data, as well as pesticide leaching data through laboratory columns, for theestimation of parameters needed for the description of both pesticide sorption and degradation.
Pesticide degradation in field or catchment scale models is commonly described using afirst-order approach (Roulier et al., 2006). However, this first-order approach is in contrast withwhat is often seen in laboratory mineralization experiments where a lag phase (i.e., a period
during which no degradation occurs) is commonly observed (Barriuso & Houot, 1996, MunierLamy et al., 2002, Mertens et al., 2008). Since in both batch mineralization and column
experiments of this study a lag phase is clearly noticeable, a first-order approximation of theMonod kinetics (Monod, 1949) is used to describe the kinetics of degradation. With respect tothe sorption kinetics of the pesticide to the soil matrix, linear equilibrium sorption is assumed.
This study uses the Shuffled Complex Evolution Metropolis (SCEM) algorithm (Vrugt etal., 2003) to identify the model parameters and their associated uncertainty. The objectives ofthe study are to compare the distributions of the model parameters that were optimized against(a) only batch mineralization data or (b) only column leachate data, and then (c) to investigate
the suitability of the optimized batch parameters to predict the observed column leachate data.
2 MATERIALS AND METHODS
Steady state pesticide transport through soil columns
The soil was collected between depths of 30 and 45 cm inside of an agricultural field (Beverst,
Belgium), which had been yearly treated with atrazine since 1973. Sixteen columns (length 6.5
cm, internal diameter 3 cm) were each packed with 80 g of disturbed soil (sieved at 2 mm) at an
average bulk density of 1.5 g cm-³. At the bottom of the column, a glass filter plate (Robuglassfilter) was installed and a suction of -100 cm H2O was applied. Columns were irrigated at a
rate of 1.7 mm day-1 using a peristaltic pump.
Before applying the atrazine, the columns were pre-leached with a salt solution, whichmimicked the soil solution, of 8 pore volumes. In the next step, atrazine dissolved in the salt
solution (final concentration of 12 mg l-1) was added as a step application for 30 days (20 porevolumes). Leachates were collected in glass bottles (250 ml), which were periodically emptied
to determine leachate volume and atrazine concentration. Atrazine was concentrated from theleachates using Oasis HLB Solid Phase Extraction (Waters corporation) cartridges andsubsequent methanol extraction from the cartridges. The methanol extract was analyzed withreverse-phase HPLC (Merck, Hitachi, acetonitrile/H2O: 60/40).

41

Batch 14C-Atrazine mineralization by the indigenous soil community
To examine the atrazine mineralization capacity of the indigenous soil community, two columns
were leached with two pore volumes (36 ml) of a 14C-atrazine containing solution (20 mg l-1
atrazine and 8 µg l-1 [U-ring-14C]-atrazine (Izotop, 1391 MBq mmol-1) dissolved in the salt
solution). The total amount of atrazine in the soil and in the pore water in the column was: 0.26
mg atrazine including 0.17 µg 14C-atrazine. These columns were incubated in a dark airtightchamber with an alkali trap (2 ml of 1 M NaOH) to trap the 14CO2 produced from 14C-atrazine
mineralization. Periodically, the NaOH solution was removed and replaced by fresh alkali. The
removed NaOH solution was mixed with 10 ml of liquid scintillation cocktail (Ultima Gold,
PerkinElmer) and analyzed in a liquid scintillation counter (Wallac 1410, Perkin Elmer).
Cumulative mineralization curves were established for both soil samples. The finalmineralization extent was evaluated from the total amount of 14C-atrazine mineralized after 45
days.

Model and inverse model development
To describe the atrazine mineralization, a first-order approximation of the Monod kinetics wasapplied, which simultaneously describes the microbial growth and pesticide degradation as
follows:dX= (μ m ,mod Cl )X − k decay X
dt

(1)

dCl
1
= − (μ m ,mod Cl )X
dt
Y

(2)

where X is the biomass concentration [mg cm-3], Cl is the liquid atrazine concentration [mg cm], μm,mod is the modified mass growth rate [cm3 mg-1 day-1], kdecay is the decay rate [day-1] and Y
is the yield [-]. This simplification of the full Monod kinetics was implemented in MATLAB for
the simulation of the batch mineralization data and in the HYDRUS-1D model (Simunek et al.
2005) for the description of pesticide leaching and degradation through the soil columns. The
inverse modeling SCEM algorithm was coupled to both the batch and the HYDRUS-1D model.
From a sensitivity analysis (data not shown), it was concluded that the sensitivity of the batchand column model to the value of Y is low. Therefore the value of Y is not incorporated in the
inverse modelling but fixed to 0.4. This value corresponds to 1 minus the final mineralized
percentage atrazine after 45 days (60 %, see Figure 1), assuming that all atrazine is either
respired to CO2 or incorporated in the biomass. The biomass decay rate, kdecay is assumed to be 0
for two reasons: (i) neither batch nor column experimental data contain any information about
the decay of the biomass and (ii) a limited time span (maximum of 40 days) is considered inboth experiments.

3

3 RESULTS AND DISCUSSION
Model fits and optimized model parameter distributions
Figure 1 (left) shows the observed and simulated mineralization rates expressed as percentage
against time. It is clear that the model is able to simulate the observed lag time (around 15 days)
and the subsequent mineralization well. Figure 1 (right) presents the observed and simulated
atrazine leaching through one (randomly chosen) of the sixteen columns. Again, the HYDRUS1D model succeeds in describing the observed atrazine leachate concentrations well. Both thebatch and column model fit was obtained using the SCEM algorithm as described above. It isinteresting to note the presence of a lag time in the column leachate data. Optimized cumulative
distributions of both batch and column model parameters are plotted in Figure 2. First of all,Figure 2 shows narrow optimized parameter distributions indicating small parameter
42

uncertainties for most model parameters (except for the xini parameter in the column model).
Although both the column and batch data are fitted well, the corresponding optimized
parameters do not coincide. Growth rates in the column are shown to be significantly higher (>5times) than for the batch data. Contrary, the initial biomass concentrations are optimized to besignificantly lower (> 1000 times) for the column case as compared to the batch case. Since bothparameters are strongly correlated, it was investigated whether the optimal batch parameter set
was able to simulate the observed column data.
0.015

70

Atrazine concentration (mg cm-3)% Atrazine broken down

60
50
40
30
20

0.01

0.005

10
0

0

10

20

30

40

0

50

0

5

10

15

20

25

30

Time (days)

Time (days)

Figure 1. Observed (x) and simulated (-) mineralization data in a batch experiment (left) and observed (x)and simulated (-) atrazine concentrations in the column leachates (right).

1
0.5

Cumulative distribution

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

3 -1

Kf (cm g )
1
0.5
0

1

1.5

2

2.5
3

3
-1

3.5

4

-9

-8

-1

log(qmmod) (cm mg day )
1
0.5
0
-13

-12

-11

-10
-3

log(xini) (mg cm )Figure 3. Optimized parameter distributions. Full lines represent the optimized parameter distributions forthe batch experiment. Dashed lines represent the optimized parameter distributions from the column
experiment.

43

Predicting column data using parameters from the batch mineralization data.
Figure 3 shows the atrazine breakthrough curve predicted using the growth rate parameter andthe initial biomass concentration optimized against the batch experiment. Despite the higher
initial biomass, the lower growth rate causes a longer lag time than the lag time observed in thecolumn experiment. Similarly, atrazine breakdown is predicted to be slower as compared to theobserved column data. This suggests that the pesticide degradation kinetics in batch and column
experiments are different. This discrepancy could possibly be explained by: (i) the difference
between the mineralization process (batch) and the degradation process (column): i.e., the firststep in the degradation path of atrazine is not the rate limiting step in the degradation, (ii) the
kinetics of the CO2 release from the soil is slower than the mineralization kinetics, (iii) liquidflow in the soil columns creates other conditions for the bacterial community than those createdby the single atrazine application in the batch experiment (e.g. biofilm formation). Currently, a
new batch experiment is being set-up, in which both mineralization (CO2 formation) anddegradation (disappearance of atrazine) are simultaneously measured to investigate the possible
difference in kinetics between these two processes.

-3

Atrazine concentration (mg cm )

0.015

0.01

0.005

0

0

50

100

150

Time (days)

Figure 4. Observed (x) and simulated (-) atrazine breakthrough in the column leachates. The full line
shows the fit to the column data, the dashed line shows the simulated breakthrough using the optimized
parameters from the batch mineralization experiment.

4 CONCLUSIONS
It is clear that the model used in this experiment was able to describe the mineralization ofatrazine in the batch experiment well. Similarly, Hydrus-1D simulated the observed
breakthrough curve of atrazine through the soil columns well. Both approaches use the Monod
kinetics to describe the degradation of atrazine. However, it is shown that predicting atrazine
breakthrough through columns using optimized batch parameters for growth rate and initialbiomass concentration is not possible! The study demonstrates that care should be taken when
upscaling pesticide degradation kinetics measured in batch experiments to larger scales (i.e., tothe column or field scale).

44

5 REFERENCES
Barriuso, E. & Houot, S., 1996. Rapid mineralization of the s-triazin ring of atrazine in soils in relation tosoil management. Soil Biology and Biochemistry. 28 (11), 1341-1348.
Monod, J. 1949. The growth of bacterial cultures. Annu. Rev. Microbiol., 3: 371-394.
Mertens, J., Sniegowski, K., Diels, J., Smolders, E., & Springael, D. 2008. Inverse modeling of pesticide
mineralization and biomass dynamics in a bioremediation system: parameterizing the Monod model,
submitted to Chemosphere
Munier-Lamy, C., Feuvrier, M.P. & Choné, T. 2002. Degradation of 14C-atrazine bound residues in brownsoil and rendzina fractions. Journal of Environmental Quality. 31, 241-247.
Roulier, S.; Baran, N.; Mouvet, C.; Stenemo, F.; Morvan, X.; Albrechtsen, H.-J.; Clausen, L., & Jarvis, N.
2006 Controls on atrazine leaching through a soil-unsaturated fractured limestone sequence at Brevilles,
France. Journal of Contaminant Hydrology. 84 (1-2), 81-105.
Šimůnek, J., van Genuchten, M. Th. & Šejna, M. 2005. The HYDRUS-1D software package for simulating
the one-dimensional movement of water, heat, and multiple solutes in variably-saturated media. Version
3.0, HYDRUS Software Series 1, Department of Environmental Sciences, University of California
Riverside, Riverside, CA, 270 pp.
Vrugt, J.A, Gupta, H.V., Bouten, W. & S. Sorooshian. 2003. A shuffled complex evolution Metropolis
algorithm for optimization and uncertainty assessment of hydrologic model parameters, Water Resour.
Res., 39 (8), 1201.

45

46

A HYDRUS based Approach for Coupled Modeling of Vadose
Zone and Ground Water Flow at Different Scales
Navin K. C. Twarakavi & J. Šimůnek
Department of Environmental Sciences, University of California Riverside, Riverside, CA 92521, USA

S. SeoColorado School of Mines, 1500 Illinois St., Golden, CO 80401

ABSTRACT: In the past, ground water flow modeling has been done with little or no consideration of the effects of vadose zone processes. This dramatic simplification and/or neglect of
the effects of vadose zone processes have been utilized owing to constraints on computational
resources. Another reason for such a simplification is to work around the relatively high nonlinearity of vadose zone flow compared to ground water. While such simplified representationof vadose zone flow processes leads to faster modeling simulations, it also lead to a poor characterization of groundwater fluxes in the area under consideration. Ideally, one would prefer tostrike a balance between the modeling efficiency and the computational demand of the adopted
approach. In this paper, a one-dimensional unsaturated flow package HYDRUS, recently developed for the three-dimensional modular finite-difference ground water model MODFLOW, isevaluated and compared to other modeling approaches accounting for fluxes through the vadosezone. The HYDRUS package for MODFLOW consists of two sub-models that interact in space
and time- (a) HYDRUS sub-model (vadose zone) and (b) MODFLOW sub-model (ground water). The HYDRUS package is based on the HYDRUS-1D program that solves the onedimensional Richards equation simulating water movement in the vadose zone. The HYDRUS
package considers the main processes and factors affecting fluxes in the vadose zone, such as
precipitation, infiltration, evaporation, redistribution, capillary rise, plant water uptake, water
accumulation at the ground surface, surface runoff, and soil moisture storage. Being fully incorporated into the MODFLOW program, the HYDRUS package provides MODFLOW with recharge fluxes into groundwater, while MODFLOW provides HYDRUS with the position of thegroundwater table that is used as the bottom boundary condition in the package. The performance of the HYDRUS package is analyzed for various case studies that involve different spatial
and temporal scales.
1 INTRODUCTION
Ground water is one of the nation's most important natural resources. Ground water provides
drinking water for more than one-half of the nation's population (Solley et al., 1998) and is thesole source of drinking water for many rural communities and some large cities. The transfer
rates between the land surface and the groundwater table is controlled by the water flow throughthe variably saturated zone. Depending upon hydrological, geological and soil characteristics,
rain and snowmelt is partitioned at the land surface into runoff, infiltration, evapotranspiration(ET), groundwater recharge, and vadose zone storage (Twarakavi et al., 2008). A proper scientific evaluation of any hydrological problem related to ground water, therefore, needs to take into account the processes affecting flow through the vadose zone. Even though the importance ofthe vadose zone has been often highlighted in ground water research, its representation has been
often simplified during modeling (Keese et al., 2005). One of the major reasons for this simplification is that modeling of vadose zone flow processes is a complex and computationally de47

manding process and that it is often handicapped by the lack of data necessary for characterizing
the hydraulic properties of the subsurface environment. There is an urgent need to circumvent
this frequent simplification and effectively simulate flow through the vadose zone in large scalehydrological models (Winter et al., 1998)
Traditional ground water models incorporate vadose zone effects on ground water using water budget or ‘residual’ approach (Scanlon, 2002). Water budget methods are based on the water
budget equation (Twarakavi et al., 2008). The water budget equation is a mathematical representation of the fact that all water arriving at the water table either leaves the system as groundwater
flow, is discharged through sinks such as surface water bodies, is evapotranspirated, or retained
as storage (Scanlon, 2002). One may indirectly estimate the water table recharge from the water
budget equation by measuring or estimating all other components in the water budget. MODFLOW (Harbaugh et al., 2000), which is a widely used ground water modeling tool, provides
this option through the REC and ET packages. MODFLOW was developed by U.S. Geological
Survey and is one of the most widely used groundwater flow models. It is obvious that the RECET package tends to oversimplify and underestimate the effects of vadose zone flow on
groundwater. In spite of its widespread use, the REC-ET approach suffers from the following
major limitations: (a) the accuracy of recharge estimates depends on the accuracy of measurements of all other water budget components (for the REC-ET package, these include precipitation and evapotranspiration), (b) the method is not reliable for deeper water tables, and (c) theapplicability of this approach is questionable for arid and semi-arid regions where soil capillarypressures play a dominant role for vadose zone flow (Gee and Hillel, 1988; Lerner et al., 1990;
Twarakavi et al., 2008).
Recently, more promising approaches have been suggested that involve coupling vadose zone
models to MODFLOW. A coupled model simulates the effects of near-surface hydrologicalprocesses on groundwater flow by linking a groundwater model with a selected vadose zone
model in space and time. The majority of currently available vadose zone models is based on either the Richards (Richards, 1931) or kinematic wave (Colbeck, 1972; Smith, 1983) equations.
The kinematic wave equation is a simplification of the Richards equation where the capillary effects on flow are neglected. Unlike the water balance approach, coupled models estimate the interaction between the ground water and vadose zone by solving the corresponding governing
equations simultaneously in space and time. Simultaneously solving the two governing equations is achieved through a numerical strategy such as finite difference and/or finite element approaches. Evaluation of interactions between the near-surface and groundwater flow processes
using coupled models has been a desirable but difficult goal. Even though the idea of a coupledapproach has been floating for several years, it has been realized only in the recent decade.
Twarakavi et al. (2008) outlined criteria that need to be addressed for a high usability of thecoupled approach: (a) accuracy in representation of physical processes that drive the vadosezone flow, (b) usability for different groundwater modeling scenarios, (c) applicability to different spatial and temporal scales, i.e., from lab or field to regional spatial scales and from hourly
to decadal temporal scales, and (d) applicability to different meteorological and climactic conditions, such as humid, arid, and semi-arid regions. Three MODFLOW packages accounting forprocesses in the vadose zone have been recently developed: (a) the Variably Saturated Flow
(VSF) package (Thoms et al., 2006), (b) the Unsaturated Zone Flow (UZF1) package (Niswonger et al., 2006), and (c) the HYDRUS package (Seo et al., 2007). Twarakavi et al. (2008) also
provide a detailed comparison of the three packages and readers are referred to the aforementioned references for further details. It was noted that the HYDRUS package provides an ideal
option for solving coupled vadose zone and ground water flow models at different scales. The
objective of this paper is to present the HYDRUS package and briefly demonstrate its applicability using a few examples.
2 MODELING APPROACH
2.1 Governing Equations
The HYDRUS package couples the one-dimensional Richards equation based vadose zone flow
model into the MODFLOW suite. The necessity to couple the one-dimensional vadose zone48

flow model to a three-dimensional groundwater flow model was driven by computational andaccuracy considerations (Twarakavi et al., 2008). The HYDRUS package for MODFLOW wasdeveloped in order to consider the effects of precipitation, infiltration, evaporation, plant water
uptake, soil moisture storage, and water accumulation at the ground surface and in the vadosezone. One must note that the developed package was based on a radical simplification of the
HYDRUS-1D program (Šimůnek et al., 2005) that can simulate the one-dimensional water
movement in the variably-saturated zone in all its complexity. The simplification involved removal of subroutines simulating solute and heat transport, hysteresis in the soil hydraulic functions, and boundary conditions that were irrelevant for this purpose.
As mentioned earlier, a coupled approach involves simultaneously solving two governing
equations, one for vadose zone and one for ground water, in space and time using a numericalstrategy. The ground water flow model in MODFLOW is represented using a finite difference
approximation of the mass-conservation equation:∂h∂ ⎛ ∂h ⎞ ∂ ⎛ ∂h ⎞ ∂ ⎛ ∂h ⎞⎜ K x ⎟ + ⎜⎜ K y ⎟⎟ + ⎜ K z ⎟ − W = S s∂t∂x ⎝ ∂x ⎠ ∂y ⎝ ∂y ⎠ ∂z ⎝ ∂z ⎠(1)

where Kx, Ky, and Kz are hydraulic conductivities (LT-1) in the direction of x, y, and z coordinates, respectively; h represents the pressure head (L), W is the volumetric flux per unit volume
through sources or sinks (T-1), Ss is the specific storage of the porous material (L-1), and t is time(T).
The vadose zone flow in the HYDRUS package is represented using the one-dimensional
Richards Equation:∂θ ∂ ⎡∂h⎤=
K ( h) − K ( h ) ⎥ − S∂t ∂z ⎢⎣∂z⎦(2)

where θ is the volumetric water content (-), h is the soil water pressure head (L), t is time (T), z
is the distance from reference datum (L), S is the sink term usually accounting for root water uptake (T-1), and K(h) is the unsaturated hydraulic conductivity (LT-1) as a function of h or θ. One
may note that the Richards equation is highly nonlinear due to the dependence of the unsaturated hydraulic conductivity, K(h), and the water content, θ(h), on the capillary pressure head, h.
The Brooks and Corey (Brooks and Corey, 1964) and van Genuchten-Mualem (van Genuchten,1980) models available in the HYDRUS package are among the most widely used approaches torepresent these nonlinear relationships. The Galerkin-type linear finite element scheme is usedin HYDRUS to numerically solve the Richards equation. In the HYDRUS package for MODFLOW, the one-dimensional Richards equation (2), represented by a finite-element approximation, is coupled in space and time to the ground water flow equation (2), as represented by thefinite-difference approximation.
2.2 Spatial and time discretization
The HYDRUS package for MODFLOW comprises of two sub-models interacting with eachother in space and time: (a) HYDRUS sub-model and (b) MODFLOW sub-model. Based on
similarities in soil hydrology and topographical characteristics, the discretized MODFLOW domain can be divided into zones, which comprise one or more grids of the MODFLOW model.
One HYDRUS soil profile is then assigned to each of these zones. It is assumed that the HYDRUS soil profile adequately represents vadose zone flow for the entire zone. The HYDRUS
vertical soil profiles are then discretized vertically into finite elements which are constructed by
splitting the soil profile into one-dimensional elements that are connected to each other at nodal
points. Once the finite elements are constructed, they may not be changed during the simulation.
Care should therefore be taken to ensure that the depth of the soil profile extends from the soil
surface to the deepest possible water table that may be expected during the simulation. To ensure convergence of the numerical solution, finite-element dimensions should be relatively small
at locations where sharp pressure head gradients are expected. Such smaller elements are usuallyneeded close to the soil surface where meteorological factors can cause rapid changes in water49

MODFLOW
Sub-model

HYDRUS
Sub-model

Average watertable depths

Bottom fluxes as recharge at
the water table for the next
MODFLOW time step

Solve for bottom fluxes
in each profile using theatmospheric data and 1D
Richards Equation

Water

Zone 1

Zone 2

content and pressure
head gradients, and atinterfaces between different soil horizons.
One also needs to consider the soil texture
during the discretization
process. For example,
coarse textured soils
generally require finer
discretization than finetextured soils due to thehigher nonlinearity of
their soil hydraulicproperties. Once thespatial discretization of
the soil profile is carried out, the distributionof different soil materials in the profile needs
to be described (Seo etal., 2007).

Figure 1. Schematic description of the coupling procedure for water flow in
HYDRUS package for MODFLOW.

The computational efficiency of the coupled HYDRUS-MODFLOW system is enhanced bysimulating vadose zone and groundwater flows at their own, often different, time steps. This isneeded because a proper treatment of the Richards equation requires smaller time steps thanthose usually used in MODFLOW simulations. The two sub-models (HYDRUS and MODFLOW) interact, i.e., exchange information about the groundwater recharge and the groundwater level, only at the end of each MODFLOW time step during which HYDRUS may perform
multiple of its own time steps. MODFLOW receives the recharge flux from HYDRUS and calculates a new water table depth for a particular time step (Figure 1). A new water table depth isthen assigned as the pressure head bottom boundary condition in the HYDRUS package for the
next MODFLOW time step. The iteration procedure in the HYDRUS package is similar to that
described in the HYDRUS-1D manual (Simunek et al., 2005). The readers are referred to theabove mentioned references and Seo et al. (2007) for more details.
The HYDRUS package for MODFLOW is analyzed for two case studies: (a) one-dimensional
Las Cruces infiltration experiment of Wierenga et al. (1991), (b) water table recharge experiment of Vauclin et al. (1979).
3 CASE STUDIES
3.1 One-dimensional infiltration
The first case study is the one-dimensional infiltration experiment at the Las Cruces trench site
(Wierenga et al., 1991). The experiment was a field study, conducted in southern New Mexico,
the primary purpose of which was to develop a dataset for validating and testing numerical models. For this purpose, the study site was heavily instrumented with neutron probes, tensiometers,
and solute samplers for measuring water contents, pressure heads and solute concentrations(Wierenga et al., 1991). More than 500 soil samples (undisturbed and disturbed) were taken at
the experimental site and analyzed in the laboratory for bulk density, and to find the saturated
hydraulic conductivity and the soil water retention curve. The experimental infiltration study involved application of water to a 4 m wide area using closely spaced drips with an average sur50

face flux of 1.82 cm/day for the 86 days of the experiment. In order to reduce the disruption of
the experimental conditions by rain and evaporation, the irrigated area and its surroundings werecovered by a pond liner.
A uniform soil profile with an equivalent saturated hydraulic conductivity (Ks=270.1 cm/d)
was considered for modeling purposes. The RETC code was used to analyze the retention curve
data for undisturbed and disturbed soils (>500 soil samples), resulting in the following retention
curve parameter values (van Genuchten, 1980): θs= 0.321, θr= 0.083, α= 0.0550 cm-1 and
n=1.509. Initial pressure heads (hi = -100 cm) in the soil profile were the same as those used by
Wierenga et al. (1991). While a constant water flux was used as the upper boundary condition(q0 = 1.82 cm/d), free drainage was considered at the lower boundary. The initial and boundary
conditions in terms of the water content, θ(z, t), are described in Equation (5)θ ( z, 0) = θinit ( z ) = θ (hinit ( z ))(5a)q( z = 0, t )=1.82 cm/d(5b)

θ (b, t ) = θinit (b)(5c)

where b is the depth of the soil profile, which must be large enough so that the wetting front
does not affect the water content at the bottom of the soil profile during the simulation. A soil
profile depth of 600 cm was used and the simulation was run for 35 days.
Experimental results of the Las Cruces trench infiltration experiment (Wierenga et al., 1991)
are compared with results simulated using the HYDRUS package in Figure 2. It may be notedthat the HYDRUS is expected to perform well as it solves the Richards equation for onedimensional problems. Figure 2 shows the soil water content profiles for different days of theexperiment and compares the model predictions of the HYDRUS package with the experimental
data. The HYDRUS package adequately predicted the one-dimensional infiltration experimentobservations.

0

0.05

Water content
0.1 0.15 0.2

0.00

0.25

(a)

-100

Depth (cm)D ep th (cm )

-300
-400

-600
-700

(b)

-200

-200

-500

0.30

0

0
-100

Water content
0.10
0.20

HYDRUS
Observed

-300
-400
-500

HYDRUS
Observed

-600
-700

Figure 2. A comparison of measured water contents and corresponding estimatescalculated using HYDRUS package for (a) day 19 and (b) day 35 for the Las Cruces
trench experiment (data from Wierenga et al., 1991).

3.2 Two-dimensional water table recharge experiment
The HYDRUS package for MODFLOW was used to model the two-dimensional transient water-table experiment of Vauclin et al. (1979). The experimental setup consisted of a 6.0-m by 20 m box containing sandy soil. The initial water table elevation was 0.65 m from the bottom. A
51

constant flux of q = 3.55 m/d was applied across the center 1.0 m of the soil surface, while therest of the surface was covered to prevent evaporation. Due to the symmetry of the experiment,
only one half of the experiment was modeled and the model domain was thus 3.0 m by 2.0 m.
The initial total head was set equal to 0.65 and the right boundary cells were constrained to theinitial water table position throughout the 8-hour simulation. The grid was discretized into uniform cells of 0.1-m width and 0.05-m depth. Only two soil profiles representing the soil directly
below the recharge zone and the rest of the transport domain were used in calculations with the
HYDRUS package. A saturated hydraulic conductivity of 840 cm/d was used. The initial andboundary conditions are described in Equation (6)Domain0 ≤ x ≤ 3m, 0 ≤ z ≤ 2 m 0 ≤ t ≤ 8 hr
Δx = 0.1m, Δz = 0.05 m , Δt = 1min
Initial condition
h( x, z , 0) = 0.65 − z [m]Boundary conditionq ( 0 ≤ x ≤ 0.5, 2, t ) = 3.55 m/dq ( x > 0.5, 0, t )= 0.00 m/dq (0, z, t )= 0.00 m/dh(3, z , t )= 0.65 − z [m](6)where q(x,z,t) is the flux at spatial coordinates x and z at time t.
Figure 3 shows the modeling results from the HYDRUS package and the observations fromthe experiment. Water tables calculated using the HYDRUS package is very comparable to theobservations even though the numerical solution of the Richards equation in the HYDRUS
package is limited only to the vertical direction. It was observed that the one-dimensional natureof the vadose zone modeling used in the HYDRUS package did not significantly affect the correspondence of simulated results with experimental data.
Recharge= 3.55 m/d

2

t=8 hours
t=4 hours
t=3 hours
t=2 hours

1.8
1.6
Elevation (m)

1.4
1.2
1
0.8
0.6

Initial water table

0.4

Constrained at 0.65 m

0.2
0
0

0.5

1

1.5
Distance (m)

2

2.5

3

Figure 3. A comparison of the simulation of water table positions calculatedusing the HYDRUS packages with the experimental data of Vauclin et al. (1979).

52

4 SUMMARY AND CONCLUSIONSThe HYDRUS package for MODFLOW provides a viable alternative for incorporating vadosezone flow processes into ground water flow models. Based on the HYDRUS-1D software, the
HYDRUS package considers the effects of infiltration, soil moisture storage, evaporation, plantwater uptake, precipitation, runoff, and water accumulation at the ground surface. The HYDRUS package also offers a variety of models for describing soil hydraulic properties. Necessary parameter values required for these models in HYDRUS may be estimated experimentallyor using pedo-transfer functions.
5 REFERENCES
Brooks, R. H. and A. T. Corey. 1964. Hydraulic properties of porous media. Colorado State University,
Fort Collins.
Colbeck, S. C. 1972. A theory of water percolation in snow. Journal of Glaciology, 2(63):369-385.
Gee G.W., and D. Hillel. 1988. Groundwater recharge in arid regions: review and critique of estimation
methods. Hydrological Processes, 2:255-266.
Harbaugh, A. W., E. R. Banta, M. C. Hill, and M. G. McDonald. 2000. MODFLOW-2000, the U.S. Geological Survey modular ground-water model user guide to modularization concepts and the groundwater flow process. Denver, CO, Reston, VA: U.S. Geological Survey.
Keese, K. E., B. R. Scanlon, and R. C. Reedy. 2005. Assessing controls on diffuse groundwater recharge
using
unsaturated
flow
modeling.
Water
Resources
Research,
41:W06010,
doi:10.1029/2004WR003841.
Lerner, D. N., A. S. Issar and I. Simmers. 1990. Groundwater recharge, a guide to understanding and estimating natural recharge. International Association of Hydrogeologists, Kenilworth, Rep 8, 345 pp.
Niswonger, R. G., D. E. Prudic, and R. S. Regan. 2006. Documentation of the Unsaturated-Zone Flow(UZF1) package for modeling unsaturated flow between the land surface and the water table with
MODFLOW-2005[Online].
Available
byU.S.
Geological
Surveyhttp://pubs.usgs.gov/tm/2006/tm6a19/.
Richards, L. A. 1931. Capillary conduction of liquids through porous media. Physics, 1: 318–333.
Scanlon, B. R. 2002. Choosing appropriate techniques for quantifying groundwater recharge. Hydrogeology Journal, 10:18-39.
Seo, H. S., J. Šimůnek, and E. P. Poeter. 2007. Documentation of the HYDRUS Package for MODFLOW-2000, the U.S. Geological Survey Modular Ground-Water Model, GWMI 2007-01, Int.
Ground Water Modeling Center, Colorado School of Mines, Golden, CO, 96 p.
Šimůnek, J., M. Th. van Genuchten, and M. Šejna. 2005. The HYDRUS-1D software package for simulating
the one-dimensional movement of water, heat, and multiple solutes in variably-saturated media. Version
3.0, HYDRUS Software Series 1, Department of Environmental Sciences, University of California Riverside, Riverside, CA, 270 pp.
Smith, R. E. 1983. Approximate sediment water movement by kinematic characteristics. Soil Science Society of America Journal, 47: 3-8.
Solley, W. B. 1998. Estimates of water use in the western United States in 1990, and water-use trends,
1960–90. U.S. Geological Survey Open-File Report 97–176, 15 p.
Thoms, R. B., R. L. Johnson, and R. W. Healy. 2006. User's guide to the Variably Saturated Flow (VSF)
processfor
MODFLOW[Online].
Available
byU.S.
Geological
Surveyhttp://pubs.usgs.gov/tm/2006/tm6a18/.
Twarakavi, N. K. C., J. Šimůnek, and S. Seo. 2008. Evaluating Interactions between groundwater and vadose zone using the HYDRUS-based flow package for MODFLOW. Vadose Zone Journal, in press.
van Genuchten, M. Th. 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal, 44(5):892-898.
Winter, T. C., J. W. Harvey, O. L. Franke, and W. M. Alley. 1998. Ground water and surface water a single resource. U.S. Geological Survey Circular, 1139, 77 p.

53

54

Modeling of Coupled Water and Heat Fluxes in both Unfrozen
and Frozen Soils
Y. Zhao, S. Peth & R. Horn
Institute of Plant Nutrition and Soil Science, Christian-Albrechts-University zu Kiel, Olshausenstr. 40,
24118 Kiel, Germany

ABSTRACT: Accurate simulation of soil freezing and thawing behavior is critical to understandhydraulic processes in the vadose zone under cold and arid climatic conditions. Using an
extended freezing code incorporated in the HYDRUS-1D model, this study was conducted 1) toverify the freezing model using field soil water and temperature data collected in Inner
Mongolia grassland, and 2) to investigate the contribution of snowmelt or soil thawing to theseasonal water balance. The results showed that both the freezing model and the snow routine
matched well the measured soil water and temperature under unfrozen conditions. However,under frozen conditions, the freezing model reflected the phase change of soil water andsubstantially improved the simulation results than the snow routine. The freezing model did
obviously not produce surface runoff generated by snowmelt or soil thawing from frozen soil
layers. Instead, it overestimated water content and thus underestimated surface runoff after
spring snowmelt. We suggest that a detailed knowledge of the soil-atmosphere processes isneeded to improve the surface runoff algorithm in the frozen soil module.
INTRODUCTION
Coupled water and heat movement in the vadose zone is a central process in many agriculturaland engineering issues. Especially, in cold and arid regions, snowmelt or lateral water
movement on frozen soil layers have a non-negligible influence on seasonal water balance.
However, although the former processes are recognized widely, the mutual interactions of waterand heat flow in frozen soil are limited in laboratory observation and theoretical analysis, andrarely considered in field applications (Flerchinger & Saxton, 1989; Smirnova et al. 2000;
Hansson et al. 2004). This study addresses the field application of the hydrodynamic model
HYDRUS-1D (Šimůnek et al. 1998). In the current newest program, an extended freezing code
is incorporated, which numerically solves coupled equations governing phase change between
water and ice and heat transport with a mass- and energy-conservative method (Hansson et al.
2004). Specifically, we will focus on following questions: 1) How well does HYDRUS-1D
simulate soil water and temperature with and without “frozen soil module”? 2) How does thefrozen soil module affect soil temperature, soil moisture and runoff simulations?
MATERIALS AND METHODS
The study was performed on a long-term experimental site of the Inner Mongolia Grassland
Ecosystem Research Station (IMGERS, 43o37′50′′N, 116o42′18′′E). Detailed descriptions of theexperimental area were given in Zhao et al. (2008). The site was protected from grazing since
1979 (24 ha). Since June 2004, soil moisture and temperature were recorded by a data-logger at
30-min intervals in summer and at 1-h intervals in winter, respectively. Soil water content wasmeasured at 5, 20 and 40 cm depth by theta-probes (Type ML2x). In case the soil is frozen, it
refers to the volumetric unfrozen water content. Soil temperature was measured in five depths at
2, 8, 20, 40 and 100 cm using Pt-100 probes. Precipitation and other weather variables wererecorded by a micrometeorological station (50 m away from the monitoring plot). To determine
root length density, root samples were taken up to 100 cm depth. In addition, undisturbed soil
samples were taken at four depths of 4-8, 18-22, 30-34 and 40-44 cm to determine the water
retention characteristics and hydraulic conductivities.
The model of coupled water and heat fluxes was performed with HYDRUS-1D. Variably
55

saturated water flow for above- and sub-zero temperatures is described using the modified
Richards equation (e.g. Hansson et al. 2004):∂θu ( h ) pi ∂θi (T ) ∂ ⎡+= ⎢ K Lh ∂h + K Lh + K LT ∂T + K vh ∂h + K vT ∂T ⎤⎥ − S
∂t
pw ∂t
∂z ⎣
∂z
∂z
∂z
∂z ⎦

(1)

where θu [L3 L-3] is the volumetric unfrozen water content (=θ+θv; θ and θv are the volumetricliquid and vapor water content, respectively), θi is the volumetric ice content, pw and pi is thedensity of liquid and ice water, respectively, t is time, z is the soil depth, h is the pressure head,
T is the temperature, and S is a sink/source term. In Eq. 1, the first five terms on the right-handside represent liquid flows due to gradient in pressure head (KLh, [L T-1]), gravity, andtemperature (KLT, [L2 T-1 K-1]), and vapor flows due to gradient in pressure head (Kvh) andtemperature (KvT), respectively.
The hydraulic conductivity of frozen soil is significantly reduced by ice lenses, which isaccounted for by an impedance factor, Ω (Lundin 1990), multiplied by Q, as follows:
KfLh =10-ΩQ KLh(2)the parameter Q is the ratio of the ice content to the total water content, which accounts for the
more significant blocking effects with the increase in ice content.
The governing equation for the movement of energy in a variably saturated rigid porous
medium is given by the following conduction–convection heat flow equation (e.g. Nassar &Horton 1992):∂C p T∂θ∂ θ (T )− L f ρ i i + L 0 (T ) v
=
∂t
∂t
∂t
∂q T
∂q T
∂q
= ∂ ⎡⎢ λ (θ ) ∂T ⎤⎥ − C w l − C v v + L0 (T ) v − C w ST
∂z ⎣
∂z ⎦
∂z
∂z
∂z

(3)

where L0 and Lf are the volumetric latent heat of water vaporization and freezing, respectively.
Cp is the volumetric heat capacity of the bulk soil, which is determined as the sum of thevolumetric heat capacities including solid, organic, liquid (Cw), ice, and vapor (Cv) phasemultiplied by their respective volumetric fractions (De Vries 1963). The symbol λ(θ) denotes theapparent thermal conductivity and q the water flux density while T is the soil temperature.
The phase change between water and ice is controlled by the generalized Clapeyron equation,
which defines a relationship between the liquid pressure head and temperature when ice ispresent in the porous material. Hence, the unfrozen water content can be derived from the liquidpressure head as a function of temperature when ice and pure water co-exist in the soil.
The initial condition was set based on measured water content and temperature. An
atmospheric boundary condition and free drainage condition was imposed at the soil surface andbottom boundary of the flow domain, respectively. The soil profile was considered to be 100 cm
deep. Root water uptake was simulated using the model of Feddes et al. (1978). The
HYDRUS-1D model was parameterized by calibrating hydraulic parameters (laboratory-derived
hydraulic properties), estimating of the intercepted water and ET partitioning based on themeasurements (Zhao et al. 2008). The current HYDRUS-1D version including soil frozen
module (denoted as freezing model) was modified to consider the subsurface soil freezing andthawing processes, as well as the surface energy and water balances. To verify the performanceof freezing model, the “normal” version including snow hydrology only without soil frozen
module (snow routine) was also run as a reference.
RESULTS AND DISCUSSION
Soil water and heat fluxes are numerically simulated for the whole hydrological year in 2006
(Fig. 1). In general, the simulated and measured soil water contents (SWC) are comparable
during the studied period in terms of root mean square error (0.02-0.07 cm3 cm-3). Particularly,
the freezing model simulates well the diurnal water dynamics which coincides with soil freezing
and thawing processes, i.e. soil moisture increases with increasing soil temperature and vice
56

-10
-20
S

F

7-Oct
7-Oct

16-Nov
16-Nov

9-Feb
9-Feb

Rainfall 40

20

Soil temperature ( C)0
9-Feb 21-Mar 30-Apr 9-Jun 19-Jul 28-Aug 7-Oct 16-Nov 26-Dec
Time (days)

b

10
0
-10
-20
S

F

21-Apr

0.
10

0.
05

0.
00

Soil temperature ( C)

0
-10
-20
14

S

16

18

20

22

0.
10

0.
05

2

4

6

8

10 12 14 16
Time (hours)

18

20

22

24

7-Oct
7-Oct

16-Nov
16-Nov

26-Dec26-DecRainfall40

20

0
9-Feb 21-Mar 30-Apr 9-Jun 19-Jul 28-Aug 7-Oct 16-Nov 26-Dec
Time (days)

40 cm

21-Mar21-Mar30-Apr30-Apr

9-Jun
9-Jun

M

19-Jul19-Jul28-Aug28-Aug

S

F

7-Oct
7-Oct

16-Nov
16-Nov

26-Dec26-DecRainfall40

20

0
9-Feb 21-Mar 30-Apr 9-Jun 19-Jul 28-Aug 7-Oct 16-Nov 26-Dec
Time (days)

30
20

100 cm

10
0
-10
-20
9-Feb
9-Feb

F

F

-20

24

Soil moisture(cm cm )M

1228-Aug28-Aug00. 0. 0. 0. 0. 0. 0. 0.
00 05 10 15 20 25 30 35

10

-3

8

3

-3

6

0.
15

Soil moisture(cm cm )

3

4

0.
20

2

19-Jul19-Jul

S

-10

o

o

Soil temperature ( C)

c

9-Jun
9-Jun

M

9-Feb
9-Feb

20-Mar 25-Mar 31-Mar 5-Apr 11-Apr 16-Apr 21-Apr 27-Apr
Time (hours)

30-Apr30-Apr10

27-Apr

10

0.
00

20

Soil moisture(cm cm )M

16-Apr21-Mar21-Mar30

0. 0. 0. 0. 0. 0. 0. 0.
00 05 10 15 20 25 30 35

11-Apr

-3

5-Apr

3

-3

31-Mar0.
15

Soil moisture(cm cm )325-Mar0.
20

20-Mar20-20

26-Dec26-Dec-3

28-Aug28-Aug3

19-Jul19-Jul0-10

Soil moisture(cm cm )

9-Jun
9-Jun

Rainfall (mm)30-Apr30-Apr10o

3oSoil temperature ( C)21-Mar21-MarM

0. 0. 0. 0. 0. 0. 0. 0.
00 05 10 15 20 25 30 35-3

Soil moisture(cm cm )

9-Feb
9-Feb

20 cm

Rainfall (mm)020

Rainfall (mm)

10

20

30

o

a

0. 0. 0. 0. 0. 0. 0. 0.
00 05 10 15 20 25 30 35

20

21-Mar21-Mar30-Apr30-Apr

9-Jun
9-Jun

19-Jul19-JulS

28-Aug28-Aug

F

7-Oct
7-Oct

16-Nov
16-Nov

26-Dec26-DecRainfall40

20

Rainfall (mm)

Soil temperature ( C)30oSoil temperature ( C)versa (Fig. 1b). However, the snow routine only fits soil water contents under unfrozen
condition. Under frozen condition, there is a clear disparity between the liquid water content line
simulated by the freezing model and the total water content line simulated by the snow routine(Fig. 1a). This discrepancy is apparently caused by the function of two models, and may be used
to approximate the ice content in the soil. For instance, the SWC simulated by the freezing
model drops shortly after 6th November associated with soil freezing (Fig. 1a). However, the
SWC simulated by snow routine keeps constant. From which, an ice content of 0.07 cm3 cm-3
can be estimated. In late March, coinciding with soil temperature increasing above 0°C, both
simulated (freezing model) and measured SWC raise, while simulated SWC by snow routine
keep constant. This again proofs that the freezing model can predict well the increase in SWC
by soil thawing. A increase in SWC in the deeper soil layers (20 and 40 cm) accompanied bysoil thawing is also clearly predicted (Fig. 1). However, there is no indication of vertical water
movement since soil water content is low and it can be held by soil.

09-Feb 21-Mar 30-Apr 9-Jun 19-Jul 28-Aug 7-Oct 16-Nov 26-Dec
Time (days)

Figure 1. Measured and simulated soil moisture and temperature at 5 cm (a, b, c), 20, 40, and 100 cm
depth during the whole year of 2006 (M: Measured liquid water content; S: Simulated total water content
running snow routine; and F: Simulated liquid water content running freezing model).

In contrast to the soil water simulations, soil temperature is simulated well by both models
with minor differences (Fig. 1). However, regarding the detailed diurnal variations of soil
temperature (Fig. 1c), the freezing model matches measured values better than the snow routine.
57

16
14
4

a

2
0
30
20
10
0
-10
-20
-30

9-Feb

Runoff (mm)30-Apr9-Jun

19-Jul28-Aug7-Octb16-Nov

26-DecRainfall

Air
Soil

60
40
20
0

9-Feb

0.24
0.20
0.16
0.12
0.08
0.04
0.00

21-MarRainfall (mm)oTemperature ( C)Snow depth (mm)This gives an evidence of the impact of the function of the frozen soil module on soil
temperature simulations. In a real world, when soil becomes freezing (i.e. the phase change fromwater to ice), soil temperature decreases and energy is released which is used to warm up soil
from a cold. However, given the same total water content, the thermal conductivity of frozen
soil is higher than that of unfrozen soil due to present of ice. Consequently, the upward soil heat
flux becomes higher when the soil gets frozen, and thus it tends to cool the soil (Smirnova et al.
2000). But the effect of thermal conductivity is probably smaller than that of water phase change.
Consequently, the freezing model, which considers the both processes, provides a morereasonably and realistic simulations of soil temperature than the snow routine.

c21-Mar30-Apr9-Jun

Snow routine

19-Jul28-Aug7-Oct16-Nov

26-DecFreezing model

9-Feb 21-Mar30-Apr 9-Jun 19-Jul28-Aug 7-Oct 16-Nov26-Dec

Time (days)
Figure 2. Rainfall, measured air and soil temperature, simulated snow depth and runoff during the wholeyear of 2006.

The temperature threshold at which precipitation is considered to fall as snow is specified as-2°C in the model, suggesting that precipitation in winter is falling as snow. The snow routine
predicts up to 15 mm snow depth (Fig. 2a), however, the simulated runoff after air temperature
increasing above 0°C is negligible (Fig. 2c). This might relate with that the snow routine does
not account for surface runoff from the frozen soil layer. In fact, it is likely that surface runoff isgenerated during snowmelt while soil is fully or at least partially frozen (Fig. 2b). Unexpectedly,
the freezing model, which can account for the subsurface freezing and thawing processes, also
does not produce surface runoff during winter (Fig. 2c). Instead, we found that the simulated
SWC by freezing model is higher than the measured values in the transition time when soil
begins to thaw (Fig. 1a). This implies that the freezing model might overestimate water content
and underestimate surface runoff after spring snowmelt. Therefore, the freezing model seems
still not sensitive enough to estimate surface runoff accompanied with snowmelt from the soil
58

frozen layer. This might relate to the fact that the freezing model we applied adopts soil surface
temperature as the atmospheric boundary condition instead of air temperature, which neglects
the lag-effects of energy transfer (Fig. 2b). Consequently, the freezing model may incorrectly
partition all the snowmelt into infiltration as both soil thawing and snow melting happen
simultaneously. Therefore, to solve this, a transferable and double-layered boundary condition(e.g. one accounting for air temperature and other accounting for considering soil temperature)should be considered. Furthermore, the reduction in infiltration capacity owing to the blocking
effects of ice is possibly underestimated by the current freezing model. Additionally, the gradualrelease of water from the frozen soil profile also might reduce the maximum rate of runoff.
Frozen soil was expected to have a great influence on the runoff simulation (Pitman et al. 1999).
Although the current frozen soil module has little effect on the simulations of surface runoff, wesuggest a detailed study on the soil-atmosphere processes and effects of boundary conditions toimprove the surface runoff algorithm in the freezing code.
CONCLUSION
We used an extended frozen soil module of HYDRUS-1D to govern coupled flow equation
which solves water and heat transport under both frozen and unfrozen conditions simultaneously.
The model was evaluated using field data of soil water and temperature at a long-term
experimental site in Inner Mongolia grassland (North China). The results showed that bothfreezing model and snow routine reflect well the measured soil water and temperature underunfrozen condition, whereas the freezing model substantially improved the simulation results
under frozen condition. In addition, the freezing model did obviously not produce surface runoff
generated by snowmelt or soil thawing from frozen soil layer. We suggest that seasonal waterbalance, especially considering rainfall water stored as snow, snow drift and the lateral waterflow on frozen soil layers need to be investigated further because of the complicated interactions
at the soil-atmosphere interface and thus effects of boundary conditions on the simulation.
ACKNOWLEDGEMENTS
This work was done with the financial support of the German Research Council (DFG) under#536 MAGIM. We thank Prof. J. Šimůnek for his helps on the HYDRUS code.
REFERENCES
de Vries, D.A. 1963. The thermal properties of soils. p. 210–235. In R.W. van Wijk (ed.) Physics of plant
environment. North Holland, Amsterdam.
Feddes, R.A., Kowalik, P.J. & Zaradny, H. Simulation of Field Water Use and Crop Yield, John Wiley andSons, New York, NY, 1978.
Flerchinger, G.N., & Saxton, K.E. 1989. Simultaneous heat and water model of a freezingsnow-residue-soil system I. Theory and development. Trans. ASAE 32:565–571.
Hansson, K., Šimůnek, J., Mizoguchi, M., Lundin, L.-C., & van Genuchten, M.T. 2004. Water flow andheat transport in frozen soil: numerical solution and freeze-thaw applications. Vadose Zone J.
3:693–704.
Lundin, L.-C. 1990. Hydraulic properties in an operational model of frozen soil. J. Hydrol. 118:289–310.
Nassar, I.N., & Horton, R. 1992. Simultaneous transfer of heat, water, and solute in porous media: I.
Theoretical development. Soil Sci. Soc. Am. J. 56:1350–1356.
Pitman, A.J., Slater, A.G., Desborough, C.E. & Zhao, M. 1999: Uncertainty in the simulation of runoff
due to the parameterization of frozen soil moisture using the Global Soil Wetness Project methodology.
J. Geophys. Res. 104:16879–16888.
Smirnova, T.G., Brown, J.M. & Benjamin, S.G. 1997. Performance of different soil model configurations
in simulating ground surface temperature and surface fluxes. Mon. Wea. Rev. 125:1870–1884.
Šimůnek, J., Sejna, M. & van Genuchten, M.Th. 1998. The HYDRUS-1D software package forsimulating the one dimensional movement of water, heat, and multiple solutes in variably-saturated
media. Version 2.0. IGWMC-TPS-70. Int. GroundWater Modeling Center, Colorado School of Mines,
Golden.
Zhao, Y., Peth, S., Krümmelbein, J., Ketzer, B., Gao, Y.Z., Peng, X.H., Bernhofer, C. & Horn, R. 2008.
Modelling grazing effects on coupled water and heat fluxes in Inner Mongolia grassland. (Submitted to
Vadose Zone Journal).
59

60

HYDRUS-1D Modeling of an Irrigated Agricultural Plot with
Application to Aquifer Recharge Estimation
J. Jiménez-Martínez
Department of Geotechnical Engineering and Geosciences, Technical University of Catalonia, UPC,
Barcelona, Spain

T.H. Skaggs
U.S. Salinity Laboratory, USDA-ARS, Riverside, CA, USA

M.Th. van Genuchten
U.S. Salinity Laboratory, USDA-ARS, Riverside, CA, USA

L. Candela
Department of Geotechnical Engineering and Geosciences, Technical University of Catalonia, UPC,
Barcelona, Spain

ABSTRACT: A variety of methods are available for estimating aquifer recharge in semi-arid
regions, each with advantages and disadvantages. We are investigating a procedure for estimating recharge in an irrigated basin. The method involves computing irrigation return flows basedon HYDRUS-1D modeling of root zone and surface boundary water dynamics. In this paper, wedemonstrate the modeling approach and the calibration of the model to experimental data.
1 INTRODUCTION
Estimating aquifer recharge is important for determining water resource availability and assessing aquifer vulnerability to pollutants (Scanlon et al. 2002). Recharge estimation can be difficult, particularly in arid and semi-arid regions where water tables are typically deep and recharge is predominately focused recharge that emanates from topological depressions such as
streams and lakes. The recharge rate is limited by the availability of water at the land surface,which is controlled by temporally and spatially variable climatic factors such as precipitationand evapotranspiration (Scanlon et al. 2002). In some basins recharge is additionally complicated by irrigation, which may simultaneously remove water from focused recharge sources
while creating new sources of diffuse recharge. In irrigated regions, accurate knowledge of recharge, evaporation, and transpiration is especially important for the sustainable management ofscarce water resources (e.g., Garatuza-Payan et al. 1998).
Several methodologies have been used to estimate groundwater recharge with varying degrees of success (Scanlon et al. 2002). The methodologies can be loosely grouped into three
categories depending on whether the focus of the method is surface-water, the unsaturated zone,or the saturated zone. In each of the three cases, physical and tracers techniques are possible, as
are numerical modeling approaches. The best choice for a particular situation depends of thespatial and temporal scales being considered and the intended application of the recharge estimate (Scanlon et al. 2002).
In this paper, we report progress on a study of recharge in the Campo de Cartagena in southeast Spain, a semi-arid region where irrigated agriculture is prevalent. The goal of the study is toinvestigate various methods for estimating the contribution of irrigation to aquifer recharge, including a modeling approach that is based on HYDRUS-1D (Šimůnek et al. 2005) simulationsof root-zone water dynamics. Here, we present preliminary results demonstrating the modelingapproach and the calibration of the model to experimental data; a future publication will present
results related to recharge estimation.

61

2 MATERIALS AND METHODS
2.1 Study location: Campo de Cartagena, Spain
The Campo de Cartagena is a region of 1440 km2 located in southeast Spain. The climate is
Mediterranean with an average annual rainfall of 300 mm and a mean annual temperature of
18 ºC. Estimates of annual potential evapotranspiration (ETP) range from 800 to 1200 mm y-1,
depending of the estimation method (Sánchez et al. 1989). The primary land use is agriculture,with 128.1 km2 of row crops (principally lettuce and melon), 34.1 km2 of perennial vegetables(principally artichoke), and 136.8 km2 of fruit trees (principally citrus). Drip irrigation is usedwidely in the region due to the scarcity of water resources and the need for water conservation.
2.2 Field site and experiment
A study of root zone soil moisture was conducted on an experimental plot at the Tomas Ferro
Agricultural Science Center, a research facility operated by the Technical University of Cartagena. The 3 × 8 m plot was managed according to agricultural practices that are standard inCampo de Cartagena, including cultivating crops (melon and lettuce) that are common in thearea and utilizing drip irrigation. The drip system featured 16 mm inside diameter tubing, 4 L h-1
emitters, and an emitter spacing of 30 cm. In total, 36 emitters were installed. The plot was instrumented to monitor water dynamics in the root zone. Instrumentation consisted of: two tensiometers (Soilmoisture Equipment Corp, Goleta, CA, USA) installed vertically at each of thedepths 30, 45, 60, 90 and 120 cm (10 tensiometers total); and two 44 mm diameter, 2 m deep
TRIME-FM TDR access tubes for measuring soil moisture at various depths (Imko,Germany).
Meteorological data for the site were available from the Servicio de Información Agraria deMurcia (http://siam.imida.es). Cropping and irrigation for 2007 are summarized in Figure 1.
2.3 HYDRUS-1D modeling
We simulated root zone water dynamics for the experimental site using HYDRUS-1D (Šimůnek
et al. 2005). The simulations included irrigation, evaporation, and root water uptake processes.
Initial estimates for the soil hydraulic parameters θr, θs, α, n, and Ks for different soil layers were
estimated with the ROSETTA pedotransfer function using measured data for bulk density andthe sand, silt, and clay contents. Refinements to these parameters estimates were made subsequently based on model fitting to a subset of the measured water content and pressure head data
(details given below).

Figure 1. Summary of cropping, irrigation, and atmospheric conditions for the experimental plot, 2007.

62

Implementing the atmospheric boundary condition required the daily irrigation and precipitation values shown in Fig. 1, as well as specifying daily potential evaporation and transpiration
values. To determine these latter values, we first obtained a daily reference evapotranspiration,ET0(t), using weather station data and the Penman-Monteith method. We then estimated the potential evapotranspiration ETp(t) for specific crops according to:ETp (t ) = Kc (t )·ET 0(t )(1)where Kc(t) is a coefficient that incorporates specific crop water uptake characteristics and theaveraged effects of evaporation from the soil. Allen et al. (1998) describe a method of specifying the time variation of Kc(t) in terms of annual crop growth stage (initial, crop development,
mid-season, and late season stage), and also provide data on the length of growth stages andmagnitude of Kc for various crops. We used the Allen et al. (1998) method and data for single
annual crops to specify daily values for Kc(t). Then, with ETp given by Eq. (1), we calculated thepotential evaporation Ep(t) according to (e.g. Kroes and Van Dam, 2003)

Ep (t ) = ETp (t )·e−β· LAI (t )(2)where β (≈ 0.4) is the radiation extinction coefficient and LAI(t) is the leaf area index. The potential transpiration Tp(t) was then given by:

Tp (t ) = ETp (t ) − Ep (t )(3)

The daily values for Tp and Ep are shown in Fig. 1. Root water uptake reduction due to drought
stress was computed using the Feddes et al. (1978) model, with model parameter values for specific crops obtained from the database included in HYDRUS-1D.
3 RESULTS
The most intensive data collection occurred during the cultivation of melon, 17 May to 10 September, 2007 (Fig. 1, day of year 137-253). We used the root zone water content and pressure
head data from this period and the inverse parameter optimization routines in HYDRUS-1D totest and adjust the calibration of the soil hydraulic property model. Several possible parameterizations were considered which differed according to the number of soil layers (1 to 4) and thenumber of hydraulic parameters that were fitted for each layer (different combinations of 1 to 4
parameters). The initial estimates for the parameters when fitted, or their fixed value when not
fitted, were the Rosetta estimates discussed above. The best overall parameterization was determined informally based on diagnostic information provided by the HYDRUS-1D routines aboutmodel fit and algorithm convergence, visual inspection of the model fit to the data, and the principle of parsimony. It was determined that the best parameterization involved four soil layerswith two parameters, α and n, fitted for each layer. Fitting more than two parameters per layer
tended to cause the inverse algorithm to fail to converge, possibly due to the lack of a unique solution. Four soil layers produced a better fit to the data than was possible with fewer layers in
the profile. Overall, the chosen solution with four layers and α and n fitted for each layer provided the best correlation between measured and simulated water content and pressure head values (R2 = 0.90).
The final soil hydraulic parameters, as well as soil textural data, are given in Table 1, whereas
the measured and simulated (water contents and pressure heads at the depths 30, 45, 60, 90, and120 cm are shown in Figure 2.
Table1. Soil textural and hydraulic properties.
Textural fractions (%)
Bulk
Depth
density
Sand
Silt
Clay
g·cm-3
0-30
18.7
76.0
3.5
1.45
30-60
13.8
80.2
6.0
1.52
60-90
19.5
77.2
3.3
1.58
90-150
10.8
82.0
6.6
1.70

θs

θr
-3

cm cm
0.04
0.05
0.04
0.04
63

-3

cm cm
0.38
0.38
0.35
0.36

α (Fit)n (Fit)

Ks

-1

1.16
1.23
1.27
1.46

cm d-1
46.4
29.9
30.1
13.2

cm
0.078
0.046
0.014
0.020

Figure 2. Measured and simulated (fitted) root zone pressure heads and water contents for melon crop.

64

Using the final soil parameterization, HYDRUS-1D was subsequently used to predict rootzone soil moisture dynamics for the growth period of the other (lettuce) crops cultivated on theexperimental plot, the first grown between 7 February and 16 May 2007 (DOY 38-136) and thesecond between 25 September and 24 December 2007 (DOY 268-358). Though not shown here,
our preliminary findings indicate that the HYDRUS predictions are in reasonably good agreement with field measurements. The correlation between measured and predicted water contents
and pressure heads was 0.82 for both lettuce crops.
4 CONCLUSIONS
Overall, good agreement was obtained between the calibrated HYDRUS simulation model andthe root zone measurements made in a cultivated field plot. With this calibrated model, it will bepossible to compute for various cropping and irrigation scenarios the fraction of irrigation waterthat passes below the root zone and contributes to aquifer recharge. Combined with data on land
use and cropping patterns for the Campo de Cartagena, it will then be possible to estimate theamount of recharge associated with irrigated agriculture. A study by the Spanish Geological
Survey (IGME, 1994), using a different methodology, estimated that total recharge to the top
unconfined aquifer was 69 hm3 y-1, with 46 hm3 y-1 due to natural recharge and 23 hm3 y-1 due toirrigation return. Our preliminary calculations indicate a slightly larger value for recharge due toirrigation; this aspect of our work will be the subject of a future publication.
REFERENCES
Allen, R.G., Pereira, L.S., Raes, D., & Smith, M. 1998. Crop evapotranspiration. Guidelines for computing crop water requirements. Irrigation and Drainage. Paper No. 56, FAO, Rome, Italy.
Feddes, R. A., Kowalik, P.J., & Zaradny, H. 1978. Simulation of Field Water Use and Crop Yield. John
Wiley and Sons, NY.
Gartuza-Payán, J., Shuttleworth, W.J., Encinas, D., McNeil, D.D., Stewart, J.B., DeBruin, H. & Watts, C.
1998. Measurenent and Modelling evaporation for irrigated crops in Northwest Mexico. Hydrol. Proc.,
12: 1397-1418.
IGME. 1994. Las aguas subterráneas del Campo de Cartagena (Murcia). Spanish Geological Survey,
Spain 62 pp.
Kroes, J.G., & J.C. Van Dam. 2003. Reference manual SWAP: Version 3.0.3. Rep. 773. Alterra Green
World Res., Wageningen, the Netherlands.
Sánchez, M.I., López, F., Del Amor, F. & Torrecillas, A. 1989. La evaporación y evapotranspiración enel Campo de Cartagena y Vega Media del Segura. Primeros resultados. Anales de Edafología y Agrobiología 1239–1251.
Scanlon, B.R., Healy, R.W., & Cook, P.G. 2002. Choosing appropriate techniques for quantifying
groundwater recharge. Hydrogeology Journal 10: 18–39.
Šimůnek, J., van Genuchten, M. Th. & Šejna, M. 2005. The HYDRUS-1D software package for simulating the one-dimensional movement of water, heat, and multiple solutes in variably-saturated media.
Version 3.0, HYDRUS Software Series 1, Department of Environmental Sciences, University of California Riverside, Riverside, CA, 270 pp.

65

66

Downward Migration of Chernobyl-derived Radionuclides in Soils in
Poland and Sweden
Gerald Matisoff1, Lauren F. Vitko1, Peter J. Whiting1, Michael E. Ketterer2, Jerzy W.
Mietelski3, Edyta Lokas3, Klas Rosen4, and Henning Persson4
1

Case Western Reserve University, Cleveland OH 44106-7216 USA, 2Northern Arizona University,
Flagstaff AZ 86011-5698 USA, 3The Henryk Niewodniczanski Institute of Nuclear Physics, Krakow
Poland, 4Swedish University of Agricultural Sciences Uppsala Sweden

INTRODUCTION
The Chernobyl accident on April 26, 1986 resulted in significant fallout of radionuclides
such as 137Cs and 239,240Pu on surface soils throughout northern Europe. Knowledge ofthe dynamics and mechanisms of the migration of these radionuclides in soils is
important for determining animal and human dose exposure rates, determining exposure
from food-chain transfer, and in planning environmental remediation and clean-up. In
addition, both natural (7Be, 210Pbxs) and anthropogenic (137Cs, 239,240Pu) fallout
radionuclides have been extensively employed to determine short-term soil erosion
rates, to trace sediment source regions, to characterize and quantify erosion
mechanisms, to constrain sediment budgets, and to better understand the delivery ratios,transit distance and transit time of fine sediment and adsorbed pollutants. Despite theuse of these radionuclides and their activity-depth profiles, little is known about how the
profile shapes develop or why they differ with location. There is a substantial amount of
information that is embedded in the distribution of radionuclides with depth that could
be extracted to refine our ability to understand significant radiochemical behavior, topredict dose exposure rates, to better plan environmental remediation, to use these andother radionuclides to understand soil erosion mechanisms, to identify sediment sourceareas and to calculate watershed inventories and residence times for better
understanding of watershed retention and erosion processes. Moreover, ignoring thisevidence may lead to a real risk of drawing incorrect conclusions from an incomplete
understanding of the radionuclide profiles.
The primary goal of the research is to characterize how radionuclides migrate throughsoils and attach to fine sediment and are eroded from the landscape and transported
through stream networks. At the same time the findings have implications for questionsof human and ecosystem health, landscape evolution, geochemical cycling and water
resource degradation. The downcore distributions of natural and anthropogenic
radionuclides in soils are controlled by a variety of processes, including depositional
flux, precipitation, infiltration, porosity/permeability, water content, dispersion
coefficient, organic matter content, clay content, bulk density, particle size, adsorption,
cation exchange capacity, colloidal transport, bypass flows through cracks, bioturbation
and radioactive decay. The objective of this research is to develop a better quantitative
understanding of the processes that lead to the development of radionuclide profiles insoil columns. To help constrain the results, the work was conducted at a series of field
locations that cover a range in values for several key parameters. Presented here are soil
profiles of 137Cs, 210Pbxs and 239,240Pu from the Lazy, Tulowice and Bor za Lasem sites in
Poland and the Hille and Skogsvallen sites in Sweden so that the profiles of these
nuclides may be compared with those obtained previously (Chemicki et al., 1996-1997;
Boro et al., 2001; Mietelski et al., 2004; Rosen et al., 1999) to better understand the
downward migration.
67

MATERIALS AND METHODS
Study area
The Lazy, Tulowice and Bor za Lasem sites in Poland and the Hille and Skogsvallensites in Sweden were selected for sampling because of previous work at the locations.
The study sites exhibit a range of properties that affect radionuclide retardation or
migration and at which there are prior radionuclide profiles to provide a time series ofdata that can be used to constrain the modeling. The study sites span a range of
geography, elevation, precipitation, Chernobyl fallout flux and soil orders. In particular,
Bor za Lasem is a peat and is a relatively unique situation compared to mineral soils – it
might be expected to have very unusual characteristics with poor Cs binding but
possibly good Pu binding and immobilization.
Soil Sampling
Vertical sections of the soil profiles were collected using both tube corers and our highresolution corer (Wilson et al., 2003). Samples were collected from all sites in 2007,
although six vertical profiles from each of the two sites in Sweden were collected
previously over the 1987-2006 time period. We took samples from undisturbed,
relatively flat sites that appeared to be neither the site of erosion nor deposition. In theupper part of the soil, the sampling interval was 0.4-3.5 cm while deeper in the cores the
sampling interval ranged from 0.8-5 cm depending upon depth and type of corer. Most
soil cores were collected to a depth of 25-30 cm although a few cores were collected toa depth of 60 cm. Soil cores were returned to the laboratory for processing and analysis.
Soil samples were dried at 80-105º C, ground, placed in containers and sealed in
preparation for radionuclide analysis.
Gamma Spectroscopy and ICPMSRadionuclide activities were determined by gamma spectroscopy. Gamma spectroscopy
for 210Pb (46.52 keV), 7Be (477.6 keV), 137Cs (661.65 keV), and 214Bi (609.3 keV) is
reasonably well-established and straightforward, but the analytical details are detector
and sample specific. Among the three labs in the US, Poland, and Sweden, we usedeight HPGe gamma detectors for the radionuclide analyses. Activities measured
between the different detectors varied by less than 10%. Samples were counted for 2-24
hours to decrease the counting errors (~10%) associated with small samples and low
radionuclide activities. For the 210Pb analyses, the soil samples were counted after three
weeks to allow secular equilibrium ingrowth of gaseous 222Rn (t1/2 = 3.82 d) from thedecay of its 226Ra (t1/2 = 1600 yr) parent. The self-attenuation correction of Cutshall etal. (1983) was used. Standards were prepared using the same geometries as the
samples.
The instrumentation and counting geometries were calibrated using
commercially available mixed energies gamma-emitting standards. All measured counts
were corrected for background levels, detector and geometry efficiencies, branchingratios and decay.
Plutonium activities (239+240Pu) and atom ratios were determined in HNO3/HF digested
soils following the methods described in Cizdziel et al. (2008). 242Pu (0.007 Bq) wasused as a spike isotope; Pu fractions were processed using TEVA resin, and Pu isotopes
were measured by sector field ICPMS.
Computer Models
Preliminary calculations were performed using the model CHAIN (van Genuchen,1985) and HYDRUS-1D (Simunek, Van Genuchten, and Sejna, 2005) to illustrate thenature of solute migration in a soil profile and how model simulations can provide
insight to the key processes governing transport.

68

RESULTS AND CONCLUSIONS
The 137Cs fallout by Chernobyl is about 30 times greater than that from bomb fallout at
the three sites in Poland and about 300 times greater than at the two sites in Sweden.
Vertical profiles of 7Be, 210Pb, and 137Cs activities in soils indicate that 7Be activity
decreases rapidly downcore due to its short (53 d) half-life. Excess 210Pb fromatmospheric fallout decreases downcore more slowly due to its longer (20.4 yr) half-life.
Some 210Pb activity persists at depth because of in situ production from 226Ra decay.
137
Cs activity displays a subsurface peak reflecting Chernobyl fallout and subsequent
downcore migration. A time series of the 137Cs profiles indicates that the peak in 137Cs
activity occurred at the soil surface in 1987, but migrated downcore about 4 cm since,
and that some 137Cs activity can now be found as deep as 20 cm. The amount of
downward migration is different at the different locations. This is likely due to thedifferent types of soils, where some soils are dominated by mineral matter and others byorganic matter. At Hille at a depth of 32 cm there is a second, small peak that might be
the 1963 fallout peak. Work currently in progress is using Pu atom ratios to distinguishbetween Chernobyl and stratospheric fallout. The 240Pu/239Pu atom ratio is about 0.38
from Chernobyl and 0.18 from stratospheric fallout. There are also even larger relative
differences in 241Pu/239Pu and 242Pu/239Pu which can be measured if the activities arehigh enough. It is also expected that differences exist in micromineralogy and phaseassociation between stratospheric fallout and Chernobyl Pu; the latter is primarily
associated with refractory U-rich fuel particles.
CHAIN model simulations capture some of the major features of the 7Be, 137Cs, and210
Pbxs profiles. The simulations generate profiles that decrease downcore for both 7Be
and 210Pbxs while also accounting for the differences in half-life of these two
radionuclides. The 7Be migrates only about 1 cm before its activity approaches 0
whereas 210Pb decreases throughout the 25 cm that were simulated. The continuous
input of 210Pbxs results in a profile that remains fixed in time whereas a pulsed input
for137Cs results in a downward migration of the 137Cs peak with time. Dispersion resultsin tails on both the leading edge and the trailing side of the 137Cs pulse and theimportance of sorption and retardation is easily seen in both the 137Cs and 210Pb profiles.
However, the CHAIN model fails to effectively describe three aspects of the data: 1)
the short term transport as seen in the 7Be data is deeper than the model simulations
would suggest; 2) some of the 137Cs and 239,240Pu have been retained in the top portion of
the core whereas the model predicts that these radionuclides will be transported
downcore as a pulse; and 3) 210Pb activities decrease more rapidly near the surface thanthe model would predict.
CHAIN and other advection-dispersion models of solute transport fail to accurately
describe some aspects of the downward migration of radionuclides in the unsaturated
soil zone. HYDRUS-1D has the potential to better describe the processes by
incorporation of nonlinear adsorption, inclusion of soil properties into multiple zones,
and accounting for time-dependent processes such as water content. Future work will
involve additional modeling using the program HYDRUS-1D which should permit
exploring the causes for these differences to refine our ability to describe the transportof radionuclides in soils, for example, by permitting the use of nonlinear adsorption or
by defining soil properties into multiple depth zones.
ACKNOWLEDGEMENT
This work was supported by the Office of International Programs at NSF, grant number
OISE-0745066.

69

REFERENCES
K. Boron, J. W. Mietelski, K. Lipka, P. Gaca, M. Jasinska, 2001. Radionuclides in raised bogs on
example of „Bór za Lasem” in Orawsko-Nowotarska Valley in Tatra Mountains Foothill, Poland.
Journal of Environmental Monitoring, 3(3:) 324-329.
Chelmicki, W., M.Klimek, K.Krzemien, M.Jasinska, K.Kozak, J.W.Mietelski, 1996-1997. Spatial
distribution of Cs-137 on Pogórze Wielickie near Bochnia. Folia Geographica, Ser. Geogr. Phys.,
XXVIII: 35-47 (in Polish).
Cizdziel, J.V., M.E. Ketterer, D. Farmer, S.H. Faller, and V.F. Hodge, 2008. 239Pu, 240Pu, 241Pu
fingerprinting of plutonium in western US soils using ICPMS: solution and laser ablation
measurements. Analytical and Bioanalytical Chemistry, 390:521-530.
Cutshall, N. H., I.L. Larsen and C. R. Olsen, 1983. Direct analysis of 210Pb in sediment samples: selfadsorption correction. Nuclear Instruments and Methods, 206:309-312.
Mietelski, J.W., P.Szwalko, E.Tomankiewicz, P.Gaca, S. Malek, J.Barszcz, S.Grabowska, 2004. 137Cs,
40
K, 90Sr, 238, 239+240Pu, 241Am and 243+244Cm in forest litter and their transfer to some species of insects
and plants in boreal forests - three cases study. Journal of Radioanalytical and Nuclear Chemistry,62(3): 645-660.
Rosén, K., I. Öborn and H. Lönsjö, 1999. Migration of radiocesium in Swedish soil profiles after the
Chernobyl accident, 1987-1995. J. Environ. Radioactivity, 46:45-66.
Simunek, J., M. Th. Van Genuchten, and M. Sejna, 2005. The Hydrus-1D software package forsimulating the one-dimensional movement of water, heat, and multiple solutes in variable saturated
media (version 3.0).
Van Genuchten, M. Th., 1985. Convective-dispersive transport of solutes involved in sequential firstorder decay reactions. Computers & Geosciences, 11: 129-147.
Wilson, C.G., G. Matisoff, and P.J. Whiting, 2003. Short-term erosion rates from a 7Be inventory
balance. Earth Surface Processes and Landforms, 28(9): 967-977.

70

Overview of applications of the HYDRUS constructed wetland
moduleG. LangergraberInstitute of Sanitary Engineering and Water Pollution Control, University of Natural Resources andApplied Life Sciences (BOKU), Muthgasse 18, A-1190 Vienna, Austria.

ABSTRACT: The multi-component reactive transport module CW2D was developed to modeltransport and reactions of the main constituents of municipal wastewater in subsurface flow
constructed wetlands and is incorporated into the HYDRUS variably-saturated water flow andsolute transport program. CW2D is able to describe the biochemical elimination and the transformation processes for organic matter, nitrogen and phosphorus and considers 12 components(solutes) and 9 processes (reactions). This paper intends to give an overview of applications of
the CW2D module that differ in terms of the type of the constructed wetland investigated, thesize of the system as well as in the type of wastewater treated.
1 INSTRODUCTION
Constructed wetlands (CWs) provide a natural way for simple, inexpensive, and robust wastewater treatment. CWs can be subdivided into two main types, surface flow and subsurface flow
CWs. Surface flow, or free water surface CWs, are densely vegetated and typically have water
depths of less than 0.4 m. In subsurface flow (SSF) CWs, no free water level is visible. SSF
CWs are subdivided into horizontal flow (HF) and vertical flow (VF) systems depending on the
direction of water flow through the porous medium (sand or gravel).
A large number of physical, chemical, and biological processes contribute to the treatment.
Because these processes occur in parallel and influence each other, detailed understanding of
CW functioning is difficult. Because of their complexity, CWs have long been seen as "blackboxes" where wastewater enters and treated water leaves the system. Numerical models that describe the transformation and elimination processes in CWs are a promising tool to get a better
understanding of the processes in CWs. Only a few numerical models are available to describe
treatment processes in SSF CWs. Available models can be divided into three categories i) models describing the hydraulic behavior and single-solute transport only, ii) models describing reactive transport in saturated conditions such as in HF CWs, and iii) reactive transport models forvariably saturated conditions that can be used for modeling VF CWs as well as HF CWs, andhave the highest complexity (Langergraber, 2008). CW2D/HYDRUS belongs to the last category, both HF and VF systems can be modeled.
This paper presents the multi-component reactive transport module CW2D (Langergraber,
2001; Langergraber and Šimůnek, 2005), the constructed wetland module of HYDRUS, andshows examples for its application.
2 THE MULTI-COMPONENT REACTIVE TRANSPORT MODULE CW2D
The multi-component reactive transport module CW2D (Langergraber, 2001) was developed todescribe the biochemical transformation and degradation processes in SSF CWs. CW2D was incorporated into the HYDRUS variably-saturated water flow and solute transport program(Langergraber and Šimůnek, 2006; Šimůnek et al., 2006). The CW2D module considers71

12 components and 9 processes. The components include dissolved oxygen, organic matter(three fractions of different degradability, i.e. readily- and slowly-biodegradable, and inert),
ammonium, nitrite, nitrate, and nitrogen gas, inorganic phosphorus, and heterotrophic and two
species of autotrophic micro-organisms. Organic nitrogen and organic phosphorus are modeled
as nutrient contents of the organic matter, i.e. they are calculated as a percentage of COD. The
biochemical elimination and transformation processes are based on Monod-type expressions
used to describe the process rates. All process rates and diffusion coefficients are temperaturedependent. The processes considered are hydrolysis, mineralization of organic matter, nitrification (modeled as a two-step process), denitrification, and a lysis process (as the sum of all decay
and loss processes) for the micro-organisms. CW2D assumes a constant concentration of microorganisms (and other compounds) in each finite element. The thickness of the biofilm is notconsidered. The mathematical formulation of CW2D is based on the mathematical formulationof the Activated Sludge Models (ASMs; Henze et al., 2000). Langergraber (2005) investigatedthe plant uptake models provided by HYDRUS that describe nutrient uptake coupled with water
uptake, and concluded that it was possible to simulate plant uptake in high loaded systems, e.g.
systems treating mechanically pre-treated municipal wastewater. For low-strength wastewater,
the simulation results indicate that potential nutrient uptake is overestimated by using thesemodels. Oxygen release via roots can be modeled in a way similar to nutrient uptake (Toscano etal., 2006). For a detailed discussion of the CW2D module see Langergraber and Šimůnek
(2005).
3 APPLICATIONS
3.1 Overview
In the last years, CW2D has been applied to several types of SSF CWs and for different applications. Different types of CWs include VF beds (Langergraber and Šimůnek, 2005; Korkusuz etal., 2007; Langergraber, 2007), HF beds (Mena, 2008), hybrid systems, i.e. HF and VF beds inseries (Toscano et al., 2006, 2007), two-stage down- and upflow systems (Langergraber andŠimůnek, 2005). Different applications beside treating municipal wastewater include e.g. CWs
for polishing the effluent of a wastewater treatment plant for reuse purposes (Toscano et al.,2006, 2007) and treatment of combined sewer overflow (Dittmer et al., 2005; Meyer et al.,2006; Henrichs et al., 2007).The following chapters present examples of results presented in theabove mentioned literature.
3.2 Horizontal flow CWs
The first example presents simulation results of tracer experiments carried out for a pilot-scale
experimental installation in southern Spain (Mena, 2008). The installation consisted of threeparallel operated HF beds (each bed had a size of 2.5 m × 0.65 m with a depth of 0.6 m). Different species of macrophytes were planted. All beds were filled with gravel with a diameter of 69 mm (porosity 0.38), apart from the first and last 10 cm, for which the diameter was 9-12 mm(porosity 0.40) to improve the distribution of the wastewater in the inlet zone and to collect thetreated wastewater in the outlet zone, respectively. The CWs were continuously fed with artificial urban wastewater with mean influent concentrations of 127 mg COD.L-1 and 10 mg NH4N.L-1. Tracer experiments were carried out by adding one litre solution of sodium bromide witha concentration of 5000 mg.L-1 by means of an impulse in each CW. After the pulse, wastewater
feeding was kept contant and effluent samples were collected. Figure 1 shows the experimental
tracer concentration and the fit of different models. Mena (2008) showed that for all experiments
simulation with HYDRUS fitted better than detention time gamma distribution (DTGD) anddispersed plug-flow (DPF) models. Especially the DPF model did not fit very well to the experimental data mainly because the long skewed "tail" increases the mean HRT and this fact
makes the DPF to be delayed in relation to the experimental data.

72

Concentration Br- (mg/L)

12

10

Cexp
DPF
DTGD
HYDRUS-2D

8

6

4

2

0
0

5

10

15

20

25

30

Time (d)Figure 1: Comparison of experimental tracer concentrations to simulation results for CW1 (planted withPhragmites australis) using different models (according to Mena, 2008).

3.3 Vertical flow CWs for treatment of urban wastewater
Simulations of water flow in VF beds have been shown by Langergraber and Šimůnek (2005).
The experience showed that simulation results match the measured data when the hydraulic behavior of the system can be described well. The good match of experimental data to reactive
transport simulations can then be obtained using literature values for the CW2D model parameters, e.g. the values of the parameters given by Langergraber and Šimůnek (2005).
Korkusuz et al. (2007) investigated the hydraulic properties of 10 different natural and artificial substrates to be used for the main layer of VF CWs. Simulation results have been presented
for lab-scale studies carried out at the technical lab of the Institute of Sanitary Engineering at
BOKU. Each lab-scale column had a diameter of 20 cm and was filled with a 50 cm main layerusing different substrates that have been selected according their potential higher nitrogen andphosphorous adsorption capacity. The following parameters have been determined for all materials: particle and bulk density, porosity, saturated hydraulic conductivity and grain size distribution. Figure 2 shows as an example results for effluent flow rates from a single loading andfrom tracer experiments for Turkish zeolite. One main result of the study was the generation of
input data sets of the hydraulic properties for different substrates for CW2D.
0.7

1.6

0.6

Measured
Simulated

EC (mS/cm)

Q (L/h)

0.5
0.4
0.3
0.2

1.4

Measured

1.2

Simulated

1.0
0.8
0.6
0.4

0.1

0.2

0.0

0.0
0

1

2

3

4

5

6

0

Time (h)

20

40

60

80

100

120

Time (h)Figure 2: Simulated water flow (left) and tracer experiments (right) for Turkish zeolite (Korkusuz et al.,2007).

Langergraber (2007) presented simulation results for outdoor experiments that had the main
goal to verify the temperature model incorporated into CW2D. The experimental plant consistedof three parallel operating VF beds with a surface area of about 20 m² each, operated with intermittent loading, and has been in operation since 2003. The organic load has been
27 g COD.m-2.d-1, which correspond to hydraulic loading rate of 43 mm/d. Results for waterflow comparing measured and simulated effluent flow rates and the cumulative effluent of theexperimental CW are shown in Langergraber (2007). Simulated effluent flow rate and cumulative effluent flow match the measured data well. Using the calibrated flow model, the effluent
concentrations during summer could be simulated using the standard CW2D parameter set73

(Langergraber and Šimůnek, 2005). However, measured COD and NH4-N effluent concentrations at low temperatures could not be simulated because hydrolysis and nitrification at low
temperatures were over predicted. Figure 3 shows measured and simulated COD and NH4-N effluent concentrations using a modified parameter set that includes temperature dependencies forthe half-saturation constants of hydrolysis and nitrification. Using this modified parameter set, it
was possible to simulate the COD and NH4-N effluent concentrations at low temperatures(Langergraber, 2007).
70

COD, NH4-N (mg/l)

60

NH4-N simulated
COD simulated
NH4-N measured
COD measured

50

40

30

20

10

0
0

2

4

6

8

10

12

14

16

18

20

Temperature (°C)

Figure 3: Measured and simulated COD and NH4-N effluent concentrations using the modified parameterset (Langergraber, 2007).

To validate the biomass simulated growth Langergraber et al. (2007) used measured data of C
and N content of the biomass (Tietz et al., 2007) and converted them to microbial COD by using
factors based on stoichiometry. Figure 4 shows calculated and simulated microbial biomass
COD in different depths of the main layer for a heterotrophic lysis rate of 0.30 d-1. Simulated
microbial biomass COD in the first cm of the main layer was between 5600 and
3400 mg COD/g DW (the range of the measured values) when using heterotrophic lysis rates ofbetween 0.25 and 0.35 d-1.

Microbial biomass COD (mg/gDW)

7000
6000
Calculated from biomass C
5000
Simulated
4000
3000
2000
1000
0
0

5

10

15

20

25

30

Filter depth (cm)Figure 4: Calculated and simulated microbial biomass COD in different depths of the main layer (Langergraber et al., 2007).

3.4 2-stage horizontal and vertical flow CWs
Simulation results for a two-stage pilot-scale CW for treating surface water have been shown by
Langergraber and Šimůnek (2005) and are included here to show an application for a more
complicated flow domain geometry. The system was designed for the treatment of heavily polluted surface water. The total surface area of the two-stage system was 2 m², divided into downflow and upflow chambers (surface area is 1 m² each). The inlet was situated on the top of thedownflow chamber; the effluent was collected on the top of the upflow chamber by means of
perforated pipes (Perfler et al., 1999, Figure 5, left). An unstructured two-dimensional finite
74

element mesh was used for simulations consisting of 1135 nodes and 2057 elements (Figure 5,
right). Figure 6 shows simulation already presented in Langergraber and Šimůnek (2005).

Δ10

OUTLETsand - gravel 0/4 mm

drainage layer - gravel 16/32 m m

CHAMBER ADOWNFLOW

45

15

sand - gravel 0/4 mm

15

top layer - gravel 4/8 mm

55

INLET

geotextile

CHAMBER B
UPFLOW

Figure 5: Schematic representation of the two-stage vertical flow CW (left, values are in cm) and 2D
mesh of the two-stage vertical flow CW (right); Langergraber and Šimůnek (2005).
8
Top of the main layer

NH4N, NO3N (mg/l)

765NH4N measured

4

NO3N measured
3

NH4N simulated

2

NO3N simulated

1

0

10

20

30

40

50

60

70

0
0

10

20

30

40

50

60

70

Vertical flow distance (cm)Figure 6: Calculated steady-state distribution of Nitrosomonas (in mg/l, left), and simulated and measured
concentrations of NH4-N and NO3-N in a cross section of the downflow chamber (right); Langergraberand Šimůnek (2005).

3.5 Two-stage pilot-scale down- and upflow CW
Toscano et al. (2007) showed simulation results for a pilot-scale two-stage SSF CW for treatment of municipal wastewater from an experimental plant located in Eastern Sicily. The system
was used for secondary or tertiary treatment and consists of four parallel lines each with two
SSF beds in series. Each bed has a rectangular shape with a surface area of 4.5 m2 (1.5 m ×
3.0 m). For their study Toscano et al. (2007) used data from line 1 (planted with Phragmites sp.)and line 2 (unplanted), each comprising a HF and a VF bed in series. Figure 7 shows measured
and simulated data for a tracer experiment of line 2 (unplanted). The simulated breakthroughcurves matched the measured data well.
2,0Electrical conductivity (mS/cm)d)Measured effluent HF2

1,8

Measured effluent VF6
Simulated effluent HF2
1,6

Simulated effluent VF6

1,4

1,2

1,0
0

12

24

36

48

60

72

84

96

Time after start of the tracer experiment (h)Figure 7: Measured and simulated data for a tracer experiment of line 2 (Toscano et al., 2007).

Since no significant changes of the influent concentrations over time occurred, a further experiment was conducted by adding ammonium chloride (NH4Cl) as pulse tracer in the influent
75

of lines 1 and 2 functioning as secondary treatment. In Figure 8 measured and simulated NH4-N
concentrations are shown. The values for the parameters describing plant uptake and oxygen release obtained in the previous simulations have been used for the simulations. CW2D seems tobe able to predict the dynamics in the effluent of the HF bed for both lines (Toscano et al.,2007).
90

90
Measured influent

80

Measured influent

80

Measured effluent HF2

Measured effluent HF1

70

Measured effluent VF5

Measured effluent VF6
Simulated effluent HF2Simulated effluent HF1

60

NH4-N (mg/L)

NH4-N (mg/L)

70

Simulated effluent VF5

50
40
30

60

40
30

20

20

10

10

0

Simulated effluent VF650

006.06.07 09:00

06.06.07 15:00

06.06.07 21:00

07.06.07 03:00

07.06.07 09:00

06.06.07 09:00

06.06.07 15:00

06.06.07 21:00

07.06.07 03:00

07.06.07 09:00

Figure 8: Measured influent concentrations and measured and simulated effluent concentrations of NH4-N
after the addition of NH4Cl solution in lines 1 (left, planted) and 2 (right, unplanted); Toscano et al.
(2007)

3.6 Vertical flow CWs for treatment of combined sewer overflow
CW2D was used to model a VF CWs for the treatment of combined sewer overflow (CSO) by
Dittmer et al. (2005), Meyer et al. (2006) and Henrichs et al. (2007). Some of the differences between using wetlands for wastewater and for CSO treatments are the loading regime and thequality parameters of the inflow. For CSO treatment, the succession of loading events and dry
periods is characterized by the stochastic nature of rainfall and the runoff behavior of the catchment area. Extreme cases involve, on the one hand, a permanent loading for weeks, and several
months without any loading event on the other. To reach the main treatment objective, the detention and reduction of peak flows (Dittmer et al., 2005), throttle valves are applied to limit the
maximum effluent flow rate and therefore the flow velocity in the filter itself. For simulation
purposes, the maximum allowed effluent flow rate had to be implemented (Langergraber andŠimůnek, 2006).
The first simulations have been performed for lab-scale columns (diameter 19 cm) by Dittmer
et al. (2005). The columns were loaded once a week with 15.7 liter of synthetic sewer, whichcorresponds to a loading rate of 0.5 m/event. The filtration rate of columns 1 to 4 was controlled
manually by a throttle valve, whereas columns 5 and 6 had free drainage. Figure 9 compares
simulated effluent flow rates for a column with free drainage and a limited seepage face flux
boundary condition. The maximum flow in the controlled effluent rate case (right) corresponds
to 1 liter/hour. While in the unrestricted flow case, the bulk of water passed through the filter
bed in 1 hour, in the flow restricted case this happened in about 15 hours.
40

Effluent flow rate (l/h)

Effluent flow rate (l/s)

45

35
30
25
20
15
10

1
0.8
0.6
0.4
0.2

5

0

0
0

0.5

1

1.5

2

2.5

0

3

5

10

15

20

25

30

Time after loading (h)Time after loading (h)Figure 9: Comparison of simulations with free drainage BC (left) and controlled effluent rate (right);
Langergraber (2008).

Figure 10 (left) shows measured and simulated breakthrough curves of a tracer experimentcarried out for column 6 with free drainage, Figure 10 (right) for column 4 and column 2, with76

controlled effluent rates of 5 liters/hour and 1 liter/hour, respectively. Using the transport parameters obtained from the free drainage simulations (column 6) the tracer experiments with
controlled effluent rate could be simulated well.
10

12

Tracer concentration (mg/l)

Tracer concentration (mg/l)

987Measured6

Simulated
5
4
3
2

10

8

6
Column 4 (5 litre/h) - measured
4

Column 4 (5 litre/h) - simulated
Column 2 (1 litre/h) - measured

2

Column 2 (1 litre/h) - simulated

1
0

0
0

5

10

15

20

25

30

0

1

Cumulated effluent volume (litre)

2

3

4

5

6

7

Time after start of the tracer experiment (h)Figure 10: Measured and simulated breakthrough curves of column 6 with free drainage (left) for a controlled effluent rate of 5 and 1 liter/hour, respectively (right); Langergraber (2008).

For multi-component reactive transport simulations, good results could be achieved when
considering the different composition of CSO compared to domestic wastewater (Dittmer et al.,2005), and adsorption for COD fractions (Henrichs et al., 2007). Figure 11 shows, as an example, the simulated effluent concentrations of COD and its fractions for a lab-scale experimentwith a controlled effluent rate. Using the parameters obtained from lab-scale experiments, it wasalso possible to successfully simulate single events of field experiments.
30

0.E+00
Outflow
COD measured

2.E-05
3.E-05

COD simulated

2

20

Outflow (m /(m .s)

4.E-053

Effluent concentrations (mg/l)

1.E-0525

CI

15

5.E-05
6.E-05O2

10

7.E-05CS

8.E-055

9.E-05CR

1.E-04

0
0

5

10

15

20

25

30

Time after loading (h)Figure 11: Effluent concentrations of COD and its fractions with COD adsorption and controlled effluent
rate (O2: dissolved oxygen, CR, CS, CI: readily, and slowly biodegradable, and inert organic matter, respectively); Langergraber (2008).

4 SUMMARY
This paper presented an overview of applications of the multi-component reactive transportmodule CW2D for HYDRUS (Langergraber and Šimůnek, 2005). Different types of CWs havebeen investigated: VF and HF beds, hybrid systems and a two-stage down- and upflow sCWs
for different applications have been investigated using CW2D: treatment of municipal wastewater, polishing the effluent of a wastewater treatment plant for reuse purposes and treatment of
combined sewer overflow. CW2D was been applied for CWs with different sizes from lab-scalecolumns over pilot-scale systems to bigger outdoor systems.
In general, results using CW2D showed that water flow, single-solute transport and reactive
transport simulations resulted in good agreements with measured data for HF and VF CWs as
well as for hybrid systems. The practical applications have shown that simulation results matchthe measured data when the hydraulic behavior of the system can be described well. A good
match of experimental data to reactive transport simulations can then be obtained for CWs treating municipal wastewater using the values for the CW2D model parameters as given by Langer77

graber and Šimůnek (2005). Therefore it is advisable to measure at least the porosity and saturated hydraulic conductivity of the filter material to obtain reasonable simulation results for water flow (Langergraber, 2008).
REFERENCES
Dittmer, U., Meyer, D. & Langergraber, G. 2005. Simulation of a subsurface vertical flow constructedwetland for CSO treatment. Water Sci Technol 51(9), 225-232.
Henrichs, M., Langergraber, G. & Uhl, M. 2007. Modelling of organic matter degradation in constructedwetlands for treatment of combined sewer overflow. Sci Total Environ 380(1-3), 196-209.
Henze, M., Gujer, W., Mino, T. & van Loosdrecht, M.C.M. 2000. Activated sludge models ASM1,
ASM2, ASM2D and ASM3. IWA Scientific and Technical Report No.9, IWA Publishing, London,
UK.
Korkusuz, E.A., Meyer, D. & Langergraber, G. 2007. CW2D simulation results of lab-scale vertical flow
filters filled with special media and loaded with municipal wastewater. In: Mander Ü, Kóiv M, Vohla
C, editors. 2nd International Symposium on "Wetland Pollutant Dynamics and Control WETPOL
2007" – Extended Abstracts II, 16-20 September 2007, Tartu, Estonia, pp.448-450.
Langergraber, G. 2001. Development of a simulation tool for subsurface flow constructed wetlands.
Wiener Mitteilungen 169, Vienna, Austria, 207p, ISBN 3-85234-060-8.
Langergraber, G. 2005. The role of plant uptake on the removal of organic matter and nutrients in subsurface flow constructed wetlands – A simulation study. Water Sci Technol 51(9):213-223.
Langergraber, G. 2007. Simulation of the treatment performance of outdoor subsurface flow constructedwetlands in temperate climates. Sci Total Environ 380(1-3), 210-219.
Langergraber, G. 2008. Modeling of processes in subsurface flow constructed wetlands – A review. Vadoze Zone J (accepted).
Langergraber, G. & Šimůnek, J. 2005. Modeling variably saturated water flow and multicomponent reactive transport in constructed wetlands. Vadose Zone J 4(4):924-938.
Langergraber, G. & Šimůnek, J. 2006. The multi-component reactive transport module CW2D for constructed wetlands for the HYDRUS Software Package. Hydrus Software Series 2, Department of Environmental Sciences, University of California Riverside, Riverside, CA, USA, 72p.
Langergraber, G., Tietz, A. & Haberl, R. 2007. Comparison of measured and simulated distribution of
microbial biomass in subsurface vertical flow constructed wetlands. Water Sci Technol 56(3), 233240.
Mena, J. 2008. Tratamiento de aguas residuales urbanas y vínicas mediante humedales artificiales de flujo
subsuperficial. PhD thesis, University of Castilla-La Mancha, Ciudad Real, Spain;. [in Spanish].
Meyer, D., Langergraber, G. & Dittmer, U. 2006. Simulation of sorption processes in subsurface verticalflow constructed wetlands for CSO treatment. In: Proceedings of the 10th IWA Specialized Group
Conference on "Wetland Systems for Water Pollution Control" – Vol.1, 23-29 September 2006, Lisbon, Portugal, pp.599-609.
Perfler, R., Laber, J., Langergraber, G. & Haberl, R. 1999. Constructed wetlands for rehabilitation andreuse of surface waters in tropical and subtropical areas – First results from Small-Scale Plots using
vertical flow beds. Water Sci Technol 40(3), 155-162.
Šimůnek, J., Šejna, M. & van Genuchten, M. Th. 2006. The HYDRUS Software Package for Simulating Twoand Three-Dimensional Movement of Water, Heat, and Multiple Solutes in Variably-Saturated Media,
Technical Manual, Version 1.0, PC Progress, Prague, Czech Republic, 241p.
Tietz, A., Langergraber, G., Sleytr, K., Kirschner, A. & Haberl, R. 2007. Characterization of microbial
biocoenosis in vertical subsurface flow constructed wetlands. Sci Total Environ 380(1-3), 163-172.
Toscano, A., Langergraber, G. & Cirelli, G.L. 2006. Simulation of hydraulics and pollutant removal of apilot-scale two-stage constructed wetlands functioning as secondary or tertiary treatment. In: Proceedings of the 10th IWA Specialized Group Conference on "Wetland Systems for Water Pollution Control" – Vol.2, 23-29 September 2006, Lisbon, Portugal, pp.1303-1311.
Toscano, A., Langergraber, G., Consoli, S. & Cirelli, G.L. 2007. Modelling pollutant removal in a pilotscale two-stage subsurface flow constructed wetlands. Ecol Eng (submitted).

78

Water flow in the different types of unsaturated soils in Lithuania
Jakimavičiūtė-Maselienė V.1, 2, Mažeika J.1, Skuratovič Ž.1
1

Institute of Geology and Geography, T. Ševčenkos str. 13, LT-03223 Vilnius, Lithuania

2

Vilnius Gediminas Technical University, Saulėtekio al. 11, LT-10223 Vilnius, Lithuania

ABSTRACT: Two in previous monitoring programs experimentally characterized sites were selected for simulations of water flow in the unsaturated zone in Lithuania. Simulations were carried out using the finite element code HYDRUS-1D. Infiltration results were compared in terms
of pressure heads and water contents versus depth.

1 INTRODUCTION

RussiaSources of potential environmental risk located at the soil surface may impact groundwater via unsaturated zone processes. In risk assessment studies analysis of water flow in unsaturated zone canbe based on experimental observation and/or numerical simulation (Mallants et al., 1994). Experimental studies yield valuable insight into the physical processes, but their implementation is often
complicated. In Lithuania, experimental studies of unsaturated zone processes in different landscapes (four sites) were temporarily carried out in 1980-1995 (Žemaitis et al., 1995). Two sites having more complex data were briefly considered in this study (Fig. 1).
As an alternative to experimental studies,
numerical simulations often provide the necessary information to analyze the flow phenomSiauliaiena at a very high spatio-temporal resolution(Jury and Roth, 1990).
In particular, this study is one of the first
Lithuanian case studies dealing with unsatuKaunas
Tallinnrated flow simulations. There are plans to exDidziasalis site
Estonia
tend this computational approach during theVilniusongoing radioactive waste management proRiga
Latviajects which will be actively implemented durLithuania
Rykantai site
Vilniusing the decommissioning stage of Ignalina
PolandNPP. Furthermore, simulations may be used toplan and interpret experiments and to predict
future flow behavior for complex conditions
Figure 1. Location of the experimental sites
including natural and engineered barriers andtheir interactions in case of near surface repository for low and intermediate level radioactive waste
disposal.
2 SITES AND EXPERIMENTAL DATA
Two sites – Rykantai and Didziasalis – carefully experimentally characterized in previous monitoring programs (Žemaitis et al., 1995) were selected for the simulation (Fig. 1). Both sites are located79

in the Baltic Upland crossing in the eastern part of Lithuania, in the Nemunas River catchment andin sub-catchments of the Neris (Rykantai site) and Sventoji Rivers (Didziasalis site), respectively.
The region of East Lithuania as a whole is located in the eastern part of the Baltic artesian basin.
This region is part of the recharge area for the basin. The subsurface aquifer system, containing
fresh groundwater and taking part in active water exchange, is attributed to Quaternary intertill deposits. The unconfined aquifer is spread in the uppermost part of the Quaternary succession composed of inequigranular material. From the geomorphological point of view, the Rykantai site is located in the plain valley (prevailing glaciofluvial deposits); the Didziasalis site is attributed to thehummocky morainic upland (till deposits). The Lithuanian climate is tyrical of the northern part ofthe midlatitudal climate zone. According to the data of the Lithuanian Hydrometeorological Service, the annual sum of precipitation in the Dukstas station (close to Didziasalis site) is 594 mm for
1971-2000 and in the Vilnius station (close to Rykantai site) is 670 mm for the same period.
In the Rykantai site, the typical soil profile itself is composed of glaciofluvial inequigranularsand with groundwater table at a depth of 7 m. The soil profile of the Didziasalis site is represented
by marginal till deposits (clayey loam) in upper part and sand in lower part with groundwater table
at a depth of 10 m. (Fig. 2).

Figure 2. Soil profiles of the Didžiasalis site (A – main lithology and A’ – model geometry and parameters)
and the Rykantai site (B – main lithology and B’ – model geometry and parameters): 1 – topsoil formed oftill; 2 – loam; 3 – various sand; 4 – technogenic soil; 5 – clayey sand; 6 – sandstone; 7 – coarse sand; 8 – aperture; 9 – manometric tensiometer and its number; 10 – pore water catcher and its number; 11 – groundwaterlevel

The two most important soil physical characteristics are the water content, θ, and the soil water
potential head or hydraulic head, h. The hydraulic conductivity of the saturated soil (Ks, m/s), which
can be determined in the laboratory on undisturbed core samples, is another very important characteristic. Based on primary parameters, the soil hydraulic functions (θ(h) and K(h)) can be set. Description of water flow through unsaturated soil requires parameterization of θ(h) and K(h). Different soil types exhibit quite different hydraulic functions. When the soil becomes drier due tointernal drainage and/or evapotranspiration, air replaces water first in the coarse parts of the pore80

space and at lower (negative) values of the water potential also in the finer pores. In the unsaturated
zone, water is held in the soil pores by the so-called capillary forces (due the to surface tension of aliquid and the contact angle between liquid and solid phases). This leads to the existence of a negative pressure or matric head.
Undisturbed core samples were taken from the soil profiles at Rykantai and Didziasalis sites prior
to the installation of experimental equipments in lysimeters (Žemaitis et al., 1995). The K values
were estimated depending on soil type by constant head and falling head methods. The hydraulicfunctions (the water content and unsaturated hydraulic conductivity dependence on the matric head)
were estimated for different soils (Fig. 3).
1.0E+00

0.28
D
Water content

1.0E-03
1.0E-04
1.0E-05
1.0E-06-2.5

Hydraulic conductivity.
m/dHydraulic conductivity,
m/d

D

0.26

0.24
-1.5

-0.5

Pressure head, m

-3

-2

-1

0

0.3
R

R
Water content

1.0E-02

1.0E-02

1.0E-04

1.0E-06
-2.5

0.18
0.12
0.06

-1.5

-0.5

Pressure head, m

Pressure head, m

0.24

-3

-2
-1
0
Pressure head, mFigure 3. Estimated hydraulic functions for the Didziasalis site (clayey loam) profile (parametrization gives
approximately α=0.7 1/m, n=1.1, when measured values are as follows θr=0.08; θs=0.28, Ks= 0.0028 m/d)
and for the Rykantai site (sand) profile (parametrization gives approximately α=21 1/m, n=1.3, when measured values are as follows θr=0.01, θs=0.30, Ks=3.1 m/d) (D – Didziasalis site; R –Rykantai site)

The hydrometeorological observations (precipitation rate, air and soil temperatures) were carriedout continuously in the period of 1980-1995. The equipment installed in the lysimeters consisted ofmanometric tensiometers, porous porcelain cups for water sampling and thermometers. Next tolysimeters, the piezometers for groundwater level observations were installed.
3 MAIN FEATURES OF UNSATURATED FLOW SIMULATIONThe unsaturated flow simulation is based on soil hydraulic functions. To describe the soil water
characteristic and the unsaturated hydraulic conductivity function, the model defined by van
Genuchten (van Genuchten, 1980) was used. This model provides a predictive closed-form equation
for the unsaturated hydraulic conductivity function by using the statistical pore-size distribution
model of Mualem (Mualem, 1976). The degree of water saturation is defined:

Se =

θ −θr
θs −θrwhere θ is soil water content, θr and θs are the residual and saturated water contents, respectively.
The soil water characteristic and the unsaturated hydraulic conductivity relationship adjusted toexperimental data are shown in Fig. 3. For the soil groups, when initial estimates of the soil water
characteristic were not available, a set of average van Genuchten parameters estimated by Carsel
and Parish (Carsel and Parish, 1988) was used.
Under natural climatic conditions, water flow in soils varies in time and in space. The transient
behavior of water flux in the unsaturated soil, q (m/s), usually affects the water content θ and thematric pressure head h. Simulations for both soil profiles were carried out using the van GenuchtenMualem hydraulic functions without hysteresis. The simulations were performed using the finite
element code HYDRUS-1D (Simunek et al., 1998).

81

The top boundary condition corresponds to the long-term average net precipitation (precipitation
rate minus potential evapotranspiration). Potential evapotranspiration was calculated using theThornthwaite model (Ward and Elliot, 1995) using local hydrometeorological data. After that, thetop boundary flux for the Rykantai soil profile is set as a constant q of 0.00048 m/d and for theDidziasalis soil profile it is 0.00035 m/d. These values are maximal assuming zero surface run-off.
The groundwater table 7 m below the soil surface at the Rykantai site and 10 m below the surface atthe Didziasalis site was used as the bottom boundary condition. This is mathematically implemented by putting the pressure head = 0 at the water table. The Rykantai soil profile of 3 materials
is discretized in 0.07-m-thick layers, which results in 101 nodes for the one-dimensional simulation.
The Didziasalis soil profile of 3 materials is discretized in 0.1-m-thick layers, resulting in 101 nodes
for the one-dimensional simulation.
The top water flux is only one of many critical parameters in the natural environment. Infiltration
through the unsaturated zone occurs via the porous matrix and/or through open fissures and other“fast” pathways. The geometry of this porous network and the degree of water saturation provide an
opportunity for water to flow through both slow and fast pathways. In our case, we considered only
the porous matrix flow. In temperate climate, generally less than 5 to 25% of precipitation infiltrates
to the water table (Zimmermann et al., 1967). More details on soil profiles model and main input
parameters are shown in Fig. 2.
4 SIMULATION RESULTS AND WATER BALANCE
The numerically calculated infiltration in both soil profiles is shown in Figs. 4 and 5 in terms ofpressure heads versus depth (h – z) and water contents versus depth (θ - z). At t < 0, the initial pressure head is equal to the height above the water table, i.e. h = 0 at the groundwater table and h = -10
m and h = -7 mm at the soil surface respectively. The corresponding water content versus depth profile has the same shape as the soil water characteristic curve. At t > 0, water starts to infiltrate into
the soil and the infiltration front moves downward. After 870 days water content at the Didziasalis
profile reaches its steady-state value and a constant water flux is in the entire profile (0.00035
m/day). A less negative pressure head develops and reaches -0.024 m in the upper part and -0.34 m
in the lower part of the profile. After 760 days water content in the Rykantai profile reaches its
steady-state value and a constant water flux is in the entire profile (0.00048 m/day). A less negative
pressure head develops and reaches -0.4 m in the upper part and -0.7 m in the lower part of the profile. Under steady-state conditions the constant water flux is driven only by gravity (i.e., gravityflow).
0

0

1
2
3
4

1
-2

2
3

-4

Depth, m

Depth, m

-2

4
5

-6
-8

-4

5

-6

6

-8

6

-10

-10
-10

-8

-6

-4

-2

0

0.1

Pressure head, m

0.15

0.2

0.25

0.3

Water content

Figure 4. Calculated pressure heads (left) and water contents (right) versus depth at the Didziasalis profile.
Solid line is for t<=0 (initial condition), while other times are as follows: 1 – 100 days, 2 – 200 days, 3 – 300
days, 4 – 400 days, 5 – 500 days, 6 – 2000 days.

82

0

0

1

-1

2
3

-3

Depth, m

Depth, m

-2

1

-1

4

-4

5

-5

-2

3

-3

4

-4
-5

6

-6

2

5
6

-6

-7

-7
-7

-6

-5

-4

-3

-2

-1

0

0.05

Pressure head, m

0.1

0.15

0.2

0.25

Water content

Figure 5. Calculated pressure heads (left) and water contents (right) versus depth of the Rykantai profile.
Solid line is for t<=0 (initial condition), while other times are as follows: 1 – 100 days, 2 – 200 days, 3 – 300
days, 4 – 400 days, 5 – 500 days, 6 – 2000 days.

5 CONCLUSIONS
Experimental data and results of numerical simulations of unsaturated water flow for two Lithuanian soil profiles were preliminarily considered. Infiltration results were compared in terms of pressure heads and water contents versus depth. The time to reach steady-state water flux was calculated
to be between 760-870 days from the beginning of infiltration.
REFERENCES
Carsel R. F., Parish R. S. 1988. Developing joint probability distribution of soil water retention
characteristics. Water Resour. Res., 24(5):755769.
Environmental Hydrology. Ward A. D., Elliot W. J. (Eds.). Lewis publishers, Boca Raton, New
York, 1995.
Jury W.A., Roth K. 1990. Transfer functions and solute movement through soil. Theory and applications. Birkhäuser Verlag, Basel.
Mallants D., Vanclooster M., Meddahi M., Feyen J. 1994. Estimating solute transport in undisturbed soil columns using time-domain reflectometry. Journal of Contaminant Hydrology,
17:91-109.
Mualem Y. 1976. A new model for predicting the hydraulic conductivity of unsaturated porous
media. Water Resour. Res., 12:513-522.
Simunek J., Huang K, van Genuchten M. Th. 1998. The HYDRUS-1D code for simulating the onedimensional movement of water, heat, and multiple solutes in variably-saturated media. Version
6.0, Research Report No. 144, USSL-ARS-USDA, Riverside, California.
Van Genuchten M. TH. 1980. A closed-form equation for predicting the hydraulic condutivity ofunsaturated soils. Soil Sci. Soc. Am. J, 44:892-898.
Žemaitis V., Chomčenko R., Kučenko L., Kilda L., Dmukauskas A. 1995. Ataskaita už projektą“Aeracijos zonos monitoringas”, Vilnius.
Zimmermann U., Münnich K.O, Roether W. 1967. Downward movement of soil moisture traced bymeans of hydrogen isotopes. In: Isotope Techniques in the Hydrologic Cycle, Geophysical
Monograph Series 11, American Geophysical Union.

83

84

Simulation of the effect of time-variable water flow on thespeciation and mobility of reactive elements in the soil
D. Jacques & D. Mallants
Institute for Environment, Health and Safety, Belgian Nuclear Research Centre (SCK•CEN), Boeretang
200, B-2400 Mol, Belgium

J. Šimůnek
Department of Environmental Sciences, University of California Riverside, Riverside, CA 92521, USA

M. Th. van Genuchten
U.S. Salinity Laboratory, USDA, ARS, 450 W. Big Spring Rd., Riverside, CA 92507, USA

ABSTRACT: Speciation determines element mobility. It is commonly known that speciation
depends on the geochemical conditions in the soil. However, it is less well known that speciation is also influenced by the soil water content and the soil water fluxes. These dependencies
are numerically investigated using the reactive coupled simulator HP1 that was developed bycoupling HYDRUS-1D to the geochemical code PHREEQC-2.15. Two examples illustrate theeffect of time-variable soil water fluxes on the speciation and mobility of heavy elements, bothon the short term (seasonal variations) and on the long term leaching to the groundwater.

1 INTRODUCTION
The mobility and bio-availability of major and trace elements are the result of a complex interplay between physical, chemical and biological processes in the soil. The mobility of an element
depends on the chemical nature of the element itself and on the other elements in the aqueous
phase and the surface of the solid phase. Geochemical processes such as aqueous complexation,
precipitation and dissolution, ion exchange and surface complexation all play a role in the speciation, and consequently on the mobility, of elements in the soil system. Numerous studies exist
on speciation and mobility of elements. For example, Jansen (2003) investigated the mobility of
Al- and Fe- complexes with dissolved organic matter as a function of the pH, metal/ligand ratio
and redox state. At higher pH and a high metal/ligand ratio, the mobilities of Al and Fe decrease
due to the formation of insoluble metal – organic matter complexes. The water phase is the major (if not the only) carrier of many elements in the soil. Physical flow processes also affect themobility. The water flow velocity and direction is mainly dictated by the conditions at the soil
surface (precipitation and evapotranspiration) and the groundwater.
Apart from determining the average solute transport velocity and dispersion, transient flow
conditions substantially influence the geochemical conditions in the soil. This is supported bysome experimental evidence during laboratory experiments (e.g., Öztürk & Özkan, 2004), seasonal variations in the field (e.g., Berner et al., 1998), and artificial recharge of groundwater(Greskowiak et al., 2005). The interaction of discharging groundwater and shallow groundwater
dynamics with temporal variations in atmospheric conditions also affects the soil water composition (e.g., de Mars et al., 1997; de Mars & Garritsen, 1997; Joris, 2005). Slattery & Ronnfeldt(1992) showed that seasonal variations in pH were partly related to ionic strength. Gonçalves etal. (2006) reported an experimental study to quantify salinization and alkalization risks inlysimeters exposed to natural atmospheric conditions (rainfall and evapotranspiration) and85

manually irrigated with waters of different quality. Overall salinity and Na concentrations
changed significantly during the year. Na concentrations near the soil surface peaked during dry
periods after irrigation events, but subsequently decreased again to initial values during rainfall
periods.
To fully explore such coupled effects, a reactive coupled model is needed. In such tool, thespeciation of elements and solid surfaces are calculated together with water flow and solutetransport equations. It enables to account for the effects of spatial and temporal variations in water flow conditions on the speciation and mobility of elements in open systems such as soils. A
number of coupled reactive codes based on the HYDRUS software are available. HYDRUS-1D
coupled with UNSATCHEM model (Šimůnek & Suarez, 1994, 1997; Šimůnek et al., 1996,2005) and HYDRUS-2D coupled with CW2D (Langergraber & Šimůnek, 2005) are examples of
specific geochemistry codes coupled with HYDRUS. This paper describes HP1 (Jacques &Šimůnek, in prep.) which couples HYDRUS-1D with the generic geochemical code PHREEQC(Parkhurst & Apello, 1998). After a short introduction, two examples illustrate the effect of transient water flow on the speciation and mobility of elements in the soil at a seasonal and a longterm time scale.
2 DESCRIPTION OF THE HP1 CODE
The comprehensive simulation tool HP1 (version 2.1) has recently been developed by coupling
HYDRUS-1D with PHREEQC (Parkhurst & Apello, 1998). Both codes represent currently
state-of-the-science models in their domains of application. HYDRUS-1D version 4.0 is used tomodel water flow, solute transport and heat transport processes in the vadose zone. It is developed for soil systems by incorporating many features specific for soils such as atmosphericboundary conditions, several analytical functions describing the soil hydraulic properties for
granular or structured soils, models for water flow and solute transport in structured soils (dualporosity and dual-permeability models), and root water uptake processes, amongst others.
PHREEQC-2.15 is a geochemical code for low temperature aquatic geochemical processes
including aqueous complexation, precipitation and dissolution reactions, ion exchange, surface
complexation and solid-solutions. PHREEQC-2 provides a large flexibility in defining the geochemical system. It offers a choice of ion exchange models (Gaines & Thomas, Vanselow, orGapon convention), surface complexation models (non-electrostatic, electrostatic, diffuse double
layer models based on Dzombak and Morel (Dzombak & Morel 1990) or CD-MUSIC (Hiemstra& van Riemsdijk 1996) with explicit calculation of the diffuse layer), and multiple surfaces
models such as IE-2SNE (ion exchange–two non-electrostatic surface sites model of Bradbury& Baeyens, 2002), or multiple ion exchange models (e.g., Jacques et al., 2008b). Kinetic reactions and networks (both in parallel and in sequence) can be incorporated by using the BASICinterpreter in PHREEQC. This powerful feature allows the calculation of aqueous(bio)degradation and redox reactions, mineral dissolution and precipitation, kinetic adsorption
and even biomass growth and decay.
Both codes allow incorporating heterogeneity, both at the larger scale (e.g., different soil horizons) or at the smaller scale (dual-porosity, heterogeneous surfaces, different surfaces or processes in inter-aggregate and intra-aggregate domains). Consequently, HP1 can simulate also a
broad range of low-temperature biogeochemical reactions in water, the vadose zone and inground water systems, including interactions with minerals, gases, exchangers and sorption surfaces based on thermodynamic equilibrium, kinetic, or mixed equilibrium-kinetic reactions.
HP1 uses an operator-splitting approach for solving the governing equations, i.e., water flow,heat transport, solute transport and geochemical equations are solved in sequence within a single
time step (see Jacques et al., 2006, for a flowchart). Jacques et al. (2006) evaluated the accuracyof the operator-splitting approach in HP1 for a kinetic reaction network, and for mixed equilibrium – kinetic reactions for different flow conditions (steady-state and transient). To minimize
operator-splitting errors, two time step controlling parameters can be used. The maximum allowed time step is the most effective measure to reduce operator-splitting errors during time periods with a high evapotranspirative atmospheric demand. The performance index (i.e., theproduct of the grid Peclet and Courant numbers) minimizes operator-splitting errors during pre86

cipitation or infiltration events. Guidance for optimal spatial and temporal discretization for reactive transport modeling during transient flow conditions was given in Jacques et al. (2006).
A number of studies demonstrated the versatility of HP1 (Jacques & Šimůnek, in prep.) including heavy metal transport, phosphate and uranium transport and leaching (Jacques et al.
2008a,b) and transport of TNT and its daughter products (Šimůnek et al., 2006). This paper willshow selected results presented in these studies to illustrate the effects of water flow on speciation and mobility of elements. Details of the model assumptions are given in the cited papers.
3 EXAMPLE SIMULATIONS ILLUSTRATING THE EFFECT OF WATER FLOW ON
SPECIATION AND MOBILITY

3.1 Heavy metal migration in a podzol soil
Jacques et al. (2008a) simulated the transport of major cations (Na, K, Mg, and Ca) and heavy
metals (Cd and Zn) in a soil during transient water flow over a period of 30 years. The focus
was on the leaching of Cd and Zn in a dry Spodosol in the sandy region of Northern Belgium.
Soils in the region were contaminated by atmospheric deposition of Cd and Zn from non-ferrous
industry (Seuntjens, 2000). The various soil parameters and initial conditions were taken from
Seuntjens (2000). The soil consisted of six horizons for which the bulk density, organic matter
content, the parameters of the water retention characteristic and the saturated hydraulic conductivity were available. Geochemical processes taken into account were aqueous complexation andion exchange. The atmospheric boundary conditions between 1969 to 1998 were obtained fromclimatological data from two weather stations. Element concentrations in the rain water were obtained from Stolk (2001).
Figure 1 shows time series between 1972 and 1982 for the water content, pH, and the amountof Cl and Cd in the aqueous phase at two depths. The figure clearly illustrates a relation between
water flow conditions and changes in geochemistry in the soil. Alternating wet (precipitation)and dry (evaporation) conditions as dictated by the atmospheric boundary conditions influenced
the geochemistry near the soil surface. The mobile elements - such as the anion Cl, but alsomonovalent cations as Na – moved upwards and accumulate near the soil surface during periods
with high evaporative demand. At the same time, the pH decreased during the dryer soil conditions. Two factors contributed to this: a decreasing water content with increasing concentrationsof the elements, and the greater mobility of anions compared to cations and their accumulationnear the surface. The total amount of heavy metals near the soil surface remained nearly constant between seasons within a single year. However, the aqueous concentrations of the metals
varied significantly during the seasons. The changes in water content resulted in changes in concentrations which in turn caused changes in the cation exchange equilibrium. The increased
concentrations promoted monovalent cations to sorb on the exchange complex and, consequently, bivalent cations transferred to the aqueous phase. The competition between monovalentand bivalent cations for the exchange complex was further intensified by the increased amountof monovalent cations near the soil surface due to upward flow and transport. The increased
amounts and concentrations of Cl near the surface also resulted in a larger amount of Cd in the
aqueous phase due to the formation of relatively stable aqueous Cd-Cl complexes.
Thus, the physical factors of decreasing water content and upward water and solute fluxes
during periods with high atmospheric evaporative demand caused a pH decrease near the soil
surface and an increased amounts of Cd in the aqueous phase. The latter may have a significant
effect on the transport of Cd in a soil profile, and on its bioavailability since many uptake andbiological-controlled processes are concentration dependent.

87

5 cm
35 cm

0.25

θ

0.2
0.15
0.1

TCl,aq (mol / 1000 cm³)

0.3

-5

3x10

-5

2x10

1x10-50x100

0.05
TCd,aq (mol 1000 cm³)

4.9

pH

-5

4x10

4.7

4.5

-8

1x10

-9

5x10

0

4.3

0x10
1972

1974

1976

1978

1980

1982

1972

1974

Year

1976

1978

1980

1982

Year

Figure 1 Time series of water content θ, pH, total amounts of Cl and Cd in the aqueous phase per unit soilvolume (TCl,aq and TCd,aq) at 5 cm (A horizon) and 35 cm (Bh1-horizon) (adapted from Jacques et al.
2008a).

3.2 P and U migration in a soil profile from inorganic P-fertilizers
Inorganic P-fertilizers as monocalcium phosphate or triple superphosphate contain uranium due
to its natural occurrence in the raw material apatite. Overall, the environmental radiological impact of inorganic P-fertilizers is relatively modest. However, the mobility of U in the soil profile
depends on many factors such as pH, soil moisture, mineralogy and texture (Zielinski et al.,1997). Transient flow conditions represent a relatively unexplored factor effecting U mobility.
Jacques et al. (2008b) analyzed this factor by comparing long term (200 year) transient flow
simulations with steady-state flow conditions in terms of the linear distribution coefficient Kd
and U-fluxes to the ground water.
The same type of soil as in the previous example was used, i.e. a podzol consisting of six soil
horizons. Basic soil properties (the same as mentioned above and the Fe-oxide content) weretaken from Seuntjens (2000). The 30 year time series of atmospheric boundary conditions was
repeated seven times to obtain daily values of precipitation and potential evaporation for 200
years. For the steady-state flow simulation, the cumulative flux after 200 years was equal to thelong term net infiltration calculated during the transient flow simulation. The fertilizer scenario
assumed an average annual inorganic P-application representative for Flanders (Northern Belgium) for the last decade. In the scenario, P-fertilizers are applied each year on May 1 at a rateof 1 g P/m² in the form of Ca(H2PO4)2 in 1 cm of water (completely dissolved in the irrigationwater). Ca(H2PO4)2 contains 10-3 mol U, thus 3.77×10-6 mol U / m² was added annually.
The geochemical processes considered are aqueous complexation, ion exchange and surface
complexation reactions. Ion exchange involved the cations H+, Na+, K+, Mg2+, Ca2+ and UO22+.
Species involved in the surface complexation reactions were H+, Ca2+, Mg2+, UO22+, PO33-, SO42and F-. The cation exchange capacity was related to the amount of organic matter, assuming thatthe amount of exchangeable protons on organic matter was 6 meq g-1 (Tipping, 2000). The total
capacity of the complex was divided in six functional groups (each with a different proton dissociation constant) following the approach described in Appelo et al. (1998). The capacity of thesurface site was derived from the amount of Fe2O3 in the soil profile assuming 0.875 reactive
sites per mole Fe (Waite et al., 1994). Equilibrium constants were taken from Dzombak & Morel(1990).

88

100% U(VI) adsorbed

Total

CEC

75

SC

50

25
(b)

0
2

3

4
pH

5

6

Figure 2 pH-dependency of the U(VI) sorption in the A-horizon on both the cation exchange complex(CEC) and the surface complexation complex (SC) (adapted from Jacques et al., 2008b).

0

0

10

10

10

20
30
40

20
30
40

(a)

50
0.0x10Depth (cm)0Depth (cm)Depth (cm)Prior to the reactive transport simulation, the sorption of U on the solid phase was simulated.
Figure 2 shows the ratio of the total adsorbed U over the total U in the system as a function of
pH. The distribution between U adsorbed on the cation exchange complex and U on the surface
complexation sites are also shown. At low pH, most U is in the aqueous phase since both thecation exchange complex and the surface complexation sites are protonated. With increasing
pH, more and more sites become deprotonated and the U adsorption increases drastically. The
relative contribution of U adsorbed by cation exchange compared to surface complexation decreased with increasing pH. This is consistent with experimental findings of McKinley et al.
(1995). Above pH 6, sorption of U decreases again due to aqueous complexation of U with carbonate (Davis et al., 2004; Waite et al., 1994, results not shown). P is almost completely sorbed
in the same pH range (results not shown). It is expected that pH variations due to variations inwater content and fluxes will influence the mobility of U. Note that also the presence of competing elements influence the sorption of U. For example, when concentrations of U, Ca, and P
were present in the irrigation water with dissolved fertilizer, sorption of U significantly decreased.
Figure 3 shows depth profiles of total amounts of Ca, P and U after 25, 50 and 200 years fortransient and steady-state flow simulations. Ca accumulated mainly in the first horizon during
the first 50 years. During the next 100 years, the amount of Ca in the two Bh horizons gradually
increased. It moved a little bit deeper in case of transient flow simulations. P migrated fasterthrough the A and E horizons compared to Ca due to the small amount of Fe2O3 in these horizons (low adsorption potential for P in these horizons). In the Bh-horizon with a high Fe2O3
amount, the migration of P significantly slowed down. P did not reach the BC horizon after 200
years. U migration was quite different between two simulations, especially during the first 50
years. U transport was again slightly faster for transient flow simulations than during steadystate flow.

0

8.0x10-3

1.6x10-2

Ca (mol / 1000 cm³ soil)

20
30

25 year
50 year
200 year

40
(b)

50
0.0x100-2

4.0x10-2

8.0x10P (mol / 1000 cm³ soil)

(c)

50
0

0.0x10-6

4.0x10-6

8.0x10U (mol / 1000 cm³ soil)

Figure 3 Concentration profiles versus depth of total Ca, P and U at selected times for the steady-stateflow simulation (gray lines) and the transient flow simulation (black lines) (adapted from Jacques et al.
2008b).
89

Figure 4 shows time series of the pH and the Kd of Ca, P, and U at a depth of 5 cm. For both
simulations, an observed long term trend is related to long time changes in Ca, P and U in thesoil solution due to addition of the P-fertilizer. For steady-state flow conditions, steady-state Kdvalues at a 5-cm depth were obtained after 40 to 50 years for P and U, and after 100 years for
Ca. This corresponds to the mobility of the three elements through top horizons as discussed in
Figure 3. In this first horizon, the effect on Kd of an increased competition for sorption due tohigher concentrations is larger than the effect of the pH increase. The latter would result in a Kd
increase whereas the former resulted in a Kd decrease. A different picture was observed fordeeper horizons (see Jacques et al., 2008b).
The figure also illustrates the effect of transient flow conditions on short term variations in
pH and Kd. Seasonal variations in pH were induced by temporal variations in water contents andwater fluxes as discussed in the previous example. These pH variations then influenced the mobility and Kd values of mainly Ca and U. The effect on the Kd of P was small since the Padsorption is almost pH-independent. The amplitude of short-time variations of Kd at the sametime decreased with time. The adopted sorption isotherm is typically nonlinear, having a slope(the first derivative) that gradually decreases with increasing aqueous concentration. The Kd, being the slope of the isotherm at a particular concentration, then also decreases with in creasing
aqueous concentration. In addition, variations in Kd with changing concentrations (the second
derivative of the isotherm) go to zero with increasing aqueous U concentration.

3

U (Bq/y/m²))

U (Bq/y/m²)

3

1x10
2
1x10
1x101
0
1x10
-1
1x10
-2
1x10
-3
1x10
-4
1x10
-5
1x10
-6
1x10
-7
1x10
-8
1x10
-9
1x10
-10
1x10

7 cm
0

25

U (Bq/y/m²)

-3

1x10
-4
1x10
-5
1x10
-6
1x10
-7
1x10
-8
1x10
-9
1x10
-10
1x10
1x10-11
-12
1x10
-13
1x10
-14
1x10

50
Time (year)

75

100

80
120
Time (year)

160

200

1x10
2
1x10
1x101
0
1x10
-1
1x10
-2
1x10
-3
1x10
-4
1x10
-5
1x10
-6
1x10
-7
1x10
-8
1x10
-9
1x10
-10
1x10

28 cm
0

40

80
120
Time (year)

160

200

100 cm040

Figure 4 Time series of pH and the Kd of Ca, P and U at a depth of 5 cm for the steady-state flow simulation (dashed line) and the transient flow simulation (full line) (adapted from Jacques et al. 2008b).

Figure 5 shows the U-fluxes at three depths for the bottom of the A (7 cm), Bh2 (28 cm) andC2 (100 cm) horizons for both transient and steady-state flow conditions. Steady-state fluxes
were obtained after 40 and 100 years for the A and Bh2 horizon, respectively. For the transient
flow simulation, solute fluxes continued to vary by almost three orders of magnitude after
steady-state solute transport conditions were reached for the steady-state flow simulation. While
U reached the bottom of the soil profile after 80 years, a steady-state flux was not yet obtained
within the 200-year simulation period. U fluxes increased earlier for the transient flow simulation, which reflects faster leaching due to interactions caused by short-time variations in watercontents, water fluxes and geochemical conditions.

90

4

4.8

1x10

Kd of Ca (l kg -1)

4.4

pH

4
3.6
3.2
(a)

2.8
0

25

50
Time (years)

75

1x102
1

1x10

(b)

1x100

100

0

25

1x104

Kd of P (l kg -1)

Kd of U (l kg -1)

1x103

1x103

50
Time (years)

75

100

1x103
2

1x10

1

1x10
(c)

1x102
0

25

50
Time (years)

75

(d)

1x100
0

100

25

50
Time (years)

75

100

Figure 5 Downward U fluxes (Bq/y/m²) at the bottom of the A (7 cm, 100 years), Bh2 (28 cm, 200 years),and C2 (100 cm, 200 years) horizons for the transient (dots) and steady-state (dashed line) flow simulations (adapted from Jacques et al. 2008b).

4 CONCLUSIONSThe two presented examples illustrate the effect of water flow processes on the speciation andmobility of elements in the soil profile. Variable water contents and fluid fluxes during transient
flow in the vadose zone significantly affect prevailing geochemical conditions such as soil pH.
Consequently, the aqueous speciation and the adsorption potential and thus the element mobility
may change with time. Also the higher concentration of more mobile elements near the soil surface originating from upward water flow may increase the aqueous concentrations of, for example, heavy metals as Cd. Both cations (by increased competition for adsorption) and anions(forming aqueous complexes) play a role in that process.
REFERENCES
Appelo, C.A.J., Verweij, E., & H. Schäfer. 1998. A hydrogeochemical transport model for an oxidation
experiment with pyrite/calcite/exchangers/organic matter containing sand. Applied Geochemistry 13:257-268.
Berner, R.A., Rao, J.-L., Chang, S., O’Brien, R., & Keller K,. 1998. Seasonal variability of adsorption
and exchange equilibria in soil waters. Aquatic Geochemistry 4: 273-290.
Bradbury, M.H., & Baeyens, B. 2002. Sorption of Eu on Na- and Ca-montmorillonites: Experimental investigations and modelling with cation exchange and surface complexation. Geochimica et Cosmochimica Acta 66: 2325-2334
Davis, J.A., Meece, D.E., Kohler, M., & Curtis, G.P. 2004. Approaches to surface complexation modeling of Uranium(VI) adsorption on aquifer sediments. Geochimica et Cosmochimica Acta 68: 36213641.
de Mars, H., & Garritsen, G. 1997. Interrelationship between water quality and groundwater flow dynamics in a small wetland system along a sandy hill ridge. Hydrological Processes. 11: 335-351.
de Mars, H., Wassen, M.J., & Olde Venterink, H. 1997. Flooding and groundwater dynamics in fens ineastern Poland. Journal of Vegetation Science 8: 319-328.
Dzombak, D.A., & Morel, F.M.M. 1990. Surface complexation modeling – Hydrous ferric oxide. John
Wiley, New York.
91

Greskowiak, J., Prommer, H., Massmann, G., Johnston, C.D., Nützmann, G., & Pekdeger, A. 2005. The
impact of variably saturated conditions on hydrogeochemical changes during artificial recharge ofgroundwater. Applied Geochemistry 20: 1409-1426.
Gonçalves, M.C., Šimůnek, J., Ramos, T.B., Martins, J.C., Neves, M.J., and Pires, F.P. 2006. Multicomponent solute transport in soil lysimeters irrigated with waters of different quality. Water Resources
Research 42, W08401, doi:10.1029/2005WR004802.
Hiemstra, T. & Van Riemsdijk, W.H. 1996. A surface structural approach to ion adsorption: the charge
distribution (CD) model. Journal of Colloid and Interface Science 179: 488-508.
Jacques, D., Šimůnek, J., Mallants, D., & van Genuchten, M.Th. 2006. Operator-splitting errors in coupled reactive transport codes for transient variably saturated flow and contaminant transport in layeredsoil profiles. Journal of Contaminant Hydrology 88: 197-218.
Jacques, D., & Šimůnek, J. 2008. User Manual of the Multicomponent Variably-Saturated Transport
Model HP1 (Version 2.0): Description, Verification and Examples. SCK•CEN, Mol, Belgium (inpreparation).
Jacques, D., Šimůnek, J., Mallants, D., & van Genuchten, M.Th. 2008a. Detailled modeling of coupledwater flow, solute transport and geochemical reactions: Migration of heavy metals in a podzol soil profile. Geoderma (in press).
Jacques, D., Šimůnek, J., Mallants, D., & van Genuchten, M.Th. 2008b. Coupling hydrological andchemical proceses in the vadose zone: A case study on long term uranium migration following mineral
P-fertilization. Vadose Zone Journal (in press).
Jansen, B., 2003. The mobility of aluminium, iron and organic matter in acidic sandy soils. Ph.D., University of Amsterdam, Amsterdam, The Netherlands.
Joris, I. 2005. Soil water dynamics in alluvial wetlands: field and modelling study. Ph.D. no. 646
Faculteit Bio-ingenieurswetenschappen, University of Leuven, Belgium.
Langergraber, G. & Šimůnek, J. 2005. Modeling Variably-Saturated Water Flow and Multi-Component
Reactive Transport in Constructed Wetlands. Vadose Zone Journal 4: 924-938.
McKinley, J.P., Zachara, J.M., Smith, S.C., & Turner, G.D. 1995. The influence of uranyl hydrolysis andmultiple site-binding reactions on adsorption of U(VI) to montmorillonite. Clays and Clay Mineralogy
43: 586-598
Öztürk, H.S., & Özkan, L. 2004. Effects of evaporation and different flow regimes on solute distribution
in soil. Transport in Porous Media 56: 245-255.
Parkhurst, D.L., & Appelo, C.A.J. 1999. User's guide to PHREEQC (Version 2) – A computer programfor speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations. Water-Resources Investigations, Report 99-4259, Denver, Co, USA.
Seuntjens, P., 2000. Reactive solute transport in heterogeneous porous medium: Cadmium leaching inacid sandy soils. Ph.D., University of Antwerp, Belgium.
Šimůnek, J., & Suarez, D. L. 1994. Two-dimensional transport model for variably saturated porous media withmajor ion chemistry. Water Resources Research. 30: 1115-1133.
Šimůnek, J., Suarez, D. L., & Šejna, M. 1996. The UNSATCHEM software package for simulating onedimensional variably saturated water flow, heat transport, carbon dioxide production and transport,
and multicomponent solute transport with major ion equilibrium and kinetic chemistry, Version 2.0.
Research Report No. 141, U.S. Salinity Laboratory, USDA, ARS, Riverside, California, 186 pp.
Šimůnek, J., & Suarez, D.L. 1997. Sodic soil reclamation using multicomponent transport modeling.
ASCE Journal of Irrigation and Drainage Engineering 123: 367-376.
Šimůnek, J., van Genuchten, M.Th., & Šejna, M. 2005. The HYDRUS-1D software package for simulating the one-dimensional movement of water, heat, and multiple solutes in variably-saturated media.
Version 3.0. HYDRUS Software Series 1, Department of Environmental Sciences, University of California Riverside, Riverside, CA, 270 pp.
Šimůnek, J., Jacques, D., van Genuchten, M.Th., & Mallants, D. 2006. Multicomponent geochemicaltransport modelling using HYDRUS-1D and HP1. Journal of American Water Resources Association
46: 1537-1547.
Slattery, W.J., & Ronnfeldt, G.R. 1992. Seasonal variation of pH, aluminium and manganese in acid soils
from north-eastern Victoria. Australian Journal of Experimental Agriculture 32: 1105-1112.
Stolk, A.P., 2001. Landelijk Meetnet Regenwatersamenstelling, Meetresultaten 1999. RIVM, The Netherlands, RIVM Rapport 723101 056.
Tipping, E., 2002. Cation binding by humic substances. Cambridge University Press, Cambridge, UK.
Waite, T.D., Davis, J.A., Payne, T.E., Waychunas, G.A., & Xu, N. 1994. Uranium(VI) adsorption to ferrihydrite: Application of a surface complexation model. Geochimica et Cosmochimca Acta 58: 54655478.
Zielinski, R.A., Asher-Bolinder, S., Meier, A.L., Johnson, C.A., & Szabo, B.J. 1997. Natural or fertilizerderived uranium in irrigation drainage: a case study in southeastern Colorado, U.S.A. Applied Geochemistry 12: 9-21.
92

Evaluating Evapotranspiration from Meteorological Data Using
HYDRUS-1D
M. Sakai & J. Šimůnek
Department of Environmental Sciences, University of California Riverside, Riverside, CA 92521, USA

H. Saito
Department of Ecoregion Science, Tokyo University of Agriculture and Technology, Fuchu, Tokyo 1838509, Japan

ABSTRACT: Evapotranspiration rates strongly affect movement of water and solute in the vadose zone. We have implemented into the HYDRUS-1D code Penman-Monteith and Hargreaves equations to estimate potential evapotranspiration rates. While Penman-Monteith equation requires information about net radiation, wind speed, minimum and maximum temperaturesand relative humidity, Hargreaves equation requires only minimum and maximum temperatures.
The developed code provides users with a lot of flexibility by allowing an input of various combinations of meteorological information. The performance of the code is demonstrated using
meteorological data for Tokyo, Japan (humid climate) and Riverside, CA (semi-arid climate).
Liquid water and water vapor flow, heat transport, root water uptake, and evapotranspiration
rates were quantitatively evaluated using meteorological data publicly available from local
weather stations.
1 INTRODUCTION
Evaluating evapotranspiration rates from soils is needed for many research applications involving lysimeter and field studies, as well as for water management of agricultural fields. Potential
evapotranspiration rates from vegetated soils can be calculated from meteorological data using
either the Penman-Monteith combination equation (Monteith and Unsworth, 1990; FAO, 1990)
or Hargreaves equation (Jensen et al., 1997). Since actual evaporation from soils and transpiration from crops are additionally also affected by water contents in the root zone, their evaluation
requires simultaneous calculations of water flow and root water uptake in soils.
Penman-Monteith and Hargreaves equations were implemented into the HYDRUS-1D software package (Šimůnek et al., 2008). Various combinations of meteorological and crop data canbe entered in the graphical user interface (GUI) of HYDRUS-1D to evaluate potential
evapotranspiration rates. Example simulations are presented here to demonstrate this new feature of HYDRUS-1D and to evaluate evaporation rates, transpiration rates, water contents, androot water uptake using readily available data during the year 2007 from weather stations at Tokyo in Japan and Riverside in California.
2 MODELS
2.1 Potential evapotranspiration
The FAO-recommended Penman-Monteith combination equation is as follows (FAO, 1990):

ρ c p (e a - e d ) / r a ⎤1 ⎡ Δ ( R n - G)(1)+⎢
λ ⎣ Δ + γ (1 + r c / r a ) Δ + γ (1 + r c / r a ) ⎥⎦where ETp is the potential evapotranspiration rate [mm d-1], λ is the latent heat of vaporization [MJ
ET p =

93

kg-1], Rn is the net radiation at surface [MJ m-2d-1], G ( ≈ 0 in this study) is the soil heat flux [MJ m2 -1d ], ρ is the atmospheric density [kg m-3], cp is the specific heat of moist air [kJ kg-1oC-1], ea is thesaturation vapor pressure [kPa], ed is the actual vapor pressure [kPa], Δ is the slope of the vapor
pressure curve [kPaoC-1], γ is the psychrometric constant [kPaoC-1], rc is the crop canopy resistance[s m-1], and ra is the aerodynamic resistance [s m-1].
The potential evapotranspiration can also be evaluated using the much simpler Hargreaves equation (Jensen et al., 1997; Droogers and Allen, 2002):

ETp = 0.0023 ⋅ 0.408 ⋅ Ra (Tm + 17.8 ) TR(2)

where Ra is the extraterrestrial radiation [MJ m-2s-1], Tm is the daily mean air temperature [oC], andTR is the temperature range between mean daily maximum and minimum air temperatures [oC].
Potential evaporation, Ep, and transpiration, Tp, fluxes can be calculated from potential evapotranspiration using Beer’s law (Ritchie, 1972):

Tp = ETp ⋅ SCF = ETp (1 − e − k ⋅LAI )(3)

E p = ETp ⋅ (1 − SCF ) = ETp e − k ⋅LAI

where SCF is the soil cover fraction [-], LAI is the leaf area index [-], k ( = 0.463) is the constantgoverning radiation extinction by the canopy [-].
2.2 Root water uptake
Root water uptake rate can be defined according to Feddes et al. (1978) as follows:

S = α ( h ) S p = α ( h ) b ( z ) Tp

(4)

where S is the sink term in Richards’ equation [s-1], α (h) is the water stress response function [], Sp is the potential root water uptake rate [s-1], and b(z) is the normalized water uptake distribution [m-1].
3 HYDRUS-1D GUI FOR METEOROLOGICAL INFORMATION
Different meteorological data are available for different applications. While a complete set of meteorological variables (net radiation, wind speed, temperature, relative humidity, vapor pressure,etc) may be available for one application, limited information (or no information at all) may beavailable for other applications. FAO suggested various equations to compensate for the lack ofdata (FAO, 1990). Some of these, and others, were implemented in the version 4 of HYDRUS-1D
to calculate potential evapotranspiration
rates from several combinations of data.
Figure 1 shows the “Meteorological Parameters” dialog window of the HYDRUS1D GUI where users select what type of information they have available. While parameters in the upper part of the window relate mainly to net radiation calculations,parameters in the lower part concern crop
information (crop height, root depth, and
LAI). Figure 2 shows the “Meteorological
Conditions” dialog window of the HYDRUS-1D GUI where one inputs time variable meteorological data (e.g. radiation, airtemperature, and relative humidity, etc) andtime variable crop data (crop height, root Figure 1. The Meteorological Parameters dialog
window of the HYDRUS-1D GUI.
depth, and LAI).
94

Figure 2. The Meteorological Conditions dialog window of the HYDRUS-1D GUI.

4 EXAMPLE SIMULATIONS
5

30
LAI
Root depth

LAI (-)

4

20

3
2

10

1
0

0
0

100

200

300

DOY

Figure 3. LAI and root depth as a function oftime for example simulations.

1

ETp (cmd-1)

Penman-Monteith
Hargreaves

Penman-Monteith
Hargreaves

0.8
0.6
0.4
0.2
0

RiversideTokyo

ETp (cmd-1)

0.8

Potential evaporation
Potential transpiration

Potential evaporation
Potential transpiration

0.6
0.4
0.2
0
0

100

200

300

0

100

200

300

DOY

DOY

Figure 4. Estimated potential evapotranspiration rates (top), and potential evaporationand transpiration rates (bottom) for the Tokyo (left) and Riverside (right) sites during
the year 2007.
95

Root depth (cm)Model performance is demonstrated using
meteorological data for 2007 for two sites,
i.e., Tokyo in Japan (35º41’N, 139º46’E, Altitude 6 m) and Riverside in California(33º58’N, 117º21’W, Altitude 311 m). Daily
maximum and minimum air temperatures,
averaged relative humidity, wind speed, sunshine hour, and rainfall data for Tokyo were
downloaded from the Japan Meteorological
Agency. Daily maximum and minimum airtemperatures, averaged relative humidity,
wind speed, solar radiation, and rainfall data
for Riverside were downloaded from the
University of California, Agriculture andNatural Resources. Time variable LAI androot depth shown in Figure 3 were used as
crop data for both sites.

Potential evapotranspiration rates, ETp, calculated using Penman-Monteith equation (eq. (1)),
and precipitations were used as a surface boundary condition. Free drainage was used as the lower
boundary condition. Loam and Sandy loam from the HYDRUS catalog were used for soil hydraulic properties of Tokyo and Riverside, respectively. Coupled movement of liquid water, water vapor, and heat (Saito et al., 2006) was simulated using HYDRUS-1D.
Figure 4 (top) shows ETp calculated using both Penman-Monteith equation (eq. (1)) and Hargreaves equation (eq. (2)). Although ETp calculated using Penman-Monteith equation fluctuatedmore than the one calculated using Hargreaves equation, reflecting more pronounced changes ofradiation, relative humidity, and wind speed, the two estimated values agreed well. ETp for Riverside was larger than for Tokyo, because of its semi-arid climate. Figure 4 (bottom) shows potential
evaporation and transpiration rates obtained by dividing ETp calculated with Penman-Monteithequation using Beer’s law (eq. (3)). The potential transpiration rate increased from 90th day to250th day depending on the growth of the crop (LAI and root depth).
Figure 5 shows simulated actual evaporation rates, transpiration rates, and water contents forTokyo and Riverside. Both in Tokyo and Riverside, large evaporation or transpiration rates occurred immediately after precipitations. Since Tokyo has large amounts of precipitation all year
round, soil water evaporation was large throughout the year, while transpiration reflected thegrowth of the crop. Water contents in the soil also fluctuated throughout the year due to precipitation events and evapotranspiration. Water content changes at a depth of 50 cm were smaller than atother depths, because of a decreasing impact of evaporation and root water uptake. On the other
hand, it was raining only during the first 110 days of the year in the semi-arid climate of Riverside.
Therefore actual evaporation rates were much smaller in Riverside than in Tokyo in spite of its larger values of potential evaporation rate. The transpiration rate was negligibly small in Riverside
because of the shortage of soil water in the root zone.

0.4

4
6

0.2

8
10

Tokyo
0.4
0.3
0.2
5cm
15cm
50cm

0.1
0
0

100

200
DOY

0.2

Evap.
Trans.
Precip.

0.4

0.1
0.6
0

Water content (cm3cm-3)

Water content (cm3cm-3)00.2
ET (cmd-1)

ET (cmd-1)

Evap.
Trans.
Precip.

Precipitation (cmd-1)

2

0.6

0

0.3

0

Precipitation (cmd-1)

0.8

0.8

Riverside
5cm
15cm
50cm

0.12

0.08

0.04

300

0

100

200
DOY

300

Figure 5. Precipitations and simulated actual evaporation and transpiration rates (top),
and water contents (bottom) for the Tokyo (left) and Riverside (right) sites during theyear 2007.

SUMMARY
This study demonstrates capabilities of the new version (4.0) of the HYDRUS-1D program that
implements Penman-Monteith equation for calculating potential evapotranspiration rates. Easy estimations of evapotranspiration rates and water contents in the root zone using readily available
meteorological information can be useful for water management decisions involving irrigation anddrainage systems.
96

REFERENCES
Droogers, P., and R. G. Allen, Estimating reference evapotranspiration under inaccurate data conditions,
Irrigation and Drainage Systems, 16, 33-45, 2002.
Feddes, R. A., E. Bresler, and S. P. Neuman, Field test of a modified numerical model for water uptake
by root systems, Water Resour. Res., 10(6), 1199-1206, 1974.
Food And Agriculture Organization of the United Nations, Expert consultation on revision of FAO methodologies for crop water requirements, ANNEX V, FAO Penman-Monteith Formula, Rome Italy, 1990.
Jensen, D. T., G. H. Hargreaves, B. Temesgen, and R. G. Allen, Computation of Eto under nonideal conditions, J. of Irrig. Drainage, 123(5), 394-400, 1997
Monteith, J. L. and M. H. Unsworth, Principles of Environmental Physics, Edward Arnold, London, 1990.
Ritchie, J. T., Model for predicting evaporation from a row crop with incomplete cover, Water Resour.
Res., 8(5), 1204-1213, 1972.
Saito, H., J., J. Šimůnek, and B. Mohanty, Numerical analyses of coupled water, vapor and heat transport in
the vadose zone, Vadose Zone Journal, 5, 784–800, 2006.
Šimůnek, J., M. Šejna, H. Saito, M. Sakai, and M. Th. van Genuchten, The HYDRUS-1D Software Package for Simulating the Movement of Water, Heat, and Multiple Solutes in Variably Saturated Media,
Version 4.0, HYDRUS Software Series 3, Department of Environmental Sciences, University of California Riverside, Riverside, California, USA, pp. 315, 2008.

97

98

Impact of soil micro-morphological features on character of waterflow and solute transport
R. Kodešová, M. Kočárek, J. Kozák
Department of Soil Science and Geology, University of Life Sciences, Kamýcká 129, 16521 Prague 6,
Czech Republic

V. KodešDepartment of Water Quality, Czech Hydrometeorological Institute, Na Šabatce 17, 14306 Prague 6,
Czech Republic

J. Šimůnek
Department of Environmental Sciences, University of California Riverside, Riverside, CA 92521, USA

N. Vignozzi, M. Pagliai
CRA-ABP, Piazza M Azeglio 30, I-50121 Florence, Italy

ABSTRACT:
Soil pore structure and water flow and chlorotoluron transport were studied in humic
horizons of various soil types: Haplic Luvisol, Greyic Phaeozem, Haplic Cambisol.
Single-porosity, and either dual-porosity or dual-permeability, flow and transportmodels in HYDRUS-1D were used to estimate the soil hydraulic parameters fromlaboratory multi-step outflow experiments, laboratory and field ponding infiltration
experiments via numerical inversion. Appropriate models were selected based on thesoil micromorphological study. Herbicide (chlorotoluronu) transport parameters were
measured using a standard laboratory procedure and some of them were estimated using
numerical inversion from to transient data obtained from laboratory ponding infiltration
experiments. Resulting parameters were used to simulate the herbicide transport
experimentally studied in the field.

99

100

What comes before upscaling? - Modeling transport andspeciation in the rhizosphere with the new tool RhizoMath
K. Szegedi, D. Vetterlein & H.-U. Neue
Department of Soil Physics, UFZ – Helmholtz Centre for Environmental Research, Halle/Saale, Germany
Institute of Agricultural and Nutritional Sciences, Martin Luther Univ. Halle-Wittenberg, Halle, Germany

R . Jahn
Institute of Agricultural and Nutritional Sciences, Martin Luther Univ. Halle-Wittenberg, Halle, Germany

ABSTRACT:
Plants are one of the "principal performers" in the vadose zone as their activity induces
or influences water flow and solute transport in the soil. The most striking feature of therhizosphere is the formation of solute gradients extending from the root surface into thesurrounding soil that occur as depletion and accumulation of solutes at the root surface,
as well as gradients in the soil pH-value. The root uptake rate of a given nutrient isprimarily a function of its concentration in the soil solution at the root surface that
depends on the combined effects of many plant and soil processes and differs from theconcentration in the bulk soil. Understanding the processes that occur in the vicinity ofroots is also essential for modeling transport at larger scales. The RhizoMath (Szegedi
et al, 2008) approach for modeling co-occurring processes in the rhizosphere, including
speciation in the soil solution, is based on coupling the mathematical package
MATLAB with the geochemical code PHREEQC (Parkhurst & Appelo, 1999).
Different geochemical models (with and without charge balance) and geometries(planar and radial) are included and a built-in initialization module performs calibration
against experimental data. RhizoMath was applied to describe a compartment system
experiment involving a complex speciation with a special emphasis on the interactionsof arsenate and phosphate with the soil substrate. RhizoMath is able to describe
concentration gradients found in experiments (Vetterlein et al. 2007). However, a moresophisticated plant response mode still has to be included.

REFERENCES
Szegedi, K., Vetterlein, D., Nietfeld, H., Jahn, R. & Neue, H.-U. 2008. The New Tool RhizoMath forModeling Coupled Transport and Speciation in the Rhizosphere, Vadose Zone Journal, Special Issue“Vadose Zone Modeling”. In press. DOI: 10.2136/vzj2007.0064.
Parkhurst, D.L. & Appelo, C. A. J. 1999. User's guide to PHREEQC (version 2) - A computer program forspeciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations. USGS WaterResources Investigations Report 99-4259, Denver, Colorado.
Vetterlein, D., Szegedi, K., Ackermann, J., Mattusch, J., Neue, H.-U., Tanneberg H. & Jahn, R. 2007.
Competitive mobilisation of phosphate and arsenate associated with goethite by root activity. Journal ofEnvironmental Quality 36, 1811-1820.

101

102

Field data and modeling of salt transport in a coastal irrigatedplain.

F. SLAMA1,2, R. Bouhlila1 & J. Tarhouni3
1 Laboratoire de Modélisation en Hydraulique et Environnement, Ecole Nationale d’Ingénieurs de Tunis, BP 37 Le
Belvédère 1002, Tunis, Tunisia
fairouz.slama@unine.ch
2 Centre of Hydrogeology (CHY), University of Neuchâtel,. 11, rue Emile-Argand. CH – 2007, Neuchâtel, Switzerland
3 Institut National Agronomique de Tunisie (INAT), 43, Avenue Charles Nicolle 1082 -Tunis- Mahrajène Tunisia

ABSTRACT
Groundwater and soil resources are being threatened by important degradation due tosalinisation processes. This is truer in the case of regions lacking for water resources andcharacterized by arid and semi-arid climate with annual precipitation that does not exceed
500 mm and with important evapotranspiration reaching annual value of 1100 mm. Thissituation, also aggravated by irrigation expansion, leads to overexploiting water
resources. Earlier researches have been performed to study soil and groundwatersalinisation processes. However, in most studies, the unsaturated zone is assimilated to
the root zone not exceeding 2 meters. Moreover groundwater and soil are in most cases
treated separately. In this study, we propose to model both groundwater and soil salinity
in the purpose of understanding and predicting these processes.
Soil samples were taken from three drip irrigated plots in Korba coastal aquifer in
Tunisia at two periods of the year: summer (after irrigation has ended) and winterseasons. Irrigation water comes from groundwater pumping and presents a high electricalconductivity reaching 9 mS/cm in summer. Experiment involved soil sampling of
unsaturated zone reaching 15 m deep. Soil samples were used to determine vertical soil
salinity profiles on extractions from saturated soil paste. Electrical conductivity andmajor ions were also determined for groundwater samples.
HYDRUS1D was used to carry out preliminary simulations of vertical soil salinity andimpacts on groundwater salinisation. Simulations were first performed for a cycle
corresponding to the irrigation season and second for some years.
Experimental data showed salt accumulation in the surface reaching 16 mS/cm.
HYDRUS1D also simulated surface salinisation during the irrigation season. Impacts on
groundwater quality are rather visible after some years of simulation
Key words Irrigation, HYDRUS1D, salinity, evaporation; precipitation, unsaturated zone

103

104

Simulation of water dynamics in two irrigated soils
Lúcia Barão, Pedro Chambel Leitão, Ramiro J.Neves, Maria Conceição Gonçalves, Tiago
B. Ramos, Nádia Castanheira
Department of Soil Science, Estação Agronómica Nacional, Oeiras, Portugal

ABSTRACT:In this present work three models were used to simulate the water content distribution in
the several soil horizons during seven years (2001-2007): Hydrus, RZWQM (Root Zone
Water Quality Model) and MohidLand. The aim of the work is to show that different
models with different modelling approaches reach the same results and that thecomparisons with the data field also validate the results.
The Hydrus and the MohidLand model used the van Genuchten parameters for thehydraulic curves, while in the RZWQM model the parameters used were the modified
Brooks and Corey. The main difference between Hydrus and MohidLand is the use of
the finite elements instead of finite volumes, respectively. Potential evapotranspiration is
calculated by the RZWQM model through the atmospheric data and was used as an input
for Hydrus and MohidLand. To evaluate the partition between transpiration andevaporation during the simulation, the results of LAI (Leaf Area Index) from RZWQM
previously validated with data fields, were also used.
The simulations were done for two different soil types: Hortic Anthrosol and Eutric
Fluviosol which have different hydraulic properties and textures, in order to observe thedifferent behavior of water distribution in the soil. These soils are located in AlentejoPortugal where irrigation has a very important contribution for the water management inagriculture.

105

The Second HYDRUS Workshop
Edited by Jirka Šimůnek and Radka Kodešová
Czech University of Life Sciences Prague
Faculty of Agrobiology, Food and Natural Resources
Czech Republic
Workshop proceedings
CD-ROM
50
Full papers are also available at the address http://www.pc-progress.cz/
110 pages
First Printing
2008
ISBN: 978-80-213-1783-3
106