Proceedings of Workshop on
HYDRUS Applications

October 19, 2005
Department of Earth Sciences
Utrecht University, The Netherlands

Preface:
Soil and water resources are of great importance to the mankind. As the world
population grows and increasing incidences of pollution of soil and groundwater are
discovered, prudent management and protection of these resources becomesincreasingly important. Most subsurface pollution problems stem from activities
involving the unsaturated (vadose) zone between the soil surface and the groundwatertable. As natural systems are extraordinarily complicated and the understanding andmanagement of these resources cannot be accomplished effectively using heuristicmeans, mathematical modeling has played a critical role in the investigation andmanagement of soil and water resources for the last three decades. The software
packages HYDRUS-1D and 2D are a typical example of commonly used state-of-theart codes involving variably saturated flow and transport.
This proceedings contains the collection of papers presented at the User’s Workshop
on Hydrus that was held in Utrecht, The Netherlands, following a short course on
Hydrus. The purpose of the workshop was to bring together the users and developersof the HYDRUS software packages, to present the latest innovations in the model
applications, and to discuss the capabilities and limitations of HYDRUS. This was aunique opportunity for all HYDRUS users to give a brief overview of their work anddiscuss future needs and directions for obtaining more reliable modelling results. This
collective work includes contributions by many of the users of the HYDRUS software
packages ranging from the very fundamental to the most compelling and important ofapplications.
The editors:
Saeed Torkzaban
S. Majid Hassanizadeh

Proceedings of Workshop on HYDRUS Applications
October 19, 2005
Department of Earth Sciences
Utrecht University, The Netherlands

Address for correspondence:
S. Majid Hassanizadeh
Department of Earth Sciences; Utrecht University
Budapestlaan 4
P.O. Box 80021; 3508 TA UTRECHT; The Netherlands
Tel: +31 30 253 7464; Fax: +31 30 2534900
E-mail: hassanizadeh@geo.uu.nl

Contents
Modeling of Colloid Transport and Deposition in Porous Media
Scott A. Bradford, Martinus Th. van Genuchten, and Jiří Šimůnek…………………..1
Effects of Canopy Shading and Irrigation on Soil Water Content and Temperature
Fares A., J. Šimůnek, L.R. Parsons, Th. M. van Genuchten, and K.T. Morgan……....6Transport of sulfadiazine in soil
Anne Wehrhan, Roy Kasteel, Jirka Simunek, Joost Groeneweg, Harry Vereecken…11
Coupled reactive transport modelling based on the new biogeochemical code HP1
D. Jacques, J. Šimůnek, D. Mallants and M.Th. van Genuchten…………………….15
The multi-component reactive transport module CW2D for subsurface flow
constructed wetlands
Guenter Langergraber, Jirka Simunek ……………………………………………….18
Numerical Analysis of One-Dimensional Coupled Water, Vapor and Heat Transport
in the Vadose Zone using HYDRUS
Hirotaka Saito, Jiri Šimůnek, and Binayak P. Mohanty ……………………………..22
Imprecise simulation of salt dynamic and balance – a Hydrus-1D flaw?
Irina Forkutsa and Rolf Sommer …………………………………………………….26
Application of Hydrus2D to derive a water balance in a karst dominated catchment
Manfred Fink & Marcel Wetzel ……………………………………………………..30
Parameter estimation of soil hydraulic properties from tension infiltrometer data in
PortugalT. B. Ramos, M. Th. van Genuchten, M. C. Gonçalves, J. C. Martins & F. P. Pires .34
Using HYDRUS to simulate water and solute transports in soil lysimeters
M. C. Gonçalves, J. Šimůnek, T. B. Ramos, J. C. Martins, M. J. Neves& F. P. Pires.38
Multi-step outflow experiments to link soil hydraulic properties with electricalcharacteristics
H.-M. Münch, A. Kemna, M. Herbst, E. Zimmermann, H. Vereecken ......................42
PARSWMS: a parallelized version of SWMS_3D
Javaux, M., Hardelauf, H., Gottschalk S., Herbst, M., Vanderborght, J. and H.
Vereecken ……………………………………………………………………………46
Non-equilibrium Isoproturon Transport in Structured Soil Columns: Experiments andModel Analysis
J. Maximilian Köhne, Sigrid Köhne, and Jirka Šimůnek ……………………………49

Numerical Analysis of Soil Limiting Flow from Subsurface Sources
Naftali Lazarovitch, Jiri Šimůnek, Martinus Th. van Genuchten and Uri Shani ……53
Modelling groundwater - surface water interactions in Dutch fens
Paul P. Schot and Stefan C. Dekker …………………………………………………57
Chlorotoluron transport in the soil profile affected by non-equilibrium flow
Radka Kodešová, Martin Kočárek, Jiří Šimůnek, Josef Kozák ……………………..61
Simulating simultaneous nitrification and denitrification with a modified 2D-mobile
immobile model.
Sigrid Köhne1, Jiri Šimůnek, J. Maximilian Köhne and Bernd Lennartz …………..65
Recent Applications of HYDRUS in Irrigated Agriculture
T. H. Skaggs, T. J. Trout, T. J. Kelleners, P.J. Shouse, J.E. Ayars ………………….69
Seasonal Furrow Irrigation Modelling with HYDRUS2
Thomas Wöhling …………………………………………………………………….74
Modelling the effects of soil tillage on water and solute transport: the HYDRUS2D/SISOL coupling
Coquet Y., Šimůnek J., Roger-Estrade J. ……………………………………………78
Inverse model simulation of Br and herbicide transport using multi-process transportmodels with HYDRUS-1D
V. Pot, J. Šimůnek, P. Benoit, Y. Coquet ...................................................................82
New Features of HYDRUS-1D, Version 3.0Jiří Šimůnek, Rien van Genuchten, S. A. Bradford, and R. J. Lenhard .....................86Modeling of Colloid Transport and Deposition in Porous Media
Scott A. Bradford1, Martinus Th. van Genuchten1, and Jiří Šimůnek2
1

George E. Brown Jr. Salinity Laboratory, USDA-ARS, Riverside, CA 92507

2

University of California Riverside, Department of Environmental Sciences,
Riverside, CA 92521

Abstract
A conceptual model for colloid transport is presented that accounts for colloid
attachment, straining, and exclusion. Colloid attachment and detachment are modeled using
conventional first-order rate expressions, whereas straining employs a first-order term that is
depth-dependent. Time dependent colloid deposition as a result of blocking of favorable
attachment sites or filling of straining sites is simulated using a Langmuirian type expression.
The influence of colloid-colloid interactions on straining is also included using at term that is
proportional to the concentration of colloids in suspension and strained on the solid phase.
Finally, size exclusion is modeled by adjusting transport parameters for the colloid-accessible
pore space. Numerical simulations are presented to highlight the various processes
controlling colloid transport and deposition. The sensitivity of simulation results to selectivemodel parameters is also demonstrated.
Introduction
Colloid transport models have typically been derived from a combination of
convection-dispersion-retardation theory for solute transport and colloid attachment theory.
Clean-bed attachment theory employs an irreversible first-order kinetic attachment term andpredicts an exponential decrease in colloid concentration with transport distance. Over the
past decade, a growing body of literature indicates that attachment theory frequently does notprovide an accurate characterization of experimental deposit profiles [Camesano and Logan,
1998; Bolster et al., 1999; Redman et al., 2001; Bradford et al., 2002; Tufenkji et al., 2003; Liet al., 2004; Tufenkji and Elimelech, 2005ab]. This is especially true for larger colloids andfiner textured porous media, and when there are repulsive electrostatic interactions between
colloids and grain surface.
Some of the discrepancies between colloid transport data and attachment theory may
be due to the fact that colloid attachment theory does not account for straining. Straining, isthe trapping of colloid particles in down-gradient pore throats that are too small to allow
particle passage. Straining occurs when colloids are physically deposited in pores that aresmaller than some critical size, although colloid transport may still occur in larger pores. Sizeexclusion is a second related factor that is generally neglected in colloid transport modeling.
Size exclusion influences the transport behavior of colloids by physically limiting theirmobility to the larger pore networks and the more conductive ranges of the pore watervelocity distribution. Hence, colloids may be transported faster than a conservative solute
tracer.
This extended abstract reports on modifications to the HYDRUS 1D code [Simunek etal., 1998; Simunek et al., 2005] to account for colloid attachment, straining, and size
exclusion [Bradford et al., 2003, 2004, 2005]. Selective model simulations are presented to1

illustrate the roles of these processes. Special attention has been given to the roles ofstraining and size exclusion in these examples.
Numerical Modeling
The aqueous colloid mass balance equation is written as:∂ (θ w C )∂J= − T − E sw
(1)
∂t
∂z
where C [N L-3; N denotes number and L length] is the aqueous phase colloid concentration, t[T] denotes time, JT [N L-2 T-1] is the total colloid flux (sum of advective and dispersive
fluxes), and Esw[N L-3 T-1] is the colloid mass transfer terms between the aqueous and solid
phases due to attachment and straining.
A generalized term for Eswatt is given as:∂ (ρ b (S att + S str ))E sw == θ w k attψ att C − ρ b k det S att + θ w k strψ str C ± ρ b k ci S str C(2)∂twith Satt [N M-1] as the solid phase concentration of attached colloids, katt [T-1] is the firstorder colloid attachment coefficient, kdet [T-1] is the first order colloid detachment coefficient,ψatt [-] is a dimensionless colloid attachment function, kstr [T-1] is the straining coefficient, ψstr[-] is a dimensionless straining function, Sstr [N M-1] is the solid phase concentrations ofstrained colloids, and kci [–1 M T-1] is the colloid interaction coefficient.
The first and second terms on the right hand side of (2) accounts for colloid
attachment and detachment, respectively. When ψatt equals 1 and kdet = 0, clean-bed
attachment theory is incorporated into the katt term. Time dependency in colloid attachment
rates as a result of blocking can be accounted for using the ψatt term in (2). During blocking
the value of ψatt decreases with increasing Satt as a result of filling of favorable attachment
sites such as metal oxide surfaces or clay edges. According to the Langmuirian approach, ψatt
= 1- Satt/Sattmax; where Sattmax [Nc M -1] is the maximum solid phase concentration of attachedcolloids. In contrast, the random sequential adsorption (RSA) approach employs a ψatt that is
a nonlinear function of Satt.
Attachment theory has been used to predict colloid transport behavior for a wide
variety of colloid and porous media characteristics. Little attention has been paid to thecolloid and porous media parameter ranges that may limit attachment model applicability,
because straining has been assumed a priori to be insignificant. The third and fourth terms on
the right hand side of (2) accounts for colloid straining and colloid-colloid interaction,
respectively. The influence of straining is modeled using a value of ψstr that is a function ofboth depth and Sstr as:−β⎛
S str ⎞⎟ ⎡ d 50 + z − z o ⎤⎜ψ str = H (z − z o ) 1 −(3)⎥⎜ S max ⎟ ⎢⎣
d
50
⎦
str ⎠
⎝
Here β [-] is a parameter that controls the shape of the colloid spatial distribution, H(z-zo) is
the Heaviside function (0 for z < zo and 1 for z ≥ zo), z [L] is depth, zo [L] is the depth fromthe column inlet or textural interface (layered systems), and Sstrmax [N M-1]is the maximum
solid phase concentration of strained colloids. The second term on the right hand side of (3)accounts for filling and accessibility of straining sites in a manner similar to the Langmuirian
blocking approach.
The fourth term on the right hand side of (2) is used to account for the interaction ofcolloids in solution and those retained by straining. Retention of strained colloids ishypothesized to be hindered (term is negative) at greater collision frequencies (higher
concentrations) when there are repulsive colloid (aqueous phase) - colloid (strained)2interactions. Conversely, attractive interactions of colloids in solution and those retained on
the solid phase produce greater straining and can be simulated using a positive value for the
fourth term of (2). Colloid-colloid interactions (repulsive or attractive) will dependent on thecolloid surface chemistry and the aqueous chemistry. Our approach assumes that thefrequency of collisions is proportional to the colloid concentration in solution and the amountretained on the solid phase.
To account for size and anion exclusion, values of θw and the Darcy water velocity, qw[L T-1] are adjusted to account for accessibility by colloids according to:θ mw = θ w - εγ

(4)q mw = q wk rmw / k rw

(5)

Here ε [-] is the porosity, γ [-] is the water saturation that is not accessible to mobile colloids,
krw [-] is the water relative permeability, and krmw [-] is the colloid accessible water relative
permeability. The ratio of qmw and θmw yields the colloid accessible pore-water velocity. The
value of krw is typically determined from capillary pressure data using a pore size distribution
model. The value of krmw is determined using the model of Burdine [1953] by adjusting its
limits of integration to reflect pore space accessible to colloids:
k rmw

Sw

2

1

2

∫ R( x ) dx
(S )= S γ
∫ R( x ) dx
w

2
w

Sw>γ

(6)

0

Here Sw [-] is the water saturation, R [L] is the pore size, and x [-] is a dummy saturation
variable of integration. When Sw is less than γ then the value of krmw is set equal to zero. The
value of γ can be theoretically related to a critical pore radius for exclusion using the Laplace
equation and a capillary pressure curve. The critical pore radius for exclusion is likely afunction of many variables, including: colloid size and concentration, diffuse double layer
thickness, surface roughness, and the presence of retained colloids on the grain and poresurfaces. Alternatively, the value of γ can be fitted to transport data and a critical pore radius
for exclusion can be calculated according to the outlined procedure.
Example Simulations
This section presents example simulations and experimental data that illustrate theoutlined model capabilities and processes. Figure 1 presents o served and simulated
breakthrough curves (Figure 1a) and deposition profiles (Figure 1b) for 1 µm carboxyl latex
microspheres in 150 µm Ottawa (quartz) sand. Simulations considered attachment anddetachment (ψatt=1, kstr = 0; and kci = 0), attachment and Langmuirian blocking (kdet = 0; kstr =0; and kci = 0), straining and Langmuirian filling (katt = 0; kdet = 0; and kci = 0), and straining
and colloid interactions (katt = 0; kdet = 0; and Sstr/Sstrmax = 1). Notice that the straining
simulations provide a reasonable description of both transport and deposition data, whereas
attachment simulations only accurately charactize the breakthrough curves. Bradford et al.
[2003] present more physically realistic simulations that consider both attachment andstraining.

3

Figure 1.
1.0
0.0

Data
katt/kdet

0.6

C/C i

katt/smax
kstr/smax

0.4

kstr/kci

0.2

Dimensionless Depth

0.8

0.2

Data
katt/kdet
katt/Smax
kstr/Smax
kstr/kci

0.4
0.6
0.8
1.0

0.0
0

50

0

100 150 200 250 300 350

2

4

6

8

10

-1

N N i / g dry sandTime (min)Figure 2 presents simulated breakthrough curves (Figure 2a) and deposition profiles(Figure 2b) when the straining parameters are fixed and the colloid interaction coefficient isvaried. Observe that transport and deposition behavior is highly sensitive to this parameter,suggesting that colloid-colloid interactions will have an important role on straining behavior.
Figure 2.
1.0
0.0

kci=0
kci=0.03(-)
kci=0.06(-)
kci=0.03(+)
kci=0.06(+)

0.6
0.4

Dimensionless Depth

C/C i

0.8

0.2
0.0
0

0.2
kci=0.0

0.4

kci=0.03(-)

0.6

kci=0.06(-)
kci=0.03(+)
kci=0.06(+)

0.8
1.0

50 100 150 200 250 300 350

0

5

10

15

20

25

30

-1

Time (min)N N i / g dry sandFigure 3 presents the predicted velocity enhancement (the ratio of the pore watervelocity for mobile colloids compared to a conservative solute) for mobile colloids in ahypothetical soil for various values of Sw and γ. Observe that for Sw>γ the velocity
enhancement tends to increase with increasing values of γ and decreasing values of Sw. For
Sw<γ no colloids are mobile.
Figure 3.
54v cw /v w

Sw=0.2
Sw=0.4

3

Sw=0.6
Sw=0.8Sw=1.0

2

1
0

0.2

0.4

0.6

0.8

1

γ

4

References
Bolster, C. H., A. L. Mills, G. M. Hornberger, and J. S. Herman. 1999. Spatial distribution ofbacteria experiments in intact cores. Water Resour. Res. 35:1797-1807.
Bradford, S. A., S. R. Yates, M. Bettahar, and J. Simunek. 2002. Physical factors affecting thetransport and fate of colloids in saturated porous media. Water Resour. Res. 38:1327,
doi:10.1029/2002WR001340.
Bradford, S. A., J. Simunek, M. Bettahar, M. Th. van Genuchten, and S. R. Yates. 2003.
Modeling colloid attachment, straining, and exclusion in saturated porous media.
Environ. Sci. Technol. 37:2242-2250.
Bradford, S. A., M. Bettahar, J. Simunek, and M. Th. Van Genuchten. 2004. Straining andattachment of colloids in physically heterogeneous porous media. Vadose Zone J.
3:384-394.
Bradford, S. A., J. Simunek, M. Bettahar, Y. F. Tadassa, M. Th. Van Genuchten, and S. R.
Yates. 2005. Straining of colloids at textural interfaces. Water Resour. Res., In press.
Burdine, N. T. 1953. Relative permeability calculations from pore size distribution data.
Trans. Am. Inst. Min. Metall. Pet. Eng. 198:71-78.
Camesano, T. A., and B. E. Logan. 1998. Influence of fluid velocity and cell concentrationon the transport of motile and nonmotile bacteria in porous media. Environ. Sci.
Technol. 32 (11), 1699-1708.
Li, X., T. D. Scheibe, and W. P. Johnson. 2004. Apparent decreases in colloid deposition
rate coefficient with distance of transport under unfavorable deposition conditions: ageneral phenomenon. Environ. Sci. Technol. 38 (21), 5616-5625.
Redman, J. A., S. B. Grant, T. M. Olson, and M. K. Estes. 2001. Pathogen filtration,
heterogeneity, and the potable reuse of wastewater. Environ. Sci. Technol. 35 (9):1798-1805.
Šimůnek, J., K. Huang, M. Šejna, and M. Th. van Genuchten, The HYDRUS-1D software
package for simulating the one-dimensional movement of water, heat, and multiple solutes
in variably-saturated media. Version 1.0, IGWMC - TPS - 70, International Ground Water
Modeling Center, Colorado School of Mines, Golden, Colorado, 186pp., 1998.
Šimůnek, J., M. Th. van Genuchten, and M. Šejna, The HYDRUS-1D software package forsimulating the one-dimensional movement of water, heat, and multiple solutes in variablysaturated media. Version 3.0, HYDRUS Software Series 1, Department of Environmental
Sciences, University of California Riverside, Riverside, CA, 270pp., 2005.
Tufenkji, N., J. A. Redman, and M. Elimelech. 2003. Interpreting deposition patterns of
microbial particles in laboratory-scale column experiments. Environ. Sci. Technol.
37(3):616-623.
Tufenkji, N. and M. Elimelech. 2005a. Spatial distributions of Crypotsporidium oocysts inporous media: Evidence for dual model deposition. Environ. Sci. Technol., 39:36203629.
Tufenkji, N., and M. Elimelech. 2005b. Breakdown of colloid filtration theory: Role of thesecondary energy minimum and surface charge heterogeneities. Langmuir, 21:841852.

5Effects of Canopy Shading and Irrigation on Soil Water Content andTemperature

Fares1 A., J. Šimůnek2, L.R. Parsons3, Th. M. van Genuchten4, and K.T. Morgan3
1
University of Hawaii-Manoa, Email AFares@Hawaii.Edu
2
University of California Riverside, 3University of Florida, 4USSL-USDA-ARS, Riverside, CA.
Abstract
Citrus root systems are exposed to different thermal and hydrologic conditions as a result of tree
canopy shading and under-tree micro-irrigation. Because micro-sprinklers wet only part of the soil
surface and are located under the tree, roots under the canopy usually receive more water than thoseoutside the tree canopy. The combined effects of different soil temperature and water input onwater redistribution under field conditions have not been fully studied in Florida sandy soils. The
objective of this study was to investigate shading and irrigation effects on spatial distribution ofwater content and soil temperature at different soil depths. Real-time capacitance soil water sensors
and thermocouples were used to monitor soil water content and temperature at depths of 0, 10, 20,
40, 80 and 150 cm. Weather parameters were monitored simultaneously at the same location.
HYDRUS-2D, a two dimensional computer package for simulating movement of water, heat, andmultiple solutes in variably saturated media, was used to simulate water flow and heat transport
under such conditions. The predicted water contents and soil temperatures compared favorablywith their corresponding observed parameters. In addition to its accuracy in simulating this system,HYDRUS-2D helped to improve the analysis of this research project.
Introduction
Citrus-tree root system is exposed to different thermal and hydrologic conditions as a resultof tree canopy and under-tree microirrigation system. Thus, roots under the tree canopy are shadedby the canopy and receive irrigation water. In contrast, roots outside tree canopy are uncovered,
receive full solar radiation, and only receive water from rain. Assuming isothermal conditions in
studying water movement and root uptake could result in large errors.
Temperature plays an important role in the chemical, physical and biological functioning ofa soil. Temperature gradients can cause water to move in the vapor phase from a hot site to a coldsite. Water moves as a result of combined effects of temperature, water potential and soluteconcentration gradients. The movement of water under combined temperature and matric potential
head gradient has been studied under wide ranges of soil and weather conditions over last four
decades and have been reported (Philip and De Vries, 1957; Horton and Chung, 1991).
Numerical models, with varying complexities were developed to describe coupled heat,water and solute transfer (Nassar and Horton, 1992; Bear and Gilman, 1995 and Šimůnek et al.,1999). The combined effects of different soil temperature and water input under Florida sandy soils
on water redistribution under field conditions has not been fully studied.
The objectives of this study were: i) to investigate, in the field, shading and irrigation effects
on spatial distribution of water content and soil temperature at different soil depths (0, 10, 20, 40, 80
and 150 cm) and six field locations, ii) to use HYDRUS to simulated coupled water and heattransport.
Materials and Methods
The field study was conducted on a citrus grove of Hamlin orange trees (Citrus sinensis (L.) Osb.)
on Swingle citrumelo (Citrus paradisi Macf. X Poncirus trifoliata (L.) Raf.) rootstock at the Water
Conserv II, Avalon, Orange County, Florida. The soil is a Candler fine sand (hyperthermic,
uncoated, Typic Quartzipsamments). Two replicates of three main locations within the tree canopy,
at the dripline, and in the open between two tree rows were selected for measuring soil water6

content and temperature. At each location, soil water contents were monitored at soil surface,and10, 20, 40, and 80 cm below soil surface at 30 min intervals. Air temperature, surface, and soil
temperature at previously mentioned depths were monitored at the same time interval (30 minutes).
Simulations
HYDRUS-2D software package was used to simulate the movement of water and heat transport. A
1-D soil transect, 0.8 m deep, was used in this study. The Candler soil profile consists of very deep,
excessively drained, rapidly permeable soils formed on thick beds of eolian or marine deposits ofcoarse-textured materials.
Soil water release curves and their corresponding hydraulicconductivities are presented in Fares et al. (2000).
Results and Discussion
Figure 1 shows the rainfall data (A) that was measured on the research site and that was used as
input for the simulations. There were several rainfall events during the first 10 days of thesimulation period; however less rainfall was registered in the second part of this period. As a result,the soil water content showed at all depths fluctuated more frequently in the first part of thesimulation as compared to the second part. The simulated and measured water contents at 10, 20
and 40 cm are shown in Fig. 1A, Fig. 1B and Fig. 1C, respectively. There is an excellent agreement
between the measured simulated data at the 20 and 40 cm.

7

12

A

10cm
Rain
10cm

B

20cm
20cm

10
8

Water Content at Different Depths (%)

6
4
2
0
10

8

6

4
11

C

40cm
40cm

10
9
8
7
6

250

260

270

280

290

300

Day of the Year (2001)
Fig. 1. Daily rainfall rates (solid bars, A), and measured (open circles) and simulated (lines) soil watercontents at 10 cm (A), 20 cm (B), and 40 cm (C) depths.

8

Simulated and measured data of the soil temperatures at 40 and 80 cm depth are shown in Fig. 2A
and Fig. 2B, respectively. Daily rainfall and daily net radiation measured outside the tree canopy inthe open are shown in Fig. 2 A and 2B, respectively. The first 10 days of the simulation showed
substantial decrease in the soil temperature at the two shown depths. This period was a rainy periodwith low net radiation values; thus, the lower the net radiation the sharper was the decrease in thesoil temperatures. During this time of the year, rainfall in Florida occurs as a result of cold front
spills that continue for days. During these cold fronts, there were almost overcast conditions that
result in a substantial decrease of the net radiation.

Sim-40cm
Obs-40cm
Rain (open)A

Temperature OC

30

10

28

8

26

6

24

4

22

2

20
28

0

B

Obs-80cm
Sim-80cm
Net Radiation (open)

27
Temperature OC

12

Rainfall (mm)Simulated and Observed Soil Temperatures

300

26
200

25
24

100

23
22

0

250

255

260

265

270

275

Day of the Year (2001)
Fig. 2. Daily rainfall (A) and net radiation (B) rates (bar graphs), and measured (lines) and simulated(dots) soil temperatures at 40 cm (A) and 80 cm (B) depths.

9

Solar Radiation (W m2)

32

The lowest soil temperatures were during the rainy and low net radiation days (days 256-258). This
wet and overcast period was followed by a week of dry and clear sky conditions; during this periodsoil temperatures at both depths 40 and 80 cm increased substantially to over 26 and 25 0C,
respectively. This period was followed by a cooling period with 11 mm rainfall event andsubstantial decrease in net radiation.
Despite the periodical decreases and increases in the soil temperature, there was an overall soil
temperature decrease over this simulation period at both soil depths (Fig. 2 A and B). As expected,
the amplitude of soil temperature decreased with depth as a result of the effect of damping depth.
The predicted soil temperatures compared favorably with their corresponding observed data.
Incoming solar radiation is the driving force of soil temperature fluctuations. Differences between
simulated and measured soil temperatures as shown in (Fig. 2 A&B) are more pronounced during
clear sky and warmer days. However, during relatively cooler days (overcast), differences between
measured and simulated temperatures were the lowest.
References
1.
2.
3.
4.
5.
6.

7.
8.

Bear, J., and A. Gilman. 1995. Migration of salts in unsaturated zone caused by heating. Transp. Porous Media
19:139-156.
Fares, A., A.K. Alva., P. Nkedi-Kizza and M. A. Elrashidi. 2000. Determination of soil water physicalproperties under field conditions using capacitance probe and Guelph permeameter. Soil Sci. 165: 768-777.
Horton, R., and S.O. Chung. 1991. Soil heat flow. P. 397-438. In J. Hanks and J. Ritchie (ed.) Modeling plantand soil systems. Agron. Monogr. 31. ASA, CSSA, and SSSA, Madison, WI.
Nassar, I.N., and R. Horton. 1992. Simultaneous transfer of heat, water and solute in porous media: I.
Theoretical development. Soil Sci. Soc. Am. J. 56:1350-1356.
Philip, J.R., and D.A. de Vries. 1957. Moisture movement in porous materials under temperature gradients.
Trans. Am. Geophys. Union 38:222-232.
Šimůnek, J., and M. Th. van Genuchten. 1999. Using the HYDRUS-1D and HYDRUS-2D codes for
estimating unsaturated soil hydraulic and solute transport parameters, in van Genuchten, M. Th., F. J. Leij, andL. Wu (eds.) Characterization and Measurement of the Hydraulic Properties of Unsaturated Porous Media,
University of California, Riverside, CA, 1523-1536.
Šimůnek, J., M. Šejna, and M. Th. van Genuchten. 1999. The HYDRUS-2D software package for simulating
to-dimensional movement of water, heat, and multiple solutes in variable saturated media. Version 2.0,
IGWMC-TPS-53, International Ground Water Modeling Center, Colorado School of Mines, Golden, CO.
Vrugt, J. A., J. W. Hopmans, and J. Šimůnek. 2001. Calibration of a two-dimensional root water uptake model
for a sprinkler-irrigated almond tree. Soil Sci. Soc. Am. J., 65(4):1027-1037.

10

Transport of sulfadiazine in soil
Anne Wehrhana,c , Roy Kasteela , Jirka Simunekb , Joost Groenewega , Harry Vereeckenaa Agrosphere Institute, ICG IV, Forschungszentrum Jülich GmbH, Germany;b Department of Environmental Sciences, University of California, Riverside, US;c e-mail: a.wehrhan@fz-juelich.deIntroduction
Among other veterinary pharmaceutical Sulfadiazine (SDZ) is a widely used antimicrobial substancein intensive livestock production to prevent and treat diseases (Boxall et al., 2004; Thiele-Bruhn,
2003). Up to 40 % of the administered sulfonamides are eliminated as microbial active substanceswith the excretions (Kroker , 1983). Antibiotics such as sulfadiazine reach agricultural soils directly
through grazing livestock or indirectly through the spreading of manure or sewage sludge on the field(Jørgensen and Halling-Sørensen, 2000). Knowledge about the fate of antibiotics in soil is crucial
to assess the environmental risk of those compounds, including possible transport to groundwater.
Despite the assumed high mobility of SDZ in soil from its physicochemical properties (Tolls,
2001), SDZ was scarcely found in the leachate of column and plot studies and most of the applied
SDZ was retained in the upper part of the soil (Kreuzig and Höltge, 2005). However, the fate of
SDZ is not yet thoroughly investigated, partly because of analytical difficulties at the relevant low
environmental concentration levels.
The objective of this study was to investigate the transport behavior of sulfadiazine in disturbed
soil columns at a constant flow rate near saturation. Sulfadiazine was applied at different concentrations and for different pulse lengths. To identify relevant sorption processes, measured breakthroughcurves and soil concentration profiles of sulfadiazine were fitted with a convective-dispersive transport model considering different sorption concepts. One, two and three sites sorption models were
tested and their assumptions and limitations are discussed.

Materials and Methods
SDZ (Figure 1) belongs to the chemical group of sulfonamides and has a water solubility of77 mg L−1 , a KOW of 0.76 and two pKa values at 1.6 and 6.5. The soil material (silty loam)
was a mixed, sieved sample (2 mm) of the upper 30 cm of an Eutric Cambisol used as grassland
near lake Greifensee, Switzerland, having an organic matter content of 3.3 % (mass) and a pH valueof 6.1.

Figure 1: Chemical structure of SDZ. The 14 C-labelling is in the phenyl-ring.

11
1

The transport of 14 C labelled SDZ was investigated in disturbed soil columns (8.5 cm in diameter, 9 cm in height) at a constant flow rate near saturation. SDZ was applied in different
concentrations (5.70 vs. 0.57 mg L−1 ) for either a short or a long pulse duration (7 vs. 68 h).
Breakthrough curves (BTC) of SDZ and the non-reactive tracer Cl− where measured during 500 h.
At the end of each experiment the concentration profile in the soil column was determined. Trace
analysis of 14 C-SDZ in the water samples was done directly by liquid scintillation counting (LSC).
For the soil samples total combustion was followed by LSC of the evolved and trapped 14 CO2 .
The BTCs were fitted with a convective-dispersive transport model considering different sorption concepts to identify relevant sorption processes. These concepts were chosen according to
the previously observed sorption behavior: In batch experiments sorption of SDZ is non-linear(Freundlich) and rate-limited. An apparent sorption equilibrium was reached only after 20 days,
the desorption seems to be much slower. Time-dependent formation of non-extractable residues isfrequently reported (Kreuzig et al., 2003). Thus, the model needs to account for the concentrationdependent, rate-limited and possibly irreversible sorption processes.
An adapted version of HYDRUS-1D (Simunek et al., 1998) was used to describe the transport experiments with SDZ. HYDRUS-1D is a finite element code which numerically solves theconvection-dispersion equation. The possible solute-soil-water interaction models include one, twoor three sorption sites, showing instantaneous sorption (linear or non-linear) or rate-limited reversible or irreversible sorption (Figure 2). Isotherm-based concepts were compared to the kinetic
attachment/detachment concept, where all interactions are described by first-order rate processes.
Attachment/detachment
model

Isotherm-based models
Linear sorption
S2

α2- C Kd- S1
β3

Freundlich sorption
S2K
α2- C  f- S1
m

S2

β1β2
- C 
S1
γ2

γ1

β3

β3

?

?

?

S3

S3

S3

Figure 2: Applied solute-soil-water distribution concepts. The boxes labelled with C represent theliquid phase with the concentration C [M L−3 ], the boxes S1 , S2 and S3 represent the three sorption
sites with the respective concentrations S [M M−1 ]. The arrows indicate the sorption process, where
Kd is the distribution coefficient, Kf and m are the Freundlich coefficient and exponent, α is thereversible ad- and desorption rate [T−1 ], β and γ are the one-way attachment and detachment
rates [T−1 ], respectively. Less complex versions of each model were derived by omitting one or two
sorption sites.

Results
The BTCs are characterized by an early arrival of SDZ and an extended tailing (Figure 3(a)). The
peak maxima of the different treatments were delayed relative to chloride by a factor of 2 to 5 anddiffered in the maximum concentration as well as the eluted mass fraction (18 % to 83 %) within
the leaching period of experiments. Due to the 14 C-analysis the error in mass balances was as small2

as 1 to 3 %. The way of application determined the fate of SDZ: (i) the longer the pulse durationand the higher the applied concentration, the more mass was leached, and (ii) the more mass wasleached the later the peak concentration arrived.
The concentration profiles of SDZ for the columns with the long pulse applications show thehighest concentrations at the top of the column, steadily decreasing towards the bottom (Figure3(b)). In the column with the short pulse application the highest concentrations were found between
Soil concentration profiles for column experiments
2 and 6 cm
ofA,thecolumn.
BTCdepthof column
B, C and
3S2R and 3S5R fitted curves
experimental and modeled data

05.70 mg L-1, 68 h
3 site kinetic model
0.57 mg L-1, 68 h
3 site kinetic model
5.70 mg L-1, 7 h
3 site kinetic model

10-1

-2

Depth [cm]Relative concentration C C0-1

100

10-2

-4
5.70 mg L-1, 68 h
3 site kinetic model
0.57 mg L-1, 68 h

-6

3 site kinetic model
5.70 mg L-1, 7 h

10-3
-8

3 site kinetic model
10-4

-10
0

20

40

60

80

100

120

140

0

1

2

3

4

5

-1

Cummulative outflow [cm]Sorbed concentration [mg kg ]

(a)

(b)

Figure 3: Breakthrough curves (a) and soil concentration profiles of the residual 14 C concentrations(b) in the three columns. Measured data are represented by the symbols, modelled curves by thesolid and dashed lines. Input concentration and pulse duration of the application scenarios are
given in the legend. Note the logarithmic scale of the concentration axis for the BTCs.
The BTC of chloride were used to determine the water flow parameters volumetric water content and dispersivity. They were kept constant during the modelling of the reactive tracer. Simple
equilibrium models overestimate the mobility of SDZ in soil, especially for lower application concentrations or short pulse applications. The extended tailing cannot be modelled with an one
site equilibrium concept. Kinetic sorption processes allow both, an early breakthrough and a pronounced tailing. Irreversible sorption processes can account for the large mass fractions found in
the soil column at the end of the experiment. Only the three site kinetic model with two reversible
attachment/detachment sites and one irreversible site fitted the BTCs of all three experiments(Figure 3(a)). However, none of the models was able to describe the observed profiles (Figure3(b)).

Conclusion
The major difficulties to fit the BTCs were the tailing and the mass recovery in the leachate.
Only the complex three site kinetic sorption model with reversible and irreversible sorption sites is
flexible enough to describe the observed BTCs under the variable application scenarios. Despite thegood agreements for the BTCs, the observed and calculated concentration profiles in soil differed
substantially. Therefore, common approaches for process identification on the basis of the main peak
breakthrough without the observation of the tailing and the concentration profile are precarious.
3

As none of the models satisfactorily described the observed transport behavior of SDZ, thesoil-water-SDZ system seems to be more complex than represented by the tested models. Dueto the restricted analysis of 14 C only in this study, transformation of the applied SDZ cannot be
ruled out. Transformation of SDZ in soil was observed by Kreuzig and Höltge (2005). In manure,
one of the identified transformation products is even able to transform back into the parent compound (Kroker , 1983). In case transformation occurred in our experiments, the radio-chemical
and mathematical analysis presented here described the effective behavior of the parent and possible transformation products. Because the active, re-transferable or unidentified transformation
products can put the environment at an additional risk, such an effective description is still ofinterest.

References
Boxall, A. B. A., L. A. Fogg, P. A. Blackwell, P. Kay, E. J. Pemberton, and A. Croxford, Veterinary
medicines in the environment, Reviews of Environmental Contamination and Toxicology, 180 ,
1–91, 2004.
Jørgensen, S. E., and B. Halling-Sørensen, Drugs in the environment, Chemosphere, 40 , 691–699,
2000, Editorial.
Kreuzig, R., and S. Höltge, Investigations on the fate of sulfadiazine in manured soil: Laboratory
experiments and test plot studies, Environmental Toxicology and Chemistry, 24 (4), 771–776,
2005.
Kreuzig, R., C. Kullmer, B. Matthies, S. Höltge, and H. Diekmann, Fate and behavior of pharmaceutical residues in soils, Fresenius Environmental Bulletin, 12 (6), 550–558, 2003.
Kroker, R., Aspects of the elimination of antibiotics after treatment of domestic animals, Wissenschaft und Umwelt, 4 , 305–308, 1983.
Simunek, J., M. Senja, and M. T. van Genuchten, The HYDRUS-1D software package for simulating
the one-dimensional movement of water, heat, and multiple solutes in variably saturated media.
Version 2.0 , ICGWMC-TPS-70. International Ground Water Modeling Center, Colorado Schoolof Mines, Golden, Colorado, 1998.
Thiele-Bruhn, S., Pharmaceutical antibiotic compounds in soil - A review, Journal of Plant Nutrition and Soil Science, 166 , 145–167, 2003.
Tolls, J., Sorption of veterinary pharmaceuticals in soils: A review, Environmental Science andTechnology, 35 (17), 3397–3406, 2001.

14
4

Coupled reactive transport modelling based on the new biogeochemical code
HP1
D. Jacques(1)*, J. Šimůnek(2), D. Mallants(1) and M.Th. van Genuchten(3)(1)

Waste and Disposal Department, SCK•CEN, Boeretang 200, B-2400 Mol, Belgium
Department of Environmental Sciences, A135 Bourns Hall, University of California
Riverside, 900 University Avenue, Riverside, CA USA 92521(3)
George E. Brown, Jr. Salinity Laboratory, USDA, ARS, 450 W. Big Spring Rd,
Riverside, CA USA 92507* Corresponding author (djacques@sckcen.be)(2)

The migration of many naturally occurring elements and contaminants in thesubsurface is affected by a multitude of complex, interactive physical, chemical,
mineralogical, geological, and biological processes. Cycles of precipitation andevaporation largely determine if contaminants remain near the soil surface. Changes
in the chemical composition or pH of the soil solution may impact the retention ofheavy metals on organic matter or iron oxides. Dissolution and precipitation ofminerals generally buffer the transport of a solution with a different pH through thesoil profile. Simulation of these and related processes requires a coupled reactive
transport code that integrates the physical processes of water flow and advectivedispersive solute transport with a range of biogeochemical processes.
Recently, a new comprehensive simulation tool HP1 (HYDRUS1D-PHREEQC,
Jacques and Šimůnek, 2005) was developed by coupling the HYDRUS-1D onedimensional variably-saturated water flow and solute transport model (Šimůnek et al.,1998) with the PHREEQC geochemical code (Parkhurst and Appelo, 1999). The HP1
code incorporates modules simulating (1) transient water flow in variably-saturated
media, (2) transport of multiple components, and (3) mixed equilibrium/kinetic
geochemical reactions. The program numerically solves the Richards equation forvariably-saturated water flow and advection-dispersion type equations for heat andsolute transport. The flow equation incorporates a sink term to account for water uptake
by plant roots. The heat transport equation considers transport due to conduction andconvection with flowing water. The solute transport equations consider advectivedispersive transport in the liquid phase. The program can simulate a broad range of lowtemperature biogeochemical reactions in water, soil and ground water systems
including interactions with minerals, gases, exchangers, and sorption surfaces, basedon thermodynamic equilibrium, kinetics, or mixed equilibrium-kinetic reactions. The
program may be used to analyze water and solute movement in unsaturated, partiallysaturated, or fully saturated porous media. The flow region may be composed ofnonuniform soils or sediments. Flow and transport can occur in the vertical, horizontal,
or a generally inclined direction. The water flow part of the model can deal withprescribed head and flux boundaries, boundaries controlled by atmospheric conditions,as well as free drainage boundary conditions. The governing flow and transport
equations were solved numerically using Galerkin-type linear finite element schemes.
The coupled transport and reaction equations are solved by the non-iterative sequential
approach. Its numerical accuracy was evaluated for a number of simulations undertransient flow conditions (Jacques et al., 2005a). Graphical interfaces for both original15

codes can be used with HP1 as well, enabling an easy generation of the input data anddisplay of computational results.
HP1 is specifically developed for simulating water flow, transport and geochemical
reactions in environmental soil quality problems. HP1 has been tested for a number ofapplications with other (less general) models (reported in Jacques and Šimůnek,2005). Several studies (Jacques and Šimůnek, 2005, Jacques et al., 2003, Šimůnek etal., 2005) dealing with typical coupled transport phenomena demonstrated HP1’scapabilities for environmental problems such as transport of heavy metals (Zn2+, Pb2+,
and Cd2+) subject to multiple cation exchange reactions, cadmium leaching in acidsandy soils, and transport of organic contaminants such as TNT or PCE (?)and its
degradation products through soil.
We present here two examples: (1) Long term transient flow and transport of major
cations and heavy metals in a soil profile (based on Jacques et al., 2005a), and (2)long term leaching of uranium from agricultural field soils following mineral fertilizer
application (based on Jacques et al., 2005b,c) .
In the first example, transport of major cations and heavy metals in a soil profile
under atmospheric boundary conditions is simulated. Two types of geochemical
reactions are taken into account: (i) equilibrium aqueous speciation reactions (e.g.
complexation reactions between Cd and Cl, (ii) and cation exchange reactions of boththe major cations and the heavy metals. As will be illustrated, there is a strong
interaction between water flow, water content, geochemical conditions (e.g. pH),
speciation and mobility. For example, aqueous concentrations of the heavy metals
increase significantly during summer near the soil surface. In other words, transient
flow conditions also influence the bioavailability of heavy metals, since uptake
processes by plants and soil micro-organisms are often concentration dependent (i.e.,
passive uptake of solutes together with uptake of water by roots, as well as active
uptake which can be simulated with HP1 using Monod or Michealis-Menten kinetics).
Moreover, the high heavy metal concentrations occur during periods with the highest(micro)biological activity during the year (the summer months). HP1 hence is thusa
very useful tool for evaluating the effects of atmospheric boundary conditions on thelong-term bioavailability and leaching of heavy metals in and from soils. While notfurther addressed here, such evaluations possibly could also account for the transportof colloids such as micro-organisms, humics or fulvics.
The second examples considers the transport of uranium (U) present in many mineral
fertilizers, particularly (super)phosphates. Field soils that receive P-fertilizers
accumulate U and thorium (Th) and their daughter nuclides, which eventually may
leach to groundwater. Our objective was to numerically assess U migration in soils.
Calculations were based on the new reactive transport model, HP1. In our simplified
model, we account for interactions between U and organic matter, phosphate, andcarbonate. Solid phase interactions were simulated by the combination of a cationexchange complex module (adsorption of major cations and U) and a surface
complexation module (adsorption of cations and P). Furthermore, all geochemical
processes were coupled with a model accounting for dynamic changes in the soil
water content and the water flux. The capabilities of the code in calculating natural U
fluxes to groundwater were illustrated using a semi-synthetic 200-year-long timeseries of climatological data for Belgium. Based on an average fertilizer application,

16

the input of phosphate and uranium in the soil was defined. In this example, the effectof the upper boundary condition (transient versus steady-state flow) on the migrationand leaching of U was specifically addressed. Results showed that due to changing
water content, water flow and, consequently, geochemical conditions (mainly pHchanges), migration of U is faster under atmospheric boundary conditions than understeady-state flow.
HP1 is free software. Information on how to obtain HP1 can be found on
www.sckcen.be\hp1 or by contacting one of the authors of this paper. A list of key
publications is given in the references.
References
Jacques, D., Šimůnek, J., 2005. User manual of the multicomponent variablysaturated flow and transport model HP1: Description, Verification, and Examples(Version 1.0). SCK•CEN, Mol, Belgium, BLG-998, 79 p.
Jacques, D., J. Simunek, D. Mallants, and M.Th. van Genuchten, 2003. The
HYDRUS-PHREEQC multicomponent transport model for variably-saturated
porous media: Code verification and application. In ‘Proceedings of MODFLOW
and MORE 2003: Understanding through modeling’, (Eds. E. Poeter, C. Zheng, M.
Hill, and J. Doherty) MODFLOW and MORE 2003, September 16-19, 2003,
International ground water modeling center, Colorado, Volume 1: 23-27
Jacques, D., Šimůnek, J., D. Mallants, and M.Th. van Genuchten, 2005a. HP1: A
coupled numerical code for variably saturated water flow, solute transport andbiogeochemistry in soil systems. Submitted to J. of Cont. Hydrol.
Jacques, D., Šimůnek, J., D. Mallants, and M.Th. van Genuchten, 2005b. Long term
uranium migration in agricultural field soils following mineral P-fertilization.
Proceedings of ICEM ’05: The 10th International Conference on Radioactive Waste
Management and Environmental Remediation, September 4 – 9, Glasgow,
Schotland.
Jacques, D., Šimůnek, J., D. Mallants, and M.Th. van Genuchten, 2005c. Modellinguranium leaching from agricultural soils to groundwater as a criterion forcomparison with complementary safety indicators. Proceedings of MRS05,
Material Research Society 2005, Symposium on the Scientific Basis for Nuclear
Waste Management, September 12-15, Ghent, Belgium.
Parkhurst, D.L., and C.A.J. Appelo, 1999. User's guide to PHREEQC (Version 2) – Acomputer program for speciation, batch-reaction, one-dimensional transport, andinverse geochemical calculations. Water-Resources Investigations, Report 99-4259,
Denver, Co, USA, 312 pp.
Šimůnek, J., D. Jacques, M.Th. van Genuchten, and D. Mallants, 2005.
Multicomponent geochemical transport modeling using the HYDRUS-1D computer
software package. 2005 Summer Speciality conference, Institutions for sustainable
watershed management: Reconciling physical and management ecology in the
Asia-Pacific, American Water Resources Association, Hyatt Regency Waikiki
Resort & Spa, Honolulu, HI, June 27-29 2005.
Šimůnek, J., M. Šejna, and M.Th. van Genuchten, 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, IGWMC-TPS-70, International
Ground Water Modeling Center, Colorado School of Mines, Golden, CO, 202 p.

17

The multi-component reactive transport module CW2D for subsurface flow
constructed wetlands
Guenter Langergraber1, Jirka Simunek2
1

Institute of Sanitary Engineering and Water Pollution Control, BOKU - University of Natural
Resources and Applied Life Sciences, Vienna, Muthgasse 18, A-1190 Vienna, Austria.
Guenter.Langergraber@boku.ac.at2
Department of Environmental Sciences, University of California Riverside, Riverside, CA
92521, USA. Jiri.Simunek@ucr.edu
Abstract
Constructed wetlands (CWs) are becoming increasingly popular worldwide for removing organic
matter, nutrients, trace elements, pathogens, or other pollutants from wastewater and/or runoff
water. CWs represent a very complex mixture of water, substrate, plants, litter, and a variety of
micro-organisms, in which a large number of simultaneous processes takes place. To describe
this complex system of mutually dependent physical, chemical and biological reactions andchemical and biological compounds is almost impossible without the help of a numerical modelthat can consider the majority of these interactions. Water flow regime in subsurface vertical flow
constructed wetlands is also highly dynamic, adding to the complexity of the overall system byrequiring a transient variably-saturated flow model.
The multi-component reactive transport module CW2D, which is an extension of the variablysaturated water flow and solute transport program HYDRUS-2D (Simunek et al., 1999), includes
biochemical transformation and degradation processes in subsurface flow constructed wetlands.
The biochemical components defined in CW2D include dissolved oxygen, three fractions oforganic matter (readily- and slowly-biodegradable, and inert), four nitrogen compounds(ammonium, nitrite, nitrate, and dinitrogen), inorganic phosphorus, and heterotrophic andautotrophic micro-organisms. Organic nitrogen and organic phosphorus were modelled as part ofthe organic matter. The biochemical degradation and transformation processes were based on
Monod-type rate expressions. All process rates and diffusion coefficients were assumed to betemperature dependent. Heterotrophic bacteria were assumed to be responsible for hydrolysis,
mineralization of organic matter (aerobic growth) and denitrification (anoxic growth).
Autotrophic bacteria were assumed to be responsible for nitrification, which is modelled as a
two-step process. Lysis was considered to be the sum of all decay and sink processes.
CW2D/HYDRUS-2D was developed to model processes in subsurface flow constructed wetlands
intended for wastewater treatment. However, the model has already been successful applied toconstructed wetlands for the treatment of combined sewer overflows (Dittmer et al., 2004;
Henrichs et al., 2005). Possible future applications could also include a) constructed wetlands forother types of wastewater, such as industrial and agricultural wastewaters having compositions
similar to municipal wastewater, b) infiltration of treated (waste)water into soil for groundwater
recharge or other applications, c) carbon and nitrogen cycling in soils, d) processes in natural
wetlands, and e) processes in riparian areas.

18

Implementation of CW2D into the HYDRUS-2D user interface
Pre-processing
To activate the CW2D module in the HYDRUS-2D graphical user interface the heading "CW2D"
has to be written in the "Main Processes" window (Figure 1). The main changes are concentrated
in the solute transport part. In the "Solute Transport" the number of solutes has to be set to 13 (12CW2D module compounds and one non-reactive tracer compound, that is independent from theother 12 compounds). Figure 2 and Figure 3 show the two user interface windows for the CW2D
model parameters.

Figure 1: "Main Processes" window.

Figure 2: "Solute Transport – CW2D Model Parameters I" window (time unit: d).

19

Figure 3: "Solute Transport – CW2D Model Parameters II" window (time unit: d).
Figure 4 shows the definition of the pressure head initial conditions in the "Boundary Conditions"window. Names of all compounds appear in the listbox in the left sidebar, with the "L" letter
denoting the concentration in the liquid phase and "S" the concentration in the solid phase.

Figure 4: "Boundary Conditions" window – initial conditions.
Post-processing
Listboxes in the "Graphical Display of Results" window (not shown) and the "Observation
Nodes" window (Figure 5) list names of all variables, including all biochemical compounds, thatcan be displayed.

20

Figure 5: "Observation Nodes" window.
Publications:
Dittmer, U., Meyer, D., Langergraber, G. (2005): Simulation of a Subsurface Vertical Flow ConstructedWetland for CSO Treatment. Water Science and Technology 51(9), 225-232.
Henrichs M., Langergraber G., Uhl M. (2005): Modelling of organic matter degradation in constructedwetlands for treatment of combined sewer overflow. In: Proceedings International Symposium on"Wetland Pollutant Dynamics and Control – WETPOL", 4-8 September 2005, Gent, Belgium.
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. (2003): Simulation of subsurface flow constructed wetlands - Results and further
research needs. Water Science and Technology 48(5), 157-166.
Langergraber, G. (2005a): The role of plant uptake on the removal of organic matter and nutrients in
subsurface flow constructed wetlands – A simulation study. Water Science and Technology 51(9), 213223.
Langergraber, G. (2005b): Simulation of the treatment performance of outdoor subsurface flow
constructed wetlands in temperate climates. In: Proceedings International Symposium on "Wetland
Pollutant Dynamics and Control – WETPOL", 4-8 September 2005, Gent, Belgium.
Langergraber, G., Šimůnek, J. (2005): Modeling variably-saturated water flow and multi-component
reactive transport in constructed wetlands. Vadose Zone Journal (accepted).
Šimůnek, J., Šejna, M., van Genuchten, M.Th. (1999): The HYDRUS-2D software package for simulating
two-dimensional 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, 251p.
Toscano, A., Langergraber, G., Consoli, S., Cirelli, G.L., Haberl, R. (2005): A simulation study for pilotscale 2-stage subsurface flow constructed wetlands for wastewater treatment. In: ProceedingsInternational Symposium on "Wetland Pollutant Dynamics and Control – WETPOL", 4-8 September
2005, Gent, Belgium.

21

Numerical Analysis of One-Dimensional Coupled Water,
Vapor and Heat Transport in the Vadose Zone using
HYDRUS
1

2

Hirotaka Saito1, Jiri Šimůnek1, and Binayak P. Mohanty2

Department of Environmental Sciences, University of California, Riverside, CA
Biological and Agricultural Engineering, Texas A&M University, College Station, TX
ABSTRACT

Although it is well known that water and/or vapor flow and heat transport processes areclosely coupled and strongly affect each other in the vadose zone of arid or semi-arid regions,
their simultaneous interactions are rarely considered. In this study, we implemented in the
HYDRUS-1D code the coupled movement of water, vapor, and heat in the subsurface, as well asinteractions of these subsurface processes with the mass and energy balance at the soil surface.
The code considers the movement of liquid water and water vapor in the subsurface to be driven
by both pressure head (isothermal transport) and temperature (thermal transport) gradients. The
heat transport module considers movement of soil heat by conduction, convection of sensible heatby water flow, transfer of latent heat by diffusion of water vapor, and transfer of sensible heat bydiffusion of water vapor. The developed code allows a very flexible way of using various types of
meteorological information at the soil-atmosphere interface for evaluating surface water andenergy balance.

INTRODUCTION
The simultaneous movement of water, vapor and heat in the vadose zone of arid
or semi-arid regions is of great interest in evaluating water and energy balance ofsubsurface environments for both agriculture and engineering applications. Although it is
widely recognized that the movement of water, vapor, and heat are closely coupled andstrongly affect each other, their mutual interactions are rarely considered in practicalapplications. The effect of heat transport on water flow is often neglected, arguably
because of model complexity and a lack of data to fully parameterize the model, butpartly also because of the unavailability of a user friendly simulation code.
The main objective of this study thus was to develop a numerical model forcoupled water, vapor, and heat transport that, together with a surface water and energybalance, can be used to predict volumetric soil water contents and soil temperatures underfield conditions. The goal was to develop a model that would provide much flexibility inaccommodating various types of meteorological data that can be collected at various timeintervals (daily, hourly, or other time intervals). The coupled movement of water, vapor,and heat in the subsurface, as well as an energy and water balance at the soil surface, wasimplemented in the HYDRUS-1D code (Šimůnek et al., 1998).
Model performance was subsequently evaluated by simulating continuouschanges in water contents, temperatures, and fluxes in a bare field soil. The coupledmodel was compared against field soil temperature and water content data collected at
different depths at a field site near the University of California Agricultural Experiment
Station in Riverside, California, during the fall of 1995.

22

THEORY
Water and Vapor Flow
The governing equation for one-dimensional flow of water and vapor in a variably
saturated rigid porous medium is given by the following mass conservation equation:
∂θ
∂ ⎡
∂h
∂T
∂h
∂T ⎤
= ⎢ K Lh ( h )+ K Lh ( h ) + K LT ( h )+ K vh+ K vT−S∂t ∂z ⎣
∂z
∂z
∂z
∂z ⎥⎦
∂ ⎡
∂h
∂T ⎤
= ⎢ KTh+ K Lh + KTT
−S
∂z ⎣
∂z
∂z ⎥⎦
where θ is the total volumetric water content (-), h is the pressure head (L), T is thetemperature (K), t is time (T), z is the spatial coordinate positive upward, S is a sink term,
KLh (LT-1) and KLT (L2K-1T-1) are the hydraulic conductivities for liquid phase fluxes due
to gradient of h and T, respectively, Kvh (LT-1) and KvT (L2K-1T-1) are the isothermal andthermal vapor hydraulic conductivities (Philip and de Vries, 1957), respectively, and KTh(LT-1) and KTT (L2K-1T-1) are the isothermal and thermal total hydraulic conductivities,
respectively.

[1]

Heat Transport
The governing equation for the movement of energy in a variably saturated
porous medium is given by the following energy conservation equation:∂C p T
∂θ
∂q
∂q T
∂q T
∂ ⎡
∂T ⎤
+ L0 v = ⎢λ (θ ) ⎥ − Cw L − L0 v − Cv v − Cw ST
∂t
∂t
∂z ⎣
∂z ⎦
∂z
∂z
∂z
where T is the temperature (K), θn is a volumetric fraction of solid phase (-), Cn, Cw, Cv
and Cp are volumetric heat capacities (Jm-3K-1, ML-1T-2K-1) of solid phase, liquid water,water vapor, and moist soil (de Vries, 1963) , respectively, L0 is the volumetric latent heatof vaporization of water (Jm-3, ML-1T-2), λ(θ) is the apparent thermal conductivity of soil
(Jm-1s-1K-1, MLT-3K-1), and qL and qv are flux densities of liquid water and water vapor(LT-1), respectively.

[2]

Surface Energy Balance
Surface precipitation-irrigation-evaporation and heat fluxes are used as boundary
conditions for water and vapor flow and heat transport, respectively. Evaporation andheat fluxes can be calculated precisely only from the surface mass and energy balance(e.g., van Bavel and Hillel, 1976):
Rn − H − LE − G = 0-2
where Rn is net radiation (Wm , MT-3), H is the sensible heat flux density (Wm-2, MT-3),
LE is the latent heat flux density (Wm-2, MT-3), L is the latent heat (J kg-1, L2T-2), E is theevaporation rate (ML-2T-1), and G is the surface heat flux density (Wm-2, MT-3). While Rn
and G are positive downward, H and LE are positive upward. The modified version of
HYDRUS-1D can either directly use measured continuous values of net radiation, Rn, andother meteorological variables, or it can calculate continuous vales from daily
meteorological data (e.g., Saito et al., submitted).

23

[3]

Numerical Simulations
Study Site
Soil temperature and water content variations were measured near the soil-air
interface at 20 and 40 minutes intervals, respectively, at depths of 2, 7, and 12 cm near
the University of California Agricultural Experiment Station in Riverside, California,
during the fall of 1995 (November 23, “day of year” or DOY 327, through December 5,
DOY 339; Mohanty et al. 1998). The experimental field was irrigated twice during themeasurement period on DOY 334 and DOY 335 with a sprinkler system. Daily standard
meteorological information obtained from the weather station in Riverside was used. The
soil at the experimental site was characterized as a Arlington fine sandy loam (coarseloamy, mixed, thermic, Haplic Durixeralf).

ResultsoTemperature C

Volumetric Water Content

Soil water contents and soil
temperatures at three different depths (2,0.3depth=2cm
7, and 12 cm) at the experimental sitewere numerically simulated from DOY
328 to DOY 339 using the coupled water,
0.2
vapor, and heat transport features
implemented in HYDRUS 1D. Initialconditions were determined from0.1
measured values on DOY 328.
Fig. 1 shows measured andsimulated water contents at the 2 cm
328
330
332
334
336
338
340
depth. Predicted water contents follow
DOY
the measured values fairly well during
Fig. 1: Simulated (line) and measured (openthe entire simulation period. Rapid
circles) water content at 2 cm depth.
increases of the water content after
irrigation on days 334 and 335 are
35
depth=7cm
predicted well. Fig. 2 depicts simulated
30
and measured soil temperatures at the 725
cm depth. Both simulated and measured
temperatures show typical sinusoidal
20
diurnal behavior.
15
Depth profiles of simulated vapor
10
fluxes and liquid water fluxes after5
irrigation are depicted in Fig. 3. Liquidwater fluxes in the top layer are upward
0
328
330
332
334
336
338
340
due to evaporation, while they are
DOY
downward (i.e., negative fluxes) in theFig. 2: Simulated (line) and measured (openrest of the soil profile due to further
circles) soil temperatures at 7 cm depth.
infiltration and redistribution of irrigatedwater. Water vapor flows downward
during the day across the soil profile due to temperature gradients, while it moves upward
during the night. Vapor fluxes become much smaller after irrigation due to increased
24

Simulated
MeasuredDepth [cm]Depth [cm]gaseous tortuosity
in the reduced soil00
air content. During-4-4
the day, downward
vapor fluxes thus
-8-8transport latent heat-12
-12
to deeper depths,
while during night-16-16
the latent heat flux
-20
-20
-0.2
-0.1
0
0.1
0.2
-0.1
-0.05
0
0.05
moves
in
the
Water Flux [cm/day]
Vapor Flux [cm/day]upward direction.
Fig. 3: Simulated vapor flux, and liquid water flux at noon of DOY
Liquid water fluxes
334, and 6, 12, 18, 24 hours after. Light irrigation was applied atreflect mainly thenoon of DOY 334.
water
content
distribution in the soil and not the temperature distribution. This coupled physical processthus indeed requires a mathematical/computer model that describes the coupled
movement of water, vapor, and heat.
This study demonstrates capabilities of the modified HYDRUS-1D program tosimulate coupled water, vapor, and heat transport. Limited daily meteorologicalinformation can be used to predict soil temperatures and water contents reasonably wellusing the model.
DOY 334.5
DOY 334.75
DOY 335
DOY 335.25
DOY 335.4

Acknowledgements
This work was supported in part by SAHRA (Sustainability of semi-Arid Hydrology and Riparian
Areas) under the STC Program of the National Science Foundation, Agreement No. EAR-9876800
and the Terrestrial Sciences Program of the Army Research Office (Terrestrial Processes andLandscape Dynamics and Terrestrial System Modeling and Model Integration). Mohanty
acknowledges the support by NASA-Global Water and Energy Cycle grant NAG5-11702.

REFERENCES[1] 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.
[2] Mohanty, B.P., P.J. Shouse, and M.T. van Genuchten. 1998. Spatio-temporal
dynamics of water and heat in field soil. Soil Tillage Res. 47:133-143.
[3] Philip, J.R., and D.A. de Vries. 1957. Moisture movement in porous materials undertemperature gradient. Trans. Am. Geophys. Union. 38:222-232.
[4] Saito, H., J. Šimůnek, and B.P. Mohanty. Using Surface Energy Balance and Limited
Meteorological Information for Numerical Analysis of Coupled Water, Vapor and Heat
Transport in the Vadose Zone, Vadose Zone J. Submitted.
[5] Šimůnek, J., M. Sejna, and M.T. van Genuchten. 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. IGWMC-TPS-70. International
Groundwater Modeling Center, Colorado School of Mines, Golden, CO.
[6] van Bavel, C.H.M., and D.I. Hillel. 1976. Calculating potential and actual evaporationfrom a bare soil surface by simulation of concurrent flow of water and heat. Agri.
Meteorol. 17:453-476.

25

Imprecise simulation of salt dynamic and balance – a Hydrus-1D flaw?
Irina Forkutsa and Rolf Sommer
Center for Development Research (ZEF), University of Bonn
Walter-Flex-Str.3, 53113 Bonn
Tel: +49 228 73 4636
e-mail: forkutsa@yahoo.com
web-page: http://www.zef.de
Introduction
The arid and semi-arid region of Khorezm in Uzbekistan is characterized bylong and dry summers and short rainy/snowy winters; the agricultural production fully
relies on irrigation. The development of huge-scale and massive irrigation schemes
under conditions of improperly working drainage systems led to rising groundwatertables, drastically changing the local and regional hydrology. Salt accumulation in theroot zone of the crops is an almost inevitable consequence adversely affecting crop
growth. However over the past century farmers gained knowledge and experience how
to manage the situation and use saline soil and partly saline water to produce crops.
First of all, by applying high water amounts of water before the cropping period in
March ("leaching") salts are washed out of the soil profile. Secondly, high amount ofirrigation water might also temporarily decrease salt concentration in the root zone("irrigation-leaching").
The objectives of this study was to estimate the water and salt balances as wellas the capillary rise of saline shallow groundwater at irrigated field in Khiva district of
the Khorezm region. In this paper the 2003-results of the salt balance estimation for one
location at a sandy loam field are presented. To understand the processes of water flow
and salt transport in soils of Khorezm, we applied Hydrus-1D.
Methods
A field experiment was conducted to monitor the water and salt distribution andto assess small-scale spatial differences. The research was focused on three points in asandy loam field during the vegetation seasons of 2002 and 2003. The groundwatertable and its salinity, the electrical conductivity of soil saturation paste (ECp), and thegravimetric soil water content and soil tension were repeatedly determined at different
depths along the field. A factor of 3.5 was used to convert ECp into electricalconductivity of the saturated extract (ECe). Furthermore, based on own data andregression analyses, a factor 0.042 was applied to convert ECe into salt concentrations
26

(TDS in g 100 g-1 soil). Tensiometers were installed prior to the first irrigation event (16July) and were kept in the sandy loam field until the end of the vegetation season (30September). Meteorological data were monitored with a weather-station in the vicinityof the fields; applied water during three leaching and five irrigation events wasquantified with Cipoletti weirs. Soil hydraulic properties were estimated from bulk
density and soil texture with Rosetta DLL v.1.1. The single porosity model of Van
Genuchten-Mualem with an air entry value of -2 cm was selected as the hydraulicmodel. The Van Genuchten parameters were inversely optimized based on tensiometers
readings and volumetric soil moisture. The atmospheric boundary conditions with asurface layer and the groundwater table and its salinity were selected as upper and low
boundary conditions. Solute transport and reaction parameters as the longitudinal
dispersivity (Ldis) and the adsorption isotherm coefficient (Kd) were inversely optimized
based on salt concentration along the 1m-depth soil profile and of groundwater. They
were constrained to 1<Ldis≤1000 and 0<kd<0.35.
The salt balance of the rooting zone can be calculated using two Hydrus output
files: solute.out and nod_inf.out. The nod_inf.out file provides information about soil
moisture, water fluxes and salt concentration of the simulated profile for any selectedday. In the study two ways to calculate the salt balance were used:(a)

Method 1
Multiplying the soil moisture (θ; cm3 cm-3) with salt concentration (C; mg cm-3)and the length of each particular soil increment between two consecutive nodes (z; cm)
at a considered point in time (t) and adding it up over the rooting zone results in theamount of salt stored in the soil profile at that date. The difference between two
consecutive dates reveals the amount of salt that has been added to, or was leached outof the rooting zone (∆Q; mg cm-2):
N

N

N=1

N=1

∆Q = ∑ [θ i ⋅ Ci ⋅ (zi+1 − zi )]t1 − ∑ [θ i ⋅ Ci ⋅ (zi+1 − zi )]t 2

(b)

Eq. 1

Method 2
Adding cumulative inputs (via irrigation [file solute.out] and/or upward

movement of salts [file nod_inf.out]) and outputs of salts (via leaching [file nod_inf.out]
at the bottom of the rooting zone) for a given period:∆Q = ∑ {[Ctop ⋅ Fluxtop ⋅ (( j + 1) − j)] − [Cbot ⋅ Fluxbot ⋅ (( j + 1) − j)]}
t2

t1

27

Eq. 2,

whereas Flux is the flow of water (cm d-1) and j and j+1 represent the current
and next time levels of the simulation run.
Results and discussion
The model only poorly reproduced the rapid salt concentration fluctuations(Figure 1) during the simulation period from 13 February till 18 November (RMSE
equal to 1.86 to 3.22 mg cm-3).

80 cm

16
-1

-1

20/Dec

5/Nov

21/Sep-1

RMSE=1.86 dS m

20/Dec

7/Aug

23/Jun9/May

25/Mar

Simulated

8/Feb

0

20/Dec

0

5/Nov

4

21/Sep

4

7/Aug

8

23/Jun

8

9/May

12

25/Mar

7/Aug

Measured

12

8/Feb

-1

23/Jun105 cm16-1

RMSE=1.97 dS m

ECe [dS m ]

9/May

8/Feb

50 cm

16

25/Mar020/Dec

0

5/Nov

4

21/Sep

4

7/Aug

8

23/Jun

8

9/May

12

25/Mar

12

8/Feb

-1

ECe [dS m ]RMSE=2.09 dS m

5/Nov

20 cm
RMSE=3.22 dS m

21/Sep16

Figure 1 Comparison of measured and simulated electrical conductivity of saturated
extract at different depths of sandy loam field in 2003 (middle location)
This might be due to wrong estimates of the solute transport and reaction
parameters (Ldis and Kd), which are essential for the simulation of solute transport.
Furthermore, there was strong influence from the saline groundwater (bottom boundary
condition), which fluctuated comparatively less. Nevertheless, the simulation results ofsolute transport provided reasonable insight into seasonal salt balance of the monitored
fields over the cotton root zone, because the simulated and observed salt concentrationin the beginning and in the end of the simulation period were in good agreement.
Determining the salt balance of the rooting zones with the two different methods
revealed some discrepancies (Table 1). Considering the whole season, the calculated28

salt storage by method 2 was two times lower than with method 1, constituting 35 t ha-1
versus 18 t ha-1 of leached salts. The reasons for discrepancies between the two methods
remain vague. A six-hour time step was used to compute the salt quantities according toequation 2. Thus a small rounding error in these calculations was tolerated(extrapolating momentary concentrations and fluxes). Nevertheless it should not amountto ~100 % deviation between these two applied methods. Possibly, discrepancies are theresult of an accumulating – and thus emphasizing – mass-balance error when using offrounded water fluxes or salt concentrations provided by the model.
Table 1 Salt balance [t ha-1 90cm-1] calculated by the two methods (eq. 1 and 2) of the
sandy loamy field (middle location), 2003
Method 2
Salt outflow
Salt input via(drainage) at theleaching or
bottom of profile
irrigation
0.0
0.0
3.9
-42.2
7.3
-46.9
7.3
-42.2

Method 1
Date

Salts stored

13-Feb
25-Apr
04-Oct31-Dec

46.3
29.1
25.8
28.5

Differences stored
13 Feb - 25 Apr
25 Apr - 4 Oct
4 Oct - 31 Dec
13 Feb - 31 Dec

-17.2
-3.3
2.7
-17.8

-38.3
-1.3
4.7
-34.9

Conclusion
The solute transport simulation failed to mimic the rapid change in the observed
soil salinity. Further research is necessary to detect the underlying causes. Additionally,
there remains the need for checking the causes of mismatching salt balances calculated
based on either salt-store or salt-in/outflow.

29

Application of Hydrus2D to derive a water balance in akarst dominated catchment
Manfred Fink & Marcel Wetzel
Institute of Geography, Department of Geoinformatics,
Friedrich-Schiller-University Jena, Loebdergraben 32, 07743 Jena

1.

Introduction

The forest of the national park Hainich is one of the largest deciduous forest ecosystems in
middle Europe. In the centre of the protected zone the forest is unmanaged since more then 50
years. The age of the trees goes up to 250 years. Therefore this site is used to determine,whether old woods are carbon sinks or they are in equilibrium.
Hence the Hainich was instrumented by a variety of sensors to measure carbon fluxes e.g. an
eddy flux tower for the measurement of carbon dioxide exchange between forest andatmosphere. Resulted balances of carbon indicates a sink term of about 5 t carbon/ha*a for thewhole system. This can not be explained by the measurements of storages and fluxes.
This result suggests to quantify the amount of carbon leaving the system as dissolved carbon
related with the movement of water. While the hydrology in the study area is characterized bya karst system it is almost impossible to measure the discharge directly. We measure soil
moisture with Frequency Domain Reflectrometry (FDR) probes to estimate the hydrologic
fluxes of the Hainich site. To get a reasonable water balance we used the Hydru2D model.

2.

Study area

The study area is located in the Hainich site in the western part of Thuringia, Germany. The
altitude varies between 480 and 400 m above sea level. The annual average precipitation isapprox. 700 mm. Annual average temperature ranges around 7° C. The geology is dominatedby limestones with karstic phenomen. Soils developed from this bedrock are cambisols with ahigh content of clay (40 – 65 %) increased with soil depth. During events of high precipitationor snow melt the clay rich subsoil with its low hydraulic conductivity which leads to lateralsubsurface flow. The forest is dominated by beech (Fagus sylvatica) combined with ash(Fraxsinus excelsior) and sycamore maple (Acer pseudoplatanus).

3.

Measurement program

To study the water dynamics we instrumented hill slope catena in the study area withequipment for hydrological measurements. The measurement period was from August 2001 to
November 2003. In particular we installed 13 FDR-Probes stations with 5 sensors in depths
varying between 10 and 100 cm. At the valley bottom we installed one groundwater gauges
which measures only periodic saturation. Additionally one runoff gauge was set up tomeasure the periodic stream (Figure 1). A rain gauge at every station was completed tomeasure trough fall. Observed data was recorded by data loggers with telemetry unit in a 5
minute interval. To characterize soils we accomplished field campaigns. Therefore we tooksamples for soil texture, organic carbon, bulk density, saturated hydraulic conductivity andsoil porosity. These have been analysed in the laboratory. The observed data are used formodel parameterization and validation.
30

475.0
463.8

H2

452.5

#

H3

#

441.3
430.1
407.6

418.8

Überhöhung
X
Hight
factor9.7
9.7

AltitudeHöhe
above
NN inlevelübersea
m [m]

# H1

#

H6

H7

#

#

H8
H9

C-Probe#

Runoff gauge

T$

Groundwater gauge
240

H5

#

#

#

120

H4

#

H10

#

H11

#

360

H12

H13

# $
T
#
480

600

720

840

960

1080

1200

Distanz
von derin mDistancefrom
theHangkuppe
hill top [m]Figure 1: Profile of the instrumented hill slope catena.

4.

Modelling

To model the soil water dynamics we used the vadose zone model Hydrus2D (Šimůnek & van
Genuchten 1999). The model domain for the instrumented hill slope was 1137m long and thealtitude difference was 64.5m which leads to a very gentle slope of 5.7%. The parameterizedsoil depth was 1.4 m at the upper hill, and 2 m in the valley. The domain consists of 246950
nodes (Figure 2).

Figure 2: Hill slope representation of the test Site in Hydrus2D (Meshgen).

For parameterising the soil we used the van Genuchten (1980) model to adjustment themeasured pF-Curves for three different soil horizons. The boundary conditions were defined
as atmospheric boundary at the surface, no flux left at the hill top, seepage face right at the31valley bottom and free drainage as lower boundary condition. As climatic input we used datafrom an adjacent climate station. Evapotranspiration rates were calculated by using Penman
Monteith method (Monteith 1975). We used the measured through fall and the measured
precipitation from that climate station to derive the interception. Since Hydrus2D requires
both evaporation and transpiration separately we divided the potential evapotranspiration
through values from the literature (Brechtel 1973). To calculate the actual transpiration the SShape method (van Genuchten 1987) was used. The modification of precipitation due to snow
dynamics was estimated according to a simple snow model (Bergström 1992). The model run
included the time period from 1/2001 to 11/2003 on a hourly time step.

5.

Results and discussion

0

80

5

60

10

40

15

kt

ct
O

O

l
Ju

r
Ap

n
Ja

kt

O

O

Ju

l

25
ct

100

0

80

5

60

10

40

15

20
0

20

R = 0,78 R² = 0,61
Ja

n

A

pr

mm

%

Ja

n

0

20

R = 0,87 R² = 0,75
r

20

mm

100

Ap

%

The data of the measured soil water dynamics were used to validate the model. Figure 3
shows an example for the station H4 in the middle of the catena in the soild depth of 10 and
20 cm. The predicted curve is damped compared to the measured one. In contrast to this thecoefficient of determination indicates a rather good fit. The damping in the predicted curve iscaused by the fact that Hydrus2D only considers matrix flux. In these soils the macro pores
are important because of the high clay content. This leads to soil cracks during the summer
period. Additionally soil cracks were build due to a high biological activity. Given thesenatural conditions high soil moisture dynamic is observed. This can not be described by themodel. Another problem was detected in spring 2003. Here an early drop of the predicted soil
moisture is indicated. It can be marked that the model Hydrus2D is not able to follow up theeffects of the frozen soil. Thereby actual soil moisture remains within the profile if comparing
to modelling results. A further problem was that interflow could not has been simulated byHydrus2D which was observed in the field. The reason for this might be the very gentle slopeof the catena. Additionally the predictions for the stations in the valley are comparatively
poor. It originates from laterally water flowing of the side slopes which can not be describedby a 2 dimensional model.

25
J

ul

ckt

O

Ja

n

predicted

A

pr

J

ul

ckt

O

measured

Figure 3: Measured and predicted soil moisture of station H4 in 10 and 20cm depth in the years 2002 and2003. On the second X-axis hourly precipitation.

In spite of the shown problems we gain a reasonable water balance for the measured catena.
Figure 4 shows typical characteristics for the components of the water balance. For example
the discharge started in winter and reached a maximum in March which fits well to general32expectations. The characteristic of the transpiration during the vegetation period shows atypical dynamic as well. Consequently the modelled discharge is valid to be used for coupling
solute carbon measurements for deriving carbon export rates (Fink et al. 2004).
120

80

40

Oct 02

Sept 02

Aug 02

July 02

June 02

May 02

Apr 02Mar 02

Feb 02

Jan 02

Dec 01-40

Nov 010-80Discharge

Evaporation

Transpiration

net. Precipitation

Storage

Figure 4: Monthly water balance in mm for the whole catena for the year 2002.

6.

References

Bergström, S. (1992): The HBV model – its structure and applications.– SMHI, RH, 4,
Norrköping, Sweden.
Brechtel, H.M. (1973): Ein methodischer Beitrag zur Quantifizierung des Einflusses von
Waldbeständen
verschiedenerBaumartenundAltersklassen
auf
dieGrundwasserneubildung in der Rhein-Main-Ebene.– Zeitschrift der Deutschen
Geologischen Gesellschaft 124, 593–605.
Fink, M., Krause, P., Gleixner, G. & Flügel, W.A. (2004): Verknüpfung von Bodenfeuchteund DOC – Messungen zur Abschätzung des Kohlenstofftransportes durch Wasser ineinem Waldökosystem – Eine Fallstudie.– In: Bronstert, A. et al.: Wasser und
Stofftransport in heterogenen Einzugsgebieten.– Forum für Hydrologie undWasserbewirtschaftung, Heft 05.04, Band 2, 11–15.
Genuchten, M.T. van (1980): A closed-form equation for predicting the hydraulicconductivity of unsaturated soils.– SSSAJ 44, 892–898.
Genuchten, M.T. van (1987): A numerical model for water and solute movement in and below
the root zone.– Research Report 121, USDA-ARS, US Salinity Laboratory, Riverside,
California.
Šimůnek, J., Sejna, M. & van Genuchten, M.T. (1999): The HYDRUS-2D software package
for simulating two-dimensional movement of water, heat, and multiple solutes invariably saturated media. Version 2.0.– IGWMC – TPS – 53, International Ground
Water Modeling Center, Colorado School of Mines, Golden, Colorado.
Monteith, J.L. (1975): Vegetation and atmosphere, Vol. 1 Principles.– Academic Press,
London.
33

Parameter estimation of soil hydraulic properties from tension
infiltrometer data in Portugal
T. B. Ramos1, M. Th. van Genuchten2, M. C. Gonçalves1, J. C. Martins1 & F. P. Pires1
1

Estação Agronómica Nacional, Quinta do Marquês, 2784-505 Oeiras, Portugal – Tel: (+351) 214
403 500 – Fax: (+351) 214 416 011 – E-mail: Tiago_Ramos@netcabo.pt
2
George E. Brown, Jr. Salinity Laboratory, USDA-ARS, 450 W. Big Springs Road, Riverside, CA
92507-4617, USA.

INTRODUCTION
Mathematical models are increasingly used to address a broad range of variablysaturated flow and contaminant transport problems. Such simulations are generally based on
numerical solutions of the Richards equation, which in turn require knowledge of theunsaturated soil hydraulic properties.
One of HYDRUS applications is the evaluation of the hydraulic properties by numericalinversion from tension disc infiltrometer data since Šimunek and van Genuchten (1996)
proposed an alternative methodology to the traditional Wooding’s (1968) analytical solution
analysis for unconfined steady state infiltration from a disc. Šimunek and van Genuchten(1996) new approach consists in the formulation of the inverse problem and is based onminimization of a objective function, Φ, during the parameter estimation process.
Minimization of the objective function Ф in HYDRUS-2D is accomplished using theLevenberg-Marquardt non-linear minimization method (Marquardt, 1963).
This brief document summarizes a recently submitted paper (Ramos et al., 2005) where
we further test the inverse modelling approach of Šimunek and van Genuchten (1996) byusing the method to characterize the hydraulic properties of 4 field sites in Portugal. We also
compare the resulting hydraulic properties with independent estimates using Wooding’s
analysis, and further interpret the results in terms of soil macroporosity measurements.
MATERIALS AND METHODS
The field tension infiltration experiments were carried out in Aljustrel and Alvalade(Alentejo), Portugal, using two experimental areas cropped with maize and irrigated with acenter-pivot irrigation system. The infiltration experiments were performed on three different
Gleyic Luvisols (LVgl), and on an Haplic Fluvisol (FLha) (soil classification according to
ISSS-ISRIC-FAO, 1998).
Soil hydraulic properties were determined from undisturbed samples, collected from thesuperficial layers. The soil water retention curve was determined using suction tables withsand or kaolin for h> -500 cm and a pressure plate apparatus for h< -1000. The evaporation
method was further used to simultaneously estimate both water retention and hydraulicconductivity data for h >-1000 cm.
The field tension infiltrometer measurements were performed twice at each location,
designated as Run A and Run B. Both runs were carried out at a distance of about 1 m fromeach other using tension infiltrometers with the discs (having a radius of 20 cm) detachedfrom the supply and tension control tubes. All four infiltration experiments were conductedwith consecutive supply pressure heads of -15, -6, -3 and 0 cm. Readings of the water supply

34

Cum. Inf. (cm3)tube were done visually. Figure 1 exemplifies a measurement cumulative infiltration rateversus time on LVgl 1. Using the capillary rise equation, the applied pressure heads
correspond to the pore radii limiting the pore classes in the classification suggested by Watsonand Luxmoore (1986).

300

Cumulative Infiltration
LVgl 1

250

0 cm

200

- 3 cm150
100- 15 cm- 6 cm

50
0
0
Run A

50

100
time (min)Adjust. A

150
Run B

200
Adjust. B

Figure 1. Measured cumulative infiltration rates at consecutive supply pressure heads of -15,
-6, -3 and 0 cm for 4 superficial layers of one of the Gleyic Luvisols used in this study.

Disturbed gravimetric samples were taken to determine the initial and final watercontents of the soils. The initial water content was determined at a different location from
where the infiltration took place to avoid disturbance of the soil during the measurements.
However, the final water content was determined under the disc membrane after reaching a
steady flux for the last pressure head.
Following Šimunek et al. (1998a; 1998b), the objective function Φ(Q, θi, θf) was
defined in terms of the measured cumulative infiltration data (Q) at multiple pressure heads,
and the initial and final soil water contents (θi and θf, respectively), while the calculations andoptimisation were carried out using HYDRUS-2D (Šimunek and van Genuchten, 1996). The
weighting coefficients wi in Eq. (1) for the different infiltration data points were all assumed
to be 1 since the observation errors of each measurement are unknown. Although being a
isolated entry in the optimisation, the final water content was given a weight of 10 toguarantee a reasonable effect on the final results relative to the cumulative infiltration data.
RESULTS AND CONCLUSIONS
Our study suggests that numerical inversion of tension infiltrometer data provides a
relatively simple and reliable alternative method for determining the water retention andconductivity curves of unsaturated soils.
The water retention and conductivity curves obtained by numerical inversion closely
matched the laboratory measured curves for the 4 superficial layers where the infiltration
experiments were carried out (Figure 2). Hydraulic conductivities obtained with the inversely
estimated Mualem-van Genuchten parameters also corresponded well with results obtainedusing Wooding’s traditional approach of disc infiltrometer data following the methodology ofAnkeny et al. (1991). This correspondence was further reflected by the very similar estimatesof the macroporosity, mesoporosity 1, and mesoporosity 2 pore classes we calculated from theestimated K(h) curves using the method of Watson and Luxmoore (1986).

35

Water retention curve
LVgl 1 (0-33 cm)

0.5

1E+3

0.4

RET C

0.3

Run A

Lab. Met.
RET C
Run A

1E+1

K (cm d-1 )Wt. cont. (cm 3 cm -3 )Lab. Met.

Run B

0.2

Hydraulic conductivity
LVgl 1 (0-33 cm)

0.1

1E-1

Run B
Ankeny et al.

1E-3
1E-5
1E-7
1E-9

0
0

1

2
3
log h (cm)

4

0

5

1

2

3
4
log h (cm)5

6

Figure 2 – Water retention θ(h) and hydraulic conductivity K(h) curves obtained through numericalinversion of the measured tension disc infiltrometer data (Runs A and B), and by means of separate
laboratory and field measurements after analysis with RETC. Results are for the surface layer of oneof the Gleyic Luvisols (LVgl 1) used in the experiments.

One major limitation of the numerical inversion method is its extreme dependence ofthe field-measured water content values. The saturated water content θs estimated bynumerical inversion in particular was very close to the final water content θf measured at theend of the infiltration tests, which indicates the extreme dependence of θs identification on the
final water content value. While this situation was somewhat expected since a weighting
coefficient of 10 was used for the final water content it does illustrate how poor estimates ofθs can have a negative effect on the water retention curve. Some problems were foundespecially with the initial soil water content, which (unlike the final water content) cannot be
determined at exactly the same location as where the tension infiltrometer measurements arecarried out. Problems of soil spatial variably can become very important. This was shown inthis study for the LVgl2 and LVgl3 field measurements which produced curves that deviatedconsiderable from the laboratory derived curves (Figure 3).
Water retention curve
LVgl 3 (0-25 cm)

1E+3

Lab. Met.
RET C
Run (θi)
Run (pF)
T
A ( F)

0.4
0.3

Lab. Met.
RET C
Run (θi)
Run (pF)
Anken y et al.

1E+1
K (cm d-1 )Wt. Cont. (cm 3 cm -3 )0.5

Hydraulic conductivity
LVgl 3 (0-25 cm)

0.2
0.1

1E-1
1E-3
1E-5
1E-7

0

1E-9

0

1

2

3

log h (cm)

4

0

5

1

2
3
log h (cm)

4

5

6

Figure 3 – Influence of a water content underestimation in the soil hydraulic curves
determination [Run (θi)] and it’s correction through the addition of water contents at 10 kPa and at
1500 kPa [Run (pF)] in the objective function.

Addressing essentially the same problem (i.e., to obtain a better description of the water
characteristic in the dry region), Schwartz and Evett (2003) recently suggested that at least
one independent measurement of θ(h) at a pressure head sufficiently less than the lowest h0 in
the objective function should be included in the objective function. To further test this, we reran the LVgl3 simulations, but now included measurements of θ(h) at pressure heads of-100 and –15,000, corresponding to pF values of 2 and 4.2, respectively. We obtained morereliable results when this independently measured water contents at -100 cm and -15000 cm

36

were added to the objective function. While this will require more time and effort (andadditional equipment), and as such negate some of the advantages of numerical inversion
methods (i.e., speed and ease of use), the results will be more reliable in the absence of having
more complete laboratory data for comparison purposes.
REFERENCES
ISSS-ISRIC-FAO. 1998. World Reference Base For Soil Resources. World Soil Resources Report 84.
FAO, Rome, Italy.
Marquardt, D. W. 1963. An algorithm for least-squares estimation of non-linear parameters. SIAM J.
Appl. Math., 11: 431-441.
Ramos, T. B.; M. Th. van Genuchten, M. C. Gonçalves; J. C. Martins and F. P. Pires, 2005. Hydraulicproperties evaluated by numerical inversion from tension disc infiltrometer data and using
laboratory methods. Vadose Zone Journal (submitted June, 2005).
Schwartz, R. C. and S. R. Evett. 2003. Conjunctive use of tension infiltrometry and Time Domain
Reflectometry for inverse estimation of soil hydraulic properties. Vadose Zone Journal, 2: 530538.
Šimunek, J. and M. Th. van Genuchten. 1996. Estimating unsaturated soil hydraulic properties fromtension disc infiltrometer data by numerical inversion. Water Resour. Res., 32(9): 2683-2696.
Šimunek, J., R. Angulo-Jaramillo, M. G. Schaap, J. P. Vandervaere and M. Th. van Genuchten.
1998a. Using an inverse method to estimate the hydraulic properties of crusted soils fromtension disc infiltrometer data. Geoderma, 86(1-2): 61-81.
Šimunek, J., D. Wang, P. J. Shouse and M. Th. van Genuchten. 1998b. Analysis of a field tension
disc infiltrometer data by parameter estimation. Int. Agrophysics, 12: 167-180.
Van Genuchten, M. Th., F. J. Leij and S. R. Yates. 1991. The RETC code for quantifying thehydraulic functions of unsaturated soils. Environmental Protection Agency, United States.
Watson , K. W. and R. J. Luxmoore. 1986. Estimating macroporosity in a forest watershed by use ofa tension infiltrometer Soil Sci. Soc. Am. J. 50: 578-582.
Wooding, R. A. 1968. Steady infiltration from large shallow circular pond. Water Resour. Res., 4:1259-1273.

37

Using HYDRUS to simulate water and solute transports in soil lysimeters
M. C. Gonçalves1, J. Šimůnek2, T. B. Ramos1, J. C. Martins1, M. J. Neves1 & F. P. Pires1
1

Department of Soil Science, Estação Agronómica Nacional, Av.República 2784-505 Oeiras, Portugal
Department of Environmental Sciences, University of California, Riverside, CA 92521, USA

2

INTRODUCTION
A variety of analytical and numerical models have been developed in the past decades topredict water and/or solute transfer processes between the soil surface and the groundwatertable. While many models quantifying solute transport in soils usually consider only one soluteand severely simplify various chemical interactions, models such as HYDRUS consider
multiple solutes and their mutual interactions.
HYDRUS has been recently updated to consider major ion chemistry [Šimůnek et al.,2005]. This new module considers transport of major ions (Ca2+, Mg2+, Na+, K+, SO42-, Cl-)and various temperature dependent equilibrium chemical reactions between these ions, suchas aqueous complexation, precipitation and dissolution of several solid phases (calcite,
gypsum, dolomite, nesquohonite, hydromagnesite, and sepiolite), and cation exchange.
The two most important chemical reactions in our application were aqueous complexation andcation exchange. Equations for the aqueous complexation reaction can be obtained using thelaw of mass action.
Partitioning between the solid exchange phase and the solution phase (cation exchange
reaction) is described in HYDRUS by the Gapon exchange equation [White and Zelazny,1986] assuming that the cation exchange capacity (CEC) is constant and independent of pH.
The first objective of this study was to carry our field experiments to quantify
salinisation and alkalisation risks of an Eutric Fluvisol irrigated with different quality watersin Alvalade-Sado (Alentejo), Portugal. The second objective was to evaluate the effectivenessof the HYDRUS-1D software package [Šimůnek et al., 2005] to predict water contents andfluxes, and concentrations of individual ions, overall salinity, as well as SAR (sodium
adsorption ratio) and ESP (exchangeable sodium percentage) indices under field conditions
where salinisation may occur.
MATERIAL AND METHODS
The quantification of salinisation and alkalisation hazards in a soil irrigated with four
different water qualities (Table 1) was carried out in Alvalade do Sado, Alentejo, Portugal.
Three soil monoliths (1.2 m2 x 1.0 m deep) were constructed in a Eutric Fluvisol (ISSSISRIC-FAO, 1998). The monoliths were laterally isolated with plastic to prevent lateral waterand solute fluxes and subjected to atmospheric conditions at the top and free drainage at thebottom. The soil monoliths were covered by annual spontaneous vegetation. TDR probes andceramic cups were installed in each soil monolith in two replicates at depths of 10, 30, 50 and
70 cm.
The monoliths were manually irrigated in years from 2001 to 2004 during the regular
irrigation period between May and September. Application depths were 10 mm each, for a
total of 500 mm per year. The monoliths were exposed to natural atmospheric conditions(rainfall and evapotranspiration) during the rest of the year. Due to variable Mediterranean
conditions, there were large differences in cumulative rainfall during the time period between
38

September and May, which corresponds with the soil’s leaching cycle. Irrigation watercompositions, named A, B and C according to the monolith to which they were applied, aredescribed in Table 1. The waters used in the experiment were produced by adding increased
concentrations of NaCl, CaCl2 and MgCl2 to the water A available in the region, with a ratioof Ca:Mg = 1:2.
Table 1. Ionic composition of irrigation waters applied to the soil monoliths and their classification according to
US Salinity Laboratory (Richards, 1954).
EC
SAR
Ca2+
Mg2+
Na+
ClUSSL
Monolith-1-1 0.5-1classification(meq L )
meq L
dS m
A
0.3
1.0
1.00
1.00
1.00
5.00
C2-S1
Waters
B
0.8
3.0
1.28
2.56
4.16
11.82
C3-S1
I
C
1.6
6.0
1.93
3.86
10.21
21.79
C3-S2
A’
0.8
1.5
1.85
3.65
2.5
13.5
C3-S1
Waters
B’
1.6
3.0
3.16
6.32
6.52
25.48
C3-S1
II
C’
3.2
6.0
5.1
10.2
16.7
47.3
C4-S2

Soil solutions were collected from ceramic cups and soil water contents were monitored
using TDRs at four depths twice a week during the irrigation periods and once a week during
remaining months. At the beginning of the experiment, at the end of each irrigation period,
and after the following winters, soil samples were collected at 5 depths (0-20, 20-40, 40-60,60-80, 80-100 cm) to measure the exchangeable cations, and the cation exchange capacity(CEC).
Evapotranspiration, soil hydraulic properties, solute transport parameters and many
physical and chemical analysis like soil texture, dry bulk density, concentrations of soluble
cations (Na+, Ca2+, Mg2+), electrical conductivity (ECe), and cation exchange capacity (CEC)
were determined during the extend of the field experiment, in order to satisfy the considerable
input data needed to evaluate the effectiveness of the HYDRUS-1D software package[Šimůnek et al., 2005].
RESULTS AND CONCLUSIONSThe use of irrigation waters with EC up to 1.6 dS m-1 during the time period from May
2001 to September 2004 did not lead to salinisation/alkalization of the medium textured
Fluvisol due to its favourable hydraulic characteristics and the winter precipitations. The
rainfall water leached the salts accumulated during the irrigation period down to a depth ofl50 cm. Irrigation with water having EC of 3.2 dS m-1 can cause significant soil salinisation
and alkalisation. After two irrigation cycles (in 2003 and 2004) EC of the soil solution
increased to about 12 dS m-1, SAR to about 8 (meq L-1)0.5, and ESP to 17% in the surface soil
horizon (0-20 cm). Although the winter rainfall leached most salts from the surface layers,
this was not sufficient to restore the soil to its initial conditions bellow the depth of 40 cm. EC
values reached 3 dS m-1, SAR 7 (meq L-1)0.5 and ESP 9% at a depth of 50 cm in March 2004,
immediately after the third rainfall season.
HYDRUS-1D has successfully simulated water regime, as well as the effects ofdifferent irrigation waters on the geochemistry of the studied Fluvissol (Figures 1 to 3). The
correspondence between observed and simulated variables is remarkable, considering thatsimulations were carried out to predict field measurements over a considerable time period (4years) without any calibration, with all input variables (i.e., soil hydraulic properties, solutetransport parameters, atmospheric demand, Gapon constants, physical and chemical
characteristics of the soil, LAI and root depth) measured independently. The agreement39between measure and simulated values was the best for soluble sodium concentrations(R2=0.7980; n=1091 observations), SAR (R2=0.8645 n=1180) and ESP (R2=0.7338 n=105),
while water contents (R2=0.6033; n=1180), ECs (R2=0.6505 n=1043), and soluble magnesium(R2=0.6323 n=1056) and calcium concentrations (R2=0.6219 n=1015) were predicted slightly
less well.
Our study shows that HYDRUS-1D can become a very powerful tool in irrigation
management, for predicting the effects of irrigation water quality on the soil and groundwater
quality, and for implementing better irrigation and fertilization practices [e.g., Gärdenäs etal., 2005], which is particularly important when designing new irrigation areas. Models suchas HYDRUS-1D, after their proper calibration and validation, should be used to establish
sound irrigation policies and to mitigate environmental risks.
A complete description of this study can be found in Gonçalves et al. (2005).

-1

EC (dS m )

20

10 cm

15
10
5
0
0

200

400

600

800

1000

1200

1400

Na Concentration [meq/L]

Figure 1 – Measured and simulated soil solution electrical conductivities at 10 cm depth in the threemonoliths A, B, and C.
50

50 cm

40
30
20
10
0
0

200

400

600

800

1000

1200

1400

Figure 2 – Measured and simulated soluble sodium at 50 cm depth in the three monoliths A, B, and C.

40

10

Measured A
Simulated A
Measured B
Simulated B
Measured C
Simulated C

70 cm

SAR (meq L )-1 0.5

8
I

6

R

I

200

400

R

I

R

I

1000

1200

4
2
0
0

600
800
Time (d)1400

Figure 3 – Sodium adsorption ratios measured and simulated at 70 cm depth in the three monoliths A,
B, and C. I and R correspond to the irrigation and rainfall periods, respectively.

REFERENCES
Gärdenäs, A., Hopmans, J. W., B. R. Hanson, and J. Šimůnek, 2005. Two-dimensionalmodeling of nitrate leaching for various fertigation scenarios under micro-irrigation,
Agric. Water Management, 74, 219-242.
Gonçalves, M. C., J. Šimůnek; T. B. Ramos, J. C. Martins, M. J. Neves and F. P. Pires, 2005.
Multicomponent solute transport in soil lysimeters irrigated with different quality waters
Water Resour. Res., (submitted July, 2005).
Richards, L. A., 1954. Diagnosis and Improvement of Saline and Alkali Soils, USDA
Handbook 60, Washington, USA.
Š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 invariably-saturated media. Version 3.0, HYDRUS Software Series 1, Department ofEnvironmental Sciences, University of California Riverside, Riverside, CA, 270pp.
White, N. L., and L. M. Zelazny, 1986. Charge properties in soil colloids. in Soil Physical
Chemistry, edited by D. L. Sparks , CRC Press, Boca Raton, Florida.

41                                                          

 ! " #  $      " % & ' (   )  * +# # " % # $   ,    - " % " "  . "  

 / 0 1 2 3 4 5 6 1 6 7 8 3 9 :9 ; 9 6 < 7 = > ? 7 @ A B C 2 1 3 D 5 ; 8 0 3 E 6 8 9 1 ; F G H I :D 5 > F J K B L M N M L G H I : D 5 B > 6 1 F O 8 P, = 6 8 9 1 O I 7 8 3 9 :9 ; 9 6 Q 2 1 R I6 D 9 1 2 8 : D 3 < S R T A B C 2 1 3 D 5 ; 8 0 3 E 6 8 9 1 ; F G H I :D 5 > F J K B L M N M L G H I :D 5 B > 6 1 F O 8 P
U V W X Y Z [ \ W ]Y V
^  " " _"  ( % +  $ _  `  a b  ( +`  $  a c ` _$ % +' $ ( +`  c % ` c " % ( +" ' ` d ' ` +_' +  _ b a " ' ( % b  ( b % " $  a' ( $ ( " + d` % # $ ( +`  a " c "  a "  ( `  e ` f $  a ( % $  ' c ` % ( c $ % $ # " ( " % '  ) _"  ( % + $ _ c ` _$ % +' $ ( +` c  "  ` # "  $  # " $ ' b % " a f +(  (  " g c "  ( % $ _ h  a b  " a i ` _$ % +' $ ( +`  j g h i k # " (  ` a  $ % "  $ b ' " a& l (  " " _ "  ( % `  a +m b ' + n " +  ( " % $  ( +`  ` d (  " ' ` + _ # $ ( % +o f + (  (  " c ` % " e b + a  ^  " c ` _ $ % + ' $ ( +` a " c "  a ' " p  `  (  " c ` % " ' + q "  `  (  " " _ "  ( % +  $ _  `  a b  ( + n + ( l ` d (  " c ` % " e b + a $ ' f " _ _$ ' `  ' $ ( b % $ ( + `  j " p  r s t u v w x y z {  | } } ~    t z w  w x y z {  | } } ~ k   " % " f " c % " ' "  (% " ' b _( ' d% ` # g h i # " $ ' b % " # "  ( ' f  +  f " % "  `  a b  ( " a +  `   b   ( +`  f +(  # b _( + ' ( " c` b ( e ` f j g  k " o c " % +# "  ( ' ( ` " o $ # + " % " _$ ( +`  '  +c ' & " ( f " "  ' ` +_  l a % $ b _+ c % ` c " % ( +" '$  a " _"  ( % + $ _   $ % $  ( " % +' ( + '    X ]  V W     W [ h  ` % a " % ( `  ` + ( _l # " $ ' b % " g h i   $ % $  ( " % +' ( + ' $  a  l a % $ b _+ c % ` c " % ( +" ' f " " ' ( $ & _+'  " a$  " f _$ & ` % $ ( ` % l # " $ ' b % " # "  ( ' " ( b c j +p   k  ^ f ` " _"  ( % ` a " ' jc ` % ` b ' & % `  q " a + '  ' k f " % "b ' " a ( ` +   "  ( $  " _"  ( % +  b % % "  ( + (  " d% "  b "   l % $  p "     ( `      + ( ` (  " ' c "  +# "   ^  " a +m " % "   " ` d (  " " _"  ( % + c ` ( "  ( +$ _ & " ( f " "  ( f ` d b % (  " % " _"  ( % ` a " ' f $ ' # " $ ' b % " a ^  " # $ +  ( $ % p " ( n $ _ b " f $ ' (  " c  $ ' " '  + d (        & " ( f " "   b % % "  ( $  a n ` _ ( $ p " f  +  
a " c "  a ' `  (  " $ c c _+" a d% "  b "   l j$ ' + b ' ` +a $ _ f $ n " d` % # f $ ' b ' " a k   ` % (  " n $ % +$ ( +`  ` d
f $ ( " % ' $ ( b % $ ( +`  $  a db % (  " % (  " a " ( " % # +  $ ( +`  ` d (  " b  ' $ ( b % $ ( " a  l a % $ b _+ c % ` c " % ( +" '$  g  d $  +_ +( l f $ ' b ' " a  ^  " $ a n $  ( $ p " ` d (  +' $ p $ +  ' ( " n $ c ` % $ ( +n " a % l +  p ( ` $ a  b ' ( (  "f $ ( " %  `  ( "  ( +' (  $ ( $  $ _ # ` ' (  `  ' ( $  ( n $ _b " ` d (  " " _"  ( % + $ _  `  a b  ( +n +( l    ¡   ` d(  " c ` % " e b +a  $  & " # $ + ( $ + " a  c % " ' ' b % "  " $ a f $ ' $ c c _+" a $ ( (  " ( ` c ` d (  "  ` _b #  $  a (  " ` b ( e ` f + ( ` $ & b % " ( ( " f $ ' # " $ ' b % " a  d( " % % " $   + p $  "  b +_+& % +b # (  " c % " ' ' b % " +p b % "   ) o c " % +# "  ( $ _ ' " ( b c '  ` f + p # " $ ' b % " # "  ( " _"  ( % `  + '   ` _ b #   $  a & b % " ( ( " 42 +p b % " |  ) o $ # c _" ` d  g  a $ ( $ '  ` f + p $ c c _+" a c % " ' ' b % " ¢  % " ' b _( $  ( ` b ( e ` f  $  a £ ¤ ¥ ¦ g¤ § (f $ ' +  % " $ ' " a  h ( +' c ` ' ' + & _" ( ` +  n " % ( d% ` # (  "  ` % % " ' c `  a + p a $ ( $ (  "  l a % $ b _ + c % ` c " % ( +" 'f +(  (  "  £ ¤ ¥ ¦ g   ¤ c % ` p % $ # ` d ¨  © ª t w « w x y z  j| } } ~ k   ` % (  +' (  " c % " ' ' b % " ¬ ­® ¯ $ ((  " ( ` c ` d (  "  ` _ b #   $ a ( ` & " ( % $  ' d ` % # " a +  ( ` ' b  ( + `  $ ( + ( ' & ` ( ( ` #  ¢ ° ® ­ ­ ® ± ²   ³´ µ ¬ ­ ® ¯ ¶    ·  ¸  ¹ º»  ¼ ¿½ ¾À Á ´ Â Ã Ä ­ Å Æ Ç ® È É ± Ê Ë Ê ­ Ì Å ° É Æ Å ­ ­ Å ²   Í Î ¢ ° Æ ® Ê Ï Å ²   Î ¢ Ç Å Æ Ä ± Ë Ç Ð Ë Ñ Ò ²  Î ¢ Ñ ® Ë È ²    ^  " ` b ( e ` f f $ ' # " $ ' b % " a $ ' c % " ' ' b % " ` d (  " f $ ( " %  ` _ b #  +  (  " & b % " ( ( "$  a  `  n " % ( " a +  ( ` f $ ( " %  " + p  ( Â ²   $  a n ` _ b # " Ó  Ô    ` % (  " b ' " ` d  £ ¤ ¥ ¦ g   ¤Ó  $ ' ( ` & "  `  n " % ( " a +  ( ` $  ` % # $ _ + ' " a ` b ( e ` f Õ ²  Ö ²  × Ø  ³ Ó ·  ½ ¾ Ù · Û × Ü  f  " % "Û ³ Ý Þ Ø a "  ` ( " ' (  " $ % " $ $  a Þ ²   (  " % $ a + b ' ` d (  " & ` ( ( ` # ` d (  "  ` _ b # ¾ Ú   + p  | '  ` f '(  " # " $ ' b % " a ` b ( e ` f $ ' f " _ _ $ ' (  " $ c c _+" a c % " ' ' b % "  " $ a d` % $  " o " # c _$ % l ' $ # c _" 
ß  W  à V  á â ]â^  " c  $ ' " ' c "  ( % b # f $ ' $  $ _l ' " a & l § ( ( + p $  ` n " % _$ l + p ã ` _"  ã ` _" # ` a " _ j ä å z w æä å z w   ç è   é w © t y  | } } } k (` (  " +# c " a $   " a $ ( $  f  " %" +# c " a $   "ê ëì í ³ î ï ð ´ ñÇÇô õ  ´ôùú Î ëöì ÷ ò í øò óû  

$  a  ëì í ³   ² ü  ý þ¾ ÿ ÿ  ^  " # b _( + ã ` _"  ã ` _" § ( l +" _a ' (  " c $ % $ # " ( " % ' % " _$ o $ ( +`  ( +# "ÿ ÿ÷      $ % p " $ & + _ + ( l Ç Ç ´   $  a ã ` _ "  ã ` _ " " o c `  "  ( ´  d ` % " $   +  a +n + a b $ _ a + ' c " % ' +` 
( " % # $ ' f " _ _ $ ' ¤ ã %ó" ' + ' ( + n + ( l î ï û      ` % a +  p ( ` (  " ` % l  (  " % " _ $ o $ ( + `  ( + # " ÷ + '  `   "  ( " a ( ` $   $ % $  ( " % + ' ( +  ' c $ ( + $ _ a + # "  ' + `      " p  ³  ÷ j  x å w x y z {  | } } | kf +(  a +m b ' +`   ` "  +"  ( Ø × Ü    ` % " o $ # c _ " (  " " m "  ( + n "  l a % $ b _ +  % $ a + b ' Þ  + '_ + . " a ( ` (  "  l a % $ b _+ c " % # " $ & +_ +( lÞ Ø & l ( " ! ` q "  l ã $ %# $  # ` a " _ jä y © y t  ç ~ k  ^  " +  n " % ' " ` d (  " ã ` _"  ã ` _" " o c `  "  ( +' $ # " $ ' b % " ` d (  " e $ (  " ' ' ` d (  " c  $ ' "' c "  ( % b #  $  a f "  `  ' +a " % × Ü $ ' $  + a + $ ( ` % d` % (  " c ` % " ' +q " a +' ( % +& b ( +`  + a " o  ^  "ã ` _"  ã ` _" c $ % $ # " ( " % Ç Ç +' $ # " $ ' b % " ` d (  " " _"  ( % +  $ _ c ` _ $ % +' $ & + _ +( l  f  +  a " c "  a ' ` 
ó
  

 

  













 







 









|
43

 +p b % " ¥ " ( "  ( +`  d b   ( +`  ` d `  " ` d (  " b ' " a & % `  q " " _"  ( % ` a " ' ' c "  + §  ' b % d $  "   _ $ l  `  ( "  (  $  a  $ ( + `  " o   $  p "  $ c $  + ( l j " p  tw¨uv t  ç ç  k   b % (  " % + ( + ' î ï× Ü
j u v w   ç è | k  f  " %"a "  ` ( " ' (  " " m "  ( +n "' $ ( b % $ ( +`  $  a (  " ' $ ( b % $ ( +`  " o c `  "  ( ^  "  £ ¤ ¥ ¦ g  ø ¤ c % ` p % $ # $ _ _ ` f ' ( ` +  n " % ( (  "  g  a $ ( $ j ( + # "    $ c c _ + " a c % " ' ' b % " " $ a ¢  ` b ( e ` f Õ k d` % (  " c $ % $ # " ( " % '   $ % $  ( " % +q + p (  " ' ` +_  l a % $ b _ + c % ` c " % ( +" ' j  +p  | k × $  a db % (  " % +( +' j y t w t ª u v x w t ^  " f $ ( " %  `  ( "  ( +' p +n "  $ ' ë ¢ í ³ÎÿÜç } k ë ¢ í × Ü ë Î ë ¢ í í× Øë¢ í ³· ´ Î ë ¢ í 
f  "%"$aa "  ` ( " (  " % " ' + a b $ _ $  a (  " ' $ ( b % $ ( " a f $ ( " %  `  ( "  (  % " ' c "  ( + n " _ l  ¢ ²  (  " c % " ' ' b % "  " $ a  ë ¢ í ²  ¡   (  " b  ' $ ( b % $ ( " a $  a²  ¡   (  " ' $ ( b % $ ( " a  l a % $ b _ +  `  a b  ( +n +(l ²  × Ü  + ' (  " +  n " % ' " ` d (  " $ + %  "  ( % l n $ _ b " $ ' f " __ $ '$ % " d` % #c $ % $ # " ( " % ' % " _ $ ( " a ( ` (  " c ` % " ' + q " a + ' ( % + & b ( + `   $  a + ' (  " ( ` % ( b ` ' + ( l  óh  ` % a " % ( `%" a b " ( "  b # & " % ` d b  .  ` f  c $ %$ # " (" %' f " '" (³  ´ ¡  ³ ¡ $a³ ¸ !ó  ` f +  p  $aj` % k  ` f_" $ a + p ( `ë¢ í ³ ë¢ í ¡ ³ ë Î ë ¢ í í$ __` f ' ( `  $ _ b _$ ( "ë ¢ í   ` % $ _ _ # $ ( " % + $ _ ' j ' ` + _  & % `  q " a + '    " % $ # +  a + '  k ó(  " ' $ ( b % $ ( " af $("% ` ("($  a  l a % $ b _+  `  a b  ( +n +( lf " % " # " $ ' b % " a  ^  " % " # $ + + p  l a % $ b _+
c % ` c " % ( +" ' ` d (  "  " % $ # + a +'  f " % " ( $ . "  d% ` #w
x w j| } } ~ k  ` % ( " & %`  q "a +'  ' f "  ` # c b ( " a $  a d% ` # # " $ ' b % " a % " ( "  ( +`  n $ _b " ' f +(  ¥ " ( ã j y t w t ª u v x w t
w x y z {  ç ç  k d% ` # (  "  £ ¤ ¥ ¦ g c $  . $ p " j +p  k  ^  "  £ ¤ ¥ ¦ g   ¤ § ( l +" _a ' (  " c ` % "' +q " a +' ( % + & b ( +`  + a " o$  a (  " & b & & _+ p c % " ' ' b % " × Ü d% ` # f  +  (  " b  ' $ ( b % $ ( " a l a % $ b _+  `  a b  ( +n +(l  $  & "  $ _ b _$ ( " a       



 ! "











%

# $

&

' "

'

( 

( )



/

* +, # +-  .

0

4

1

1

2

5

5



4

4

' "

5

5



6

7

3

5

5



6 8

1

2

' 2

1

4









' 2

'

4

5

5



9 : -

1

2

' "







1

' 2

1

2

;

4

<

 =



  >



/





?

4

 â [ W â

^  " $  $ _l ' +' ` b ( _+ " a $ & ` n " # $ a " +( c ` ' ' +& _" ( ` " o $ # + " (  " % " _$ ( +`  & " ( f " "  (  "  l a % $ b _+  c % ` c " % ( +" ' $  a (  "   $ % $  ( " % +' ( + " _"  ( % + $ _ c $ % $ # " ( " % '   +p  è '  ` f ' (  " % " _$ ( +` 
& "(f ""   ÷  $  ad` % (  " + n " ' ( +p $ ( " a ' $ # c _" ' f +(  a +m " % "  ( p % $ +  ' +q " '  ^  " % " _$ o $ ( +`  ( + # " ' d` % (  " a +m " % "  ( p % $ + ' +q " ' a +m " % ' l ' ( " # $ ( +  $ __l ` n " % (  "  `  ' +a " % " a ' $ ( b % $ ( +` 
% $  p "   ` % $ n $ _ b " ` d " m "  ( +n " ' $ ( b % $ ( +`  +( +' c ` ' ' +& _" ( ` a " ( " % # + " $ % " _$ ( +`  '  +c & " ( f " "  l a % $ b _+  `  a b  ( +n +( l $  a % " _$ o $ ( +`  ( + # "  (  "  _` ' " % + ' c "  ( +`  ` d f  +  f + __ & " ' b &  "  (` d db ( b % " f ` % . 
1





44

 +p b % " è  ¥ " _ $ ( +`  '  +c & " ( f " "   l a % $ b _+  `  a b  ( +n +( l  % " _$ o $ ( +`  ( +# " ÷  $  a" m "  ( +n " ' $ ( b % $ ( +` d` % ' $  a ' $ # c _" ' f +(  a +m " % "  ( p % $ + ' +q " '  X V \â
u v  w    j  ç è | k  ^  " " _"  ( % + $ _ % " ' + ' ( +n +( l _` p $ ' $  $ +a + a " ( " % # + +  p ' ` # " % " ' " % n ` +%   $ % $  ( " % +' ( + '  ~è | t z w   { ¨ z y x w  { { ª « w  r {ä y  y t   { j | } } ~ k  ^  " % " _$ ( +`  '  +c& " ( f " "  ' c "  ( % $ _ +  a b  " a c ` _$ % +q $ ( +`  $  a  l a % $ b _+ c % ` c " % ( +" ' ` d ' $ ( b % $ ( " a $  a b  '$ (b %$(" a '$  a '(`  " 
 + c % " ' ' t w  { { ¨ u v t  { { j  ç ç  k  % " _$ ( +`  & " ( f " "  (  "  b $ a % $ ( b % "  ` # c `  "  ( ` d" _"  ( % + $ _  `  a b  ( +n +( l $  a (  " ' c "  +§  ' b % d $  " $ % " $ ` d ' " a + # "  ( $ % l % `  . '   |  ä y © y t   ä  j  ç ~ k   _` f ` d p $ ' " ' (  % ` b p  c ` % ` b ' # " a +$  b ( ( " % f ` % (  '  `  a `  ä å z w  !  g   æ ä å z w  ¥    j  ç è  k  ¤ +' c " % ' +`  $  a $ & ' ` % c ( +`  + a +" _"  ( % +  '   è ~é w © t y   j | } } } k  ^ ` # ` p % $ c  + +  n " % ' +`  ` d  ` # c _" o % " ' + ' ( +n +( l  (  " ` % l $  a $ c c _ + $ ( +`   i   ¤  (  " ' +'  `   b # ¦  +n " % ' +(l  h g
 ç è
ç|ç
r s t u v  { r { é w © t y  { x å  é { © © w © y t t  {w wwu«wt {j | } } ~ k  ¤ " c "  a "   " ` d (  " ' c "  ( % $ _ + a b  " a c ` _ $ % +q $ ( +`  % " ' c `  ' " ` d ' $  a ' `  ' $ _+ +( l 
p % $ + ' +q " $  a ' $ ( b % $ ( +`    }| è~w
x w  ä v  j| } } ~ k  c " % ' `  $ _  ` # # b  + $ ( +`   ' " " $ _' `w
xw  ävj| } } è k
& " % c %  db  p a " '  b $ _" #   % +" a # $    ` a " __' a b %    " ' ' b  p a " '$ ''" %p "  $ _( ' b  a a " % " _" . ( % +'   "  " +( d  +p . " +( + ` a "  ' b _"  # +( ( " _ ' ^ ¤ ¥ b  ( " % ( % $  ' +"  ( "  " a + p b  p "   ¤ + c _` # $ % & " +(  ¦  +n  `   ¨  © ª t w «    y t w t ª u v x w t  r  v   æ ¨ w t y  r  j | } } ~ k   £ ¤ ¥ ¦ g   ¤  - } } ã ` a " d` % ' +# b _$ ( + p (  " `  "  a +# "  ' +`  $ _ # ` n " # "  ( ` d f $ ( " %   " $ ( $  a # b _( +c _"' ` _b ( " ' + n $ % +$ & _l  ' $ ( b % $ ( " a c ` % ` b ' # " a +$  ¥ +n " % ' +a " j¦ g k 
y t w t ª u v x w t  r  v  jç } k _` ' " a  d` % # "  b $ ( +`  d` % c % " a + ( + p (  "  l a % $ b _+ `  a b  ( +n +( l ` d b  ' $ ( b % $ ( " a ' ` +_'   ç| ç 
y t w t ª u v x w t  r  v  ¨  © ª t w «   w     ¨ w t y  r  j  ç ç  k  ¥ " ( ã  - } 
ã ` a " d` %  b $  ( +dl + p (  "  l a % $ b _+ d b   ( +`  ' ` d b  ' $ ( b % $ ( " a ' ` +_ '  ¥ +n " % ' +a " j ¦ g k 
x å  é  é å © y å  { y y å  {w x «   { j| } } | k  ^  " ` %" ( + $ _ $  a" o c " % +# "  ( $ _ ' ( b a l ` d ( + # " a ` # $ +  + a b  " a c ` _ $ % +q $ ( +`  + f $ ( " %  ' $ ( b % $ ( " a ' $  a '  è è 
1



?

@



/

A



 



Z

Y

E D bc F K

d e

B C D E F G

H

X

Y

Z

[

D K O C

\

]

^



I

G

J E F K GL

M E GN

>

O F G

L

O K G



\

P

Q

O KC G

> >

R

E S G

T U V

W 

/

O F G

B _ O

H





W 

f



g

h

i G

Q _ c F G

k



g

m

<

 =

  >

l





n

p O a q _ c F G

;

\

O F G

H

  



r F KC G

A

s



o



;

^

<

 =   >



u

h

v

g

g

v

g

/

w

/



0

x a M b/x y M G

]



i G





t

]

j _ O I

W 

g

^

` a S



H q q bM O z

p O a q _ c F G

o

{ |



} W

>

X

o

 

è
45

x a y G

w

Z



X

H

]

I

G



>

i G

U U

w

0

W 0

0



G

PARSWMS: a parallelized version of SWMS_3D
Javaux, M., Hardelauf, H., Gottschalk S., Herbst, M., Vanderborght, J. and H. Vereecken

1. Introduction
There is currently a need for large-scale models representing the water flow and solutetransport processes at large scale. Among the numerous approaches dealing with this problem,
distributed modeling based on Richards’equation is one of the most popular. Yet, a key issueof the distributed modeling is the maximum grid resolution at which the model is valid.
Theoretically, the grid resolution should be of the same order of magnitude as the Darcy (orelementary representative volume) scale, that is between 10-2 -100 m. Thus, when large-scale
problems have to be faced, the major drawback of fine grid models is the enormous demandfor computational time and resources (Harter and Hopmans, 2004). The speed and efficiency
of current models have therefore to be improved. The paralellization of the code is one
possible way to decrease the computational time by distributing a complex large geometry
problem over multiple processors working in parallel.

2. Code principles
The PARSWMS code is based on the SWMS_3D code (Simunek et al., 1995). Mostof the subroutines, functions and variable names have been kept the same whereas somechanges were made to parallelize the code. First the programming language was switchedfrom FORTRAN to C++ , which allows for dynamic allocation of all variables, i.e. there is noneed to redefine array sizes. Also the output was slightly modified. Basically, the code gives
one series of output files per processor. Therefore, a post-processing is needed to merge all
the data together after running. Though, note that the input files are the same as in
SMWS_3D.
The new code is based on MPI (Message-Passing Interface), a library specification formessage-passing, which provides source-code portability and allows efficient implementation
across a range of computer architectures. MPI is a free software for LINUX or UNIX, which

46

1needs to be installed on the working machine, and which provides several libraries for parallel
computing.
The partitioning of the problem is automatically performed based on the PARMETIS
library. With the aid of this library the partitioning process itself is parallelised, partitioning
could be made dynamically (adapted at each time step), can be carried out for irregular grids,
and the parallel code can be run on heterogeneous cluster of processors. Currently thepartitioning is static, in the sense that this is done once at the beginning of the run. The
algorithm optimizes the partitioning between processors based on weights given by the user todifferent subdomains, by minimizing the surface of the subvolumes on each processor as wellas the connection numbers between processors.
The solution of the system of linear equations is achieved for the nodes of the
subvolume allocated on each processor. The PETSc toolbox was used since it allows the user
to choose between a large range of solvers or preconditionners able to deal with linear andnon-linear problems. Therefore, the solver implemented in PARSWMS is different to
SWMS_3D, and could be further optimized.
The PARSWMS code can be compiled and run on any network (or even uniprocessor) of LINUX or UNIX workstations or on specialized parallel hardware as long as the
MPI, PETSc and PARMETIS freewares are installed.

3. Applications and results
A first release of the code is working and it performance is currently evaluated.
Several points are particularly investigated: (1) the code is submitted to a benchmarking
developed for testing numerical codes (Vanderborght et al., 2005), (ii) the code is comparedwith the current SWMS_3D code in terms of performance and results, and (iii) the
performance of the code in terms of computational time is evaluated, i.e., the time scaling ofseveral computationally-demanding problems is studied on a massive parallel system.
First results confirm that SWMS_3D and PARSWMS give convergent simulations formost cases, under several benchmarking scenarios for water and solute. Figure 1 shows thetime scaling for a steady-state solute transport scenario through a soil cylinder with 500.000
nodes. Tests were made on 1 to 128 processors and there was a linear relationship between thelog10 number of nodes and the log10 time reduction with a slope around 0,9.

47

2

3.5

3

Log10 of the time [min]log10(t)= − 0.93*log10(np) + 3.2

2.5

2

1.5

1

0.5

0

0.5

1.5
1
Number of processors (in log10)

22.5

Fig. 1.:Ttime scaling for solute transport scenario with 492264 nodes

4. References
Harter, T and J. W. Hopmans. 2004. Role of vadose-zone flow processes in regional-scale hydrology: review,
opportunities and challenges. In Unsaturated-zone modeling. Progress, challenges and applications,edited by Feddes, R. A., G. H. De Rooij, and J. C. Van Dam Wageningen Frontis Series. Dordrecht,
The Netherlands: Kluwer Academic Publisher.
Simunek, J., K. Huang, and M. Th. Van Genuchten. 1995. The SWMS_3D code for simulating water flow andsolute transport in three-dimensional variably-saturated media. Version 1.0. Research Report n°139.
Riverside, California. U. S. Salinity Laboratory. Agricultural Research Service. U.S. Department ofAgriculture.
Vanderborght, J., R. Kasteel, M. Herbst, M. Javaux, D. Thiery, M. Vanclooster, C. Mouvet, and H. Vereecken.
2005. A Test of Numerical Models for Simulating Flow and Transport in Soils Using Analytical
Solutions. Vadose zone journal. 4:206-221.

48

3

Non-equilibrium Isoproturon Transport in Structured Soil Columns:
Experiments and Model Analysis
J. Maximilian Köhne*1, Sigrid Köhne1, and Jirka Šimůnek2
1

Institute for Land Use, Faculty for Land Use and Agricultural Sciences, University Rostock,
Justus-von-Liebig-Weg 6, 18044 Rostock, Germany, max.koehne@uni-rostock.de,sigrid.koehne@uni-rostock.de; 2 University of California, Riverside, Dept. of Environmental
Sciences, 900 University Avenue, A135 Bourns Hall, Riverside, CA
92521,*
jiri.simunek@ucr.edu; Presenter
Introduction
Model predictions of pesticide transport in structured soils are complicated by multiple
processes acting concurrently. In this study, the hydraulic, physical, and chemical
nonequilibrium (HNE, PNE, and CNE) processes governing isoproturon transport undervariably-saturated flow conditions were analyzed. Objectives were (i) to investigate theapplicability of model parameters obtained from inverse analysis of transient water flow andBr- transport, and from equilibrium or kinetic batch isoproturon sorption experiments, for theforward simulation of isoproturon transport during variably-saturated flow, and (ii) to detect
potential differences in isoproturon sorption for the PFPs and the matrix domains.
Materials and Methods
An undisturbed soil column (14.7 cm diameter, 15 cm long) was extracted from a loamy
Ap horizon (2% clay, 44% silt, 54% sand, 1.8% organic carbon, bulk density 1.4 g cm-3) withsubangular aggregate structure (Meyer-Windel, 1998). Bromide (Br-) and isoproturon (IPU,
3-(4-isoprpylphenyl)-1,1-dimethylurea) were applied to the soil column which was thenirrigated at a rate of 1 cm h-1 for 3 h. Two more irrigations at the same rate and duration
followed in weekly intervals. Nonlinear (Freundlich) equilibrium and two-site kinetic
sorption parameters were determined for IPU using batch experiments. The observed waterflow and Br - transport were inversely simulated using mobile-immobile model (MIM), dualpermeability model with equilibrium sorption (DPM) or with kinetic sorption (DP-KM), andcombined triple-porosity model (DP-MIM), all implemented in HYDRUS-1D (Šimůnek etal., 1998, 2003, Pot et al., 2005, Köhne et al., 2005).
Results
The MIM, DPM, and DP-MIM inverse simulation results for water flow and Brtransport in the Ap soil column are displayed in Figure 1. Observed bulk soil water contents
varying between 0.3 and 0.36 (Fig. 1a), pressure heads between -60 and 0 cm (Fig. 1b), andcumulative outflow (Fig. 1c) for three intermittent irrigation events were more or lessmatched by all three models. Simulated Br- BTCs were within fluctuations of measured
concentrations up to 7.5 days (Fig. 1d). Both the early but low Br- peak of 5.5 mg L-1,
corresponding to only 0.5% of the applied concentration of 1000 mgL-1, and the shape of the49Br- BTC were delineated by all model approaches. The late-time low Br- concentrations
observed after two weeks could be most accurately captured using the DP-MIM (Fig. 1d).
The assumption of the two pore domains was a suitable concept to describe Br- leaching at
HNE and PNE, at least during the first two irrigation events. The accurate match of Brconcentrations after two weeks using the DP-MIM suggests that the late-time Br- transport
was dominated by diffusive transfer into the microporosity subdomain of the matrix that
acted as a sink for Br-.
For isoproturon, Figure 2 shows the measured BTC together with simulation results
using various forward (Fig. 2a) and inverse (Fig. 2b,c,d) modeling approaches. The IPU
breakthrough displayed a somewhat similar shape as the Br- breakthrough, but at much lower
concentrations below 0.2 mg L-1 (0.08% of the applied 254 mgL-1). The simultaneous start ofthe BTCs for Br - and IPU at relatively low peak concentrations suggested (weak) preferential
flow conditions. The lower peak concentrations for IPU than for Br - can be explained withsorption of IPU. The MIM forward approach predicted zero concentrations for all times andis not shown in Fig. 2a. While none of the forward modeling approaches closely
approximated the observed IPU breakthrough, the DP-KM was the only approach thatpredicted an early concentration increase (Fig. 2a). Fitting the fraction of sorption sites inequilibrium with the mobile, fmo (MIM), or fracture, ff (DPM), domains improved results forthe two-domain equilibrium sorption approaches and enabled them to almost match the early
IPU peak (Fig. 2b). Note that we refer to DPM with optimized parameter ff as DPMf in Fig.
2b and that we use similar notation for other inverse optimizations. IPU concentrations at one
week were overestimated with DPM, but were matched well when the first-order degradation
constant, µ, was additionally fitted (Fig. 2b). MIM slightly underestimated IPU
concentrations after one week, which could not be compensated by additionally fitting µ (notshown). The match of the early IPU peak was further improved using inverse DP-KM
approaches assuming HNE, PNE, and CNE (Fig. 2c). Predictions of the measured early IPU
concentration increase (Fig. 2c) improved in the same order. The DP-MIM approach with theinverse estimation of ff failed to describe IPU concentrations, but approximations improved
when fitting ff and µ; or ff, µ and β (Fig. 2d). Furthermore we note that we obtained similarvalues for the fraction of sorption sites fmo (= 0.056) for MIM and corresponding global
values ff wf (=0.04 to 0.06) for DPM and DP-MIM, which can be explained with thecomparable mathematical structure of the two models. The simultaneous inverse estimationof ff, ωc,f, and µ in the DP-KMfωµ approach decreased ff from 0.375 (batch test value for thebulk soil) to 0.187, increased the sorption rate coefficient in the fracture pore system, ωc,f,from 0.105 (batch test) to 0.756 d-1, and increased µ from 0.046 to 0.32 d-1, corresponding toa short half-life of 2.2 days. The total IPU mass loss after 14.2 days was slightly
underestimated with the DP-KMfωµ approach. While the observed mass loss was 0.090 mg
or 7.1 % of 1.27 mg IPU applied, the simulated mass loss was 0.074 mg or 5.8 %. Although
the agreement between DP-KM and observed IPU concentrations was further improved byfitting the Freundlich isotherm exponent β together with ff, ωc, and µ, resulting parametervalues appeared to be unrealistically high. Overall, the DP-KMfωµ approach was deemed togive the best representation of the IPU breakthrough.

50

Cumulative Pressure head Water content
3-3(cm cm )(cm)
outflow (cm)Br (mgL )

8
6
4
2
0

-

-1

0.40
0.35
0.30
0.25

M IM
D P M -M IM

D ata
DPM

0

5
Data
M IM
DP M
DP M -M IM

3 cm :

20
0
-20
-40
-60

10

0

a)

15

20

Data
M IM
DP M
DP M -M IM

13 c m :

5

b)

10

15

20

8
6
4
2
0

c)

0

0.2

0.4 7.1

7.3

D ata

0

0.1

0.2

7.2

0.3

7.5

14

14.2

M IM

7.3

7.4

14.4

d)
D P M -M IM

DPM

14

14.1

Figure 1. Water flow and bromide (Br-) transport data and simulation results using the mobile-immobile model(MIM), dual-permeability model (DPM), and the dual-permeability triple-porosity model (DP-MIM), a) watercontents, b) pressure heads at 3 and 13 cm depths, c) cumulative outflow, and d) Br- effluent concentrations.
0.08

Data
DPM
DPM-KIMa) forward

Isoproturon concentration (mg L-1)

0.04
0.00
0

0.1

0.2

0.3 7.22

7.32

0.08

14

14.1
TER
MIM f
DPM f
DPM f µb) inverse HNE, PNE

0.04
0.00
0

0.1

0.2

0.3 7.22

0.08

7.32

14

14.1
Data
DPM-KIM f
DPM-KIM f ω
DPM-KIM f ω µ
DPM-KIM f ω µ βc) inverse HNE, PNE, CNE

0.04
0.00
0

0.1

0.2

0.3 7.22

0.08

7.32

14

d) inverse HNE, PNE
triple porosity

0.04

14.1
Data
DPM-MIM f
DPM-MIM f µ
DPM-MIM f µ β

0.00
0

0.1

0.2

0.3 7.22

7.32

14

14.1

Figure 2. Isoproturon (IPU) transport data and simulation results using the mobile-immobile model (MIM),
dual-permeability model (DPM), dual-permeability kinetic sorption model (DP-KM), and the dual-permeability
triple-porosity model (DP-MIM), a) forward simulations, and b), c, d) different inverse simulations.

51

Conclusions
The improved DP-KM predictions of the isoproturon concentrations using the batch twosite kinetic sorption parameters (as compared to using batch equilibrium sorption parameters)
indicated that kinetic sorption was a relevant process for preferential isoproturon leaching in
the Ap soil column. The inverse DP-KM approach revealed that kinetic sorption was
different in the matrix and fracture pore system, which cannot be derived from batch tests on
the bulk soil. A good description of the isoproturon BTC was obtained with the inverse DPKM approach when fitting three parameters: the fraction of equilibrium sorption sites in thefracture pore region, ff, the first-order kinetic sorption coefficient in the fracture pore region,ωc,f, and the first-order decay constant in the liquid phase of both pore regions, µ. A smaller
magnitude of ff as compared to values derived from batch tests suggests that sorption in thefracture pore region was time-dependent, as compared to being largely at equilibrium in thematrix. The larger magnitude of ωc,f indicated that kinetic sorption in the fracture pore system
was relatively fast. Large decay rates, µ, for the bulk soil are in accordance with results ofprevious studies in laboratory column transport experiments. However, as was suggested
already by Pot et al. (2005), high degradation rates may partly reflect irreversible sorption.
Overall, modeling our data was successful and permitted identification of different processes
involved in multi-process isoproturon transport under variably saturated flow conditions.
References
Meyer-Windel, S. (1998). Transport processes of bromide and isoproturons under steadystate and transient flow conditions. PhD-dissertation, Schriftenreihe Institut fürWasserwirtschaft und Landschaftsökologie der Christian-Albrechts-Universität zu Kiel,
28, 113pp.
Köhne, J.M., Köhne, S., and J. Šimůnek (2005). Multi-Process Herbicide Transport in
Structured Soil Columns: Experiments and Model Analysis, submitted to Journal ofContaminant Hydrology.
Pot, V., Šimůnek, J., Benoit, P., Coquet, Y., Yra, A., and M.-J. Matrinez Cordόn (2005).
Impact of rainfall intensity on the transport of two isoproturons in undisturbed grassed
filter strip soil cores, J. Contam. Hydrol., accepted.
Šimůnek, J., M. Šejna, and M. T. van Genuchten (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, IGWMC - TPS - 70, International Ground Water
Modeling Center, Colorado School of Mines, Golden, Colorado, 202pp.
Šimůnek, J., N.J. Jarvis, M.T. van Genuchten, and A. Gärdenäs (2003). Review andcomparison of models for describing nonequilibrium and preferential flow and transport in
the vadose zone. J. Hydrol. 272, 14-35.

52

Numerical Analysis of Soil Limiting Flow from Subsurface Sources
Naftali Lazarovitch1,2, Jiri Šimůnek2, Martinus Th. van Genuchten3 and Uri Shani11Dept. of Soil and Water Sciences, Faculty of Agriculture, Food and Environmental Quality Sciences,
Hebrew Univ. of Jerusalem, Rehovot, Israel. lazarovi@agri.huji.ac.il
2
Dept. of Environmental Sciences, Univ. of California, Riverside, CA USA.
3
George E. Brown, Jr. Salinity Laboratory, USDA, ARS, Riverside, CA, USA.

Abstract
The infiltration rate of water from a subsurface cavity is affected by manyfactors, including the pressure in the cavity, its size and geometry, and the hydraulicproperties of the surrounding soil. When a predetermined discharge of a subsurface
source (e.g. a subsurface emitter) is larger than the soil infiltration capacity, thepressure head in the source outlet increases and becomes positive. The increased
pressure may significantly reduce the source discharge rate. The main objective of
this study was to develop a system-dependent boundary condition that describes thisprocess while considering the source characteristics, and to implement this boundary
condition into HYDRUS-2D software package. The updated numerical model was
validated against transient experimental data. Good agreement was found between
transient cavity pressures measured in laboratory experiments and those calculatedusing the updated numerical model.
Introduction
When a predetermined discharge from a subsurface source (e.g. a subsurface
emitter) is larger than the soil infiltration capacity, the pressure head in the source
outlet increases and becomes positive (Shani et al., 1996). This pressure buildup in thesoil reduces the pressure difference across the source and, subsequently, decreases thesource discharge rate. Warrick and Shani (1996) showed that the reduction in thepressure difference is larger for soils having a lower hydraulic conductivity and thatthe positive pressure in the vicinity of the source increases rapidly at the beginning of
the infiltration event and approaches a final value after only several minutes. The
source discharge is also affected by the source characteristics and the cavity size ofthe source outlet (Shani et al., 1996). Analytical solutions for water flow from a
subsurface cavity are available for calculating spatial distributions of the pressure
head in the soil, ψ [L]. However, since these solutions were derived only for steadystate water flow conditions (Philip, 1992; Ben-Gal et al., 2004; Lazarovitch et al.,2005a,b), the models can be applied only to relatively long infiltration events. The
models are also restricted to Gardner’s (1958) hydraulic conductivity function, which
represents an important simplification of the dependence of the hydraulic conductivity
on the pressure head. Since numerical models can implement any model for the soil
hydraulic properties and are applicable to transient conditions, they can overcome
many of the restrictions imposed by analytical models. The objectives of this study
were to a) implement into a transient numerical model a system-dependent boundary
condition that describes infiltration from subsurface sources and considers the source
characteristic curves, b) evaluate the dependence of the transient source discharge rateon the soil and source hydraulic properties, c) carry out laboratory experiments
involving transient water flow from subsurface water sources, and d) compare theexperimental results with those obtained with the a modified version of the HYDRUS2D software package of Šimůnek et al. (1999).

53

Theory
Warrick and Shani (1996) suggested the following relationship between thesource discharge, Q [L3T-1], and the inlet pressure of the source (Pin [L]):c⎛ P - ψs ⎞Q = Q0 ⎜ in[1]⎟
⎝ Pin ⎠
where ψs [L] is the pressure head at the source-soil interface, often called the backpressure, Q0 [L3T-1] is the nominal discharge of the source for an inlet pressure Pin of
10 m and a back pressure equal to zero, and c [-] is an empirical constant that reflects
the flow characteristics within the emitter. Normally, c = 0.5 corresponds to a turbulent
flow emitter and c = 1 to a laminar one.
In our analysis we used a numerical solution of the Richards’ equation as
implemented in the HYDRUS-2D code for axisymmetric flow. In addition to existingboundary conditions that available in HYDRUS-2D code, Eq. [1] was implemented asa new system-dependent boundary condition. This boundary condition allows
calculation of the source discharge while considering properties of the source, and theinlet and back pressures:c⎛ P -ψ s ⎞∂ψ ∂ψ
⎞⎟⎟ = − ∫ K ⎛⎜Q = Q0 ⎜⎜ in+− δ iz ⎟ni dΓ
[2]
∂z
⎠
⎝ ∂r
⎝ Pin ⎠
Γ
where δiz is the Kronecker delta [-], Γ is the boundary of the source, and ni are thecomponents of the outward unit vector normal to the boundary. At each time stepduring the calculations, the discharge rate was adjusted while ψs was calculated byaveraging pressure heads (from the previous time step) around the subsurface source:1
ψ S = ⋅ ∫ψ S d Γ[3]Γ ΓFigure 1. Finite element grid and flow
domain. Water flows from the half
sphere on the left (zoomed). Boundary
condition [1] is specified on the spheresurface, while zero flux conditions areapplied to the other boundaries.

Application
An axisymmetrical three-dimensional flow
domain was selected as shown schematically in
Fig. 1. The flow domain (50 x 30 cm) was
discretized into 1183 triangular finite elementswith triangles significantly smaller around thesource and then smoothly increasing with distancefrom the source. Except for the free drainage
boundary condition at the bottom of the flow
domain and the source boundary condition [2] atthe source-soil interface, all boundaries were
subjected to zero flow conditions.
Laboratory drip irrigation experiments werecarried out with a non-compensated dripper packed
inside a perforated plastic sphere (the source) witha radius of 0.01 m. The source was connected to a
flexible tube and buried 0.25 m deep in a Magal
clay loam soil. Sources with nominal discharges Q0
of 4 and 8 L h-1 were used, while a value of 0.5 forc in [1] was assumed for all sources. Hereafter wedenote the source by its Q0 value.
54

The saturated hydraulic conductivity was estimated using the horizontal soilcolumn constant head method (Hillel, 1971), while the retention curve θ(ψ) was
determined using a ceramic plate suction cell apparatus (Klute, 1986). The Mualem-van
Genuchten soil hydraulic model (van Genuchten, 1980) was fitted to the data using the
RETC least-square optimization program of van Genuchten et al. (1991) to obtain thefollowing values for the soil hydraulic parameters: α =3.01 m-1, n=1.57, Ks = 3.47.10-6
ms-1, θr=0.05, and θs=0.51.
The subsurface source was connected to a Mariotte apparatus that maintained
a constant input pressure, Pin. The discharge Q was measured by weighing the
Mariotte apparatus. A pressure transducer was inserted into the source to measure thepressure head in the source cavity. The initial water content was measured at several
points and depths using TDR.
Measured and calculated backpressures ψs as a function of time fortwo sources with nominal discharge rates
Q0 of 4 and 8 L h-1 are presented in Fig.
2 for a radius r0 of 0.01 m and an inlet
pressure Pin of 10 m. The initial watercontent was considered to be constant(0.07) throughout the flow domain. For
both sources ψs increased rapidly during
the first 10 to 15 minutes, but thengradually approached a final value astime proceeded. The final value of ψs
was 2 and 3.7 m for the 4 and 8 L h-1
sources, respectively. Notice theexcellent agreement between themeasurements and the calculations.
Figure 2. Measured (triangles and circles)
The transient decrease in theand calculated (lines) soil pressure heads as asource discharge is depicted in Fig. 3 forfunction of time for two subsurface sourcesa nominal discharge Q0 of 8 L h-1, an. -1(4 and 8 L h ) in the Magal clay loam with
inlet pressure Pin of 10 m, and a radius
the source radius r0 of 0.01 m and the inlet
r0 of 0.01 m. Following an increase inpressure Pin of 10 m.
ψs from zero to 3.7 m (Fig. 2), thepressure difference between the soilback pressure and the source inlet
decreased, and the discharge rate
decreased from 2.22.10-6 m3 s-1 to
1.73.10-6 m3.s-1. HYDRUS-2D predicted
the measurements again very well. The
results presented in Fig. 3 demonstrate
the transient nature of the subsurface
source discharge and the significant
decrease in discharge due todevelopment of the back pressure.
Figure 3. Measured (triangles) and calculated(lines) discharge from a subsurface source as afunction of time for an 8 L.h-1 source and theMagal clay loam with the source radius r0 of 0.01
m and the inlet pressure Pin of 10 m.

55

Conclusions
Good agreement was found between transient back pressures measured in thelaboratory experiment and those calculated using the updated HYDRUS-2D
numerical model. The modified model allows the use of any hydraulic model for soil
and source properties, simulation of both short and long duration infiltration events,
and consideration of various source geometrical shapes. The modified numericalmodel also enables complex processes in the root zone involving short-time irrigation
and root water uptake.
Previous infiltration models that describe transient flow from subsurface
sources neglect the possibility that the soil hydraulic properties will limit the
infiltration capacity. This phenomenon can lead to an increase in the water pressure
head at the source, and hence a lower source discharge rate. If neglected, thedependence of the subsurface source discharge on the soil hydraulic properties canlead to smaller than designed irrigation volumes when the irrigation amount is presetby time, and consequently to considerable non-uniformity in water application rates(Warrick and Shani, 1996).
References
Ben Gal, A., N. Lazarovitch, and U. Shani. 2004. Subsurface drip irrigation in gravel
filled cavities. Vadose Zone J. 2004 3: 1407-1413.
Gardner, W. R. 1958. Some steady-state solutions of the unsaturated moisture flow
equation with application to evaporation from a water table. Soil Sci., V. 85,
P. 228-232.
Hillel, D. 1971. Soil and Water Physical Principles and Processes. Academic Press,
New York. pp 82–83,106–107,228–230.
Klute, A. 1986. Water retention:Laboratory methods.In Methods of Soil Analysis,Part
1 2nd edition.Ed.A Klute. pp 635–662.ASA Inc.,Madison,WI.
Lazarovitch, N. J., J. Šimůnek, and U. Shani. 2005. System dependent boundary
condition for water flow from subsurface source, Soil Sci. Soc. Am. J., 69, 4651.
Lazarovitch, N., Thompson, L.T., Warrick, A.W. and Shani, U. Soil Hydraulic
Properties Effecting Discharge Uniformity of Subsurface Drip Irrigation.
Submitted to J. Irrig. and Drain. Eng.
Philip, J. R. 1992. What happens near a quasi-linear point source? Water Resour.
Res., V. 28, P.47-52.
Šimůnek, J., M. Sejna, and M. T. van Genuchten, 1999. The HYDRUS-2D software
package for simulating two-dimensional movement of water, heat, andmultiple solutes in variably saturated media, Version 2.0, Rep. IGWMC-TPS53, 251 pp., Int. Ground Water Model. Cent., Colo. Sch. of Mines, Golden.
Shani, U., Xue, S., Gordin-Katz, R. and Warrick, A. W. 1996. Soil limiting flow fromsubsurface emitters. 1: Pressure measurements . J. Irrig. and Drain. Eng.
122(5), P. 291-295.
van Genuchten, M.Th. 1980. A closed-form equation for predicting the hydraulicconductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892-898.
van Genuchten, M.Th., F.J. Leij, and S.R. Yates. 1991. The RETC code forquantifying the hydraulic functions of unsaturated soils. USEPA Rep. 600/291-065 (IAGDW1293-3-934). U.S. Salinity Laboratory, USDA-ARS,
Riverside, CA.
Warrick, A.W. and Shani, U. 1996. Soil limiting flow from subsurface emitters. 2:
Effect on uniformity. J. Irrig. and Drain. Eng. 122(5), P. 296-300.

56

Modelling groundwater - surface water interactions in Dutch fens
Paul P. Schot and Stefan C. Dekker
Department of Environmental Sciences
COPERNICUS Institute for Sustainable Development and Innovation
Utrecht University
P.O. Box 80.115, 3508TC Utrecht
The Netherlands
p.schot@geo.uu.nl

Introduction
Wetland ecosystems in the Netherlands are heavily influenced by human interventions.
From the twelfth century onward, the wetlands have increasingly been drained for
agricultural purposes. Systems of connected drainage canals (ditches) were dug to create
polders in which the surface water levels were artificially controlled using windmills todischarge excess water. Nowadays most wetland areas are drained and predominantly
used for dairy farming.
A number of wetland remnants within the agricultural land use matrix presently aredesignated as protected nature reserves. However, they still contain the drainage canals inwhich surface water levels are artificially adjusted to the needs of the surrounding
agricultural areas. Amongst the protected nature reserves are fen wetlands which arecharacterised by a species rich vegetation and high biodiversity. Inflow of nutrient poor,
alkaline groundwater is generally seen as a prerequisite for the development of thisspecies rich fen vegetation (e.g. Wheeler and Shaw, 1995).
However, drainage and groundwater abstraction have affected the hydrological
functioning of the fens. These human interventions may lead to the development of socalled rainwater lenses in the upper groundwater of fens, which prevent upward seeping
alkaline groundwater from reaching the fen root zone (fig. 1). This threatens theconservation of the species rich vegetation.
As a basis for conservation and restoration of fen biodiversity, numerical groundwaterflow simulations were performed for a case area in the Netherlands, in order to provide
insight in rainwater lens dynamics.

Figure 1. Schematic representation of a rainwater lens on top of upwelling groundwater.
57

Method
The numerical two-dimensional unsaturated-saturated groundwater flow and transportmodel Hydrus-2D (Simunek et al., 1999) was used to simulate the form of rainwater
lenses characteristic for the case area, as a function of hydrological boundary conditions
and geology.
A two-dimensional section representative for the fen remnants area was modelled. The
section is perpendicular to the parallel drainage canals that are characteristic of thepolders in the case area. Left and right model boundaries represent flow lines locatedhalfway two drainage canals making them no flow boundaries. The atmospheric uppermodel boundary automatically switches between Dirichlet and Neumann boundary
conditions if minimum or maximum surface pressure head is reached. The lower model
boundary was defined as a constant flux boundary assuring upward seepage at all times.
The model drainage canals have a free drainage face above the fixed surface water level,
while a fixed head boundary applies below the surface water level. The model was set to
instantly remove exfiltrating groundwater at the soil surface (e.g. surface runoff) when
the phreatic water level exceeds the soil surface. Ponding was not considered.
Permeabilities of the sandy aquifer and semi-confining peat and clay layers are
representative for the sediments in the case area and the aquifer is consideredhomogenous and isotropic.
The model was used to provide insight in rainwater lens dynamics as a function ofdrainage canal water levels, upward groundwater inflow fluxes and the thickness of semiconfining layers. Table 1 presents boundary conditions for the steady state simulations.
Thicknesses of the semi-confining layer of 0, 1 and 2 m were simulated.
Table 1. Boundary conditions for the steady state simulations. (W and S = weak andstrong drainage; – and + = low and high upward groundwater flux)
Run
W+
WS+
S-

Drainage canal water level(m NAP)
-0.30
-0.30
-0.60
-0.60

Lower boundary upward flux
(mm d-1)105
10
5Results
The simulations indicate that drainage canal surface water levels mainly affect thephreatic level in the fens, while groundwater inflow fluxes mainly affect the thickness of
the rainwater lens. Increasing thickness of a semi-confining surface layer results in a risein phreatic level. As soon as the phreatic level reaches the soil surface the rainwater lens
may split in two as a result of groundwater discharge windows, which come into
existence in the centre of the parcel.
58

Figure 2. Yearly average form of the rainwater lens as a function of:- weak and strong drainage (W and S) and groundwater inflow flux (- and+) (vertical axis);- presence of a semi-confining surface layer (0, 1 and 2 meter) (horizontalaxis).
For W-, W+, S-, S+ see table 1. (source: Schot et al, 2004)

Discussion
Most simulations show rainwater lens thicknesses well exceeding rooting depths of fen
plants. This is inconsistent with the fact that fen vegetation is still found in nature
reserves for which the simulations are considered representative. Chemical buffer
processes in the soil may explain this inconsistency, masking an acidification process
under way. As buffer capacities will become exhausted over time, acidification and fen
59deterioration will occur at a certain future moment. To arrive at more definite conclusions
about this hypothesis, our simulations of steady-state convective groundwater flow need
to be supplemented by transient modelling of groundwater flow, as well as ofhydrochemical processes in rainwater lenses, followed by verification at actual field
locations. These activities are presently underway.

References
Schot, P.P., Dekker, S.C. and Poot, A.. The dynamic form of rainwater lenses in drained
fens. Journal of Hydrology, 293 (2004) pp. 74-84.
Simunek, J., Sejna, M. and Van Genuchten, M.T., 1999. The HYDRUS2D software
package for simulating the two-dimensional movement of water, heat, and multiple
solutes in variably-saturate media. U.S. Salinity Laboratory, USDA/ARS, Riverside,
California, pp. 227.
Wheeler, B.D. and Shaw, S.C., 1995. A focus on fens - Controls on the composition offen vegetation in relation to restoration. In: B.D. Wheeler, S.C. Shaw, W.J. Fojt and R.A.
Robertson (Editors), Restoration of temperate wetlands, Chichester, Wiley, pp. 33-48.

60

Chlorotoluron transport in the soil profile affected by non-equilibrium flow
Radka Kodešová1, Martin Kočárek 1, Jiří Šimůnek2, Josef Kozák1
1
Czech University of Agriculture in Prague, Department of Soil Science and Geology,Kamýcká 129, 16521 Prague, Czech Republic, kodesova@af.czu.cz, kocarek@af.czu.cz,
kozak@af.czu.cz
2
University of California Riverside, Department of Environmental Sciences, Riverside, CA
92521, USA, jiri.simunek@usr.edu
Introduction
Non-equilibrium water flow and contaminant transport in soil porous media are frequently
observed phenomena. Numerical models that assume bi-modal soil porous systems have beendeveloped to describe non-equilibrium flow and solute transport. Single-porosity, dualporosity and dual-permeability models in HYDRUS-1D (Šimůnek et al., 1998, 2003) wereapplied in this study to estimate soil hydraulic parameters from laboratory multi-step outflow
experiment via numerical inversion and to simulate chlorotoluron transport in the soil profilethat was experimentally studied in the field.
Materials and Methods
The transport of chlorotoluron in the soil profile was studied under field conditions. The
herbicide Syncuran was applied on a four square meter plot using an application rate of 2.5
kg/ha of active ingredient. Chlorotoluron dissolves in water, adsorbs on soil particles anddegrades with time. Soil samples were taken after 6, 14, 22 and 36 days to study the residual
chlorotoluron distribution in the soil profile.
The soil was defined as Greyic Phaozem. The soil profile was divided into two layers: 0-35
cm and 35-80 cm. The soil hydraulic properties were defined as follows. The multi-stepoutflow experiments (van Dam et al., 1994) were performed on the 100 cm3 undisturbed soil
samples (A,B,C,D) using the 15, 30, 70, 130 and 210 cm pressure head steps. The pressure
head steps and the soil-water contents obtained from the final soil-water contents and waterbalances in the soil samples provide data points of the soil-water retention curves. The code
HYDRUS-1D (Šimůnek et al., 1998, 2003) and the numerical inversion (Hopmans et al.,2002) were used to analyze the cumulative outflow and the soil-water retention data points toobtain hydraulic parameters characterizing different soil-water flow models: the singleporosity model, the dual-porosity model, and the dual-permeability model. In all cases limited
sets of parameters had to be optimized to ensure unique results.
The chlorotoluron transport under field conditions was simulated using the single-porosity,dual-porosity, and dual-permeability models in HYDRUS-1D. The top boundary condition
was defined using daily precipitation and potential transpiration rates. Given the root zone
depth of 20 cm, a Feddes model in HYDRUS-1D with parameters defined for wheat wasapplied to simulate the root water uptake. The bottom boundary condition was defined as afree drainage. The bulk densities were 1.5 g/cm3 and 1.55 g/cm3 for the first and second
layers, respectively. The Freundlich adsorption isotherm was measured for the first layer: kF =2.73 and nF = 1.33. Degradation rate 0.002 day-1 was estimated from the herbicide content 36
days after the herbicide application (Kočárek et al., 2005). Degradation was assumed to occur
in both water and solid phases. The reasons were previously discussed in Kodešová et al.

61(2004) for Chermozem. The longitudinal dispersivity was set to 2 cm for the single- and dualporosity models, and 2 and 1 cm in matrix and macropore domain, respectively, for the dualpermeability model. The molecular diffusion was 1 cm2day-1 in all cases.
Results and Discussion
Soil hydraulic parameters (van Genuchten, 1980) for the single-porosity model obtained using
the numerical inversion of the cumulative outflow and soil-water retention curve data areshown in Tab. 1. Values of θs were obtained from the water balance in the soil samples.
Soil hydraulic parameters for the dual-porosity model obtained using the numerical inversionof cumulative outflow and soil-water retention curve data are shown in Tab. 2. Values of
θr,t, θs,t (t = total), α m and n m (m = mobile domain) are equal to values of θr , θs , α and n forthe single-porosity model. Values of θr,im and θs,im (im = immobile domain) were estimated. ωw
is the water transfer coefficient between the immobile and mobile domains.

Table 1: Soil hydraulic parameters – single-porosity model.
θr
A
B
C
D

( 1 - 35)
( 1 - 35)
( 35 - 80)( 35 - 80)

3

θs

-3

[cm cm ]
0.293
0.306
0.126
0.000

3

α

-3

[cm cm ]0.391*
0.389*
0.413*
0.410*

n

-1

[cm ]
0.01410.01010.01400.0145

1.891
2.289
1.255
1.145

Ks
[cm/day]
1.061
0.382
20.40
10.78

* not optimized
Table 2: Soil hydraulic parameters – dual-porosity model.
θr,t
A
B
C
D

( 1 - 35)
( 1 - 35)
( 35 - 80)( 35 - 80)/ θr,im[cm3cm-3]0.293* / 0.293*
0.306* / 0.306*
0.126* / 0.126*
0.000* / 0.000*

θs,t/ θs,im[cm3cm-3]
0.391* / 0.320*
0.389* / 0.320*
0.413* / 0.250*
0.410* / 0.250*

αm[cm-1]
0.0141*0.0101*0.0140*0.0145*nm
1.891*
2.289*
1.255*
1.145*

Ks,m[cm/day]
1.020
0.334
19.92
12.48

ωw

[day-1]
5.86 10-5
2.25 10-4
4.85 10-5
9.07 10-4* not optimized
Table 3: Soil hydraulic parameters – dual-permeability model.
θr,f / θr,mA ( 1 - 35)
B ( 1 - 35)
C ( 35 - 80)
D ( 35 - 80)[cm3cm-3]
0* / 0.293*
0* / 0.126*
0* / 0.000*

θs,f/ θs,m[cm3cm-3]0.450* / 0.391*
0.450* / 0.413*
0.450* / 0.410*

αf

/ αm
[cm-1]
0.023 / 0.0037
0.030 / 0.0028
0.022 / 0.0027

nf / nm
3* / 2.990
3* / 1.689
3* / 1.237

Ks,f / Ks,m
[cm/day]
9.720* / 0.012
187.2* / 0.142
98.88* / 0.219* not optimized
Soil hydraulic parameters for the dual-permeability model obtained using the numericalinversion of cumulative outflow data are shown in Tab. 3. Following parameters defining thestructure of matrix (m) and macropore (f) domains were used in all cases: w = 0.1, β = 15, d =0.3 cm, γw = 0.4. Values of θr,m a θs,m are equal to values of θr and θs for the single-porosity
model. Values of θr,f , θs,f a nf were estimated. Values of Ks,f were determined based on the62following assumptions: 0.9 Ks,m + 0.1 Ks,f equals to Ks for the single-porosity model andKs,f values are hundred times larger then initial estimates of Ks,m values. Experimental data forthe sample B could not be successfully analyzed using this model. In other three cases thesum of soil-water retention curves of both domains multiplied by corresponding ratios of eachdomain follow the measured soil-water retention data points.
The hydraulic parameters obtained for the soil samples A and C were used to simulate
chlorotoluron transport in the two-layered soil profile. Calculated chlorotoluron
concentrations 6, 14, 22 and 36 days after the application are shown in Fig. 1a for the singleporosity model, in Fig. 1b for the dual-porosity model and in Fig. 1c for the dual-permeability
model. Measured chlorotoluron concentrations are shown in Fig. 1d.
2

Concentration [µg/g]
4

6

0

8

0

0

-2

-2

-4

-4

-6

-6

Depth [cm]Depth [cm]

0

-8
-10
-12

6 days14 days
22 days
36 days

-14
-16
-18

2

4

6

8

-8
-10
-12

6 days14 days
22 days
36 days

-14
-16
-18
-20

-20

a)

b)
Concentration [µg/g]

0

2

4

6

0

8

0

0

-2

-2

-4

-4

-6

-6

Depth [cm]Depth [cm]Concentration [µg/g]

-8
-10
-12
-14
-16
-18

6 days14 days
22 days
36 days

2

4

6

8

-8
-10

6 days14 days
22 days
36 days

-12
-14
-16
-18

-20

Concentration [µg/g]

-20

c)

d)

Figure 1. Chlorotoluron distribution in the soil profile a) single-porosity model, b) dualporosity model, c) dual-permeability model, d) measured values (averages from 3 sampling
positions).
Figs. 1a and 1b show that immobile water in the dual-porosity system had only a little impact
on the solute transport. Solute moved to a depth of 10 cm in the single-porosity systemcompared to 9 cm in the dual-porosity system. In the case of the dual-permeability systemsolute moved to a depth of 60 cm (depth not shown in Fig. 1c). Comparison of all simulated
and measured chlorotoluron distributions makes obvious that the results of the dualpermeability model correspond more to the measured values then the other two models due tothe fact that observed herbicide transport was highly affected by a preferential flow. It shouldbe mentioned that the measured total chlorotoluron contents in the soil profile 22nd day after

63

the herbicide application is higher then the total chlorotoluron concentrations measured on 6th
and 14th day. This may be due to an uneven distribution of chlorotoluron on the soil surface,herbicide transport through the preferential paths and its deviation from the vertical axes.
Observed oscillations of chlorotoluron concentrations in depths of 6 – 20 cm were probably
caused by preferential flow (quick solute penetration to deeper depths) and soluteaccumulation due to the disconnection of preferential pathways. This phenomenon was notconsidered in the model.
Conclusions
It is shown that, despite having very similar total soil hydraulic prosperities for different flow
models, the simulated chlorotoluron transport is different. Observed chlorotoluron transport
was affected by preferential flow. Therefore chlorotoluron concentrations in the soil profile
simulated using the dual-permeability model are closer to observed values then those
calculated with the other two models. The impact of different definitions of the soil porous
system as single-porosity or dual-permeability systems with the similar total soil hydraulicproperties was previously discussed for Chernozem in Kodešová et al. (2005). The dualpermeability model is a powerful tool describing preferential flow and transport.
Acknowledgement
Authors acknowledge the financial support of the Czech Grant Agency GA CR No.
103/05/2143. Thanks are due to K. Němeček for his help with the fieldwork and O. Drábek
for the chemical analysis.
References
Hopmans, J. W., Šimůnek, J., Romano, N., Durner, W. (2002): Inverse methods, in Methodsof Soil Analysis, Soil Science Society of America, Inc. Madison, Wis., 963-1008.
Kočárek, M., Kodešová R., Kozák, J., Drábek, O., Vacek, O. (2005): Chlorotoluron behaviour
in five different soil types. Plant, Soil and Environment, 51 (7), 304-309.
Kodešová, R., Kozák, J., Vacek, O. (2004): Field and numerical study of chlorotolurontransport in the soil profile, Plant, Soil and Environment, 50 (8): 333-338.
Kodešová, R., Kozák, J., Šimůnek, J., Vacek, O. (2005): Single and dual-permeability modelof chlorotoluron transport in the soil profile. Plant, Soil and Environment, 51 (7), 310315.
Šimůnek J., Šejna M., van Genuchten M. Th. (1998): The HYDRUS-1D software package forsimulating the one-dimensional movement of water, heat and multiple solutes invariably-saturated media. Version 2.0. IGWMC-TPS-53. International Ground Water
Modeling Center, Colorado School of Mines, Golden, CO.
Šimůnek J., Jarvis N. J., van Genuchten M. Th., Gärdenäs A. (2003): Review and comparisonof models for describing non-equilibrium and preferential flow and transport in thevadose zone, Journal of Hydrology, 272, 14-35.
van Dam, J. C., Stricker, N. M., and Droogers, P. (1994): Inverse method to determine soil
hydraulic functions from multistep outflow experiments. Soil Sci. Soc. Am. J., 58, 647652.
van Genuchten M. Th. (1980): A closed-form equation for predicting the hydraulicconductivity of unsaturated soils. Soil Sci. Soc. Am. J., 44: 892-898.

64Simulating simultaneous nitrification and denitrification with a modified 2D-mobile immobile
model.
Sigrid Köhne1*, Jiri Šimůnek2, J. Maximilian Köhne1 and Bernd Lennartz1
1 Introduction
Permanent grassland and cattle production characterize agriculture in many lowland regions withartificial drainage systems. In addition to mineral fertilizers, cattle manure is used as nitrogen (N)source in such areas. Whereas the manure carries about mostly ammonium (NH4) as directly plant
available N, mineral fertilizer N may contain nitrate (NO3) or NH4, or both.
Nitrification and denitrification are biological processes within the soil nitrogen cycle.
Ammonium is oxidized to NO3 by autotrophic nitrifying bacteria under aerobic conditions. Under
anaerobic conditions, NO3 may be used by heterotrophic denitrifying bacteria during respiration to
be reduced into gaseous NO2 or N2. The facultative anaerobic denitrifying bacteria need a very low
oxygen (O2) concentration to accept NO3 as terminal electron acceptor. In soils, very low O2
concentrations are usually associated either with water saturation (Wilson and Bouwer 1997) orwith spots of high decomposition (Nielsen et al. 1996, De Klein et al. 1996). Moreover, nitrificationand denitrification depend on soil temperature and soil pH.
Because nitrification and denitrification proceed
under contrasting redox conditions, they are regarded
as separate processes with respect to time or location.
Recently, there is increasing evidence that nitrificationand denitrification occur simultaneously and in close
vicinity to each other in the soil (Abassi and Adams
1998, De Klein et al. 1996, Nielsen et al. 1996, Zanner
and Bloom 1995). Vertical Redox gradients across thefirst few centimetres from the surface resulted inparallel occurrence of nitrification and denitrification
in a grass covered compacted soil (Abassi and Adams
1998). Anaerobic-aerobic interfaces may prevail atspots of high decomposition of organic matter ormanure and can induce coupled nitrificationdenitrification (Nielsen et al. 1996). The NH4 released
from such anaerobic spots to the neighbouring aerobic
Fig 1: Small scale variation of soil
soil material was found to be nitrified at the anaerobicproperties in aggregated soils.
aerobic interface, whereupon the nitrification product,
NO3, diffused back into the anaerobic spot to be denitrified (Nielsen et al, 1996). Similarly, De
Klein et al. (1996) found that mineralization-induced O2 consumption in combination with highwater contents led to the development of denitrification zones within a basically aerated soil.
Adjacent cells of aerobic and anaerobic conditions may also be caused by local variations in watercontents and lead to simultaneous nitrification and denitrification in flooded soils (Lindau et al,1988). Zanner and Bloom (1995) suspected anaerobic soil aggregate centres and aerobic interaggregate space (Fig. 1) to cause simultaneous nitrification and denitrification in soil.
________________________________________________________________________________
1

Institute for Land Use, Faculty for Land use and Agricultural Sciences, University Rostock, Justus-von-Liebig-Weg
6,18044 Rostock, Germany, sigrid.koehne@uni-rostock.de, max.koehne@uni-rostock.de, bernd.lennartz@unirostock.de , * Corresponding author
2
Dept. of Environmental Sciences, University of California, Riverside, 900 University Avenue, A135 Bourns Hall,Riverside, CA 92521, jiri.simunek@ucr.edu

65

The Infeld experimental station is situated in northern Germany and has been used to monitor,
among others, NO3 and NH4 leaching losses via tile drains over decades. In a study about N
leaching losses and N balance of the site under pasture, Spatz et al. (1989) found that the highamount of N-input by fertilizer is not recovered with the exported crop plus the N loss in drainage
water. The authors conclude that a considerable amount of the fertilizer N is lost by denitrification.
Concentration (mg/l)

10Bromide measured

2D-MIM

2D-SPM8
6
4
2
0

11.21.99
12.21.99
Bromide application

1.20.00

2.19.00
Date

3.20.00

4.19.00

Fig. 2. Bromide concentrations in tile drainage water as compared to simulations with the 2D
mobile immobile model (MIM) and the 2D single porosity model (SPM).
A tracer study performed at Infeld was analysed with different mobile immobile (MIM) andsingle porosity (SPM) models and revealed that transport through the soil towards drains was
influenced by immobile water regions (Fig. 2). While the mobile regions accelerated bromidetransport through the soil, the immobile regions acted as a sink for more than 60% of the surface
applied tracer mass (Köhne et al. 2005).
Our hypothesis is that immobile regions hold high water contents over long periods and thus are
more often anaerobic than the adjacent mobile water regions that change their water contents more
dynamically in response to variable atmospheric conditions. Physically, the immobile regions may
represent soil aggregates with rather small pores. The mobile regions may represent the interaggregate space. As with decomposition induced anaerobicity described above, NO3 can beproduced from added NH4 in the aerated mobile regions, and then be transferred into the immobile
regions (aggregates) where the reductive denitrification process is taking place. During dry andsufficiently warm periods, nitrification can also occur
within the unsaturated immobile region.
The basic idea of this study is to investigate longterm nitrate transport as affected by both
immobilization and metabolization. To follow theobjective, we extended the 2D-MIM, as implemented in
HYDRUS-2D, with a relatively simple nitrogenmetabolism model, where the nitrogen reaction rates in
each, the mobile and immobile regions, are functions ofbulk soil temperature and the mobile and immobile
water contents, respectively.
2 Extension of the 2D-MIM
Figure 3 shows exemplary water contents in mobile
and immobile regions as simulated with the HYDRUS2D MIM. The simulation is based on model parameters
as given in Köhne et al (2005) for the Infeld plot B. The water contents of the mobile region reach a66

lower saturation than those of the immobile region. Fig. 4 is an extension of Fig. 3 that shows theimplementation of simultaneous aerobic and anaerobic processes. Nitrification and denitrification
were calculated as first order processes with potential rates µ’ and µ, respectively, and response
functions considering the effects of soil moisture.

Fig. 4. Coupling the dual porosity concept with water content dependent processes rates toobtain simultaneous occurrence of aerobic and anaerobic processes.
The response functions for the aerobic and anaerobic processes are depicted in the green andyellow frames of Fig. 4, respectively. Under aerobic conditions, denitrification is reduced to aminimal rate of rmin0 times µ. The relative rate r increases linearly from rmin0 to 1 over a watercontent range between θ0 and θ1. In a similar way, relative nitrification rates reduce from 1 to rmin1
over the range specified by θ2 and θ3. The soil moisture response functions are separately definedfor each, mobile and immobile water contents. Besides the soil region moistures, temperaturedependence of the process rates was considered (not shown in Fig. 4).
3 First simulation results
Fig. 5 shows simulations of NO3 concentrations in drain outflow (symbols) and in the mobile
and immobile region at the drain (light and dark lines, respectively) as compared to the 12-year
measurements at the Infeld site. Simulation 1 and 2 are based on N transformations limited to theupper 25 cm (top soil) of the 2-m soil profil with 1 m tile drain depth and 18 m drain spacing. They
differ in the water content thresholds with nitrification (denitrification) decreasing (increasing)
between 93and 95 % saturation of both, the mobile and immobile zones for simulation 1 andbetween 90 and 93 % saturation of the regions for simulation 2. Simulation 3 allows for N
transformations in the entire profile with the saturation thresholds set as in simulation 1. Increasing
denitrification (by setting the thresholds in a lower saturation range) yields lower concentrations at
the drain (compare simulations 1 and 2).
67NO3-concentration (mg/L)

20
Simulation 1
15

Simulation 2
Simulation 3

10
50
Nov-89

Nov-90Nov-91

Nov-92

Nov-93

Nov-94 Nov-95
Date

Nov-96

Nov-97

Nov-98

Nov-99

Nov-00Fig. 5. Twelve year measurements of nitrate concentrations in tile drain outflow andsimulations of NO3 concentrations in drain outflow (symbols) and in the mobile and immobile
region at the drain (light and dark lines, respectively)
To allow for water content controlled denitrification in the entire profile (simulation 3), rather
than in the upper 25 cm only, leads to zero-concentration at drain depth at the end of each drainage
period as observed with the data. The dynamic NO3 pattern across drainage seasons, with highconcentrations at the start and low concentrations at the end of each drainage season, lies between
simulations 2 (upper horizon is involved in N-metabolisms) and 3 (entire profile involved).
4 Conclusions
The NO3concentrations measured for more than one decade in tile drainage water at the Infeld
site seem to be governed by both physical non-equilibrium and N-transformations. The 2D-MIM
model extended with water content dependent rate response functions for N-transformations in
mobile and immobile regions could simulate simultaneous aerobic and anaerobic processes andcould reproduce the principal observed NO3 concentration pattern. Due to convective transport andlimited exchange between regions, NO3 concentrations fluctuated stronger in the mobile regions
than in the immobile regions. Both measured NH4 and NO3 concentrations will be used next toimprove the model calibration and understanding of complex N fate and transport.
5 References:
Abassi M.K., Adams W.A., 1998. Loss of nitrogen in compacted grassland soil by simultaneous nitrification anddenitrification. Plant and Soil 200: 265-277.
De Klein C.A.M., Van Logtestijn R.S.P., De van Meer H.G., Geurink J.H., 1996. Nitrogen loss due to denitrificationfrom cattle slurry injected onto grassland soil with and without a nitrification inhibitor. Plant Soil183: 161-170.
Köhne S., Lennartz B, Köhne J.M., Šimůnek J., 2005. Bromide transport at a tile-drained field site: experiment, andone- and two-dimensional equilibrium and non-equilibrium numerical modeling. J. Hydrol, accepted 8/05.
Nielsen T.H., Nielsen L.P., Revsbech N.P., 1996. Nitrification and coupled nitrification-denitrification associated with asoil-manure interface. Soil Sci. Soc. Am. J. 60: 1829-1840.
Šimůnek, J., Šejna, M., van Genuchten, M.Th., 1999. The HYDRUS-2D software package for simulating twodimensional movement of water, heat, and multiply solutes in variably saturated media. Version 2.0, IGWMC-TPS53, International Ground Water Modelling Center, Colorado School of Mines, Golden, Colorado, p.202.
Spatz G., Neuhaus H., Pape A., 1992. Stickstoffdynamik einer Mähweide auf Knickmarsch. J. Agronomy and Crop
Science 168: 298-309.
Wilson L.P., Bouwer E.J. 1997. Biodegradation od aromatic compounds under mixed oxygen/denitrifying conditions: A
review. J. Ind. Microbiol. Biotech. 18: 116-130.
Zanner C.W, Bloom R.R., 1995. Mineralization, Nitrification and Denitrification in Histosols of Northern Minnesota.
Soil Sci Soc. Am. J. 59: 1505-1510.

68

Recent Applications of HYDRUS in Irrigated Agriculture
T. H. Skaggs1, T. J. Trout2, T. J. Kelleners3, P.J. Shouse1, J.E. Ayars2
1

George E. Brown, Jr. Salinity Laboratory, USDA-ARS, Riverside, CA, USA
2
Water Management Research Laboratory, USDA-ARS, Parlier, CA, USA
3
Dept. of Plants, Soils, and Biometeorology, Utah State University, Logan, UT, USA

Introduction
With its wide ranging capabilities for simulating variably saturated water flow and solutetransport, HYDRUS (Šimůnek et al., 1998; 1999) has proven to be a valuable tool for
analyzing and designing irrigated agricultural systems. We briefly review three recent
projects that illustrate the versatility of the model and software.
Drip Irrigation
Drip irrigation technology gives growers greater control over the application of water,
fertilizers, and pesticides; increased precision in water and chemical applications can leadto water savings, improved environmental quality, and increased profits. Realizing the fullpotential of drip technology requires optimizing the operational parameters that areavailable to irrigators, such as the frequency and duration of irrigation, the emitter
discharge rate and spacing, and the placement of drip tubing. Numerical simulation is anefficient approach to investigating optimal drip management practices. Unfortunately, little
work has been done to demonstrate the accuracy of numerical simulations, leading some toquestion the usefulness of simulation as a research and design tool.
Skaggs et al. (2004) compared HYDRUS-2D simulations of drip irrigation withexperimental data. A Hanford sandy loam soil was irrigated using thin-walled drip tubing
installed at a depth of 6 cm. Three trials (20, 40, and 60 L m-1 applied water) were carriedout. At the end of each irrigation and approximately 24 h later, the water content
distribution in the soil was determined by gravimetric sampling. Soil hydraulic properties
were estimated using the Rosetta pedotransfer function model (Schaap et al., 2001). The
HYDRUS-2D predictions of the water content distribution were found to be in very good
agreement with the data. Figure 1 shows the simulated and measured wetted soil profile
for the 40 L m-1 trial.
The results of Skaggs et al. (2004) provide support for the use of HYDRUS-2D as a toolfor investigating and designing drip irrigation management practices. The results areespecially noteworthy given that the simulations were pure predictions, made without any
fitting to the experimental data.

69

40 L m-1 Applied Water
Time = 10.75 h
DISTANCE (cm)
10

20

30

40

DISTANCE (cm)DISTANCE (cm)
50

60

10

20

30

40

50

60

0

-10

-10

-20

-20
DEPTH (cm)

DEPTH (cm)0Time = 31 h

-30
-40

20

30

40

DISTANCE (cm)
50

60

10

20

30

40

50

60

-30
-40
-50

-50
-60

10

OBSERVED-60

PREDICTEDOBSERVED

PREDICTED

-70

-70

Figure 1 Comparison of measured and predicted volumetric water content following a 10
h irrigation (40 L m-1 applied water) from a drip line buried 6 cm below the surface (fromSkaggs et al., 2004).

Effects of Salinity and Drought Stress on Transpiration and Drainage
The disposal of agricultural drainage waters into surface waters or onto lands leads tosalinization and degradation of soil and water quality. In some instances the impacts ofagricultural drainage disposal may be reduced if drainage waters are isolated and reused for
irrigation. Over the last 25 years, the recycling of drainage water for irrigation has been
examined in numerous experiments and demonstration projects. Scientific and economic
questions remain about the suitability of different waters, soils, and crops for reuse
operations; optimal management practices are still largely unknown.
While modeling offers a cost-effective means of developing optimal management practices,questions exist about the accuracy of simulations for the dynamic, highly saline
environments the may exist in reuse operations. A recent study by Skaggs et al. (2005a,b)compared HYDRUS-1D simulations with extensive drainage and root water uptake datacollected on forage crops grown using synthetic drainage waters. Experiments wereconducted in a facility consisting of 24 volumetric lysimeters, each measuring 81.5 cm
wide × 202.5 cm long × 85 cm deep. Twelve lysimeters were planted in ‘Salado’ alfalfa(Medicago sativa) and 12 in ‘Jose’ tall wheatgrass (Agropyron elongatum). Experimental
treatments consisted of different combinations of irrigation water salinities (ranging from
ECiw = 2.5 dS m-1 to 28 dS m-1) and irrigation rates (ranging from deficit to excessive).
To model reductions in root water uptake due to salinity and water stress, Skaggs et al.
(2005b) used a multiplicative model for the uptake reduction factor in HYDRUS (Šimůnek
et al., 1998):

α ( h, hφ ) = [1 − A( hφ − B )]⋅1
1 + ( h / h50 ) p

(1)

70

CUMULATIVE DRAINAGE

MEASUREDHYDRUS' SALADO' ALFALFA
0
-10

-30

ECiw = 8 dS/m

18 dS/m

Lysimeter #75

Lysimeter #78

10

28 dS/m
Lysimeter #62

5
0

-40
250

270

290

250

270

290

250

270

IRRIGATION (cm)

-20

290

Figure 2. Measured and simulated cumulative drainage for selected alfalfa lysimeters.
Irrigation for each lysimeter is shown along the bottom of each plot (the scale is indicatedon the right axes). Although the lysimeters received roughly the same amount of water,
drainage increases as ECiw increases because uptake is reduced (from Skaggs et al.,2005b).
where h is the matric pressure head and hφ is the osmotic pressure head. Values for theparameters A, B, h50, and p for both crops where identified through a global fit to
cumulative drainage data. The simulated and measured drainage curves for selected
lysimeters are shown in Figure 2. The agreement between the model and the data for theremaining 21 lysimeters was similar to that depicted in Figure 2. Although the model was
fitted, it is noteworthy that a single set of parameters was able to reproduce theexperimental results over a very broad range of experimental conditions: irrigation waters
with EC’s ranging from 2.5 to 28 dS/m; and irrigation rates ranging from deficit toluxurious.
Inverse Analysis of Upward Water Flow from Shallow Groundwater
It has long been known that in arid and semi-arid regions upward water flow from shallow
groundwater can salinize surface soils. In the San Joaquin Valley, California, thismechanism of salinization has been receiving renewed attention. The west side of the San
Joaquin Valley is underlain by shallow groundwater, a result of widespread irrigation. In
many areas, the groundwater must be drained from the root zone using subsurface drainage
systems. The drainage is saline and contains appreciable amounts of several toxic traceelements, including selenium. Environmental problems associated with the disposal of thisdrainage have led to efforts to reduce the amount of drainage being generated. At least two
drainage reduction strategies raise questions about the salinizing potential of shallow
groundwater: land retirement, in which formerly irrigated land is left fallow; and shallow
ground water management, in which drain flow is restricted in an effort to maintain thewater table near the bottom of the root zone and induce consumption of groundwater bycertain crops.

71Figure 3. Contour plot of the objective function Φ for the parameter combinations n-θr,
n-θs, and λ-n using θ(z,t) data in the objective function (from Kelleners et al., 2005).

Calculating rates of upward water flow and root zone salinization with HYDRUS requires
identifying appropriate parameter values for the soil hydraulic properties. In Kelleners etal. (2005), the hydraulic properties of a silty clay soil in a large weighing lysimeter were
estimated by inverse modeling of a summer fallow period during which capillary rise fromthe groundwater table replenished the depleted root zone. The analysis was performed with
the HYDRUS-1D parameter optimization routines using the objective function:p

mqk nqikk =1

i =1 j =1[

Φ ( b; q ) = ∑ wk ∑∑ qk* ( zi , t j ) − qk ( zi , t j ; b)](2)where b is the vector of parameters to be optimized, q k* is the kth measured datum of type q(which may be water content, matric pressure head, cumulative outflow, etc.), and qk is thecorresponding model calculation. In addition to estimating parameter values, Kelleners etal. (2005) investigated data requirements and problems of parameter uniqueness. Figure 3
shows plots of Φ and the location of minima for different combinations of parameters using
water content data in the objective function.
References
Kelleners, T.J., R.W.O Soppe, J.E. Ayars, T.H. Skaggs, and J. Simunek. 2004. Inverse analysis of upward
water flow in a groundwater table lysimeter, Vadose Zone J., 4:558–572.
Schaap, M.G., F.J. Leij, and M.Th. van Genuchten. (2001). ROSETTA: A computer program for estimating
soil hydraulic properties with hierarchical pedotransfer functions. J. of Hydrology, 251, 163–176.
Skaggs, T.H., J.A. Poss, and P.J. Shouse. 2005a. Irrigation of forage crops with saline-sodic drainage waters,
1. Volumetric lysimeter studies. In Review.
Skaggs, T.H., P.J. Shouse, and J.A. Poss. 2005b. Irrigation of forage crops with saline-sodic drainage waters,
2. Modeling drainage and root water uptake. In Review.
Skaggs, T.H., T.J. Trout, J. Simunek, and P.J. Shouse. 2004. Comparison of HYDRUS-2D simulations ofdrip irrigation with experimental observations, J. Irrig. Drain. Eng., 130:304-310.
Šimůnek, J., M. Šejna, and M. Th. van Genuchten. 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. IGWMC. Golden, Colorado.
Šimůnek, J., M. Šejna, and M.Th. van Genuchten. 1999. The HYDRUS-2D Software package for simulating
the two-dimensional movement of water, heat and multiple solutes in variably-saturated media.
IGWMC-TPS 53, Version 2.0. Int. Ground Water Modeling Center, Colorado School of Mines, Golden,

72

CO.

73

Seasonal Furrow Irrigation Modellingwith HYDRUS2
Thomas Wöhling1
Abstract
A seasonal furrow irrigation model (FIM) was developed by iteratively coupling a 1D
surface flow model with a series of HYDRUS2 model parameterisations along a single furrow. HYDRUS2 is further linked to a crop model, which calculates daily leafarea index and final crop yield. FIM simulations of soil water content and yield are ingood agreement with field experimental data from a corn plot in Montpellier, France.

Model Formulation
The seasonal furrow irrigation model (FIM) comprises process-based simultaneous modelling of the 1D surface flow, the quasi-3D soil water transport, and thecrop growth. The surface flow sub-model FAPS utilizes an analytical solution of
the Zero-Inertia open-channel flow equations for modelling the non-uniform flow
during both the advance and the early storage phase of an irrigation event(Schmitz et al., 1992). Process adequate simplifications are applied for the latestorage phase, the depletion and therecession phase (Wöhling, 2005).
The numerical code HYDRUS2 (Simúnek etal., 1996) portrays both the two-dimensional
infiltration from arbitrary shaped furrows andthe soil water transport in a series of verticalplanes along the furrow (Fig. 1). FAPS andHYDRUS2 are iteratively coupled by the
infiltration rate and the flow depth (Wöhlinget. al, 2004). Due to the increasing/decreasing flow depth during the simulationof an irrigation event, the boundary type offurrow nodes of the HYDRUS2 calculation
mesh changes from atmospheric boundary type to prescribed head type and viceversa. Both boundary type and value are assigned to the upper boundary nodes
accordingly to the transient conditions, which are determined by the surface flow
model and an evapotranspiration module (Wöhling, 2005).
The set of non-linear partial differential equations of the coupled furrow irrigation
advance model is solved by a space discretisation in combination with the Newton
iteration method. This new solution brought a fundamental improvement as regards convergence, numerical stability and computational time in comparison to
the common solution of the considered problem, namely a time-discretisation in
combination with the Fixpoint iteration. The new coupled sub-model compares perfectly to data from various laboratory and field irrigation events (Wöhling, 2005;
Wöhling et al., 2005).
Figure 1: FIM flow domain

1

Research Associate, Dresden University of Technology, Institute of Hydrology and Meteorology, Würzburger Str. 46, D-01187 Dresden, Germany. Email: Thomas.Woehling@mailbox.tu-dresden.de74For modelling the soil
water transport during
the whole growing period, HYDRUS2 andthe crop growth submodel LAI-SIM arecoupled by a waterstress index. LAI-SIM
includes modules forpredicting daily leafarea index, potential
transpiration / evaporation, and root growth Figure 2: Common boundaries between HYDRUS2 and the(Wöhling, 2005). The crop module LAI-SIM
HYDRUS2 boundary
conditions (during times of redistribution) are:• atmospheric boundary type at the soil surface: soil evaporation and precipitation (Fig. 2a)• lateral no-flux boundaries• seepage face type at the lower boundary• internal sink type within the rooting zone: root water uptake by the crop(Fig. 2b)
The root depth during the growing season is described by an linear root growthmodel (Wöhling, 2005) and the root activity is assumed to increase with depth asproposed by Novak (1987) (Fig. 2c).

Figure 3: FIM-GUI visualisation of simulated pressure head distribution at the furrow inlet
before (left) and after (right) an irrigation event

User Interface
The framework of the seasonal furrow irrigation system, FIM, is set up by a comprehensive event control and time management unit. A user-friendly graphical user75

interface (FIM-GUI) supports the parameterisation of the FIM sub-models. The
Software package HYDRUS-2D and MESHGEN-2D are necessary to generate
the HYDRUS2 input files selector.in, boundary.in, dimensio.in, domain.dat, andmeshtria.txt. These files are subsequently imported with the help of the FIM-GUI.
The user can also load and modify previous FIM simulation scenarios and has theoption to save both the new parameterisation and the simulation results.
A comprehensive analysis and visualization of the simulation results is also provided by the FIM-GUI. Fig. 3 shows the simulated pressure head distribution of an
example run at the furrow inlet at times before (left) and after (right) an irrigation
event.
Results
FIM is validated ondata from furrow irrigation experiments ofgrowingcornonloamy soil at the
130m long experimental plot at CEMAGREF
Montpellier, France (Nemeth,
2001, Mailhol, 2001,
Mailhol et al. 2001).
Fig. 4 shows thesimulated soil mois- Figure 4: Simulated soil moisture distribution at an upstream
ture distribution at an location of the Montpellier-experiment at 44, 90, 91 and 133
upstream location at days after sowing (das)various times of thegrowing season. The soil water storage, Ssoil, is determined by integration of thesoil moisture in the root zone.
Simulated and observed soil water
storage are in good agreement(Fig. 5). The simulated corn yieldof 12.8 t/ha matches perfectly withthe harvest of 12.7 t/ha corn.
Conclusions
FIM shows a clear superiority overwater balance models (VBM). HYDRUS2 simulates correctly the soil
water transport during a growing
Figure 5: Observed and simulated water storage
season of corn at Montpellier
in the root zone during the growing season of the(France) and provides detailed inMontpellier experimentformation about the quasi 3D water
distribution in the soil. FIM estimates irrigation performance criteria, such as irrigation efficiency and distribution
uniformity, more precisely than VBM. It can be applied for improving furrow irrigation design and management, for irrigation planning, for cost-benefit analysis, andfor estimation of sustainability of furrow irrigation systems.
76

References
Mailhol, J. C. (2001). Contribution a l´amelioration des pratiques d´irrigation a la raie par une
modelisation simplidiee a l´echelle de la parcelle et de la saison (in French. PhD Thesis.
Mailhol, J. C., Ruelle, P., and Nemeth, I. (2001). Impact of fertilisation practices on nitrogenleaching under irrigation. Irr. Sci., 20:139–147.
Nemeth, I. (2001). Devenir de l´azote sous irrigation gravitaire. Application au cas d´un perimetre
irrigue au mexique. (in French). PhD thesis, University of Montpellier II.
Novak, V. (1987). Estimation of soil water extraction patterns by roots. Agric. Water Mgt.
Schmitz, G. H. and Seus, G. J. (1992). Mathematical zero-inertia modeling of surface irrigation:
advance in furrows. Journal of Irrigation and Drainage Engineering, 118(1):1–18.
Simúnek, J., Sejna, M., and van Genuchten, M. T. (1996). Hydrus-2d. simulating water flow andsolute transport in two-dimensional variably saturated media. Version 1.0. User manual, IGWMCTPS 53 International Ground Water Modeling Center, Colorado School of Mines. IGWMC-TPS 53.
Wöhling, Th. (2005). Physically based modeling of furrow irrigation systems during a growing season. PhD Thesis, Faculty of Forest, Geo und Hydro Science, Dresden University of Technology.
ISBN 3-86005-481-3.
Wöhling, Th., Fröhner, A., Schmitz, G.H. and Liedl, R. (2005). Efficient solution of interacting 1D
surface - 2D subsurface flow during furrow irrigation advance. Submitted to the Journal of Irrigation
and Drainage Engineering.
Wöhling, Th., Singh, R., & Schmitz, G.H. (2004). Physically based modeling of interacting surface-subsurface flow during furrow irrigation. Journal of Irrigation and Drainage, 130(5), 296-303.

77

Modelling the effects of soil tillage on water and solute transport: theHYDRUS-2D/SISOL coupling
Coquet Y.1, Šimůnek J.2, Roger-Estrade J.3
1

UMR INRA/INAPG Environment and Arable Crops, Institut National de la Recherche
Agronomique/Institut National Agronomique Paris-Grignon, B.P. 01, 78850 ThivervalGrignon, France
2

Department of Environmental Sciences, University of California, Riverside, CA, 92521,
USA
3

UMR INRA/INAPG Agronomy, Institut National de la Recherche Agronomique/Institut
National Agronomique Paris-Grignon, B.P. 01, 78850 Thiverval-Grignon, France
Introduction
The topsoil properties play a central role in the dynamics and fate of pollutants in agriculturalsoils. It is in the surface soil layer that biological activity is most intense and that soil
variables controlling degradation, retention and transport of solutes are the most prone tovariation. The topsoil is time- and space-variable due to the effects of agricultural practices,
mainly tillage, but also fertilizer or manure application, crop harvesting, climate, andbiological activity. This variability must be accounted for if one aims at understanding andpredicting the impact of agricultural practices on the environment, esp. groundwater quality.
The model SISOL (Roger-Estrade et al., 2000) has been developed for predicting the variation
through time and space of the structural status of agricultural soils. SISOL has been designed
specifically for frequently mouldboard-ploughed soils. It has been tested in a variety ofagricultural situations representative of European intensive agricultural systems in loamy
soils. SISOL is a 2-D model, describing time changes and spatial variability of soil structurein the ploughed layer, taking into account fragmentation by tillage and climate, soil
displacement by ploughing and compaction due to traffic. Description of soil structure isbased on the morphological method proposed by Manichon (1987). This paper describes thecoupling of HYDRUS-2D with SISOL and presents some illustrations of the effect of soil
structure on water transfer and solute leaching.
Materials and Methods
The SISOL Model
The SISOL model predicts the structural status of cultivated soils by assessing the distributionof compacted zones or clods, called Δ zones, within the ploughed layer. These Δ zones aredefined by Manichon (1987) as soil zones having no eye-visible macroporosity. They result
from severe soil compaction under the wheel tracks of farming machinery. The extent ofcompacted Δ soil below the wheels is dependent on the soil water content and on thecharacteristics of the farming machine (tyre width, pressure and characteristics, axle load).
Basically, any agricultural operation involving machinery may produce Δ zones within thesoil. Empirical relationships are used by SISOL to predict the extent of compacted Δ zones
created by each agricultural operation (e.g. sowing, fertilizing, harvesting, etc). These Δ zones
may be destroyed within the upper part of the ploughed topsoil when it is superficial tilled(for seed bed preparation or stubble breaking) or by weathering (freezing/thawing and78

drying/wetting cycles). Finally, Δ zones may be split up and displaced by ploughing that cutsand rotates soil furrow slices (Figure 1). Additional details about the SISOL model may befound in Roger-Estrade et al. (2000).

Figure 1.
Modelling of the soil displacement
during ploughing (Roger-Estrade etal., 2000).

All these effects are modelled by SISOL that describes the time-evolution of the distributionof the Δ zones within the ploughed soil according to the sequence of farming operations,which is a characteristic of the considered cropping system. A high time-averaged proportionof Δ soil within the tilled layer is considered to have a detrimental effect on crop performance
because Δ soil zones impede crop rooting and have a low hydraulic conductivity (Coutadeur
et al., 2002) thereby limiting plant access to water and nutrients. The role of these Δ
compacted zones on soil water and solute dynamics was further explored by using HYDRUS2D and coupling it with the SISOL model (Coquet et al., 2005a, 2005b).
Coupling HYDRUS-2D with SISOL
The consequences of the presence of compacted Δ zones on water and solute transport intilled soils were first evaluated through a field experiment, involving the application of a
tracer step input at the surface of a loamy soil that had been prepared for maize sowing (ieafter an autumn ploughing and a spring harrowing), and a monitoring of water content, water
potential and tracer concentration within the tilled soil (Coquet et al., 2005a). The field
experiment was then modelled using HYDRUS-2D (Šimůnek et al., 1999) by accounting
explicitly for the locations, sizes and forms of the Δ zones within the ploughed layer. The Δ
material had different hydraulic parameters from the non-Δ material. Hydraulic properties of
all types of materials have been estimated from independent measurements. The observed
water and tracer behaviours were correctly described, including preferential flow in thevicinity of the compacted Δ zones (Coquet et al., 2005b).
To account for the time-variations of the structural status of tilled soils while modelling its
water and solute dynamics, we linked HYDRUS-2D to the SISOL model in the following
way. After each agricultural operation, HYDRUS-2D reads the output file generated by
SISOL, which describes the structure of the ploughed layer using a binary status (Δ or non-Δ)
79

for each 1 by 1 cm soil pixel, and translates it into a new HYDRUS-2D material distribution.
HYDRUS-2D then proceeds with the current simulation until the next agricultural operation.
All the information for the SISOL/HYDRUS-2D interaction (ie coordinates and size of
SISOL pixels, number and timing of agricultural operations and names of the SISOL output
files) is given in a specific HYDRUS-2D input file called SeSoil.in to be placed in thecorresponding HYDRUS-2D project.
We chose 3 experimental plots from a long term field trial where 3 cropping systems arecompared. This trial is located at Mons-en-chaussée in the northern part of the French ParisianBasin. Cropping system 1 (S1) is characterised by a high compaction risk (mainly due tofrequent harvesting in late autumn and seed bed preparation in early spring, when the soil
water content is high). In this treatment, the crop sequence is maize/winter wheat/sugar
beet/winter wheat. Cropping system 2 (S2) is characterised by the same crop sequence, but
the scheduling of cultivation operations is different, in order to minimise compaction (earlysugar beet and maize harvesting in autumn and late sowing of these crops in spring, when thesoil is dry enough). Crop sequence in cropping system 3 (S3) is rape/winter wheat/pea/winterwheat. In this region of France, harvesting of these crops is usually performed in summer,
when the soil is dry, and seed bed preparation is done in the middle of autumn (except for pea
sowing, which takes place in February) when the soil is still rather dry. Compaction risk inthis system is low. All the plots are mouldboard ploughed once every year, in September
before wheat, rape or pea, or in November before sugar beet and maize.
Simulation of the time changes in soil structure of the ploughed layer of these three plots
where performed using SISOL over a seven-year period (Roger-Estrade et al., 2000). The
simulation was stopped just after winter wheat sowing in the three plots, leading to three
SISOL output files. The corresponding material distributions are presented on the left part offigure 2.
Results and Discussion
The HYDRUS-2D/SISOL tool allows the exploration of the effects of agricultural systems on
the dynamics of water and solutes in soils. It clearly demonstrated that the abundance anddistribution of the compacted Δ zones have a strong effect on pore water velocity (Figure 2).

0-9.5

0-6.0

0-2.5

Figure 2. HYDRUS-2D material distribution (Δ soil zones in blue) and nodal fluxes for the 3
cropping systems S1 (top), S2 (middle) and S3 (bottom). Numbers below nodal flux maps
show flux range in cm/h. Size of the simulation domain is 0.5 by 3 m.

80

Figure 3. Solute concentration maps for the 3 cropping systems S1 (top), S2 (middle) and S3(bottom) and the same solute transport scenario: 1 h step-input of solute at a steady-statewater flow input rate of 5 mm/h followed by a 11 h elution at the same water flow rate withsolute-free rainfall. Size of the simulation domain is 0.5 by 3 m.
The larger the abundance and sizes of the Δ zones, the more heterogeneous is the flow field.
Such a flow dynamics has a strong impact on solute transport in soil (figure 3). Solute tends tomove faster in some soil regions, essentially in the vicinity of large Δ zones, and reaches thebottom of the soil much earlier in the S1 cropping system than in the S2 and S3 cropping
systems. That means that the S1 cropping system is likely to have a greater impact on
groundwater than the other two cropping systems, because a larger part of solute may movequickly out of the ploughed layer where the essential of degradation takes place. In thatcondition, the attenuation of pollutant concentrations through soil may be less effective.
The coupling between HYDRUS-2D and SISOL is quite loose at present. It could indeed be
improved by a bi-directional communication between the two softwares. The SISOL model
could use the water contents modelled by HYDRUS-2D to evaluate the extent of compaction
due to wheeling. However, we believe that the HYDRUS-2D/SISOL tool may already be
useful for the environmental risk assessment of cropping systems.
References
Coquet, Y., Coutadeur, C., Labat, C., Vachier, P., van Genuchten, M. Th., Roger-Estrade, J., Šimůnek, J.,
2005a. Water and solute transport in a cultivated silt loam soil; 1. Field observations, Vadose Zone J. 4:558572.
Coquet, Y., Šimůnek, J., Coutadeur, C., van Genuchten, M. Th, Pot, V., Roger-Estrade, J., 2005b. Water andsolute transport in a cultivated silt loam soil; 2. Numerical analysis, Vadose Zone J. 4:573-586.
Coutadeur, C., Coquet, Y., Roger-Estrade, J., 2002. Variation of hydraulic conductivity in a tilled soil. Europ. J.
Soil Sci. 53:619-628.
Manichon, H., 1997. Observation morphologique de l’état structural et mise en évidence d’effets de compactage
des horizons travaillés. In: Soil Compaction and Regeneration (eds G. Monnier, M.J. Goss), pp. 39-52. A.A.
Balkema, Rotterdam.
Roger-Estrade, J., Richard, G., Boizard, H., Boiffin, J., Caneill, J., Manichon, H., 2000. Modelling structuralchanges in tilled topsoil over time as a function of cropping systems. Europ. J. Soil Sci. 51:455-474.
Šimůnek, J., M. Šejna, and M.Th. van Genuchten. 1999. The HYDRUS-2D software package for simulating
two-dimensional 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, CO.

81

Inverse model simulation of Br and herbicide transport using multi-process transportmodels with HYDRUS-1D
V. Pota, J. Šimůnekb, P. Benoita, Y. Coquetaab

UMR INRA – INA P-G, Unité Environnement et Grandes Cultures, BP01, 78850 Thiverval-Grignon, France
Department of Environmental Sciences, University of California, Riverside, California, USA.

Abstract
The modified one-dimensional transport model HYDRUS-1D (Šimůnek et al., 1998, 2003) wasused to identify and quantify different physical and chemical non-equilibrium transport processes
acting in a grassed soil. The convective-dispersive model (CDE), the dual-porosity model (MIM),
the dual-permeability model (DP) and the triple-porosity, dual permeability model (DP-MIM)
were used to analyze experimental bromide and the herbicide isoproturon breakthrough curves
measured from two series of displacement experiments performed on two undisturbed grassed
filter strip soil cores, under unsaturated steady-state flow conditions. The CDE, MIM and DP
models were combined with both chemical instantaneous and kinetic sorption, while DP-MIM
model was combined with chemical instantaneous sorption. Several rainfall intensities (0.070,
0.147, 0.161, 0.308 and 0.326 cm h-1) were used. Two contrasting physical non-equilibriumtransport processes occurred. Multiple (three) porosity domains contributed to flow at the highest
rainfall intensities, including preferential flow through macropore pathways. Macropores were notactive any longer at intermediate and lowest velocities, and the observed preferential transport
was described using dual-porosity-type models with a zero or low flow in the matrix domain.
Chemical non-equilibrium transport of isoproturon was found at all rainfall intensities and thedegree of non-equilibrium sorption increased with increasing rainfall intensities. Significantly
higher estimated values of degradation rate parameters as compared to batch data were correlatedwith the degree of non-equilibrium sorption.
Material and methods
Two undisturbed soil cores (14-cm i.d., 30-cm height) were sampled from the soil surface layer ofgrassed strip (La Jaillière, Loire-Atlantique, France), planted with perennial rye-grass (Lolium
perenne). The soil is a hydromorphic silt loam (24% clay, 40% silt, 36% sand). The displacement
experiments were conducted at steady-state flow under unsaturated conditions and consisted ofinjecting a pulse of the CaCl2 solution containing actual concentrations of 66.1 ± 2.1 mg l-1 Br, 7.6± 1.5 mg l-1 isoproturon (IPU) and 8.5 ± 0.2 mg l-1 metribuzin (Met). Three and two successive
displacement experiments were performed at different rainfall intensities for column I and II,
respectively. The solute pulse duration was fixed so that an approximately equal mass of Br, IPU
and Met was injected at each velocity. Met results are not shown in this paper (for full data set,
see Pot et al., 2005).
The following transport processes were used to analyze the experimental breakthrough curves(BTCs): a) physical and chemical equilibrium transport described using the classical convectiondispersion transport model (Lapidus and Amundson, 1952); b) chemical non-equilibrium transport
using the two-site chemical model (Cameron and Klute, 1977); c) physical non-equilibriumtransport using the two region model (mobile and immobile water) (van Genuchten and Wierenga,1976) and the dual-permeability model (Gerke and van Genuchten, 1993); and finally, e)simultaneous physical and chemical non-equilibrium transport using the dual-permeability model(Gerke and van Genuchten, 1993) combined with the two-site chemical model. For all
simulations, the linear sorption isotherm was used, and degradation was considered in the liquid
phase only.
Physical transport parameters were first estimated by inverting the Br BTCs for both columns
using the CDE, MIM, DP and/or DP-MIM transport models. Then, chemical parameters for IPU
were estimated by inverting the herbicide BTCs using the physical transport model that best82

described Br transport combined with equilibrium or non-equilibrium chemical transport models.
Note that non-equilibrium chemical transport is currently not yet available in the HYDRUS-1Dmodel for the physical non-equilibrium DP-MIM transport model.
Results
Br transport
At the highest rainfall intensities, two arrival fronts and a marked shift to the left of peak positions
were observed by contrast to the other rainfall intensities suggesting contrasting physical transport
processes as a function of rainfall intensity (Figure 1). Similar features of the shape of the curves
were observed with simultaneous leaching of IPU and Br (same arrival time of the peaks) at thehighest rainfall intensities and delayed leaching of IPU was observed for the lower velocities(Figure 1). A preferential flow of a macropore type was thus suggested at the highest rainfall
intensities.
The CDE model did not describe satisfactorily the Br data for either column. For both columns,
the DP model, similarly to the MIM model, was not able to describe two arrival fronts that
occurred during the Br transport at the highest rainfall intensities, but predicted well the Br
transport at the lower velocities. Results found with MIM and DP models were quite similar sinceinstead of an immobile zone considered in the MIM model, the DP model predicted very slow
water flow in the matrix zone. However, in the case of the lowest rainfall intensity, the DP model
showed that the contribution of the matrix to the total flux can be significant.
The DP-MIM model succeeded in describing the transport of Br at the highest velocity.
Differences in pore velocities and water content in the fracture domain together with opposite
results found for the two mass transfer rate coefficients (between fracture and matrix domains on
the one hand, and between mobile and immobile regions of the matrix domain on the other hand)
between both columns explained rather well the observed differences between the two
experimental BTCs. The BTC for column I shows two solute fronts of very different
concentrations. The pore velocity in the fracture domain is very high and the mass transfer
between fracture and matrix domains is small. The long concentration tail is then explained bylarge exchange between mobile and immobile regions of the matrix domain. The BTC for column
II, on the other side, shows two solute peaks of similar magnitude and relatively smooth increases
and decreases in concentrations, which requires smaller pore velocity in the fracture domain andfaster exchange between the fracture and matrix domains.
IPU transport
Results showed that for all rainfall intensities and regardless of the used physical transport model
the chemical equilibrium transport model was not well suited to describe the transport ofherbicides in the grassed soil columns while the non-equilibrium chemical model described wellthe elution tails of the BTCs.
The fitted distribution coefficient (Kd ) values were lower than those derived from batch sorption
experiments with the same soil data (Madrigal et al., 2002) while the fitted degradation rates were
up to 31 times larger than the literature values (Agritox, 2001).
The Kd values increased when rainfall intensities decreased, with comparable values between thetwo columns while a significant decrease in the sorption kinetic rate was observed as velocity
decreased. Maraqa (2001), in a review of column displacement experiments, reported similar
trends. A possible explanation is that for lower velocities more sites seem to be at relativeequilibrium with the solution, and given the lower sorption kinetic rate the overall rate of reactionis then more or less the same. Since the flow in the matrix region is slower, there is more time toreach equilibrium in the matrix domain and thus the fraction of exchange sites assumed to be in
equilibrium with the solution phase of the matrix domain is not a relevant parameter. Equilibrium83

sorption in the matrix domain is thus a suitable model. By contrast, for the fracture domain where
there is not enough time to reach equilibrium, due to faster flow, non-equilibrium sorption is anappropriate model to explain the observed long elution tails.
The first-order degradation rate was found to increase when velocity increased. Our modeling
analysis suggested that the degree of non-equilibrium sorption increased as velocity increased.
Thus, our results suggest that the estimated degradation rates are positively related to the degreeof non-equilibrium sorption. It supports the hypothesis that chemical irreversible sorption, defined
as fast attachment and much slower detachment, occurred most probably on the sorption sites
where kinetic sorption occurs. The irreversible sorption acts as a disappearance process at the timescale of our experiments and is therefore accounted for by apparently higher degradation rates by
the models.
Conclusions
For the highest rainfall intensities, multiple porosity domains and multiple permeabilities
preferential flow were probably active: a) rapid flow through macropore pathways, b) slower flow
through a mesoporosity, and c) no-flow in the remaining micropores. For the lower rainfall
intensities, macropore flow was not active anymore, the preferential transport observed is of adual porosity type with a zero or low flow in the matrix domain.
The main conclusion drawn from the modeling analyses is that non-equilibrium chemicaltransport processes occurred for all velocities and that the degree of non-equilibrium sorption
increased with velocity. Significantly higher values of degradation rate parameters as compared tobatch data were positively related with the degree of non-equilibrium sorption. This work suggeststhat the sorption and/or degradation rates obtained from batch experiments can not be used in nonequilibrium transport models to describe the transport experiments unless the actual underlying
mechanisms measured by the two techniques are fully identified.
References

Agritox: toxixological data base on pesticides, INRA, France, 2001.
Cameron D.R., Klute A., 1977. Convective-dispersive solute transport with a combined equilibrium and kinetic
adsorption model. Water Resour. Res., 13, 183-188.
Gerke H.H., van Genuchten M. Th., 1993. A dual-porosity model for simulating the preferential movement of waterand solutes in structured porous media. Water Resour. Res., 29, 305-319.
Lapidus L., Amundson N.R., 1952. Mathematics of adsorption in beds. VI. The effect of longitudinal diffusion in ion
exchange and chromatographic columns. J. Phys. Chem., 56:, 84-988.
Madrigal I., Benoit P., Barriuso E., Etiévant V., Souiller C., Réal B., Dutertre A., 2002. Capacités de stockage etd’épuration des sols de dispositifs enherbés via-à-vis des produits phytosanitaires. Deuxième partie : Propriétés derétention de deux herbicides, l’isoproturon et le diflufénicanil dans différents sols de bandes enherbées. Etude etGestion des Sols, 9, 287-302.
Maraqa, M.A., 2001. Prediction of mass-transfer coefficient for solute transport in porous media. J. ContaminantHydrol., 50, 1-19.
Pot V., Šimůnek, J., Benoit P. Coquet Y., 2005. Impact of rainfall intensity on the transport of two herbicides inundisturbed grassed filter strip soil cores. J. Contaminant Hydrol., in press.
Šimůnek, J., Šejna M., van Genuchten M. Th., 1998. The HYDRUS-1D software package for simulating the onedimensional movement of water, heat, and multiple solutes in variably-saturated media. Version 2.0, IGWMC - 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 for describing
non-equilibrium and preferential flow and transport in the vadose zone. J. Hydrol., 272, 14-35.
van Genuchten M.Th., Wierenga P.J., 1976. Mass transfer studies in sorbing porous media I. Analytical solutions. Soil
Sci. Soc. Am. J., 40, 473-480.

84

Figure 1: Br and IPU BTCs measured and calculated for columns I and II at the highest rainfall intensity. Relative concentrations of the effluents are presented against the number of porevolumes.

85New Features of HYDRUS-1D, Version 3.0Jiří Šimůnek1, Rien van Genuchten2, S. A. Bradford2, and R. J. Lenhard3
1
University of California Riverside, Department of Environmental Sciences, Riverside, CA
92521, USA, jiri.simunek@ucr.edu
2George E. Brown, Jr. Salinity Laboratory, USDA-ARS, Riverside, CA 92507, USA
3
Idaho National Laboratory, P.O. Box 1625, Idaho Falls, ID 83415, USA
Introduction
Version 3.0 of HYDRUS-1D (Šimůnek et al., 2005), released in May 2005, implements thefollowing new features as compared to version 2.1 (Šimůnek et al., 1998): a) new approaches
to simulate preferential and nonequilibrium water flow and solute transport (Šimůnek et al.,2003), b) a new hysteresis module (Lenhard et al., 1991) that avoids the effects of pumping,c) compensated root water uptake (Jarvis, 1994), d) options to simulate the transport ofviruses, bacteria, and colloids (Bradford et al., 2003), e) carbon dioxide production andtransport (Šimůnek and Suarez, 1993), and f) geochemical transport of major ions (Šimůnek
and Suarez, 1994; Šimůnek et al., 1996, 1997). This paper briefly summarizes the newfeatures.
Preferential and Nonequilibrium Water Flow
The variably-saturated flow equation may consider dual-porosity type flow with a fraction of
the water being mobile, and a fraction immobile. The approach assumes that the soil matrix
can exchange, retain, and store water, but does not permit advective flow. This
conceptualization leads to two-region, dual-porosity type flow and transport models thatpartition the liquid phase into mobile (flowing, inter-aggregate), θm, and immobile (stagnant,intra-aggregate), θim, regions with possible exchange of water and/or solutes between the tworegions, usually calculated by means of a first-order rate equation.
The dual-porosity formulation for water flow as used in HYDRUS-1D is based on a mixed
formulation, which uses the Richards equation to describe water flow in the fractures(macropores), and a simple mass balance equation to describe moisture dynamics in thematrix as follows (Šimůnek et al., 2003):∂θ m ∂ ⎡⎛ ∂h⎞⎤= ⎢ K ( h ) ⎜ + cos α ⎟ ⎥ − Sm − Γ w
∂t
∂z ⎣
⎝ ∂x
⎠⎦
∂θim
=
∂t

(1)

− Sim + Γ w

where Sm and Sim are sink terms for the two regions, and Γw is the transfer rate for water fromthe inter- to the intra-aggregate pores.
The mass transfer rate, Γw, in (1) for water between the fracture and matrix regions in dualporosity studies (e.g., Šimůnek et al., 2003) is often assumed to be proportional to thedifference in effective saturations of the two regions using the first-order rate equation:

Γw =

∂θim
= ω ⎡⎣ Sem − Seim ⎤⎦
t
∂
86

(2)

where θim is the matrix water content, ω is a first-order rate coefficient (T-1), and Sem and Seimare effective fluid saturations of the mobile (fracture) and immobile (matrix) regions,
respectively. Equation (2) assumes that the mass transfer rate is proportional to the differencein effective saturation, rather than in the pressure head, which should provide a more realistic
description of the exchange rate between the fracture and matrix regions.
When the rate of exchange of water between the fracture and matrix regions is assumed to be
proportional to the difference in pressure heads between the two pore regions (Gerke and van
Genuchten, 1993), the coupling term, Γw, becomes:

Γ w = α w ( h f - hm )(3)

in which αw is a first-order mass transfer coefficient [L-1T-1].
Hysteretic Model Without Pumping Effects

While relatively simple to implement, the hysteretic model that was previously used in
HYDRUS-1D has been found to sometimes suffer from a so-called pumping effect that cancause the hysteresis loop to move to physically unrealistic parts of the retention function. As
an alternative, we incorporated in HYDRUS the hysteresis model of Lenhard et al. (1991) thateliminates pumping by keeping track of historical reversal points.
Compensated Root Water Uptake

The ratio of the actual, Ta [LT-1], to potential, Tp [LT-1], transpiration rates for root uptake
without compensation is defined as follows:Ta
1=
Sdz = ∫α (h,hφ , z ) b( z )dz = ω
Tp Tp L∫R
LR

(4)

where α(h) is the water stress response function [-], b(x) is the normalized water uptake
distribution [L-1], LR is the rooting depth [L], S is the root water uptake rate [T-1], and ω is adimensionless water stress index (Jarvis, 1994). Following Jarvis (1994), we introduce a
critical value of the water stress index ωc, referred to as the root adaptability factor, which
represents a threshold value above which reductions in root water uptake in stressed parts of
the root zone are fully compensated by increased uptake from other parts:∫α (h,hφ , z, t )b( z, t )dz
Tac (t ) LR
ω (t )=≤1=
Tp (t )
max [ω (t ), ωc ]max [ω (t ), ωc ]S (h,hφ ,z , t ) = α (h,hφ , z , t )b( x,z , t ) Lt

(5)

Tp (t )

max [ω (t ), ωc ]where Tac is the actual transpiration rate for compensated root uptake [LT-1]. When theparameter ωc is equal to one we hence have noncompensated root water uptake, and when ωcis equal to zero fully compensated uptake results.

87

Attachment/Detachment Transport Model

Virus, colloid, and bacteria transport and fate models commonly employ a modified form of
the advection-dispersion equation:
∂s
∂s
∂θ c
∂ ⎛
∂c ⎞ ∂qc+ρ 1 +ρ 2 =- µ wθ c - µ sρ ( s1 + s2 )⎜θ D ⎟ −
∂t
∂t
∂t
∂x ⎝
∂x ⎠ ∂x(6)

where c is the (colloid, virus, bacteria) concentration in the aqueous phase [NcL-3], s is thesolid phase (colloid, virus, bacteria) concentration [NcM-1], subscripts 1 and 2 represent two
kinetic sorption sites, Nc is the number of colloids (particles), and µw and µs represent
inactivation and degradation processes in the liquid and solid phases, respectively.
Mass transfer between the aqueous and two solid kinetic phases can be described withequations of the form (note that we dropped the subscripts 1 and 2):ρ∂s= θ k aψ c − kd ρ s∂t(7)where ka is the first-order deposition (attachment) coefficient [T-1], kd is the first-order
entrainment (detachment) coefficient [T-1], and ψ is a dimensionless colloid retention function[-]. Various functions were implemented for ψ to describe blocking, ripening, or straining
phenomena. The attachment coefficient can be calculated using filtration theory.
Carbon Dioxide Production and Transport

The HYDRUS-1D software package now also includes modules for simulating carbon
dioxide production and transport. CO2 transport is assumed to occur by diffusion in both theliquid and gaseous phases, and by advection in the liquid phase. CO2 production is considered
to be the sum of production by soil microorganisms and by plant roots. Various
environmental factors, such as temperature and the water content are used to modify the
optimum CO2 production rate.
Geochemical Transport of Major Ions

We also incoporated a new geochemical transport module in HYDRUS-1D based on the
UNSATCHEM model of Šimůnek et al. (1996, 1997). The major variables of the chemicalsystem are Ca, Mg, Na, K, SO4, Cl, NO3, H4SiO4, alkalinity, and CO2. The model accounts
for equilibrium chemical reactions between these components such as complexation, cationexchange and precipitation-dissolution. For the precipitation-dissolution of calcite anddissolution of dolomite, either equilibrium or multicomponent kinetic expressions can be usedwhich include both forward and back reactions. Other dissolution-precipitation reactions
considered include gypsum, hydromagnesite, nesquehonite, and sepiolite. Since the ionicstrength of soil solutions can vary considerably in time and space and often reach high values,
both modified the Debye-Hückel and the Pitzer expressions were incorporated into the model
as options to calculate single ion activities.

88

The new UNSATCHEM module enables quantitative predictions of processes involving
major ions, such as simulations of the effects of salinity on plant growth and estimating theamount of water and amendment required to reclaim soil profiles to desired levels of salinity
and ESP (exchangeable sodium percentage).
Table 1. Chemical species included in the UNSATCHEM major ion module of HYDRUS-1D.
Aqueous components

7Complexed species

10

Precipitated species

6

Sorbed (exchangeable)
species
CO2-H2O species
Silica species

Ca2+, Mg2+, Na+, K+, SO42-, Cl-, NO3CaCO3o, CaHCO3+, CaSO4o, MgCO3o, MgHCO3+,
MgSO4o, NaCO3-, NaHCO3o, NaSO4-, KSO4CaCO3, CaSO4⋅ 2H2O, MgCO3⋅ 3H2O,
Mg5(CO3)4(OH)2⋅ 4H2O,
Mg2Si3O7.5(OH) ⋅ 3H2O, CaMg(CO3)2

4

Ca, Mg, Na, K

7
3

PCO2, H2CO3*, CO32-, HCO3-, H+, OH-, H2O
H4SiO4, H3SiO4-, H2SiO42-

References
Bradford, S. A., J. Šimůnek, M. Bettehar, M. Th. van Genuchten, and S. R. Yates, Modeling colloid
attachment, straining, and exclusion in saturated porous media, Environ. Sci. & Technology, 37(10),2242-2250, 2003.
Gerke, H. H., and M. Th. van Genuchten, A dual-porosity model for simulating the preferential
movement of water and solutes in structured porous media, Water Resour. Res., 29, 305-319,
1993.
Jarvis, N.J., The MACRO model (Version 3.1), Technical description and sample simulations. Reports
and Dissertations 19. Dept. Soil Sci., Swedish Univ. Agric. Sci., Uppsala, Sweden, 51 pp., 1994.
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.
Šimůnek, J., and D. L. Suarez, Modeling of carbon dioxide transport and production in soil: 1. Modeldevelopment, Water Resour. Res., 29(2), 487-497, 1993.
Šimůnek, J., and D. L. Suarez, Two-dimensional transport model for variably saturated porous media withmajor ion chemistry, Water Resour. Res., 30(4), 1115-1133, 1994.
Šimůnek, J., D. L. Suarez, and M. Šejna, 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, 186pp.,
1996.
Šimůnek, J., and D. L. Suarez, Sodic soil reclamation using multicomponent transport modeling, ASCE J.
Irrig. and Drain. Engineering, 123(5), 367-376, 1997.
Šimůnek, J., K. Huang, M. Šejna, and M. Th. van Genuchten, The HYDRUS-1D software package forsimulating the one-dimensional movement of water, heat, and multiple solutes in variably-saturated
media. Version 1.0, IGWMC - TPS - 70, International Ground Water Modeling Center, Colorado
School of Mines, Golden, Colorado, 186pp., 1998.
Šimůnek J., N. J., Jarvis, M. Th. van Genuchten, and A. Gärdenäs, Review and comparison of models
for describing 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 ofCalifornia Riverside, Riverside, CA, 270pp., 2005.

89