Journal of Hydrology 407 (2011) 129–144

Contents lists available at ScienceDirect

Journal of Hydrology
journal homepage: www.elsevier.com/locate/jhydrol

Field evaluation of a multicomponent solute transport model in soils irrigatedwith saline waters
T.B. Ramos a,⇑, J. Šimůnek b, M.C. Gonçalves a,c, J.C. Martins c, A. Prazeres c, N.L. Castanheira c, L.S. Pereira aa

CEER – Biosystems Engineering, Institute of Agronomy, Technical University of Lisbon, Tapada da Ajuda, 1349-017 Lisbon, Portugal
Department of Environmental Sciences, University of California, Riverside, CA 92521, USAc
Estação Agronómica Nacional, L-INIA, Instituto Nacional de Recursos Biológicos, 2784-505 Oeiras, Portugalb

a r t i c l e

i n f oArticle history:
Received 27 December 2010
Received in revised form 13 April 2011
Accepted 13 July 2011
Available online 23 July 2011
This manuscript was handled by Philippe
Baveye, Editor-in-Chief, with the assistanceof Wilfred Otten, Associate Editor
Keywords:
HYDRUS-1D
Field experiment
Modeling
Soil salinity
Osmotic stress
Nitrogen leachings u m m a r y
Soil salinization, sodiﬁcation, and non-point source pollution are among the most important and widespread environmental problems in agricultural regions with scarce water resources. Models evaluatingthese environmental problems should therefore consider an integrated approach to avoid favoring one
problem over the other. The HYDRUS-1D software package was used to simulate water movement andsolute transport in two complex experiments carried out under ﬁeld conditions in Alvalade and Mitra,
Portugal. The experiments involved irrigating maize with synthetic saline irrigation waters blended withfresh irrigation waters and waters with different nitrogen concentrations. The major ion chemistry
module of HYDRUS-1D was used to model water contents (RMSE 6 0.04 cm3 cm3), the overall salinity
given by the electrical conductivity of the soil solution (ECsw) (RMSE 6 2.35 dS m1), the concentrationof soluble cations Na+ (RMSE 6 13.86 mmol(c) L1), Ca2+ (RMSE 6 5.66 mmol(c) L1), Mg2+ (RMSE 6
4.16 mmol(c) L1), and SAR (RMSE 6 6.27 (mmol(c) L1)0.5) in different experimental plots. RMSE werealways lower for the soil with coarse texture of Mitra. The standard HYDRUS solute transport module1
1) and N–NO
) concentrawas used to model N–NHþ
3 (RMSE 6 2.60 mmol(c) L
4 (RMSE 6 0.07 mmol(c) Ltions in the soil solution while either including or neglecting the effects of the osmotic stress on nutrient
uptake. The model was able to successfully simulate root water and nutrient uptake reductions due toosmotic stress. Consequently, modeled ﬂuxes of N–NHþ
4 and N–NO3 leached from the soil proﬁlesincreased due to the effects of the salinity stress on water and nutrient uptake. Possible causes of
disagreements between the modeling and experimental data are discussed. HYDRUS-1D proved to be a
powerful tool for analyzing solute concentrations related to overall soil salinity and nitrogen species.
Ó 2011 Elsevier B.V. All rights reserved.

1. Introduction
In recent decades, countries located in regions with arid, semiarid, and even sub-humid conditions have developed irrigation
areas to satisfy the growing demands for food and raw materials.
In regions with scarce water resources, it is not always possible
to irrigate with waters of good quality, and therefore even saline
waters are seen as an important resource (Pereira et al., 2009).
As a result, human-induced salinization and sodiﬁcation are among
the most important and widespread soil degradation processes,
generally associated with irrigation practices and poor water
management. As rainfall may not be sufﬁcient to remove the salts
accumulated during irrigation, it is often necessary to monitor soil
and water quality, and to promote efﬁcient leaching management
to counteract soil salinization. Irrigation volume and frequency are
therefore adjusted to remove salts from the root zone (e.g.,⇑ Corresponding author. Tel.: +351 214403500; fax: +351 214416011.
E-mail address: tiago_ramos@netcabo.pt (T.B. Ramos).
0022-1694/$ - see front matter Ó 2011 Elsevier B.V. All rights reserved.
doi:10.1016/j.jhydrol.2011.07.016

Gonçalves et al., 2007; Pereira et al., 2007). However, promoting
salt leaching also results in increased nitrogen leaching, which
contributes to the degradation of groundwater reservoirs. As
non-point source pollution from agricultural fertilization is considered to be a leading cause of water quality problems, especially inregions like those mentioned above, an integrated approach should
always be considered to avoid emphasizing one problem over theother.
Over the last few decades the scientiﬁc community has invested
considerable time and resources in the development of analyticaland numerical models as tools to predict the long and short-term
effects of irrigation water quality on soil properties, crop yield,
groundwater, and the environment (Wagenet and Hutson, 1987;
Jarvis, 1994; van Dam et al., 1997; Ahuja et al., 2000; van den Berget al., 2002; Šimůnek et al., 2006, 2008a). Most vadose zone models
are based on the numerical solution of the Richards equation forvariably-saturated water ﬂow, and on analytical or numericalsolutions of the Fickian-based convection–dispersion equation forsolute transport. A sink term is usually included in these equations

130

T.B. Ramos et al. / Journal of Hydrology 407 (2011) 129–144

to account for root water and nutrient uptake, and the effects ofwater and osmotic stresses (Feddes and Raats, 2004; Šimůnek
and Hopmans, 2009).
Perhaps the most signiﬁcant difference among vadose zone
models is how they describe solute partitioning between the liquidand solid phases and the complex chemical reactions involved insolute transport. Most vadose zone models consider the transportof only one solute, and severely simplify various chemical interactions. The relatively complex processes of adsorption and cationexchange are often accounted for by means of empirical adsorption
isotherms. Other processes such as precipitation/dissolution andbiodegradation are frequently ignored or simulated by invoking
simpliﬁed ﬁrst- or zero-order rate equations. Only a few models
have been developed that can consider multiple solutes and theirvarious interactions, such as precipitation/dissolution and competition for sorption sites (Šimůnek and Valocchi, 2002).
The HYDRUS-1D software package (Šimůnek et al., 2008a) uses
several modeling concepts for evaluating solute transport, two of
which will be used in our study. The standard solute transportmodule, already available with earlier versions of HYDRUS, considers the transport of one or multiple solutes, which can be either
independent or involved in sequential ﬁrst-order decay reactions.
This module has been used for a wide range of applications in research and irrigation management of poor quality waters (e.g.,
Forkutsa et al., 2009; Hanson et al., 2008; Roberts et al., 2008,2009). It has also been used to simulate the fate of nutrients in soils
by evaluating and comparing fertilization strategies for different
crops (e.g., Cote et al., 2003; Gärdenäs et al., 2005; Hanson et al.,2006; Ajdary et al., 2007; Crevoisier et al., 2008). While salinity
studies have typically been limited to estimating the electrical conductivity of the soil solution (ECsw), the studies addressing fertilization strategies have been mainly theoretical, or when experimental
data existed, limited to very short periods of time, while evaluating
nutrient leaching independently of irrigation water quality. As far
as we know, no integrated studies using HYDRUS exist where bothproblems have been considered at the same time.
The major ion chemistry module adapted from the UNSATCHEM model (Šimůnek et al., 1996) has been used much less often
than the standard solute transport module. We refer here to thecomprehensive study carried out by Gonçalves et al. (2006), whoanalyzed transient water ﬂow and solute transport in three soil
lysimeters irrigated with waters of different quality over 4 years.
The study made full use of the major ion chemistry module, whichallowed comparisons between model-simulated and experimental
soil water contents, concentrations of Na+, Ca2+, and Mg2+, ECsw, thesodium adsorption ratios (SAR), and the exchangeable sodium percentages (ESP). This study demonstrated how the major ion chemistry module of UNSATCHEM is a better tool for evaluating saline
water management and for assessing the effects of irrigation water
quality on groundwater recharge (Gonçalves et al., 2006).
While it has been shown that these two solute transport modules of HYDRUS-1D can be helpful tools for developing management irrigation strategies in regions with water scarcity, there is
still a need to conduct studies where model predictions are compared against ﬁeld experimental data in order to lend greater credibility to both the simulations and extrapolations to different soil
types, crops, climatic conditions, tillage operations, and water
management schemes. As arid and semi-arid regions are particularly vulnerable to soil salinization and non-point source pollution,
there is also a need for these studies to consider an integrated
approach, as the solution of one problem can easily aggravate theother problems.
The objective of this study was to use the HYDRUS-1D software
package to evaluate soil salinization and sodiﬁcation risks byanalyzing experimental data collected in two agricultural ﬁelds(Ramos et al., 2009) irrigated with synthetic saline waters, and toquantify the effects of the salinity stress on nitrogen leaching.
We further evaluated the effectiveness of HYDRUS-1D (i) to predictwater contents and ﬂuxes, (ii) to predict salinization and sodiﬁcation risks by estimating the overall salinity given by ECsw, individual cations Na+, Ca2+, and Mg2+, and SAR, (iii) to quantify water
uptake reductions due to the use of saline waters, and (iv) to
predict N- NHþ
4 and N–NO3 concentrations in the soil and theirleaching under ﬁeld conditions.
2. Material and methods
2.1. Field experiment
Since the experimental data used in this study was presented byRamos et al. (2009), only information relevant to our modelingstudy will be given here.
2.1.1. Experimental ﬁelds
Two ﬁeld plot experiments were conducted from June 2004 to
February 2007 in Southern Portugal. One was a soil classiﬁed as Typic Xeroﬂuvent with a medium texture in the Alentejo region at
the Alvalade Experimental Station (37° 560 4800 N and 8° 230
4000 W), and the other a Typic Haplanthrept soil with a coarse texture found in Herdade da Mitra (38° 310 5500 N and 8° 000 5900 W)(Soil Survey Staff, 2006). The climate in this region is mostly dry
sub-humid to semi-arid, with hot dry summers, and mild winterswith irregular rainfall. In both experimental ﬁelds, maize was irrigated with synthetic saline waters, obtained by adding NaCl to theavailable irrigation water in the region. The irrigation water was
also used to deliver NH4NO3 to the crop.
2.1.2. Experimental design and treatments
A triple emitter source irrigation system was used to deliver
water, salts (e.g., Na+), and fertilizer (N) to the crop. This system,
adapted from Beltrão et al. (2002), consisted of three trickle laterals placed along each maize line (Fig. 1). The ﬁrst lateral was connected to the salt stock solution, the second to the nitrogen
reservoir, and the third to the source of fresh water. This last lateral
was used to obtain a constant water application rate at each dripping point. Gradients of applied salt and nitrogen concentrations
were produced by varying discharge rates at different laterals. Table 1 lists the different discharge rates of the emitters applying
salts, nitrogen, and fresh water in each experimental plot while
maintaining an overall constant cumulative discharge of 18 L/h/m (24 mm/h).
Each experimental ﬁeld was divided into four groups (I–IV)
with three triple joint laterals each, establishing a N gradientdecreasing from group I to IV. Each group was then divided into
3 sub-groups, A–C each with a surface area of 6.75 m2 (2.25 m
wide  3 m long; 0.75 m between maize lines), and the Na+ gradient decreasing from A to C. Synthetic saline waters (with EC varying between 7.5 and 14.6 dS m1) were blended with the locallyavailable water (EC 6 1.2 dS m1). The amount of nitrogen applied
to various plots varied from 0 to 65 g m2 y1, while the amount of
Na+ varied from 0 to 2800 g m2 y1. Table 2 presents the total
amount of water applied in each year at both experimental ﬁelds,
independent of its quality. A more detailed description of thetrickle irrigation scheme used can be found in Ramos et al. (2009).
2.1.3. Observations and analysis
In the plots with the highest application of the synthetic saline
waters (sub-group A), and in the plots irrigated only with thelocally available water (sub-group C), TDR probes and ceramic cups
were installed at 20, 40, and 60 cm depths to measure soil watercontents and collect soil solutions. The soil solution was monitored

T.B. Ramos et al. / Journal of Hydrology 407 (2011) 129–144

131

3m

16.5 m3m

2.25 m3m

0.75 mFig. 1. Layout of the triple emitter source design (adopted from Ramos et al. (2009)). The salt gradient decreases from sub-group A to C and the fertilizer gradient decreases
from group I to IV.

Table 1
Discharge rates of the laterals applying salt (Na+), nitrogen (N) and fresh water (W) ineach experimental plot. There is an overall constant cumulative discharge at eachdripping point of 18 L/h/m (Ramos et al., 2009).
Treatment

Application rates (L/h/m)
Group I

A
B
C

Group IIGroup III

Group IV

Na+

N

W

Na+

N

W

Na+

N

W

Na+

N

W

12
6
0

6
6
6

0
6
12

12
6
0

4
4
4

2
8
14

12
6
0

2
2
2

4
10
16

12
6
0

0
0
0

6
12
18

Table 2
Total amount of water applied (fresh water + saline water + water with fertilizer) in
Alvalade and Mitra during the three irrigation seasons.
Experimental ﬁeld

Alvalade
Mitra

Water applied (mm)
2004

2005

2006

997
1067

1012
725

1028
729

þ
for the concentrations of soluble Na+, Ca2+, Mg2+, N–NO
3 ; N–NH4 ,ECsw, and SAR. Soil water content measurements and soil solutions
were taken twice a week during irrigation seasons, generally 24 hafter an irrigation event, and twice a month during the rest of theyear. Large periods with no data were due to low soil water contents, which made it impossible to collect soil solutions with theceramic cups used in our experiments.
Modeling of water ﬂow and solute transport was carried out for
the worst case scenario, where potential root water uptake was reduced due to water and osmotic stresses, increasing the risk of salt
and nitrogen leaching to the groundwater (group I, sub-group A,
i.e., plot I-A). Results for this worst-case scenario were then compared with results for plots where it was assumed that the potential root water uptake was affected only by the water stress, andthat osmotic stress could be neglected (groups I and IV, sub-groups
C, i.e., plots I-C and IV-C).
We used the HYDRUS-1D software package to simulate onedimensional water ﬂow and solute transport in the plots where
the observations were carried out. The major ion chemistry module
UNSATCHEM was used to simulate water contents, overall salinity
given by ECsw, concentrations of individual ions (Na+, Ca2+, andMg2+), and SAR. The standard HYDRUS solute transport module,
run for every plot in sub-groups A and C, was used to simulate
water contents, ECsw, and concentrations of the nitrogen species.
Potential and actual root water uptakes obtained using the two
HYDRUS-1D modules were compared for each experimental plot.
2.2. HYDRUS-1D simulation model
The HYDRUS-1D software package (Šimůnek et al., 2008a)
numerically simulates one-dimensional (1D) water ﬂow, solute,and heat transport in variably-saturated porous media. Althoughother HYDRUS models (Šimůnek et al., 2008b) can simulate twoand three-dimensional transport, the 1D version simpliﬁes the processes involved in our experiments by neglecting water and solute
ﬂuxes, and pressure head and concentration gradients in the horizontal direction. Below we give an overview of the main HYDRUS1D processes that were involved in our experiments.
2.2.1. Water ﬂow
Variably-saturated water ﬂow is described using the Richards
equation:



@h
@
@h
¼
KðhÞ
 KðhÞ  Sðz; tÞ
@t @z
@z

ð1Þ

where h is the volumetric soil water content (L3 L3), t is time (T), z
is the vertical space coordinate (L), h is the pressure head (L), K is
the hydraulic conductivity (L T1), and S is the sink term accounting
for water uptake by plant roots (L3 L3 T1). The unsaturated soil

132

T.B. Ramos et al. / Journal of Hydrology 407 (2011) 129–144

hydraulic properties are described using the van Genuchten–
Mualem functional relationships (van Genuchten, 1980).
2.2.2. Root water uptake
The sink term, S, is calculated using the macroscopic approach
introduced by Feddes et al. (1978). In this approach, the potential
transpiration rate, Tp (L T1), is distributed over the root zone using
the normalized root density distribution function, b(z, t) (L1), andmultiplied by the dimensionless stress response function,a(h, h/ , z, t), accounting for water and osmotic stresses (Feddes
et al., 1978; van Genuchten, 1987; Šimůnek and Hopmans, 2009):

Sðh; h/ ; z; tÞ ¼ aðh; h/ ; z; tÞSp ðz; tÞ ¼ aðh; h/ ; z; tÞbðz; tÞT p ðtÞ

ð2Þ

where Sp(z, t) and S(h, h/ , z, t) are the potential and actual volumesof water removed from unit volume of soil per unit of time
(L3 L3 T1), respectively, and a(h, h/ , z, t) is a prescribed dimensionless function of the soil water (h) and osmotic (h/ ) pressure
heads (0 6 a 6 1). The actual transpiration rate, Ta (L T1), is thenobtained by integrating Eq. (2) over the root domain LR:

Ta ¼

Z

Sðh; h/ ; z; tÞdz ¼ T p

LRZaðh; h/ ; z; tÞbðz; tÞdz

ð3Þ

LR

a1 ðhÞ ¼h2 h1>1;>>>: hh4 ;h3 h4h > h1 or h 6 h4
h2 < h 6 h1

ð4Þ

h3 < h 6 h2
h4 < h 6 h3

where h1, h2, h3, and h4 are the threshold parameters. Water uptake
is at the potential rate when the pressure head is between h2 and h3,
drops off linearly when h > h2 or h < h3, and becomes zero whenh < h4 or h > h1. Soil water pressure head parameters are available
for maize in HYDRUS-1D internal database based on the work of
Wesseling et al. (1991).
Root water uptake reduction due to salinity stress, a2(h/ ), wasdescribed using the Maas’s (1990) threshold and slope function.
The threshold-slope salinity stress model is implemented in thestandard HYDRUS and UNSATCHEM solute transport modules as:

8
1;
>
>
>
< 1  ðEC  EC Þ0:01sT
a2 ðh/ Þ ¼
> or
>
>
:
1 þ ðh/  h/T Þsh/ ¼ 3:8106EC sw þ 0:5072

ðR2 ¼ 0:997Þ

h/ ¼ 3:8143EC sw þ 0:6990

ðR2 ¼ 0:996Þ

ð6Þ

respectively. Eq. (6) were determined by ﬁtting a line to the computed values of ECsw and h/ obtained in the UNSATCHEM simulations. They are very similar to the relationship reported by the US
Salinity Laboratory Staff (1954) for estimating the osmotic pressureof soil solutions from EC measurements (h/ ¼ 3:7188EC; if converted to S.I. units). HYDRUS-1D provides a database with thethreshold-slope salinity parameters for different plants, including
maize, based on work by Maas (1990). However, these values are given for the electrical conductivity of the saturation extract (ECe). In
our applications, ECe values ﬁrst had to be converted into ECsw:

EC e  kEC ¼ EC sw

In our simulations, we assumed that the potential root water
uptake was reduced due to water stress in all experimental plots.
Water stress is a function of the adopted irrigation schedule, which
may lead to insufﬁcient (during early growth stages to allow thedevelopment of the rooting system) or excessive (due to over-irrigation during the remaining growth periods) supply of water to thecrop. We also assumed that potential root water uptake was further reduced by osmotic stress resulting from the use of saline
waters in the experimental plots irrigated with the synthetic saline
waters (sub-group A in Fig. 1). We further assumed that the effectsof the water and salinity stresses were multiplicative, i.e.,a(h, h/ ) = a1(h)a2(h/ ) (van Genuchten, 1987), so that different
stress response functions (as described below) could be used forthe water and salinity stresses. It can be shown (e.g., Oster et al.,private communication) that the combined effect of the two stresses is larger when the multiplicative (compared to additive) approach is considered.
Root water uptake reduction due to water stress, a1(h), was described using the model developed by Feddes et al. (1978):

8
0;
>
>
>
>
< hh1 ;which root water uptake occurs without a reduction, h/T is the corresponding threshold value given in terms of the osmotic head (L),and s and s are the slopes determining root water uptake decline
per unit increase in salinity (in standard HYDRUS) or osmotic head(in UNSATCHEM) above or below the threshold, respectively. For
the Alvalade and Mitra soils we found the relationship between
ECsw and h/ to be:

EC 6 EC T or h/ P h/T
EC > EC T
orh/ < h/T
ð5Þ
1

respectively, where ECT is the salinity threshold (dS m ), whichcorresponds to the value of the electrical conductivity (EC), below

ð7Þ

where kEC is the ratio of the electrical conductivity of the in situ soil
water at ﬁeld capacity and the EC of the soil water in the saturation
extract. In our study we assumed a kEC of 2, which is a common
approximation (e.g., Ayers and Westcot, 1985) used for soil watercontents near ﬁeld capacity in medium-textured soils. Further details on kEC variability can be found in Skaggs et al. (2006).
2.2.3. Solute transport
The partial differential equations governing one-dimensional
advective-dispersive chemical transport under transient ﬂow in avariably-saturated rigid porous medium are deﬁned in HYDRUS1D as:



@hck
@ck
@
@ck
@qck
hD
þq
¼

þ /k  Scr;k
@z
@t
@t
@z
@z

ð8Þ

where h is the volumetric water content (L3 L3), c, c and cr are solute concentrations in the liquid phase (M L3), solid phase (M M1),and sink term (M L3), respectively, q is the soil bulk density(M L3), q is the volumetric ﬂux density (L T1), D is the hydrodynamic dispersion coefﬁcient (L2 T1), / represents chemical reactions of solutes involved in a sequential ﬁrst-order decay chain,such as nitriﬁcation of nitrogen species (M L3 T1), and subscript
k represents chemical species present in our study (e.g., Na+, Ca2+,and other major ions in the UNSATCHEM module and EC, N–NO
3,
and N–NHþ
4 in the standard HYDRUS solute transport module).
The last term of Eq. (8) represents a passive root nutrient uptake(Šimůnek and Hopmans, 2009).
The parameter / in Eq. (8), which is only considered in the standard HYDRUS solute transport module, represents nitriﬁcation ofþ
the N–NHþ
4 species to N–NO3 , and appears in Eq. (8) for N–NH4

and N–NO3 species as follows, respectively:/N—NHþ ¼ /N—NO3 ¼ lw;N—NHþ hcN—NHþ  ls;N—NHþ qcN—NHþ
4

4

4

4

4

ð9Þ

where lw and ls are the ﬁrst-order rate constants for solutes in theliquid and solid phases (T1), respectively. In our study, we consid
ered only the nitriﬁcation process from N–NHþ
4 to N–NO3 , which resulted from the application of NH4NO3 fertilizer. Other reactions,such as the nitriﬁcation from N–NO
2 to N–NO3 , the volatilizationþof N–NHþ
andsubsequent
N–NH
transportbygaseous diffusion,
4
4
mineralization of crop residues and soil humus, and the denitriﬁcation of N–NO
3 into N–N2 or N–N2O, were neglected. Some of thesereactions, such as mineralization of crop residues and soil humus,
simply cannot be described with sequential ﬁrst-order decay chain

T.B. Ramos et al. / Journal of Hydrology 407 (2011) 129–144

reactions, while others occur at a rate so fast that they are often
lumped, such is the case of the nitriﬁcation from N–NO
2 to N–NO
(e.g.,
Hansonetal.,2006).
3
The two different HYDRUS-1D modules used here apply two
different approaches for relating solutes in the liquid and solid
phases. The major ion chemistry module takes into account the factthat the soil liquid phase always contains a mixture of many ionsthat may interact, create complex species, precipitate, dissolve,
and/or compete with each other for sorption sites on the solid phase(van Genuchten and Šimůnek, 2004). Thus, the UNSATCHEM module considers these interactions, including aqueous complexation,
precipitation/dissolution, and cation exchange, described using
the Gapon exchange equations (Šimůnek and Suarez, 1994).
The standard HYDRUS solute transport module accounts for the
relatively complex processes of adsorption and cation exchange bymeans of empirical linear or nonlinear adsorption isotherms. In ourapplication, the adsorption isotherm relating c and c in Eq. (8) is
described using the following linear equation:

ck ¼ K d;k ck

ð10Þ
3

1

where Kd,k (L M ) is the distribution coefﬁcient of a chemical
species k.
The electrical conductivity of the soil solution (ECsw) is determined in the UNSATCHEM module from individual anions andcations following the method of McNeal et al. (1970), while inthe standard HYDRUS solute transport module, ECsw was run as
an independent solute, available only in the liquid phase (i.e.,
Kd = 0 cm3 g1). SAR, which is only determined in the UNSATCHEM
module, was calculated as follows:Naþ
SAR ¼ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
ðCa2þ þMg2þ Þ2

ð11Þ

2.2.4. Root nutrient uptake
The parameter cr in the last term of Eq. (8) is the dissolved
nutrient concentration taken up by plant roots in association withroot water uptake, and is deﬁned as:

cr ðz; tÞ ¼ min½cðz; tÞ; cmax ð12Þ

where cmax is the a priori deﬁned maximum concentration of theroot uptake. We considered unlimited passive nutrient uptake fornitrogen species, which means that cmax was set to a larger concentration value than the dissolved concentrations, c, allowing all
dissolved nutrients to be taken up by plant roots, and zero uptake
for other species (EC, major ions), which means that cmax was set
to zero. Since root nitrogen uptake likely involves both passive
and active mechanisms (e.g., Šimůnek and Hopmans, 2009), considering only passive uptake will likely underestimate the total N
uptake. By integrating passive nutrient uptake over the root
domain, LR, we obtained an equation similar to Eq. (3), given as:

Pa ðtÞ ¼ T p ðtÞZaðh; h/ ; z; tÞbðz; tÞ min½cðz; tÞ; cmax dz

ð13Þ

LR

where Pa is the passive root nutrient uptake for the whole root
domain (M L2 T1) (Šimůnek and Hopmans, 2009).
2.3. Statistical analysis
In addition to a visual check, ﬁeld measured values werecompared with the results of the HYDRUS-1D simulations using
the mean absolute error and the root mean square error. The mean
absolute error (MAE) given by

MAE ¼

N
1X
jOi  P i j
N i¼1

ð14Þ

133

describes the difference between observations (Oi) and model predictions (Pi) in the units of a particular variable, with N being thenumber of observations. The root mean square error (RMSE) given
bysﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
PN2
i¼1 ðOi  P i Þ
RMSE ¼N1

ð15Þ

is the square root of the mean square error, also given in the units ofa particular variable. In general, RMSE P MAE. The degree in which
the RMSE value exceeds MAE is usually a good indicator of the presence and extent of outliers, or the variance of the differences between the modeled and observed values (Legates and McCabe,1999).
3. Input data
3.1. Initial conditions
The initial soil water content was set to a uniform value of0.25 cm3 cm3 throughout both soil proﬁles. Initial conditions for
the UNSATCHEM module were given in terms of concentrationsof Na+, Ca2+, Mg2+, and K+ in the liquid and solid phases. Initial conditions for the standard HYDRUS solute transport module were
speciﬁed in terms of ECsw, and concentrations of N–NHþ
4 and N–NO3:
To obtain these values, chemical analyses were performed on
soil samples collected in the beginning of the experiments. Concentrations of soluble cations Na+, Ca2+, Mg2+, and K+ were measured in the soil solutions collected from saturation extracts
using atomic absorption spectrophotometry. Exchangeable cations
Naþ , Ca2þ , Mg2þ , and Kþ were determined with the Bascomb method (Bascomb, 1964), using a solution of BaCl2 + Triethanolamine at
pH 8.1. The CEC was determined as deﬁned in Šimůnek and Suarez(1994). EC was determined by electrometry, and converted to ECsw.
N–NHþ
4 was determined using a modiﬁed Bertholot method(Searle, 1984). N–NO
3 was determined by an automated segmented ﬂow analyzer, using the cadmium reduction method toquantify nitrate–N (Hendrilsen and Selmer-Olsen, 1970). Table 3
presents the initial soluble and exchangeable cation concentrations
determined in representative soil proﬁles of both experimentalﬁelds.
3.2. Time-variable boundary conditions
Atmospheric and free drainage conditions were deﬁned asboundary conditions at the surface and the bottom of each ﬁeld
plot, respectively. Atmospheric boundary conditions were speciﬁed using meteorological data, from which daily values of the reference evapotranspiration rate (ET0) were calculated using thePenman-Monteith method (Allen et al., 1998). Crop evapotranspiration rates (ETc) were then calculated using the product of ET0
and Kc, where Kc is a crop coefﬁcient accounting for both soil evaporation and crop transpiration. Values of Kc were taken from Allenet al. (1998). As required by HYDRUS-1D, ETc daily values were divided into two components: crop transpiration (T) and soil evaporation (E) rates. These two components were estimated as afunction of the Leaf Area Index (LAI) and the corresponding Soil
Cover Factor (SCF), following Ritchie (1972). LAI values were measured in each plot of both experimental ﬁelds during different
stages of the maize cycle using a LI-COR area meter (Model LI3100C, LI-COR Environmental and Biotechnology Research Systems, Lincoln, Nebraska), and were linearly interpolated between
measurement dates. Although the combined effect of salinity andnitrogen fertilization is present in LAI values, in order to reduce
the number of variables in our study, and to simplify the analysis

134

T.B. Ramos et al. / Journal of Hydrology 407 (2011) 129–144

Table 3
Physical and chemical soil characteristics (initial conditions)a.
Alvalade

Mitra

Depth (cm)
Coarse sand (g kg1)
Fine sand (g kg1)
Silt (g kg1)
Clay (g kg1)
TextureBulk density (g cm3)
EC (dS m1)
SAR (mmol(c) L1)0.5
pH (H2O)

0–30
83
329
458
13
Loam
1.49
0.42
3.25
7.00

30–75
65
245
510
18
Silty–loam
1.51
1.22
3.87
7.13

75–100
58
215
481
246
Loam
1.61
0.96
2.96
7.33

0–30
461
284
176
79
Sandy–loam
1.51
0.48
0.41
6.63

30–50
431
306
178
85
Sandy–loam
1.70
0.63
0.56
6.59

50–90
423
323
160
9.4
Sandy–loam
1.69
0.25
0.57
7.11

Soluble cations (mmol(c) L1)
Ca2+
Mg2+
Na+
K+
Cl (mmol(c) L1)b

1.225
0.978
3.050
0.440
5.693

1.055
0.830
2.440
0.595
4.920

0.705
0.590
2.980
0.085
4.360

2.960
1.278
0.668
0.713
5.619

2.940
1.175
0.800
0.570
5.485

1.520
0.640
0.644
0.235
3.039

Exchangeable cations (mmol(c) kg1)
Ca2+
60.01
Mg2+
17.50
+
Na
2.53
K+
5.71
CEC (mmol(c) kg1)c
85.75

61.64
18.51
2.28
7.1
89.53

62.77
20.35
3.49
2.95
89.56

98.10
13.16
0.99
6.16
118.41

78.10
12.14
1.03
4.68
95.95

71.07
10.72
1.11
3.63
86.53

Gapon selectivity coefﬁcients (mol L1)1/2
KCa/Na2.92
KMg/Ca
0.33
KCa/K
0.11

1
N–NO3 (mmol(c) L )
0.276
1
0.001
N–NHþ
4 (mmol(c) L )

2.87
0.34
0.11
0.287
0.001

2.85
0.35
0.10
0.307
0.001

1.73
0.20
0.09
0.350
0.001

1.58
0.25
0.10
0.760
0.001

1.51
0.23
0.09
0.060
0.001

a

EC, electrical conductivity; SAR, sodium adsorption ratio; CEC, cation exchange capacity; K, Gapon selectivity coefﬁcient.
Calculated to maintain the charge balance.
cCalculated from the sum of ion exchange species.

b3.3. Soil hydraulic properties
Undisturbed soil samples (100 and 630 cm3) were collected at
the beginning of the experiment from different soil layers of eachsoil proﬁle to measure soil hydraulic properties. The soil water
retention curve, h(h), was determined in the laboratory using
suction tables with sand or kaolin for suctions below 500 cm,
and a pressure plate apparatus for suctions above 1000 cm. The
evaporation method (Wind, 1968; Halbertsma and Veerman,
1994) was further used to simultaneously estimate water retention
and hydraulic conductivity data between pressure heads ofapproximately 50 and 800 cm. The same samples had beenused previously to determine the saturated hydraulic conductivity
Ks using a constant-head method (Stolte, 1997), completing thehydraulic conductivity function, K(h). The parameters of the van

Alvalade

10

80

8

60

6

40

4

20

2

0
01-06-04

18-12-04

06-07-05

22-01-06

10-08-06

0
26-02-07

Mitra

100

10

80

8

60

6

40

4

20

2

0
01-06-04

19-12-04

07-07-05

23-01-06

Rainfall

11-08-06

Irrigation

Reference evapotranspiration (mm)100

Water applied (mm)of potential and actual root water uptake and solute leaching, weused only those LAI values measured at plot I-A (Fig. 1), where potential root water uptake was reduced due to water and osmotic
stresses thereby increasing the risk of salt and nitrogen leaching
to the groundwater, and plot IV-C (Fig. 1), where potential root
water uptake was affected only by water stress, and where we assumed that osmotic stress could be neglected since osmotic heads
in these plots were above the threshold value and thus a2(h/ ) = 1,
i.e., no reduction occurred. For the nitrogen simulations, whichwere run for every plot in sub-groups A and C, we adopted LAI values measured in plots I-A and IV-C, respectively.
Experimental ﬁelds were irrigated three times per week between June and September. In Alvalade, application amounts averaged 23 mm per irrigation event, while in Mitra the mean
application amount was 18 mm per irrigation event. Daily valuesof precipitation, irrigation, and ET0 for both experimental ﬁelds
are presented in Fig. 2.

0
27-02-07
ET0

Fig. 2. Daily values of precipitation, irrigation and reference evapotranspiration
rate in Alvalade (top) and Mitra (bottom) from 1st June 2004 to 27th February 2007.

Genuchten–Mualem equations were optimized using simultaneously retention and conductivity data determined by all threemethods with the RETC computer program (van Genuchten et al.,1991). Table 4 lists the van Genuchten–Mualem parameters forthe two ﬁeld sites. Identical soil hydraulic parameters were considered for all experimental plots at each ﬁeld, thus neglecting the135

T.B. Ramos et al. / Journal of Hydrology 407 (2011) 129–144

probable effect of the spatial variability of soil hydraulic properties
on water ﬂow and solute transport.

Table 4
Soil hydraulic parameters of the van Genuchten–Mualem functions (van Genuchten,1980) and solute transport parameters.
Alvalade

3.4. Solute transport parameters
Solute transport parameters were obtained from solute displacement experiments carried out on undisturbed 9040 cm3
cylindrical samples with a cross-sectional area of 452 cm2. A
0.05 M KCl pulse was applied during steady-state ﬂow. The sampling and preparation of the soil columns were carried out according to the method described by Mallants et al. (1994). All
experimental procedures are explained in Gonçalves et al. (2001).
The chloride breakthrough curves were expressed using thedimensionless concentration as a function of the number of porevolumes leached through the soil column. Transport parameters
were obtained using Toride et al’s. (1995) non-linear parameter
estimation code CXTFIT 2.1, implemented in the STANMOD software package (Šimůnek et al., 1999) by ﬁtting analytical solutionsof the CDE to observed breakthrough data. Dispersivity (k) values
were calculated from the ratio D/t and are presented in Table 4.
3.5. Relation between liquid and solid phases
The Gapon selectivity coefﬁcients K Ca=Na , K Mg=Ca , and K Ca=K arepresented in Table 3. They were calculated from the initial solubleand exchangeable cations concentrations presented in the same
Table, according to the system of equations described by Šimůnek
and Suarez (1994) for the case of exchange of the four cations simulated with UNSATCHEM.
In the standard HYDRUS solute transport module, nitrate (N–NO3 ) and ECsw were assumed to be present only in the dissolved
phase (Kd = 0 cm3 g1), while ammonium (N–NHþ4 ) was assumed
to adsorb to the solid phase using a distribution coefﬁcient Kd of
3.5 cm3 g1. The ﬁrst-order decay coefﬁcients lw and ls, represent
ing nitriﬁcation from N–NHþ
4 to N–NO3 in the liquid and solid
phases, were set to be 0.2 d1. The parameters Kd, lw, and ls weretaken from a review of published data presented by Hanson et al.
(2006), and represent the center of the range of reported values.
3.6. Ionic concentration of irrigation waters
Irrigation water was monitored for concentrations of Na+, Ca2+,Mg2+, electrical conductivity (ECiw), N–NHþ
4 , and N–NO3 :The average concentrations of individual ions in different irrigation waters (saline waters, waters with nitrogen, and fresh waters)
mixed during application to the crop in the two experimental ﬁelds
are presented in Table 5. Fertigation was limited to a certain number of irrigation events in the beginning of maize growing season.
In Alvalade, fertigation was applied during 8, 6, and 6 irrigation
events in the ﬁrst, second, and third years, respectively. In Mitra,
the fertilizer was applied during 4, 6, and 4 irrigation events during
the 3 years of our experiment. During remaining irrigation events,
the drip emitters that were ﬁrst used to apply water with nitrogen
applied only fresh irrigation water.
The ECiw values used in the standard HYDRUS solute transportmodule were based on relations between EC and solute concentrations derived using UNSATCHEM. This was done in order to obtain
similar osmotic stresses and actual root water uptake rates by thetwo modules. The relationships for Alvalade and Mitra were foundto be:

EC sw ¼ 0:1063csum þ 0:0915

ðR2 ¼ 0:999Þ

EC sw ¼ 0:1030csum þ 0:1585

ðR2 ¼ 0:998Þ

ð16Þ

respectively, where csum is the sum of the solute concentrations(M L3) in the liquid phase as given by UNSATCHEM. The relations

Depth (cm)hr (cm3 cm3)hs (cm3 cm3)a (cm1)
g (–)
‘ (–)
Ks (cm d1)k (cm)

0–30
0.050
0.380
0.027
1.21
4.41
16.6
25.8

Mitra
30–75
0.108
0.380
0.115
1.19
5.37
84.4
25.8

75–100
0.000
0.375
0.045
1.17
6.48
21.0
12.2

0–30
0.000
0.340
0.238
1.15
7.33
57.0
8.5

30–50
0.000
0.319
0.091
1.12
0.00
52.4
6.5

50–90
0.000
0.312
0.077
1.18
0.00
98.9
6.5

(16) are only slightly different from the common dilute-solution
approximation EC = 0.1csum (e.g., Bresler et al., 1982).
3.7. Root distribution and root-water uptake
In each plot of both experimental ﬁelds, the root depth was set
to 60 cm and the root density was assumed to decrease linearlywith depth. These estimates were based on ﬁeld observations
made with the minirhizotron technique of Machado and Oliveira(2003). Soil water pressure head parameters in the Feddes et al.
(1978) model were taken from the HYDRUS-1D internal database,which is based on Wesseling et al. (1991), i.e., h1 = 15, h2 = 30,
h3 = 325 to 600, h4 = 8000 cm. In the Maas (1990) function,
the salinity threshold (ECT) for maize corresponds to a value of
1.7 dS m1 for ECe, and a slope (s) of 12. These values were converted into ECsw using Eq. (7) and a kEC of 2. The correspondingh/T and s were obtained according to Eq. (6).
4. Results and discussion
The HYDRUS-1D simulations began on 22 April 2004 in Alvalade, and on 27 April 2004 in Mitra, at the beginning of the maize
growing season. The actual experiments began later, on 6 June
2004 in Alvalade, and on 15 June 2004 in Mitra, when irrigation
started. Beginning the HYDRUS-1D simulations long before irrigation events ensured that the initial soil water conditions were not a
factor in solute transport simulations.
The amount of applied water (irrigation and rainfall), the irrigation schedule, and daily ETc values were the same for all plots located in each experimental ﬁeld. The effect of irrigation water
quality on root water uptake (transpiration) was assumed to be
the main difference between experimental plots. Although measured and simulated data were compared at depths of 20, 40,
and 60 cm, only results for the depth of 40 cm are presented graphically, in order to limit the number of ﬁgures and to maintain consistency. The statistical analysis presented in Table 6 involves
results obtained for all three depths.
4.1. Volumetric water contents
Fig. 3 shows the water contents measured with TDRs for thetwo experimental ﬁelds in plots I-A and IV-C, and compares thesevalues with the results of the HYDRUS-1D simulations between 1stJune 2004 and 26th February 2007 (i.e., 1000 days). During
irrigation periods, the amount of water applied in both experimental ﬁelds (Table 2) was considerably higher than the usual amountof water used to irrigate maize in the Alentejo region(500–700 mm). The objective was to increase water contents
above soil ﬁeld capacity so that soil solution samples could be
easily collected using the installed ceramic cups. As a result, soil
water content rapidly increased in the beginning of each irrigationseason, and then varied between soil saturation and soil ﬁeld

136

T.B. Ramos et al. / Journal of Hydrology 407 (2011) 129–144

Table 5Weighted average ionic composition of irrigation waters applied in the experimental ﬁelds.
IrrigationwatersCa2+(mmol(c) L1)

Fresh waters
Alvalade
2.80
Mitra
1.10
Saline waters
Alvalade
2.80(2004)
Alvalade
2.80(2005)
Alvalade
2.80(2006)
Mitra
1.10(2004)
Mitra
1.10(2005)
Mitra
1.10(2006)
Waters with fertilizer
Alvalade
2.80(2004)
Alvalade
2.80(2005)
Alvalade
2.80(2006)
Mitra
1.10(2004)
Mitra
1.10(2005)
Mitra
1.10(2006)

Mg2+(mmol(c) L1)

Na+
(mmol(c) L1)b

K+
(mmol(c) L1)

Cl(mmol(c) L1)

EC
(dS m1)

SAR(mmol(c) L1)

USSL
classiﬁcationa

N–NHþ4(mmol(c) L1)

N–NO
3
(mmol(c) L1)

3.00
1.36

4.16
1.00

0.24
0.04

10.20
3.50

1.2
0.5

2.4
0.9

C3S1
C2S1

0.03
0.03

0.15
0.24

3.00

66.13

0.24

72.17

7.8

38.8

C4S4

0.03

0.15

3.00

66.13

0.24

72.17

7.8

38.8

C4S4

0.03

0.15

3.00

136.660.24

142.70

14.6

80.2

C4S4

0.03

0.15

1.36

74.23

0.04

76.73

8.1

66.9

C4S4

0.03

0.24

1.36

68.57

0.04

71.07

7.5

61.8

C4S4

0.03

0.24

1.36

68.33

0.04

70.83

7.5

61.6

C4S4

0.03

0.24

3.00

4.16

0.24

10.20

6.0

2.4

C4S1

46.1

46.1

3.00

4.16

0.24

10.20

6.0

2.4

C4S1

46.3

46.3

3.00

4.16

0.24

10.20

6.0

2.4

C4S1

45.6

45.6

1.36

1.00

0.04

3.50

8.0

0.9

C4S1

80.4

80.4

1.36

1.00

0.04

3.50

6.0

0.9

C4S1

49.5

49.5

1.36

1.00

0.04

3.50

8.0

0.9

C4S1

83.0

83.0

a
US Salinity Laboratory Staff (1954). C2, medium-salinity water (EC 0.25–0.75 dS m1); C3, high-salinity water (EC 0.75–2.25 dS m1); C4, very high salinity water(EC > 2.25 dS m1); S1, low-sodium water [SAR 0–10 (mmol(c) L1)0.5]; S4, very high sodium water [SAR > 26 (mmol(c) L1)0.5].
bCalculated to maintain the charge balance.

capacity. Between the end of the irrigation seasons (September)and the beginning of rainy seasons (October or November), soil
water contents gradually decreased, allowing maize to matureand be harvested. During rainy seasons, soil water contents weredependent on rainfall events (Fig. 2).
Some differences found between measurements and simulations during irrigation seasons can be explained by the fact that
TDR measurements were usually taken 24 h after an irrigation
event, which means that the highest observed water content values correspond roughly to values close to soil ﬁeld capacity, while
HYDRUS-1D, which calculates water ﬂow continuously, produced
higher water contents during irrigation events. Mismatch between
simulations and measurements depended on the time step used forspecifying boundary conditions (1 day) and were more signiﬁcant
for depths near the soil surface, which reﬂected these inter-daily
variations. While in the model we assume constant daily ﬂuxes,
soil evaporation varies during this time interval (has peaks during
the day and very low during the night). Soil water content at thesoil surface reﬂects these variations, which were not consideredin the model, which focused on longer-term changes. Additionally,
the distribution of water when drip emitters are used is known tobe spatially non-uniform. Soil water content is highest near the
drip line after water application and then water redistributes
throughout the soil proﬁle, as controlled by soil physical properties(Gärdenäs et al., 2005; Hanson et al., 2006). This non-uniformity
could explain deviations between simulated and measured watercontents, especially at the 20-cm depth, and a relatively high RMSE
of 0.04 cm3 cm3 obtained for both soils. Also, in Mitra the apparently chaotic TDR data during irrigation periods may have caused
larger deviations between observations and simulations. Using
two-dimensional version of HYDRUS-1D would likely resulted inbetter agreement between experimental and modeled data,
especially for shallower depths.

4.2. Root water uptake and transpiration
Fig. 4 shows the cumulative potential root water uptake in Alvalade and Mitra during the 3 years of the experiment. It also shows
the actual root water uptake during the same time period, while
either considering only water stress (plots IV-C) or both waterand salinity stress (plots I-A). In Alvalade, cumulative potential
root water uptake rates (transpiration) in plots I-A and IV-C were2074 and 2301 mm, respectively (Table 7). Since ETc was assumed
to be the same in all experimental plots, lower transpiration rates
were balanced by higher evaporation rates. Variations in transpiration reﬂected the effects of the salinity and fertilization gradientson plant growth, characterized by LAI values measured in different
plots. As an example, maximum LAI values (LAImax) varied between
3.6 and 5.4 m2 m2 in plot I-A, and between 5.3 and 5.5 m2 m2 inplot IV-C. Differences between LAImax in the different plots wereusually less than 1 m2 m2, which means that the results obtained
for the other experimental plots would not differ much from thosepresented here for plots I-A and IV-C, since transpiration was the
only factor that was different between experimental plots.
In all plots of sub-group C, we assumed that potential root
water uptake was reduced only due to water stress, since theseplots were irrigated with fresh irrigation water (Table 5). Rootwater uptake (transpiration) was reduced from 2301 to 1396 mm
when the UNSATCHEM module was used, and to 1420 mm when
the standard HYDRUS solute transport module was used. These
differences in cumulative actual transpiration were minimal (lessthan 3%). We also ran simulations with UNSATCHEM for plots
IV-C while either considering or neglecting the effect of the salinity
stress on the potential root water uptake. Osmotic heads in theseplots were above the threshold value and thus a2(h/ ) = 1. As aresult, the values for actual root water uptake were identical whileconsidering or neglecting the salinity stress, conﬁrming our initial

771
0.84
2.01
637
1.53
2.60
642
0.02
0.05
648
0.04
0.07
361
0.62
0.99
N, number of observations; MAE, mean absolute error; RMSE, root mean square error. Units for MAE and RMSE are the units of a particular variables.

244
1.60
2.04
361
0.61
0.87
244
1.91
2.35
333
2.14
3.91
216
4.68
6.27
223
3.28
4.16
N
MAE
RMSE

446
0.03
0.04

698
0.03
0.04

239
10.65
13.86

356
3.78
6.44

221
4.46
5.66

336
2.76
3.54

342
1.32
1.75

Mitra
Alvalade
Alvalade
Alvalade
Mitra
Alvalade

Mitra

Alvalade

Mitra

Mitra

Alvalade

Mitra

Alvalade

Mitra

Alvalade

Mitra

Alvalade

Mitra

N–NO
3
N–NHþ
4
ECsw (HYDRUS)
ECsw (UNSATCHEM)
SAR
Mg2+
Ca2+
Na+
Water content
Statistics

Table 6
Results of the statistical analysis between measured and simulated soil water contents, soluble Na+, Ca2+, Mg2+, electrical conductivity of the soil solution (ECsw), N–NHþ
4 and N–NO3 obtained for all three studied depths (i.e., 20, 40, and60 cm).

T.B. Ramos et al. / Journal of Hydrology 407 (2011) 129–144

137

assumption that no salinity stress reduction occurred in the plots
in sub-group C.
In all plots of sub-group A, we assumed that potential root
water uptake was reduced due to water and osmotic stresses, sincethese plots were irrigated with saline waters. As ECsw calculationsin the two codes followed different methodologies, root water uptake reductions due to osmotic stress were slightly different. In the
UNSATCHEM module, transpiration in plot I-A decreased from
2074 to 1208 mm due to water stress, and further to 932 mm
due to the osmotic stress. In the standard HYDRUS solute transportmodule, cumulative actual root water uptake in plot I-A was
949 mm, which is very close to the value obtained using UNSATCHEM. The effect of osmotic stress on water content can be observed in the simulation results (Fig. 3), in which HYDRUS-1D
produced water content values slightly lower in plot IV-C than inplot I-A, where root water uptake reduction was higher, during irrigation periods.
In Mitra, cumulative potential transpiration in plots I-A and IVC were 1809 and 2016 mm, respectively. Root water uptake was
reduced to 458 mm in plot IV-C due to water stress, and to346 mm in plot I-A due to water and osmotic stresses (values
determined by the UNSATCHEM module). The reduction due towater stress was considerably higher than in Alvalade. This dramatic reduction was caused mainly by soil pressure heads being
kept above the optimum range of the Feddes model (Feddes
et al., 1978) for maize (h > 30 cm). The irrigation schedule wasthus the main factor responsible for causing such high reductionsin actual root water uptake. Also, soil hydraulic properties determined in Mitra resulted in some atypical values (e.g., the g parameter for different layers ranged from 1.12 to 1.18 and was lower
than usual for coarse textured soils, and the ‘ parameter was 0.0
in layers under 30 cm), which may have contributed to errors in
the water content and ﬂux simulations. Underestimated values of
the K(h) function could cause HYDRUS-1D to predict longer periodsof soil proﬁle saturation than in reality, maintaining pressure
heads above the optimum range of the Feddes stress responsemodel, thus contributing to the lower root water uptake. Note thatmethods used to determine the K(h) function (i.e., the evaporationand hot air methods) produce more precise K(h) estimates for pressure heads below -50 cm and are less precise close to saturation.
Nevertheless, cumulative actual transpiration agreed with lower
maize yields as reported in Ramos et al. (2009). However, we were
not able to relate yields with the Ta/Tp ratio.
4.3. Overall salinity
Soil salinity increased considerably in the plots irrigated withsaline waters (Fig. 5). In Alvalade, measured ECsw in plot I-A
reached values higher than 6.0 dS m1 during the ﬁrst two irrigation seasons, and higher than 16.0 dS m1 during the third year,
as the applied waters became more saline. In the plots irrigatedwith the locally available water (IV-C), soil salinity remained below
5.0 dS m1 throughout all irrigation seasons. Soil salinity decreased
in all plots during rainfall periods, due to soil leaching. Only in therainy season of 2004–2005, rainfall was not sufﬁcient tocompletely remove salts from the root zone in plot I-A, where ECsw
values between 2.0 and 6.3 dS m1 were observed in the soil proﬁle. In Mitra, irrigation with saline waters led to similar ECsw peaks
as in Alvalade, but with ECsw values reaching only 9.0 dS m1 during the ﬁnal irrigation season. In plot IV-C, where the applied fresh
irrigation water was of much better quality than in Alvalade, ECsw
was kept below 2.0 dS m1 throughout irrigation seasons. Soilsalinity also decreased in Mitra during rainfall periods, due to soilleaching. Since the saturated hydraulic conductivity of the top soil
layer in Mitra (Ks = 57.0 cm d1) was higher than in Alvalade(Ks = 16.6 cm d1), salts were more easily removed from the root

Water content (cm 3cm -3 )T.B. Ramos et al. / Journal of Hydrology 407 (2011) 129–144

0.5

Watercontent (cm 3 cm-3)

138

0.5

Alvalade
40 cm

0.4
0.3
0.2
0.1

I

0.0
01-06-04

I

18-12-04

I

06-07-05

22-01-06 10-08-06

26-02-07

Mitra
40 cm

0.4
0.3
0.2
0.1

I

0.0
01-06-04

I

18-12-04

I

06-07-05

22-01-06

10-08-06

26-02-07

Time
HYDRUS-1D simulations (plot I-A)
HYDRUS-1D simulations (plot IV-C)

Measureddata (plot I-A)
Measureddata (plot IV-C)

Fig. 3. Measured (TDR) and simulated (HYDRUS-1D) volumetric water contents at a40-cm depth in Alvalade (top) and Mitra (bottom). Symbol I denotes the irrigation
periods.

Root water uptake (mm)zone, although during the rainy season of 2004–2005, ECsw values
above 2.0 dS m1 were still observed in the deeper layers.
Fig. 5 shows ECsw values measured by electrometry in samples
collected in plots I-A and IV-C, and ECsw values simulated using
the major ion chemistry module (UNSATCHEM module) and thestandard HYDRUS solute transport module. Although ECsw was calculated using different methodologies, the two modules produced
very similar results during irrigation seasons. The main differenceswere found at the end of the irrigation seasons when soil watercontents decreased signiﬁcantly below ﬁeld capacity, due to rela-

Alvalade
2500

2500
Plot: I-A

Plot: IV-C

2000

2000

1500

1500

1000

1000

500

500

0

0
0

Root water uptake (mm)tively high air temperatures and soil evaporation rates. During
these periods, the standard HYDRUS solute transport module simply increased ECsw linearly as the soil dried out, while the UNSATCHEM module produced a nonlinear increase of ECsw as a result of
cation exchange. The UNSATCHEM model also considers processesof precipitation/dissolution of solid phases, such as calcite and gypsum. However, the reported pIAP values (the negative logarithm of
the ion activity product) indicated that conditions were undersaturated with respect to calcite and gypsum. In Alvalade, minimum
calculated pIAP values for calcite and gypsum were 12.8 and
31.2, respectively. In Mitra, corresponding minimum pIAP values
were 12.5 and 30.8. Since precipitation only occurs at values of
8.37 for calcite and 4.8 for gypsum (Truesdell and Jones, 1974), differences found between the two HYDRUS-1D modules were caused
mainly by cation exchange. Consequently, during irrigation seasons the standard HYDRUS solute transport module produced
higher peaks of ECsw than the UNSATCHEM module. Similar trends
were obtained during the rainy season of 2004–2005, when rainfall
was not sufﬁcient to remove salts from the root zone. During thisseason, the standard HYDRUS solute transport module clearly
showed that high concentrations of salts remained in the soil proﬁle of Alvalade, while the UNSATCHEM module produced lower
simulated ECsw values. RMSE between measured and simulated
ECsw values (Table 6) were not sufﬁcient to conclude which module
produced better results. In addition, since it was difﬁcult to collectsoil solution samples with the ceramic cups during dry periods,
model simulations could not be compared with measured data
during these periods, which would have provided conclusive evidence regarding which module better ﬁtted measured data. On
the other hand, the statistical indicators showed that both modules
reproduced ECsw measured values equally well, with RMSE being
similar for both modules for each soil.
Model simulations were very helpful in evaluating the management of irrigation with saline waters in the two experimentalﬁelds, allowing for a few additional observations. Based on model
predictions, blending saline waters, in which ECiw varied between
7.5 and 14.6 dS m1, with irrigation waters of good quality
increased the amount of water available to meet crop water

400

800

1200

0

400

800

1200

Mitra
2500

2500

Plot: I-A

2000

2000

1500

1500

1000

1000

500

500

0

Plot: IV-C

0
0

200

400

600

800

1000

0

200

Time (days)

400

600

800

1000

Time (days)

Potential root water uptake
Actual root water uptake due to water stress
Actual root water uptake due to water stress plus salinity stressFig. 4. Cumulative potential and actual root water uptake in plots I-A (left) and IV-C (right) in Alvalade (top) and Mitra (bottom).

139

T.B. Ramos et al. / Journal of Hydrology 407 (2011) 129–144
Table 7
Cumulative potential and actual root water uptake determined in the experimental plots.a.
Experimental plots

I-A
II-A
III-A
IV-A
I-C
II-C
III-C
IV-C

Mitra

Potential root water uptake (mm)Actual root water uptake (mm)Potential root water uptake (mm)Actual root water uptake (mm)2074 (2074)
2074
2074
2074
2301
2301
2301
2301 (2301)

949 (932)
999
1005
1021
1418
1420
1418
1422 (1396)

1809 (1809)
–
1809
1809
2016
–
2016
2016 (2016)

362 (346)
–
368
374
483
–
482
484 (458)

Values obtained with the standard HYDRUS solute transport module. Values in brackets were determined with the UNSATCHEM module.

4.4. Individual cations
Measured and simulated concentrations of soluble Na+, Ca2+,and Mg2+ for plots I-A and IV-C at a depth of 40 cm are presented
in Fig. 6. As the only cation being added to the synthetic saline irrigation waters, the general behavior of sodium was similar to ECsw
The highest concentrations were reached in plots I-A during theirrigation seasons. In Alvalade, the highest measured Na+ concentrations were about 50 mmol(c) L1 during the ﬁrst 2 years and
130 mmol(c) L1 during the third year. In Mitra, as saline waters
were similar during the three irrigation seasons, measured Na+concentrations peaked always at about 50 mmol(c) L1. During
the rainy seasons, Na+ concentrations decreased considerably due
to soil leaching, similar to ECsw. RMSE calculated for Na+ concentrations (Table 6) resulted in 13.86 and 6.44 mmol(c) L1 for Alvaladeand Mitra, respectively.
The general dynamics of calcium and magnesium concentrations were similar to those of sodium and ECsw. In both experimental ﬁelds, the highest measured Ca2+ and Mg2+ concentrations wereabout 20 and 10 mmol(c) L1, respectively. During rainfall events,
these values also decreased due to soil leaching. As Ca2+ andMg2+ were not added to the synthetic saline irrigation waters
applied in plots I-A, concentrations of these cations in saline waters
were the same as in fresh waters in plot IV-C (Table 5). Since Na+
was applied in large concentrations to plot I-A, both Na+ concentrations in the soil solution and in the solid phase (not shown)increased, leading to soil sodiﬁcation. While the exchangeable
Na+ concentration in the solid phase increased, the other cations,
namely Ca2+ and Mg2+, were inevitably released to the soil solution.
For this reason concentrations of Ca2+ and Mg2+ were higher inplots I-A than in plots IV-C where low concentrations of Na+ were
present in the irrigation waters. As the UNSATCHEM module is able
to consider cation exchange, compared to simpler models based on
the adsorption isotherms, this module is more adequate forrepresenting reality than linear models with adsorption isotherms.
RMSE obtained for Ca2+ and Mg2+ in Alvalade were 5.66 andAlvalade
40 cm

25ECsw (dSm -1)

requirements while maintaining some level of maize yield (Ramos
et al., 2009). Both measured and simulated values revealed thatduring irrigation seasons, ECsw values became high enough to leadto signiﬁcant maize yield reductions (e.g., Ayers and Westcot,
1985; Mass, 1990; Steppuhn et al., 2005), but never exceeded
levels that would lead to zero yield, which for our crop would
correspond to a ECsw of 20 dS m1 (ECe = 10 dS m1; kEC = 2).
Additionally, variations in simulated ECsw (Fig. 5), which increased
and decreased regularly with irrigation events, showed the contribution of blending in providing some control on soil salinity status
for both experimental ﬁelds. These variations were caused byblending saline waters with fresh waters, while keeping soil salinity relatively constant, which otherwise would have increased
considerably.

20

I

I

I

15
10
5
0
01-06-04

18-12-04

06-07-05

20

22-01-06

10-08-06

26-02-07

Mitra
40 cm

25

ECsw (dSm -1)

aAlvalade

I

I

I

15
10
5
0
01-06-04

18-12-04

06-07-05

22-01-06

10-08-06

26-02-07

Time
Measured data (plot I-A)
Measured data (plot IV-C)HYDRUS-1D simulations (plot I-A)
HYDRUS-1D simulations (plot IV-C)

Fig. 5. Measured and simulated soil solution electrical conductivities at a 40-cm
depth in Alvalade (top) and Mitra (bottom). Symbol I denotes the irrigation periods.
Full lines represent the results of the major ion chemistry module (the UNSATCHEM
module). Dashed lines represent the results obtained with the standard HYDRUS
solute transport module.

4.16 mmol(c) L1, respectively. In Mitra, RMSE obtained for Ca2+and Mg2+ were 3.54 and 1.75 mmol(c) L1, respectively.
We believe that the inability of HYDRUS-1D to produce lower
RMSE for Mitra, and to a lesser extent at Alvalade, was caused bymultiple factors. In general, measurement errors, model input errors, and model structure errors could cause disagreement between simulated results and experimental data. First, whilemodels usually report point values, measurements with suction
cups are averaged over a sampling area of a certain volume, thesize of which depends on soil hydraulic properties, the soil watercontent, and the applied suction in the ceramic cup (Weihermüller
et al., 2005). Measured values thus do not represent point values.
Second, while measured soil hydraulic and solute transport parameters were determined in the laboratory on soil samples of acertain size, the obtained parameters may not always be representative of simulated ﬂow, transport and reaction processes at amuch larger ﬁeld scale. Third, while soil hydraulic and solutetransport properties are spatially variable at the ﬁeld scale, our
HYDRUS-1D simulations assumed homogeneous soil environment

T.B. Ramos et al. / Journal of Hydrology 407 (2011) 129–144-1-1

Mgconcentration (mmol(c) L ) Caconcentration (mmol (c) L )-1

Na concentration (mmol(c) L )140

Alvalade
40 cm

140
120

I

I

I

120

100

100

80

80

60

60

40

40

20

20

I

I

I

0

0
01-06-04

18-12-04

06-07-05

22-01-06

10-08-06

26-02-07

40 cm

60
I

01-06-04

20

20

0

0
40 cm
I

I

I

30
20

10

10

22-01-06

I

10-08-06

26-02-07

I

40 cm

40

20

I

I

I

0

0
01-06-04

06-07-05

40 cm
I

I

40

40

18-12-04

60

I

40

30

Mitra
40 cm

140

18-12-04

06-07-05

22-01-06

10-08-06

26-02-07

01-06-04

18-12-04

06-07-05

22-01-06

10-08-06

26-02-07

Time

Time
Measured data (plot I-A)
Measured data (plot IV-C)HYDRUS-1D simulations (plot I-A)
HYDRUS-1D simulations (plot IV-C)

Fig. 6. Measured and simulated soluble sodium (top), calcium (middle), and magnesium (bottom) concentrations at a 40-cm depth in Alvalade (left) and Mitra (right). Symbol
I denotes the irrigation periods.

and thus might not fully account for inherent spatial variability.
Finally, as already discussed above, while only a one-dimensionalmodel was used to evaluate experimental data, the drip irrigation
used in the ﬁeld experiments produces certain multidimensional
phenomena that were not accounted for by the model.
4.5. Sodium adsorption ratio
Fig. 7 presents the results for measured and simulated SAR values. SAR is an integral variable that characterizes salt-affected soils
and provides information on comparative concentrations of Na+,
Ca2+, and Mg2+ in soil solutions. This variable takes into consideration that the adverse effects of sodium are moderated by the presence of calcium and magnesium ions. This makes SAR an importantvariable to consider when managing saline irrigation waters.
The standard HYDRUS solute transport module is unable todetermine SAR. As a result, many studies found in the literature(e.g., Forkutsa et al., 2009; Hanson et al., 2008; Roberts et al.,2009) do not consider the effects of saline waters and management
practices on soil sodiﬁcation. The major ion chemistry module istherefore a better tool for saline water management, as it is capable of evaluating interactions between individual ions, from which
SAR can be determined. It is thus no surprise that the statisticalindicators calculated for SAR showed the same trends as obtainedfor individual ions. RMSE obtained for Alvalade and Mitra were
6.27 and 3.91 (mmol(c) L1)0.5, respectively.

The general behavior of SAR was very similar to that documented in Gonçalves et al. (2006). In plots I-A, SAR increased
rapidly in the surface layers after irrigation events, and thengradually at deeper depths. The application of saline irrigationwaters, in which SAR ranged from 38.8 to 80.2 (mmol(c) L1)0.5(Table 5), led to very high SAR values in the soil solution, reaching
about 54 and 21 (mmol(c) L1)0.5 at the end of the experiments in
Alvalade and Mitra, respectively. In those plots irrigated with fresh
waters (plots IV-C), SAR values did not vary signiﬁcantly, having
roughly the same values in the beginning as at the end of theexperiments. While HYDRUS-1D simulations with the standard
solute transport module revealed that overall salinity decreased
as a result of soil leaching after rainfall, simulations using the
UNSATCHEM module provided predictions of SAR, which is an
important component in designing better management of irrigation practices with saline waters. Consequently, soil salinizationand sodiﬁcation could be considered simultaneously.
4.6. Nitrogen concentrations
Figs. 8 and 9 present the measured and simulated concentra
tions of N–NHþ
4 and N–NO3 in the soil solution at a depth of
40 cm. Since the amount of nitrogen applied in sub-groups A–C
was the same, N–NHþ
4 and N–NO3 concentrations for sub-groups
A and C are compared to determine the effect of salinity stresson nutrient uptake. Fig. 8 compares N–NHþ
4 concentrations in plots

141

T.B. Ramos et al. / Journal of Hydrology 407 (2011) 129–144

I

30
20
10
18-12-04

06-07-05

22-01-06

10-08-06

26-02-07

Mitra
40 cm

60
50

0.75

I

I

I

-1

40

0
01-06-04

I

I

I

40
30
20

0.50
0.25
0.00
01-06-04

18-12-04

06-07-05

10-08-06

26-02-07

Mitra
40 cm

1.00
0.75

22-01-06

I

I

I

0.50

+

SAR (mmol (c) L-1)0.5

I(mmol (c) L )I

Alvalade
40 cm

1.00

N-NH 4 concentration-1(mmol (c) L )SAR (mmol (c) L-1)0.5

50

N-NH4+ concentration

Alvalade
40 cm

60

10
0
01-06-04

18-12-04

06-07-05

22-01-06

10-08-06

26-02-07

Time
Measured data (plot I-A)
Measured data (plot IV-C)HYDRUS-1D simulations (plot I-A)
HYDRUS-1D simulations (plot IV-C)

Fig. 7. Measured and simulated sodium adsorption ratios (SAR) at a 40-cm depth inAlvalade (top) and Mitra (bottom). Symbol I denotes the irrigation periods.

I-A and I-C for both experimental ﬁelds. As the results obtained forexperimental plots of the remaining three groups (II to IV) werevery similar, results obtained in group I are used to exemplify N–NHþ
4 behavior. Fig. 9 compares measured and simulated N–NO3
concentrations in experimental plots of groups I (I-A and I-C) andIV (IV-A and IV-C). The results obtained for the remaining groups
in the middle of the nitrogen gradient (groups II and III) are not given, as the results in groups I and IV represent the maximum andminimum N–NO
3 dynamics.
Since nitriﬁcation in soils is usually a fast process, measuredN–NHþ
4 concentrations in every experimental plot of both experimental ﬁelds were always lower than 0.30 mmol(c) L1. While thehighest residual values were observed during fertigation events,
most values observed during remaining months were much lower
and close to zero. Corresponding HYDRUS-1D simulations also pre1dicted the highest N–NHþ
4 concentrations of about 0.13 mmol(c) L
during fertigation events, but only in the top layer of Alvalade (notshown graphically). Most simulated values were also close to zero,but even lower than the N–NHþ
4 concentrations measured in thesoil solution. Since most measured and calculated N–NHþ
4 concentrations were close to zero, RMSE resulted in values lower than
0.07 mmol(c) L1 (Table 6), which may be quite signiﬁcant sincevalues measured in the ﬁeld were practically residual.
There are a number of reasons that can explain low levels ofagreement for N–NHþ
4 concentrations. The nitriﬁcation rate coefﬁcient of 0.2 d1 was an average value taken from Hanson et al.
(2006). Since measured N–NHþ
4 values were higher than those simulated, nitriﬁcation in HYDRUS-1D was apparently faster than that
occurring in the soil. This was more relevant for Alvalade where
measured N–NHþ
4 concentrations were higher than in Mitra. Whenconsidering the lowest nitriﬁcation value (0.02 d1) found in
Hanson et al. (2006) for the I-A plot in Alvalade, simulations
1
showed higher N–NHþ
4 peaks, reaching values of 1.0 mmol(c) L
in the top 20 cm, and lower values with increasing depth. N–NO
3
concentration peaks decreased as ammonium was converted into
nitrate at a rate 10 times slower than before (Fig. 10). However,this explanation can explain disagreement only during the days
following fertigation events. A more likely explanation is related

0.25
0.00
01-06-04

18-12-04

06-07-05

22-01-06

10-08-06

26-02-07

Time
Measureddata (plot I-A)
Measureddata (plot I-C)

HYDRUS-1D simulations (plot I-A)
HYDRUS-1D simulations (plot I-C)

Fig. 8. Measured and simulated ammonium concentrations at a 40-cm depth in theexperimental plots of group I in Alvalade (top) and Mitra (bottom). Symbol I
denotes the irrigation periods.

to the processes that were not considered in our simulations, butwhich are important for accurately describing N–NHþ
4 concentrations in the soil. Namely, mineralization of crop residues or otherorganic wastes, mineralization of the soil humus fraction, the release of N–NHþ
4 adsorbed to the solid phase into the liquid phase
due to cation exchange, or other N processes that cannot be described with sequential ﬁrst-order decay chains available in HYDRUS-1D. Since such processes were not considered in our
simulations, here we are discussing only whether the model adequately simulated residual concentrations in the soil.
High nitrate concentrations were measured in both soil proﬁlesduring irrigation seasons due to fertigation and nitriﬁcation of the
NH4NO3 fertilizer. In both experimental ﬁelds, the highest measured N–NO
3 concentrations were observed, as expected, in theplots located in group I, where the largest amount of fertilizer
was applied. Measured N–NO
3 concentrations in group I reacheda maximum value of about 15.0 and 16.5 mmol(c) L1 in Alvaladeand Mitra, respectively. Only a few isolated measurements produced values higher than these. In group IV, where no NH4NO3 fertilizer was added, the only source of nitrogen was the naturally
occurring nitrogen in the fresh water. In both soils, the highest ob1
served N–NO
.
3 concentrations were approximately 1.0 mmol(c) L
Similar values were also observed during rainy seasons. As discussed for N–NHþ
4 , some of the observed N–NO3 ; namely the values measured during the rainfall seasons, likely had origin innitrogen processes that were not considered in our simulations.
As a result, RMSE calculated for N–NO
3 concentrations were 2.60
and 2.01 mmol(c) L1 (Table 6).
The effect of the salinity stress on nutrient uptake is presented
in Table 8. According to HYDRUS-1D simulations, cumulative nitrogen ﬂuxes through the bottom of the plots located in sub-groups A
were always higher than those predicted in the plots of sub-groups
C. Such results were expected since only passive nutrient uptake
was considered in our study. Irrigation of the experimental plots
in sub-group A with saline waters led to the reduction of root
water uptake due to osmotic stress, and consequently to thedecreased mass ﬂow of nutrients into roots. As a consequence of

T.B. Ramos et al. / Journal of Hydrology 407 (2011) 129–144

N-NO3- concentration (mmol(c) L-1)

N-NO3- concentration (mmol(c) L-1)

142

Alvalade
GroupI
40 cm

30
25

I

30

I

I

25

20

20

15

15

10

10

5

5

I

I

I

0

0
01-06-04

18-12-04

06-07-05

22-01-06

10-08-06

26-02-07

GroupIV
40 cm

2.0
1.5

Mitra
GroupI
40 cm

I

01-06-04

18-12-04

I

1.5

1.0

1.0

0.5

0.5

22-01-06

10-08-06

26-02-07

GroupIV
40 cm

2.0

I

06-07-05

I

I

I

0.0

0.0
01-06-04

18-12-04

06-07-05

22-01-06

10-08-06

26-02-07

01-06-04

18-12-04

06-07-05

Time

22-01-06

10-08-06

26-02-07

Time

Measured data (plot I-A)
Measured data (plot I-C)

HYDRUS-1D simulations (plot I-A)
HYDRUS-1D simulations (plot I-C)

Alvalade
PlotI-A
40 cm

1.00
0.75

I

I

I

0.50
0.25
0.00
01-06-04

N-NO3- concentration(mmol (c) L-1)

N-NH4+ concentration(mmol (c) L-1)

Fig. 9. Measured and simulated nitrate concentrations at a 40-cm depth in the experimental plots of groups I (top) and IV (bottom) in Alvalade (left) and Mitra (right). Symbol
I denotes the irrigation periods.

18-12-04

06-07-05

30
25

22-01-06

10-08-06

26-02-07

40 cm
I

I

I

20
15
10
5
0
01-06-04

18-12-04

06-07-05

22-01-06

10-08-06

26-02-07

Time
Measured data (plot I-A)
-1
Nitrification rate 0.2 d
Nitrification rate 0.02 d-1
Fig. 10. Measured and simulated ammonium (top) and nitrate (bottom) concentrations at a 40-cm depth in the experimental plot I-A in Alvalade when considering
nitriﬁcation rates of 0.2 and 0.02 d1. Symbol I denotes the irrigation periods.

the salinity stress and reductions in nutrient uptake, the mass ﬂow
through the bottom of the soil proﬁle increased.
Although this simple approach was sufﬁcient to describe ourresults, one could pose the obvious question whether active
nutrient uptake should have been considered, and if it were,would we have obtained different results. Answers to these ques-

tions can be found in hypothetical examples discussed in Šimůnek and Hopmans (2009). First, nutrient concentrations
measured in a soil do not provide sufﬁcient evidence about therelative importance of active and passive nutrient uptake. Second,
had active uptake been considered, root nutrient uptake would
have been larger, which would have meant that N–NHþ
4 and N–NO
ﬂuxes
throughthebottomofthesoil
proﬁles
would
have
3
been smaller than those presented in Table 8. Nutrient uptake
would have always been reduced due to the effect of osmotic
stress, which would likely reduce the nutrient demand. However,since our objective was to simulate our ﬁeld experiments in a relatively simple way, considering only passive nutrient uptake
seemed to be sufﬁcient to describe nitrogen dynamics in thetwo studied soils. The values presented in Table 8 are merely
indicative and serve only as a reference. The simulated cumulative N–NO
3 ﬂuxes, across the bottom of the soil proﬁle for example, in plots I-A and I-C of Alvalade differed by only 22 g m2
during the 3 years of the experiment. It is obvious that had we
not simpliﬁed the LAI values used in our study, the LAI values thatshould have been used in I-C would have been higher than thosemeasured in plot IV-C, which were simply extrapolated to theremaining plots in sub-group C. LAI values in plot I-C would thendescribe a crop well supplied with nitrogen, while plot IV-C wasnot, and with no salinity stress. In this scenario, higher LAI values
would produce higher transpiration rates. Consequently, root
water and nutrient uptakes would also increase, resulting in lower nitrogen ﬂuxes to the groundwater. In this more realistic scenario the effect of the salinity stress on non-source pollution
would have been much more considerable.
5. Summary and conclusions
The HYDRUS-1D numerical model successfully simulated waterand solute transport in two multifactoral experiments, in whichwaters with different salinities and nitrogen concentrations were143

T.B. Ramos et al. / Journal of Hydrology 407 (2011) 129–144
Table 8Cumulative nitrogen ﬂuxes in the soil proﬁles of Alvalade and Mitra.
Alvalade

Mitra
2
N–NO
3 (g m )2N–NHþ
4 (g m )

Plots

A

C

2
N–NO
3 (g m )2N–NHþ
4 (g m )

A

C

A

C

A

C

1.30E+02
6.70E+01
3.65E+01
6.72E+00

1.33E+02
6.66E+01
3.70E+01
6.34E+00

6.37E+01–
1.27E+01
1.22E+00

6.37E+01–
1.17E+01
1.08E+00

7.07E+01–
1.99E+01
8.45E+00

6.98E+01–
1.81E+01
7.46E+00

Cumulative amount of solute removed from the ﬂow region by root water uptake
I
5.95E01
1.03E+00
4.97E+01
8.00E+01
II
3.31E01
5.26E01
2.88E+01
4.24E+01
III
1.81E01
2.87E01
1.57E+01
2.37E+01
IV
8.13E03
1.11E02
2.70E+00
3.50E+00

3.59E02–
8.76E03
1.75E03

1.45E01–
2.99E02
3.46E03

7.12E+00–
2.88E+00
1.64E+00

1.34E+01–
4.52E+00
2.10E+00

Cumulative solute ﬂux across the bottom of the soil proﬁleI
1.46 E05
1.42 E05
2.04 E+02
II
1.42 E05
1.42 E05
1.01 E+02
III
1.42 E05
1.42 E05
5.26 E+01
IV
1.42 E05
1.41 E05
6.36 E+00

7.15 E05–
7.15 E05
7.15 E05

7.14 E05–
7.14 E05
7.16 E05

1.28 E+02–
3.07 E+01
9.04 E+00

1.21 E+02–
2.62 E+01
7.47 E+00

Cumulative solute ﬂux across the soil surface
I
1.24E+02
1.29E+02
II
6.18E+01
6.17E+01
III
3.12E+01
3.20E+01
IV
1.34E+00
1.26E+00

1.82 E+02
8.63 E+01
4.67 E+01
5.12 E00

used. In these experiments, irrigation with waters blended usingsynthetic saline irrigation waters (ECiw 6 14.6 dS m1) and fresh
irrigation waters (ECiw 6 1.2 dS m1) led to the salinization andsodiﬁcation of the two studied soils.
The major ion chemistry module of HYDRUS-1D successfully
simulated the water regime (RMSEAlvalade = 0.04 cm3 cm3;
RMSEMitra = 0.04 cm3 cm3), the overall salinity characterized by
ECsw(RMSEAlvalade = 2.04–2.35 dS m1;
RMSEMitra = 0.87–0.99
1
dS m ), the concentration of soluble Na+ (RMSEAlvalade = 13.86
mmol(c) L1; RMSEMitra = 6.44 mmol(c) L1), Ca2+ (RMSEAlvalade =5.66 mmol(c) L1; RMSEMitra = 3.54 mmol(c) L1), Mg2+ (RMSEAlvalade= 4.16 mmol(c) L1;
RMSEMitra = 1.75 mmol(c) L1),and
SAR(RMSEAlvalade = 6.27 (mmol(c) L1)0.5; RMSEMitra = 3.91 (mmol(c)L1)0.5) in different plots of each experimental ﬁeld. RMSE were always lower in the soil with coarse texture of Mitra than in the soil
with medium texture of Alvalade. Possible causes of disagreements
between the modeling and experimental data were discussed in
the Results section. In addition, the less favorable hydraulic conditions of Alvalade, as compared to Mitra, namely the lower Ks valuesof the soil with medium texture, were decisive for not obtaining
higher levels of performance in the goodness-of-ﬁt tests for thissoil. Note that in Alvalade, solute concentrations in the irrigationwater were low (e.g., for Ca2+ and Mg2+), resulting only from what
was present in the available water in the region.
HYDRUS-1D further predicted root water uptake reductions due
to water and osmotic stresses. Root water uptake was reduced byabout 39% in the experimental plots at Alvalade and 77% at Mitra,
due to water stress. These water stress reductions were a consequence of the soil physical and hydrodynamic characteristics andplant type, but mostly of the irrigation schedule. Potential transpiration was further reduced by about 59% at Alvalade and 83% at
Mitra due to the effects of the osmotic stress in the experimental
plots irrigated with saline waters.
Root water uptake reductions obtained with the major ion
chemistry module were reproduced with the standard HYDRUS
solute transport module in order to study the effect ofsalinity stress on nutrient uptake. HYDRUS-1D successfully
1modeled N–NHþ
; RMSEMitra = 0.05
4 (RMSEAlvalade = 0.07 mmol(c) L
mmol (c) L 1 ) and N–NO
(RMSE=
2.60 mmol (c) L 1 ;Alvalade
3
RMSEMitra = 2.01 mmol(c) L1) concentrations while either assuming or neglecting the effect of the osmotic stress on nutrient uptake. HYDRUS-1D does not account for N processes that cannot
be described with sequential ﬁrst-order decay chains. As a result,

RMSE obtained for N–NHþ
4 and for N–NO3 reﬂect errors caused
by not considering processes other than nitriﬁcation, which wereapparently relevant for describing residual concentrations. According to HYDRUS-1D simulations, irrigation with saline waters led toroot water and nutrient uptake reductions due to osmotic stress.
Consequently, the ﬂuxes of N–NHþ
4 and N–NO3 through the bottomof the soil proﬁles increased.
In our study, most model inputs were independently measured
in the laboratory and used in simulations without any further
adjustments and/or calibration. The correspondence between measurements and model results would have obviously been better,
had the input parameters been calibrated. However, only a modelthat can be successfully run with independently measured input
parameters is sufﬁciently robust for practical applications.
In spite of the considerable demand on input data, HYDRUS-1D
proved to be an effective and versatile tool that may become very
useful for irrigation management in regions with scarce water resources where suitable waters are not always available for irrigation. HYDRUS-1D was able to analyze two of the most importantsoil processes (i.e., transport and reactions of salts and nitrogenspecies in the soil proﬁle) resulting in degradation of groundwaters
in these regions in an integrated way. HYDRUS-1D, after proper
calibration and validation, should be considered to be a useful toolfor establishing sound irrigation policies in order to mitigate soil
salinization/sodiﬁcation and non-point source pollution from agricultural applications of fertilizers in irrigated areas of countries located in regions with arid, semi-arid, and even sub-humid
conditions.
Acknowledgments
This work was possible due to the funding provided by the Project PTDC/AGR-AAM/66004/2006 of the Fundação para a Ciência ea Tecnologia (FCT) and the Project AGRO 727 of the Portuguese
Ministry of Agriculture, Fisheries, and Rural Development. T.B. Ramos was funded by the FCT grant (Contract SFRH/BD/60363/2009).
References
Ahuja, L.R., Rojas, K.W., Hanson, J.D., Shaffer, M.J., Ma, L. (Eds.), 2000. The Root Zone
Water Quality Model. Water Resources Publ., LLC, Highlands Ranch, CO, 372pp.
Ajdary, K., Singh, D.K., Singh, A.K., Khanna, M., 2007. Modelling of nitrogen leachingfrom experimental onion ﬁeld under drip irrigation. Agric. Water Manage. 89,
15–28.
Allen, R.G., Pereira, L.S., Raes, D., Smith, M., 1998. Crop Evapotranspiration –Guidelines for Computing Crop Water Requirements. Irrig. Drain. Pap. 56. FAO,
Rome, Italy.
Ayers, R., Westcot. D., 1985. Water Quality for Agriculture. Irrig. Drain. Pap., 29. FAO,
Rome, Italy.

144

T.B. Ramos et al. / Journal of Hydrology 407 (2011) 129–144

Bascomb, C.L., 1964. Rapid method for the determination of cation-exchange
capacity of calcareous and non-calcareous soils. J. Sci. Food Agric. 12, 821–823.
Beltrão, J., Jesus, S.B., Silva, V., Sousa, P.B., Carvalho, I., Trindade, D., Rodrigues, M.H.,
Machado, A., 2002. Efﬁciency of triple emitter source (TES) for irrigation
experiments of horticultural crops. Acta Hortic. 573, 183–188.
Bresler, E., McNeal, B.L., Cárter, D.L., 1982. Saline and Sodic Soils. Principles–
Dynamics–Modeling. Advanced Series in Agricultural Sciences 10. SpringerVerlag.
Cote, C.M., Bristow, K.L., Charlesworth, P.B., Cook, F.J., Thorburn, P.J., 2003. Analysis
of soil wetting and solute transport in subsurface trickle irrigation. Irrig. Sci. 22,
143–156.
Crevoisier, D., Popova, Z., Mailhol, J.C., Ruelle, P., 2008. Assessment and simulationof water and nitrogen transfer under furrow irrigation. Agric. Water Manage.
95, 354–366.
Feddes, R.A., Raats, P.A.C., 2004. Parameterizing the soil–water–plant–root system.
In: Feddes, R.A., de Rooij, G.H., van Dam, J.C. (Eds.), Proceedings of theUnsaturated Zone Modelling: Progress, Challenges and Applications.
Wageningen UR Frontis Series, vol. 6. Kluwer Academic Publishers, Dordrecht,
The Netherlands, pp. 95–141.
Feddes, R.A., Kowalik, P.J., Zaradny, H., 1978. Simulation of ﬁeld water use and crop
yield. Simulation Monographs Pudoc., Wageningen, The Netherlands.
Forkutsa, I., Sommer, R., Shirokova, Y.I., Lamers, J.P.A., Kienzler, K., Tischbein, B.,
Martius, C., Vlek, P.L.G., 2009. Modeling irrigated cotton with shallow
groundwater in the Aral Sea Basin of Uzbekistan: II. Soil salinity dynamics.
Irrig. Sci. 27, 319–330.
Gärdenäs, A., Hopmans, J.W., Hanson, B.R., Šimůnek, J., 2005. Two-dimensionalmodeling of nitrate leaching for various fertigation scenarios under microirrigation. Agric. Water Manage. 74, 219–242.
Gonçalves, M.C., Leij, F.J., Schaap, M.G., 2001. Pedotransfer functions for solutetransport parameters of Portuguese soils. Eur. J. Soil Sci. 52, 563–574.
Gonçalves, M.C., Šimůnek, J., Ramos, T.B., Martins, J.C., Neves, M.J., Pires, F.P., 2006.
Multicomponent solute transport in soil lysimeters with waters of different
quality. Water Resour. Res. 42, W08401, doi:10.1029/2005WR004802.
Gonçalves, J.M., Pereira, L.S., Fang, S.X., Dong, B., 2007. Modelling and multicriteria
analysis of water saving scenarios for an irrigation district in the upper Yellow
River Basin. Agric. Water Manage. 94, 93–108.
Halbertsma, J.M., Veerman, G.J., 1994. A New Calculation Procedure and Simple Setup for the Evaporation Method to Determine Soil Hydraulic functions. DLO
Winand Staring Centre, Report 88, Wageningen, The Netherlands.
Hanson, B.R., Šimůnek, J., Hopmans, J.W., 2006. Evaluation of urea–ammonium–
nitrate fertigation with drip irrigation using numerical modeling. Agric. Water
Manage. 86, 102–113.
Hanson, B.R., Šimůnek, J., Hopmans, J.W., 2008. Leaching with subsurface drip
irrigation under saline, shallow ground water conditions. Vadose Zone J. 7, 810–
818.
Hendrilsen, A., Selmer-Olsen, A.R., 1970. Automatic methods for determination of
nitrate and nitrite in water and soil extracts. The Analyst 95, 514–518.
Jarvis, N.J., 1994. The MACRO Model (Version 3.1). Technical Description andSample Simulations. Reports and Dissert. 19. Dept. Soil Sci., Swedish Univ. Agric.
Sci., Uppsala, Sweden, 51pp.
Legates, D.R., McCabe Jr., G.J., 1999. Evaluating the use of ‘‘goodness-of-ﬁt’’measures in hydrologic and hydroclimatic model validation. Water Resour.
Res. 35, 233–241.
Machado, R., Oliveira, M.R., 2003. Comparison of tomato root distribution byminirhizotron and destructive sampling. Plant Soil 255, 375–385.
Mallants, D., Vanclooster, M., Meddahi, M., Feyen, J., 1994. Estimating solutetransport in undisturbed soil columns using time domain reﬂectrometry. J.
Contam. Hydrol. 17, 91–109.
Mass, E.V., 1990. Crop salt tolerance. In: Tanji, K.K. (Ed.), Agricultural Salinity
Assessment and Management. Manual Eng. Pract., vol. 71. Am. Soc. of Civ. Eng.,
Reston, VA, pp. 262–304.
McNeal, B.L., Oster, J.D., Hatcher, J.T., 1970. Calculation of electrical conductivity
from solution composition data as an aid to in-situ estimation of soil salinity.
Soil Sci. 110, 405–414.
Pereira, L.S., Gonçalves, J.M., Dong, B., Mao, Z., Fang, S.X., 2007. Assessing basin
irrigation and scheduling strategies for saving irrigation water and controllingsalinity in the upper Yellow River Basin, China. Agric. Water Manage. 93, 109–
122.
Pereira, L.S., Cordey, I., Iacovides, I., 2009. Coping with Water Scarcity. Addressing
the Challenges. Springer.
Ramos, T.B., Gonçalves, M.C., Castanheira, N.L., Martins, J.C., Santos, F.L., Prazeres, A.,
Fernandes, M.L., 2009. Effect of sodium and nitrogen on yield function ofirrigated maize in southern Portugal. Agric. Water Manage. 96, 585–594.
Ritchie, J.T., 1972. Model for predicting evaporation from a row crop withincomplete cover. Water Resour. Res. 8, 1204–1213.
Roberts, T., White, S.A., Warrick, A.W., Thompson, T.L., 2008. Tape depth andgermination method inﬂuence patterns of salt accumulation with subsurface
drip irrigation. Agric. Water Manage. 95, 669–677.
Roberts, T., Lazarovitch, N., Warrick, A.W., Thompson, T.L., 2009. Modeling salt
accumulation with subsurface drip irrigation using HYDRUS-2D. Soil Sci. Soc.
Am. J. 73, 233–240.
Searle, P.L., 1984. The Berthlot or Indophenol reaction and its use in the analysis
chemistry for nitrogen. The Analyst 109, 549–565.

Šimůnek, J., Hopmans, J.W., 2009. Modeling compensated root water and nutrient
uptake. Ecol. Model. 220, 505–521.
Šimůnek, J., Suarez, D.L., 1994. Two-dimensional transport model for variably
saturated porous media with major ion chemistry. Water Resour. Res. 30, 1115–
1133.
Šimůnek, J., Valocchi, A.J., 2002. Geochemical transport. In: Dane, J.H., Topp, G.C.
(Eds.), Methods of Soil Analysis, Part 1. Physical Methods. Soil Sci. Soc. of Am.,
Madison, Wis., USA, pp. 1511–1536.
Šimůnek, J., van Genuchten, M. Th., Šejna, M., 2006. The HYDRUS Software Package
for Simulating Two- and Three-Dimensional Movement of Water, Heat, andMultiple Solutes in Variably-saturated Media, Technical Manual. Version 1.0, PC
Progress, Prague, Czech Republic, 241pp.
Šimůnek, J., Suarez, D.L., Šejna, M., 1996. The UNSATCHEM Software Package forSimulating One-dimensional 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, Res.
Rep. 141, US Salinity Lab., Agric. Res. Serv., Riverside, Calif., 186pp.
Šimůnek, J., van Genuchten, M.Th., Šejna, M., Toride, N., Leij, F.J., 1999. The
STANMOD Computer Software for Evaluating Solute Transport in Porous Media
Using Analytical Solutions of Convection–Dispersion Equation. Versions 1.0 and
2.0, IGWMC-TPS-71, International Ground Water Modeling Center, Colorado
School of Mines, Golden, Colorado, 32pp.
Šimůnek, J., Šejna, M., Saito, H., Sakai, M., van Genuchten, M.Th., 2008a. The
HYDRUS-1D Software Package for Simulating the Movement of Water, Heat,and Multiple Solutes in Variably-saturated Media. Version 4.0. HYDRUS
Software Series 3, Department of Environmental Sciences, University ofCalifornia Riverside, Riverside, CA, USA, 315pp.
Šimůnek, J., van Genuchten, M.Th., Šejna, M., 2008b. Development and Applicationsof the HYDRUS and STANMOD Software Packages, and Related Codes. Vadose
Zone J., Special Issue ‘‘Vadose Zone Modeling’’, 7, 587–600. doi:10.2136/
VZJ2007.0077.
Skaggs, T.H., Shouse, P.J., Poss, J.A., 2006. Irrigating forage crops with saline waters:
2. Modeling root uptake and drainage. Vadose Zone J. 5, 824–837.
Soil Survey Staff, 2006. Keys to Soil Taxonomy, 10th ed. Natural Resources
Conservation Service, United States Department of Agriculture, Washington,
DC.
Steppuhn, H., van Genuchten, M.Th., Grieve, C.M., 2005. Root-zone salinity: II.
Indices for tolerance in agricultural crops. Crop Sci. 45, 221–232.
Stolte, J., 1997. Determination of the saturated hydraulic conductivity using theconstant head method. In: Stolte, J. (Ed.), Manual for Soil PhysicalMeasurements. Technical Document 37, DLO Winand Staring Centre,
Wageningen, The Netherlands.
Toride, N., Leij, F.J., van Genuchten, M.Th., 1995. The CXTFIT Code for Estimating
Transport Parameters from Laboratory or Field Tracer Experiments. Version 2.0.
Research Report 137, US Salinity Laboratory, Riverside, California, USA.
Truesdell, A.H., Jones, B.F., 1974. Wateq, a computer program for calculatingchemical equilibria of natural waters. J. Res. US Geol. Surv. 2, 233–248.
US Salinity Laboratory Staff, 1954. Diagnosis and Improvement of Saline and Alcaly
Soils. USDA Handbook 60, Washington, USA.
van Dam, J.C., Huygen, J., Wesseling, J.G., Feddes, R.A., Kabat, P., van Valsum, P.E.V.,
Groenendijk, P., van Diepen, C.A., 1997. Theory of SWAP, Version 2.0. Simulationof Water Flow, Solute Transport and Plant Growth in the Soil–Water–
Atmosphere–Plant Environment. Dept. Water Resources, WAU, Report 71,
DLO Winand Staring Centre, Wageningen, Technical Document 45.
van den Berg, M., Driessen, P.M., Rabbinge, R., 2002. Water uptake in crop growthmodels for land use systems analysis. II. Comparison of three simple
approaches. Ecol. Model. 148, 233–250.
van Genuchen, M.Th., 1987. A numerical model for water and solute movement inand below the root zone. Res. Rep. 121, US Salinity Laboratory, USDA, ARS,
Riverside, California.
van Genuchen, M.Th., Šimůnek, J., 2004. Integrated modeling of vadose zone ﬂow
and transport processes. In: Feddes, R.A., de Rooij, G.H., van Dam, J.C. (Eds.),
Proc. Unsaturated Zone Modelling: Progress, Challenges and Applications,
Wageningen, The Netherlands, pp. 37–69 (October 3–5).
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., Leij, F.J., Yates, S.R., 1991. The RETC Code for Quantifying theHydraulic Functions of Unsaturated Soils. Report No. EPA/600/2-91/065 . R.S.
Kerr Environmental Research Laboratory, US Environmental Protection Agency,
Ada, OK, 85pp.
Wagenet, R.J., Hutson, J.L., 1987. LEACHM: Leaching Estimation and Chemistry
Model, A Process-based Model of Water and Solute Movement,
Transformations, Plant Uptake and Chemical Reactions in the Unsaturated
Zone, Continuum 2, Report. Water Resour. Inst., Cornell Univ., Ithaca, NY.
Weihermüller, L., Kasteel, R., Vanderborght, J., Pütz, T., Vereecken, H., 2005. Soilwater extraction with a suction cup: results of numerical simulations. Vadose
Zone J. 4, 899–907.
Wesseling, J.G., Elbers, J.A., Kabat, P., van den Broek, B.J., 1991. SWATRE: Instructions
for Input. Report, Winand Staring Cent., Wageningen, Netherlands.
Wind, G.P., 1968. Capillary conductivity data estimated by a simple method. In:
Rijtema, P.E., Wassink, H. (Eds.), Water in the Unsaturaded Zone. Proceedings ofa Symposium, vol. I. June 1966, Wagningen, The Netherlands. IASH/AIHS–UNESCO, pp. 181–191.