PUBLICATIONS
Water Resources Research
RESEARCH ARTICLE
10.1002/2016WR018930
Key Points:
 Three-dimensional multicomponent
solute transport model
parameterization
 Multicomponent solute transport in
drip-irrigated lysimeters
 Lysimeter design effects on leaching
fractions and ESP
Supporting Information:
 Supporting Information S1
Correspondence to:N. Lazarovitch,lazarovi@bgu.ac.il
Citation:
Raij, I., J. Simůnek, A. Ben-Gal, and
N. Lazarovitch (2016), Water ﬂow andmulticomponent solute transport in
drip-irrigated lysimeters, Water Resour.
Res., 52, 6557–6574, doi:10.1002/
2016WR018930.
Received 15 MAR 2016
Accepted 7 AUG 2016
Accepted article online 16 AUG 2016
Published online 25 AUG 2016

Water flow and multicomponent solute transport in dripirrigated lysimeters
2, Alon Ben-Gal3, and Naftali Lazarovitch1

Iael Raij1, Jirı Simůnek
1
Wyler Department of Dryland Agriculture, French Associates Institute for Agriculture and Biotechnology of Drylands,
Jacob Blaustein Institutes for Desert Research, Ben-Gurion University of the Negev, Sede Boqer Campus, Midreshet BenGurion, Israel, 2Department of Environmental Sciences, University of California, Riverside, Negev, California, USA, 3Instituteof Soil, Water and Environmental Sciences, Agricultural Research Organization, Gilat Research Center, Israel

Abstract Controlled experiments and modeling are crucial components in the evaluation of the fate ofwater and solutes in environmental and agricultural research. Lysimeters are commonly used to determine
water and solute balances and assist in making sustainable decisions with respect to soil reclamation, fertilization, or irrigation with low-quality water. While models are cost-effective tools for estimating and preventing environmental damage by agricultural activities, their value is highly dependent on the accuracy of
their parameterization, often determined by calibration. The main objective of this study was to use measured major ion concentrations collected from drip-irrigated lysimeters to calibrate the variably saturated
water ﬂow model HYDRUS (2D/3D) coupled with the reactive transport model UNSATCHEM. Irrigation alternated between desalinated and brackish waters. Lysimeter drainage and soil solution samples were collected for chemical analysis and used to calibrate the model. A second objective was to demonstrate thepotential use of the calibrated model to evaluate lower boundary design options of lysimeters with respect
to leaching fractions determined using drainage water ﬂuxes, chloride concentrations, and overall salinity ofdrainage water, and exchangeable sodium percentage (ESP) in the proﬁle. The model showed that, in thelong term, leaching fractions calculated with electrical conductivity values would be affected by the lower
boundary condition pressure head, while those calculated with chloride concentrations and water ﬂuxes
would not be affected. In addition, clear dissimilarities in ESP proﬁles were found between lysimeters withdifferent lower boundary conditions, suggesting a potential inﬂuence on hydraulic conductivities and ﬂow
patterns.

1. Introduction

C 2016. American Geophysical Union.
V

All Rights Reserved.

RAIJ ET AL.

A decrease in the availability of freshwater resources and an increase in food demand have led to the widespread utilization of lower-quality water for irrigation. The use of water with high salt concentrations necessitates irrigation in amounts higher than actual evapotranspiration (ET) in order to prevent soil salinizationand subsequent yield loss [Rawlins, 1973]. Traditionally, irrigation water requirements, designed to minimize
the negative effects of salinity on crop growth and yield when irrigating with brackish or high-salinity water,
are determined by multiplying potential ET by a crop/cover factor and dividing by one minus a leaching
fraction [Allen et al., 1998; Corwin et al., 2007; Dudley et al., 2008a]. This approach is problematic for real-timescheduling of crop irrigation for a number of reasons. Miscalculations and high daily variability of ET or timing and magnitude of rain events can lead to inefﬁcient use of water and fertilizers, resulting in economic
losses and to groundwater pollution from leached agrochemicals [Oren et al., 2004; Corwin et al., 2007].
Leaching fractions are commonly based on crop response functions determined for an entire season andaverage root zone salinity, which may not be relevant for all actual temporal or spatial ﬁeld conditions [Leteyet al., 2011; Groenveld et al., 2013]. Crop factors have additionally been shown to be sensitive to a reduction
in transpiration due to salinity and other types of stresses [Rhoades et al., 1992; Bhantana and Lazarovitch,2010]. Negative effects from these shortcomings can be minimized by scheduling irrigation and fertilization
based on in situ monitoring of water and solute balances and by using pressurized irrigation systems, whichallow frequent automated applications. Fertigation via drip systems has been shown to increase the efﬁciency of water and fertilizer applications [e.g., Hanson et al., 2006], to reduce salinity and drought stress in
the root zone [e.g., Hanson et al., 2008], and to prevent leaching of fertilizers and agrochemicals into thegroundwater [Rhoades et al., 1992; Hartz and Hochmuth, 1996].

WATER FLOW AND MULTICOMPONENT SOLUTE TRANSPORT

6557

Water Resources Research

10.1002/2016WR018930

In situ monitoring of salts, fertilizers, and other contaminants leaching into deep layers of the soil proﬁlecan provide information for improved decision making with respect to efﬁcient and sustainable irrigation
practices. Monitoring techniques include soil solution sampling with suction devices [Rhoades et al., 1999;
Corwin, 2002; Weiherm€uller et al., 2007; Wang et al., 2012; Turkeltaub et al., 2015], destructive soil sampling€m,and analysis, soil salinity sensors [Vanclooster et al., 1995; Rhoades et al., 1999], and lysimeters [Bergstro
1990]. Lysimeters isolate a volume of soil from its surroundings and allow monitoring of drainage amountand composition. The water and solute balances in lysimeters are obtained by monitoring exactly how
much water and solutes enter, exit, and remain in the system. Lysimeters are commonly used in both thelaboratory and the ﬁeld to estimate the fate of pesticides and herbicides [Winton and Weber, 1996; Schoen
et al., 1999]. Lysimeters are also widely used to determine crop factors (Kc) [Marek et al., 2006; Bhantana andLazarovitch, 2010; Bryla et al., 2010] for irrigation scheduling based on meteorological data.
Monitoring of the chemical composition of drainage and solute loads allows the use of lysimeters as man~alver et al., 2015] and salt leaching [Tripler et al., 2012].
agement tools for fertigation scheduling [Ruiz-Pen
However, the boundary condition at the bottom of a lysimeter may limit or otherwise inﬂuence water ﬂow
and solute transport and thus, subsequently, water drainage and salt leaching into the underlying soil proﬁle [Flury et al., 1999; Ben-Gal and Shani, 2002].
Lysimeters are generally categorized into two groups according to their lower boundary condition (LBC)€m, 1990]. First, free-drainage lysimeters are characterized by an
and drainage collection method [Bergstro
atmospheric pressure at the bottom of the lysimeter, with drainage from them occurring only when saturation is reached. Studies comparing solute transport in free-drainage lysimeters and under ﬁeld conditions
reported differences not only due to the LBC but also in relation to soil texture, water ﬂuxes, pore watervelocity, soil heterogeneity, and sorption [Flury et al., 1999; Gasser et al., 2002; Abdou and Flury, 2004]. Twodimensional numerical simulations showed that bromide moved faster in free-drainage lysimeters than in
the ﬁeld when lysimeters contained soil with vertical heterogeneities rather than with horizontal or isotropic
heterogeneities [Abdou and Flury, 2004]. Numerical experiments reported by Flury et al. [1999] also showedthat the coarser the soil and/or the smaller the water ﬂux, the more pronounced the differences in the
breakthrough of solutes between lysimeters and the ﬁeld. This effect is due to the higher water contents inthe lysimeter proﬁle and therefore lower pore water velocities.
Second, suction-controlled lysimeters prevent saturated conditions at the lower soil boundary by applying anegative pressure head. Suction lysimeters can be either active or passive. Active suction-controlled lysimeters allow application of the pressure head, equivalent to that measured by a tensiometer under ﬁeld conditions, at the bottom boundary [e.g., Hannes et al., 2015; Groh et al., 2016]. This type of lysimeter has high€m, 1990].
installation and maintenance costs and is difﬁcult to repair when malfunctioning [Bergstro
In passive suction-controlled lysimeters, highly conductive media is used to connect the soil at the bottomof the lysimeter with the lysimeter’s drainage exit and to extend it below this point [Ben-Gal and Shani,
2002]. These extensions hydraulically deepen the soil using a relatively small volume of material and movethe saturated conditions to the end of the extensions themselves. The extensions are designed to not inﬂuence water ﬂow in and from the soil above, and to maintain water content in the lysimeters similar to that
expected in the soil proﬁle under natural ﬁeld conditions. This allows the use of a shallower soil component
in lysimeters for agricultural research, simplifying installation and operation. While the effects of different
lysimeter LBCs on water ﬂow as a function of soil hydraulic properties are well understood [Ben-Gal andShani, 2002], effects on drainage water composition or on dissolved and adsorbed solute distribution in the
proﬁle are not well understood and need to be further explored. It should be taken into account that inboth passive and active suction-controlled lysimeters, materials used to apply the lower boundary suction€m, 1990; Weiherm€
may react with dissolved salts in the drainage solution [Bergstro
uller et al., 2007].
Models can be cost-effective tools for evaluating and consequently preventing environmental damage by
agricultural activities. However, model accuracy and reliability are highly inﬂuenced by the quality of theirparameterization, often determined by calibration and validation [Roberts et al., 2009; Ramos et al., 2011;
Skaggs et al., 2014]. Using replicated lysimeters for calibration of transport models is a good compromise
between the accuracy of laboratory columns and the representativeness of ﬁeld measurements [Skaggs
et al., 2012]. Models such as HYDRUS-1D [Simůnek et al., 2008] coupled with UNSATCHEM [Suarez andSimůnek, 1997] have often been used to assess the implications and risks arising from various agricultural

RAIJ ET AL.

WATER FLOW AND MULTICOMPONENT SOLUTE TRANSPORT

6558

Water Resources Research

10.1002/2016WR018930

practices such as irrigation with low-quality water in arid and semiarid regions [Gonçalves et al., 2006; Kaledhonkar and Keshari, 2006; Ramos et al., 2011; Rasouli et al., 2012; Skaggs et al., 2014]. However, the use of
one-dimensional models to represent three-dimensional patterns created under drip irrigation may lead toerroneous conclusions [Skaggs et al., 2004; Warrick and Lazarovitch, 2007; Chen et al., 2010]. Evaluation of
the reliability of lysimeters for monitoring solute leaching from agricultural lands irrigated with brackish
water is expected to beneﬁt from multidimensional solute transport and water movement modeling thatconsiders various reactive processes, such as cation exchange and/or precipitation/dissolution.
The development of the UNSATCHEM module for the HYDRUS (2D/3D) software [Simůnek et al., 2012],
which allows consideration of water ﬂow and solute transport with complex major ion chemistry in unsaturated conditions, has opened a window to explore the effects of the LBC on the quality of drainage fromlysimeters and the distribution of ions in the soil proﬁle of the lysimeter. This solution takes into consideration the complex interactions between the soil matrix and the major ions in the soil.
The objectives of this study thus were (a) to carry out water ﬂow and solute transport experiments involving
drip irrigation with brackish and desalinated water in lysimeter systems, (b) to use collected experimental
data to calibrate the three-dimensional ﬂow and transport model HYDRUS (2D/3D) [Simůnek et al., 2008]
coupled with the reactive transport model UNSATCHEM [Simůnek et al., 2012], and (c) to use the calibrated
three-dimensional reactive transport model to evaluate possible design options for the bottom boundary inlysimeters and their effects on (1) transient leaching fractions deﬁned by the daily ratio of drainage and irrigation, either by volume, electrical conductivity, or chloride concentration and (2) spatial ESP (exchangeablesodium percentage) proﬁles in the lysimeter.

2. Experiments
2.1. Laboratory Column Experiments for Estimating Soil Retention Curve Parameters
A 0.85 m long PVC column with a diameter of 0.15 m was ﬁlled with loamy sand soil in intervals of 5 cm inorder to maintain the same bulk density as in the lysimeters. The column was slowly wetted from below
and once saturated a pressure head of h 5 0 was imposed as the LBC using a Mariotte bottle. In addition,
the top of the column was covered to prevent any evaporation. Once equilibrium was achieved, the column
was opened and soil samples were taken every 5 cm and then oven dried for gravimetric water content calculation. In addition, a tensiometer-like device was ﬁlled with a known volume of oven dried loamy sandsoil. This device consisted of a plastic cylinder 0.3 m long and a diameter of 0.02 m, with a ceramic cup at
the lower boundary and a connection to pressurized nitrogen gas at the upper boundary. The device wassaturated with water from below and then a pressure of 0.5 bar was applied at the top until equilibrium
was reached, followed by another pressure application of 0.8 bar until equilibrium was reached again. Equilibrium was determined by weighing the device until its mass was stable and the gravimetric water content
of the soil in the device was calculated. During this stage, the ceramic cup was covered with a plastic ﬁlm toprevent evaporation.
2.2. Lysimeter Experiments
Water ﬂow and solute transport experiments were conducted in an automated rotating lysimeter system[Lazarovitch et al., 2006], located in a greenhouse at The Jacob Blaustein Institutes for Desert Research, Sede
Boqer Campus, Israel. The system consisted of eight lysimeters located in a circular arrangement. All lysimeters had the same surface area (0.70 m length and 0.35 m width) but four of them had a depth of 0.5 m andthe other four had a depth of 0.25 m. The shallow lysimeters were used for the calibration of the Gapon
exchange parameters and the deep lysimeters for the model validation. Each lysimeter had an independent
automated irrigation delivery system and drainage collection system, and was located on a load cell for continuous weighing. The system rotated at two revolutions per hour and automatically stopped for either ascheduled irrigation application or drainage collection when a lysimeter was positioned at the irrigation/
drainage collection station. All lysimeters were ﬁlled with loamy sand soil. Leachate ﬂowed through two
extensions (0.4 m long and 5 cm diameter) ﬁlled with rockwool [Ben-Gal and Shani, 2002] into a 0.0065 m3
collection container. The two extensions were located at X 5 17.5 cm and Y 5 25 and 45 cm, respectively.
Lysimeter geometry and drainage extension locations are described in Figure 1. Evaporation was prevented
by plastic sheeting placed on the soil surface. The plastic mulch was covered with 2 cm of dry soil to minimize temperature variations. In each lysimeter, a ceramic suction sampling device (a ceramic cup, 6 cm

RAIJ ET AL.

WATER FLOW AND MULTICOMPONENT SOLUTE TRANSPORT

6559

Water Resources Research

10.1002/2016WR018930

Figure 1. A schematic showing the simulated transport domain representing the lysimeter (0.7 m length, 0.35 m width, 0.5 m depth) in arotating lysimeter system, ﬁnite element nodes, and the boundary conditions used in HYDRUS (2D/3D). Gray nodes are those with a noﬂow boundary condition (BC), red nodes are those with a time-variable ﬂux BC (top boundary), and blue nodes are those with a seepage
face BC (bottom boundary). The yellow point represents the soil solution suction device.

long and 1.5 cm in diameter) was located in the middle of the container at depths of 0.25 and 0.125 m indeep and shallow lysimeters, respectively, to collect pore water samples every 2 or 3 days.
One liter of irrigation solution was applied twice daily to each lysimeter via four nonpressure compensateddrippers (2 L h21, Netaﬁm, Israel). To ensure similar initial conditions, low-salinity water (DW in Table 1) wasapplied ﬁrst for 60 days until a constant drainage electrical conductivity (EC) was obtained in all lysimeters.
During the second step, irrigation with brackish groundwater (GW in Table 1) was applied for 6 days. The
time period of brackish water application was followed by low-salinity water irrigation until a condition ofsteady state was again reached regarding the EC of the drainage water. Drainage mass was measured automatically twice a day. Drainage water was manually sampled every 1 or 2 days, and the soil solution fromthe ceramic suction sampling device was collected every other weekday. Samples from the drainage andsoil pore solutions were analyzed for EC by means of a Multimeter MM 401, Crison; for sodium, magnesium,
calcium, potassium, and sulfate concentrations by inductively coupled plasma emission spectrophotometry(ICP); and for chloride concentration by Chloride Analyzer (926, Sherwood).

3. Multicomponent Solute Transport: HYDRUS (2D/3D) Coupled With UNSATCHEM
Advective-dispersive chemical transport under transient water ﬂow conditions in a partially saturated
porous medium is described in the model as follows [Simůnek et al., 2012]:
Table 1. Irrigation Water Composition: GW Is Brackish Groundwater and DW Is Desalinated Watera

GW
DW
a

RAIJ ET AL.

[Ca21][Mg21]

[Na1]

[K1]

[SO22
4 ]

[Cl2]

Measured EC (dS m21)

EC From Cations (dS m21)

10.13
1.72

9.76
0.38

13.89
1.52

0.33
0.02

7.81
1.08

22.67
1.37

3.7
0.36

3.4
0.36

Units are in meq L21 unless indicated otherwise. EC is electrical conductivity.

WATER FLOW AND MULTICOMPONENT SOLUTE TRANSPORT

6560

Water Resources Research

10.1002/2016WR018930



@hck
@c k
@^c k
@
@ck
@qi ck
2
1qb
1qb
5
hDij
k51; 2; . . . :; Nc ;@t@t@t @xi
@xj
@xi

(1)

where ck is the total dissolved concentration of the component k [ML23], c k is the total sorbed concentration of the component k [MM21], ^c k is the total concentration of the component k in the minerals, which
can precipitate or dissolve [MM21], qb is the bulk density of the soil [ML23], Dij is the dispersion coefﬁcienttensor [L2 T21], i and j are spatial dimensions, qi is the volumetric water ﬂux [LT21] (qi 5 ui h,where u is thepore water velocity), h is the volumetric water content [L3 L23] and Nc is the number of aqueous components. The dispersion coefﬁcient tensor in the liquid phase, Dij , is given by Bear [1972]:
hDij 5kT jqjdij 1ðkL 2kT Þ

qj qi
1hDw sw dij ;
jqj

(2)

where kL and kT are the longitudinal and transverse dispersivities, respectively [L], dij is the Kronecker delta
function (dij 5 1 if i 5 j, and dij 5 0 if i 6¼ j), Dw is the molecular diffusion coefﬁcient [L2 T21], and sw is a tortuosity factor in the liquid phase.
Cation exchange between aqueous and exchangeable phases is described using the Gapon equation, under
the assumption that the cation exchange capacity is constant and independent of pH [Simůnek et al., 2012]: 1=x
c y1
cjx1
i
Kij 5

1=y ;
c x1
ciy1
j

(3)

where Kij is the Gapon selectivity (KG) coefﬁcient, y and x are the valences of species i and j, respectively, c isthe exchangable concentration in mmolc kg21 of soil, and c is the soluble concentration in mmolc L21.
When all four cations (Mg21, Ca21, Na1, and K1) considered by UNSATCHEM are present, three Gapon coefﬁcients for cation pairs Mg/Ca, Ca/Na, and Ca/K are needed.
Calcite and gypsum precipitation and dissolution are described by equilibrium equations (equations (4) and(5)) in the presence of CO2. The UNSATCHEM module also includes precipitation and dissolution of nesquehonite, hydromagnesite, and sepiolite which were not included in this work.
CaSO4  2H2 O ! Ca21 1SO22
4 12H2 O;(4)

CaCO3 1CO2 ðgÞ1H2 O ! Ca21 12HCO2
3 :

(5)

G
C
The solubility products KSP[meq4 L24] and KSP[meq2 L22] for gypsum and calcite, respectively, are given by
equations (6) and (7) [Simůnek et al., 2012].



G
KSP
5 Ca21 SO22
ðH2 O Þ2 ;
(6)
4




C
KSP
5 Ca21 CO22
:
3

(7)

where parentheses represent ion activities. A more thorough explanation can be found in Suarez andSimůnek [1997] and Simůnek et al. [2012].
Numerical solution of the Richards equation, which describes water ﬂow in soils, requires the knowledge ofsoil hydraulic functions that relate the pressure head, the water content, and the hydraulic conductivity ofthe soil. Soil hydraulic parameters required for the water ﬂow model when the van Genuchten-Mualem
model [van Genuchten, 1980] is used are the residual water content, hr [L3 L23], the saturated water content,
hs [L3 L23], and the shape parameters a [L21], and n for the soil water retention function, and the saturated
hydraulic conductivity Ks [L T21], and the tortuosity factor l for the hydraulic conductivity function.

4. Parameterization of HYDRUS (2D/3D) and Its UNSATCHEM Module
The HYDRUS-1D and HYDRUS (2D/3D) software packages [Simůnek et al., 2008] were used to inversely estimate selected transport and reaction parameters using different techniques. As stated previously, parameterization processes can highly inﬂuence the quality and relevance of modeled results. Each calibration step
was performed at a speciﬁc dimension (1-D or 3-D) and with a different optimization tool while taking into

RAIJ ET AL.

WATER FLOW AND MULTICOMPONENT SOLUTE TRANSPORT

6561

Water Resources Research

10.1002/2016WR018930

consideration the computational times
and the scale/dimension dependency ofeach parameter to be calibrated.

Table 2. Soil Texture, Hydraulic and Transport Parameters and Exchange
Coefﬁcients of the Loamy Sand Soila
Soil Parameter
hr
hs
a
n
l
KS
qb
Sand
Clay
Silt
KG (Mg/Ca)
KG (Ca/Na)
KG (Ca/K)
CEC
kL
kT

Units

Value

(cm3 cm23)(cm3 cm23)(cm21)

0.004
0.36
0.016
3.43
0.5
145.44
1.5
89
4
7
0.43
0.14
0.60
43.7
1.7
0.17

(cm d21)(g cm23)
(%)
(%)
(%)
(L meq20.5)(L meq20.5)(meq kg21)
(cm)
(cm)

S.E. Coefficient

Hydraulic parameters of the loamy sandsoil were either estimated using the RETC
software [van Genuchten et al., 1991] ormeasured in the laboratory. The dispersivities, which are highly scale-dependent[Neuman, 1990; Vanderborght and Vereecken, 2007], were manually calibratedusing HYDRUS (2D/3D) in threedimensional simulations. The dimensionindependent Gapon exchange coefﬁcientswere optimized using the external software UCODE and HYDRUS-1D. HYDRUS-1D
was used rather than HYDRUS (2D/3D)
because of high computational requirements when running three-dimensional
simulations. Finally, the initial adsorbed
and precipitated concentrations were estimated by trial and error.

0.007
0.0006
0.36

0.08
0.04
0.1

a
hr, residual water content; hs, saturated water content; a and n, empiricalshape parameters; l, tortuosity factor; KS, saturated hydraulic conductivity;
qb, bulk density, KG, Gapon exchange coefﬁcients for different cation pairs;
CEC, cation exchange capacity; kL and kT, longitudinal and transverse dispersivities, respectively.

Optimized parameters using either RETC, UCODE, or HYDRUS were obtained by minimizing the objectivefunction deﬁned as the sum of weighted least squares. This function compares each observed data point tothe corresponding modeled value
SðbÞ5

n
X

2

xi ½yi 2y 0 i ðbÞ ;

(8)

i51

where S(b) is the objective function, b is a vector of optimized parameters, n is the number of data points,
xi is the weight associated with each data point, yi is the observed data point, and y 0 i ðbÞ is the data point
calculated using the parameters b.
Constant weights were used in the RETC and HYDRUS (2D/3D) optimizations. In the UCODE optimization,
the coefﬁcients of variation, calculated as the ratio between the averaged and standard deviations of thefour replications, were used as weights.
4.1. Soil Hydraulic Parameters
The hs, a, and n parameters of the van Genuchten-Mualem hydraulic model [van Genuchten, 1980] were estimated using the RETC computer program [van Genuchten et al., 1991] (Table 2). Water contents at 16 depths
under hydrostatic conditions from the 0.85 m long column and water contents at pressure heads of 25 and28 m were used as input data for the objective function of the optimization. The saturated hydraulic conductivity, KS, was determined using the constant head method [Ali, 2010] in 20 cm long backﬁlled soil columns with a diameter of 5 cm. The residual water content, hr, was assumed to be equal to the measured
soil water content at a soil water pressure head of 28 m. Bulk density, qb, was measured using undisturbed
soil core samples. The shape parameter l was assumed to be 0.5 as used before for similar soils [Ityel et al.,2011]. The hr, KS, and l parameters were used as prior information for the inverse problem. A summary of
the hydraulic properties and their related statistics can be found in Table 2.
4.2. Gapon Exchange Parameters
Gapon exchange parameters can be calculated from measured values of soluble and exchangeable concentrations using equation (3) [Gonçalves et al., 2006; Ramos et al., 2011]. However, since some of these concentrations are relatively low, small measurement errors can have a great inﬂuence on calculated values [Smiles
and Smith, 2004]. The Gapon exchange coefﬁcients can also be estimated from a predeﬁned range ofknown values for a speciﬁc soil [MDH Engineered Solutions Corp., 2003; Simůnek et al., 2012] or inversely estimated from measured soil solution and adsorbed concentrations [Jacques et al., 2012]. Here the Gapon

RAIJ ET AL.

WATER FLOW AND MULTICOMPONENT SOLUTE TRANSPORT

6562

Water Resources Research

10.1002/2016WR018930

Table 3. Initial Adsorbed and Precipitated Concentrations for the 1-D and 3-D Simulations
Adsorbed Concentrations (meq kg21)

Precipitated Concentrations (meq kg21)

Ca21

Mg21

Na1

K1

CalciteGypsum

8.7

20

14

1

3197

27 (1-D)/22 (3-D)

coefﬁcients were estimated using the general optimization code UCODE_2014 [Poeter et al., 2005, 2014]
and HYDRUS-1D coupled with UNSATCHEM. As shown in Poeter et al. [2014], UCODE can be used to inversely estimate parameters that cannot be estimated using the built-in inverse module in HYDRUS-1D. The
observation data were obtained by averaging the drainage solute concentrations from the four shallow
lysimeters. Five different solute concentrations (Na1, Ca21, Mg21, K1, and SO22
4 ), and EC values measured
in the drainage solution at 19 different time points, resulted in a total of 114 data points. Modeled concentrations obtained at an observation point located at the lower boundary of a 1-D proﬁle were ﬁtted against
the measured drainage concentrations data. This approach assumes that the resident concentration at theobservation point in the seepage face equals the modeled drainage concentration. Estimated values of theGapon coefﬁcients and associated statistics are given in Table 2.
Adsorbed Ca21, Mg21, Na1, and K1, as well as precipitated calcite and gypsum, concentrations were measured in soil samples obtained in the same area as the soil used for the experiment and were used as a starting point for the simulations. They were all, except for the calcite concentrations, changed as needed inorder to match the ﬁrst drainage concentration measurement by the model (Table 3). This exercise was performed in 1-D due to its shorter computational times, the relative scale independence of exchange coefﬁcients and the assumption of homogeneous adsorbed and precipitated initial concentrations along thelysimeter proﬁle. The 1-D proﬁle had a depth of 0.25 m and a LBC deﬁned as a seepage face with a pressure
head of 230 cm. Initial pressure heads and soil solution concentrations were uniform along the proﬁle withvalues of 2100 cm and the solution composition of GW (Table 1), respectively.

Objective function S(b)4.3. Longitudinal and Transverse Dispersivities
Longitudinal dispersivity (kL) for the deep lysimeters was optimized in 3-D by manually ﬁnding the levelwithin a range between 0.5 and 5 cm best ﬁtting measured results. The transverse dispersivity (kT) was maintained at one tenth of the kL as suggested in the literature [Bear, 1972; Skaggs and Leij, 2002; Simůnek and van
Genuchten, 2006]. Results from a total of 17 simulations were compared to observed data and the objectivefunction was calculated according to equation (8). A total of 56 data points were included in the objectivefunction comparing modeled and measured chloride concentrations of drainage solution at 38 times and 18
resident concentrations representing suction cup solution. In coarse soils, resident concentrations can be used
to represent pore water solutions
0.16[Jacobsen et al., 1992]. The lowest
objective function value, asdeﬁned in equation (8), was found
0.14
for kL 5 1.7 cm and kT 5 0.17 cm(Figure 2). Such manual calibration
enabled a less tedious optimizationin 3-D because the transverse dis0.12
persivity was a function of the longitudinal dispersivity, eliminating theneed to determine them separately.
0.1

0.08
0

1

2

3

4

5

λL (cm)RAIJ ET AL.

5. Modeling ThreeDimensional Multicomponent Solute Transport

Figure 2. Objective function as a function of dispersivity. The black dot represents theminimum of the objective function.

Optimized transport and exchange
parameters were used to model
the three-dimensional water ﬂow

WATER FLOW AND MULTICOMPONENT SOLUTE TRANSPORT

6563

Water Resources Research

10.1002/2016WR018930(a) Ca2+

1
0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

Simulated
Measured
Irrigation

0
0

10

20

30

40

50

0

60

(c) Na+1

Ci/Co(b) Mg2+

1

10

20

30

40

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

60

(d) SO42-

1

0.8

50

0
0

10

20

30

40

50

60

0

10

20

30

40

50

60

2
1

(f) K+(e) Cl1.5

0.8
0.6

1
0.4
0.5
0.2
0

0
0

10

20

30

40

50

60

0

10

20

30

40

50

60

(g) ECEC (dS/m)

3

2

1

0
0

10

20

30

40

50

60

DOE (days)

Figure 3. Simulated (dashed lines) and measured (symbols) relative concentrations of major ions in the drainage water. Error bars represent one standard deviation. DOE stands for day of experiment.

and multicomponent solute transport in the lysimeter experiment. The three-dimensional transport domain
had surface dimensions of 0.35 and 0.7 m and a depth of 0.5 m (Figure 1). An observation point was deﬁned
in the center of the lysimeter at a depth of 0.25 m, at the location of the soil pore water suction device (Figure1). Twenty horizontal equally spaced planes were discretized using an unstructured ﬁnite element mesh,resulting in a total of 14,640 nodes and 26,068 three-dimensional tetrahedral elements (Figure 1).
The initial condition for water ﬂow was deﬁned as hydrostatic equilibrium with a pressure head of 230 cm
at the bottom of the soil proﬁle. The initial soil solution was assumed to have the same solute concentrations as the GW because the ﬁeld, from which the soil was obtained, was previously irrigated with thatwater quality. Adsorbed and precipitated concentrations were the same as in 1-D simulations, except for
the gypsum concentrations, which had to be reduced to 22 meq kg21 in the 3-D domain where its leaching
is not as effective as in 1-D. The initial CO2 concentration was assumed to be equal to average atmosphericvalues (400 ppm) [NOAA-ESRL, 2015].

RAIJ ET AL.

WATER FLOW AND MULTICOMPONENT SOLUTE TRANSPORT

6564

Water Resources Research

10.1002/2016WR018930

1.2(a) Ca2+(b) Mg2+

1.2
0.8

Simulated
Measured
Irrigation

0.8

0.4

0.4

0
0

10

20

30

40

50

60

(c) Na+

1

0

20

30

40

50

60

(d) SO42-

1

0.6

0.6

0.4

0.4

0.2

0.2

0

0
0

10

20

30

40

1

50

60

0

10

20

30

40

2

(e) Cl-

0.8

1.6

0.6

1.2

0.4

0.8

0.2

0.4

0

50

60

(f) K+

0
0

10

20

30

40

50

60

0

10

20

30

40

50

60

(g) EC3

EC (dS/m)

10

0.8

Ci/Co

0.8

0

2

1

0
0

10

20

30

40

50

60

DOE (days)
Figure 4. Simulated (dashed lines) and measured (symbols) relative concentrations of major ions in the soil solution collected using ceramic suction devices. Error bars represent one standard deviation. DOE stands for day of experiment.

The upper BC was deﬁned as no ﬂow, except for four separate areas of 2.4 cm2 each representing drippers,
in which a variable ﬂux BC with speciﬁed solute concentrations was assigned. Two liters of water per day ofvarying quality were applied during the entire simulation as in the actual experiment. The ﬂushing periodwith DW (ﬁrst 60 days) water is not shown in the results and the day of experiment (DOE) begins 6 days
before the pulse with GW starts. The lower boundary was deﬁned as no ﬂow, except for two circles, representing the outﬂow drainage extensions, for which a seepage face BC with the pressure head of 230 cm
was speciﬁed (Figure 1). This pressure head at the lower BC was an average of measurements in lysimeterswith the same LBCs. All remaining boundaries were deﬁned as no-ﬂow boundaries (Figure 1). A third type,
Cauchy BC, was used for solute transport at both the upper and lower boundaries.
5.1. Calibration Results
Measured and simulated relative major ion concentrations are presented in Figures 3 and 4. Relative concentrations were calculated by dividing the actual concentration of a particular ion by its concentration in

RAIJ ET AL.

WATER FLOW AND MULTICOMPONENT SOLUTE TRANSPORT

6565

Water Resources Research

10.1002/2016WR018930

Table 4. Correlations Between Modeled and Measured Dataa

Drainage
R2
Slope
N
r
Soil Water Solution
R2
Slope
N
r
All Samples
R2
Slope
N
r

[Ca21][Mg21]

[Na1]

[K1]

[SO22
4 ]

[Cl2]

All Solutes

EC

0.96
1.13
22
0.99

0.92
0.99
22
0.96

0.63
1.04
22
0.81

0.70
0.89
22
0.95

0.69
0.79
22
0.86

0.98
1.07
22
1.00

0.94
0.91
132
0.92

0.94
1.20
40
0.97

0.95
1.04
13
0.98

0.95
0.68
13
0.99

0.87
1.28
13
0.93

0.90
0.92
13
0.98

0.86
0.94
13
0.95

0.98
0.98
12
0.99

0.93
0.95
77
0.94

0.97
1.03
20
0.99

0.95
1.09
35
0.98

0.83
0.84
35
0.93

0.74
1.14
35
0.86

0.83
0.89
35
0.96

0.78
0.85
35
0.90

0.97
1.02
34
0.96

0.94
0.92
209
0.93

0.94
0.87
60
0.97

a
Data for soil water and drainage solution concentrations from the rotating lysimeter system experiment. Measured EC (electricalconductivity) was compared to simulated EC calculated from cations.

the brackish water (GW in Table 1). Both ﬁgures also show the timing of the period of saline irrigation water
application and the error bars associated with the concentration measurements in the four lysimeters. Simulated and averaged measured data were compared using linear regression with an intercept equal to zero.
Numbers of observations, slopes, coefﬁcients of determination (R2), and correlation coefﬁcients (r) are
shown in Table 4. Regressions were performed for each ion, all solutes together, and the measured EC as afunction of the EC calculated from simulated cations (as in Bresler et al. [1982]).
The modeled drainage concentrations of Ca21, Mg1, and Cl2 were in good agreement with the measured
concentrations (Figure 3 and Table 4). The modeled concentrations of sodium were over and underestimated at the beginning and end of the pulse of saline water, respectively. However, the peak relative valuesof sodium were well predicted by the model. Sulfate drainage concentrations were lower than those predicted by the model at the beginning of the experiment. Potassium modeled drainage concentrations werein the range of one standard deviation from the average measured values.
2
1
The modeled soil solution concentrations of Ca21, SO22
4 , Cl , and K were in good agreement with themeasured concentrations (Figure 4 and Table 4). Measured sodium concentrations of the soil solution samples were higher than the modeled concentrations during the breakthrough of saline water. Magnesium
soil solution concentrations were underestimated by the model (slope 5 0.68). However, the overall trendof measured concentrations is described well as indicated by a high coefﬁcient of determination (R2 5 0.93)
and a slope value close to one (0.95). Measured relative concentrations of potassium were characterized byhigh variability and values higher than one for both the drainage and soil solution relative concentrations.
This is probably due to the failure to completely leach potassium before applying the saline water and apossible presence of different initial amounts of potassium in the soil of each lysimeter.

Overall salinity, evaluated using simulated cations (as in Bresler et al. [1982]), was comparable to measured
EC. The values presented in Figure 3g for the drainage solutions and Figure 4g for the soil water solutions
showed good agreement between the simulated and measured data with an R2 of 0.94 and a slope of 0.87
(60 samples, Table 4).
The effect of uncertainty and variability in Gapon coefﬁcient values on the results was studied with a sensitivity analysis in the range of the 95% conﬁdence interval for the optimized Gapon coefﬁcients (Supporting
Information). It was found that the variability caused by using different Gapon coefﬁcients was similar tothe variability of the measured data.

6. Lysimeter Design Evaluation
A hypothetical example was analyzed, in which a crop was grown in a lysimeter with the same dimensions
as in the rotating lysimeter system while three different seepage face values, situated at the same location
as in Figure 1, were considered. Simulation time was selected so that transient drainage concentrations andRAIJ ET AL.

WATER FLOW AND MULTICOMPONENT SOLUTE TRANSPORT

6566

Water Resources Research

10.1002/2016WR018930

leaching fractions reached steady state. The seepage
face BCs were considered with pressure head values of0, 230, and 2100 cm. The ﬁrst case represented a sceWater Stress Response Function Parameters
nario in which a free-drainage lysimeter has no drainP50
cm
2800
age extensions and no suction devices [e.g., Flury et al.,
P3
3
PW
cm
21.00E110
1999; Abdou and Flury, 2004; Gonçalves et al., 2006].
Salinity Stress Response Function Parameters
This lower boundary condition is the one generally
c50
cm
22417
applied to containers used in nurseries and in detachedP3
1.65
Osmotic coefﬁcient
1
growing media. In the second case, the pressure head
was the same as during the actual experiment, representing a lysimeter with highly conductive drainage
extensions. In the third case, the pressure head of 2100 cm represented a scenario in which the bottom
suction is applied using suction plates or similar devices [e.g., Evett et al., 2012; Skaggs et al., 2012]. This
range of LBCs is applicable not only to lysimeters but to crop experiments performed in pots, or to horticulture implemented in detached soilless media with limited root zone volume [e.g., Ben-Gal et al., 2009;
Shenker et al., 2011; Caron et al., 2015b].
Table 5. Root Water Uptake Parameters for the
Hypothetical Example

Irrigation (I) was applied continuously through four drippers at a rate of 1.1 L/d (representing 4.5 mm/d
according to the lysimeter area) with the solution composition deﬁned as GW in Table 1. Potential transpiration (Tp) was 3 mm/d and evaporation was neglected, so that the ratio of I/Tp was 1.5. The initial pressure
head was set to 230 cm at the bottom of the lysimeter and 280 cm at the soil surface, with a linear interpolation between these two depths. The initial soil solution concentrations were set as equal to the irrigation water quality. The initial adsorbed and precipitated concentrations were the same as the ones used in
the previous 3-D simulation.
Contrary to the experiment, in which there were no plants present, root water uptake was considered in the
hypothetical scenarios. The vertical root distribution was deﬁned according to the model of Vrugt et al.
[2001] with the following parameters: a maximum rooting depth of 30 cm, a depth of maximum intensity of
20 cm, and a shape parameter Pz of 2. S-shaped stress response functions were used for both water andosmotic stresses, the effect of which was assumed to be multiplicative [van Genuchten, 1987]. Parametersdeﬁning the S-shaped water uptake reduction functions for water and osmotic stresses are deﬁned in Table
5 as calculated in Groenveld et al. [2013] for similar conditions. Compensated root water uptake was considered by setting the critical stress index to 0.5 [Simůnek and Hopmans, 2009].
Leaching fractions were calculated daily and independently for water ﬂuxes (LFw), chloride concentrations(LFCl2), and EC (LFEC) for each simulation according to [Ayers and Westcot, 1994]
DLFw 5 ;I
LFCl2 5(9a)

½Cl2 I;½Cl 2 D

(9b)

ECI
:
ECD

(9c)

LFEC 5

Exchangable sodium percentage (ESP) was calculated as the percentage of adsorbed sodium from the total
CEC (Table 2), which was assumed to be constatnt through the simulations. The ESP spatial proﬁles werecompared at four different times (50, 100, 150, and 200 days). The results are presented for a cross sectional
area at the dripper location X 5 20.26 cm.
6.1. Effects of the Lysimeter Design on the Transient Leaching Fraction and ESP Profile Distribution
Relative drainage concentrations for various solutes are presented in Figure 5 for three different lysimeter
designs. With a pressure head of 0 cm at the seepage face there was no drainage during the ﬁrst 29 days.
As expected, the lower the pressure head at the seepage face, the earlier the concentration peak due to theleaching of dissolved solutes from adsorbed and precipitated species. The drainage concentrations for thethree different LBCs during the season show the dynamic effects of solute leaching as a result of solute

RAIJ ET AL.

WATER FLOW AND MULTICOMPONENT SOLUTE TRANSPORT

6567

Water Resources Research

10.1002/2016WR018930

610(a) Ca2+
543(b) Mg2+

9

h = 0 cm

8h = -30 cm

7

h = -100 cm

6

Irrigation

5
4

2

3
2

1

1
0

0
0

50

100

150

200

250

2

0

50

100

150

200

250

12

(d) SO42-(c) Na+
10
1.5

Ci/Co

8
1

6
4

0.5
2
0

0
0

50

100

150

200

250

0

50

100

150

200

250

10

(e) Cl-(f) K+

3
8

6

2

4
1
2

0

0
0

50

100

150

15

200

250

0

50

100

150

200

250

(g) ECEC (dS/m)

12
9
6
3
0
0

50

100

150

200

250

Time (days)
Figure 5. Simulated relative major ion concentrations and EC in the drainage water for three different lysimeter designs represented by aseepage face bottom BC with pressure heads (h) of 0, 230, and 2100 cm.

precipitation, adsorption and dissolution processes. An impact of these processes is visible in the increaseof Na1 concentrations in drainage after day 190 and the simultaneous decrease of Ca21 concentrations.
Since a constant theoretical leaching fraction of 0.33 was set according to applied irrigation and potential transpiration rates, it was expected that steady state would be reached for both water and solutes (Figure 6) [Leteyet al., 2011; Tripler et al., 2012]. Steady state was assumed when each leaching fraction reached a value of 99%of the ﬁnal corresponding leaching fraction. Initial LFCl were high due to relative lower Cl2 concentrations at
the beginning of the simulation and subsequent accumulation during the growing season until steady state
was achieved. Initial LFEC were lower than those of the LFCl due to a combination of solute leaching at the
beginning of the simulation (mainly SO22
4 from gypsum dissolution, Figure 5) and subsequent solute accumulation with time. Identical steady state values were achieved for all lysimeters for LFW and LFCl (Table 6), except
for the LFCl for the 2100 cm seepage face lysimeter. The lower the seepage face pressure head, the higher thesteady state value for LFEC (Table 6). Initial LFEC values were almost identical for all treatments. However, these

RAIJ ET AL.

WATER FLOW AND MULTICOMPONENT SOLUTE TRANSPORT

6568

Water Resources Research

10.1002/2016WR018930

1

0.8

LF (-)

0.6

0.4

0.2

0
0

50

100

150

200

250

Time (days)
W 0 cmW -30cm
W -100cm

Cl 0 cm
Cl -30cm
Cl -100cm

EC 0 cm
EC -30cm
EC -100cm

Figure 6. Simulated leaching fractions (LF) calculated from water ﬂuxes (black), chloride concentrations (blue), and EC values (red) according to equation (9).

values changed with the progress of the simulation. This was due to a combination of higher water content in
the root zone, higher pore water velocity, and its effect on cation exchange and precipitation/dissolution processes that affect calculated EC. Higher pore water velocities result in higher hydrodynamic dispersion, therefore causing LFCl and LFEC to achieve steady state earlier for lysimeters with lower seepage face pressure heads(Table 6 and Figure 6). LFW reached steady state already during the ﬁrst day of ﬂow for the simulations with the
0 and 230 cm seepage face BCs (day 29 is the ﬁrst day of outﬂow for the 0 cm lysimeter), while it took 16 days
for the simulation with the 2100 cm seepage face BC to reach steady state. This is due to the relative difference between the initial pressure head distribution along the proﬁle, which was identical for all simulations(hydrostatic equilibrium with a pressure head of 230 cm at the bottom of the soil proﬁle), and the ﬁnal steadystate distributions, which were different, reﬂecting the different lysimeter designs. The lysimeter with the seepage face BC value of 230 cm had a pressure head proﬁle distribution at steady state very similar to theimposed initial condition, reaching steady state for the LFW already at the ﬁrst day of the simulation.
There has been wide interest in the concept of the leaching fraction and whether it should be evaluatedusing transient or steady state models [Corwin et al., 2007; Letey and Feng, 2007; Dudley et al., 2008b; Leteyet al., 2011]. By considering the complex multidimensionality of ﬂow and transport patterns in lysimeters
resulting from imposed BCs (e.g., drippers and drainage outlet designs), the modeling approach presented
in this paper offers an opportunity for a better understanding of how the spatial distributions of water contents, water ﬂuxes, and dissolved, precipitated, and adsorbed solute concentrations inﬂuence measured
leaching fractions. The model showed that, in the long term, after the system reaches steady state conditions, there was no effect on LFW and LFCl. However, different LBCs did affect salt distribution in the proﬁleand the steady state LFEC values. If transient response, and not only steady state conditions, is of interest,
the effects of the different drainage outlet designs are expected to become even more relevant.
ESP values in the simulated proﬁles for different LBCs ranged between 10 and 30 and decreased with time
(Figure 7). Differences in ESP values between the simulations at the same time step were as high as 8, withthe largest variations found comparing seepage face pressure head values of 0 and 2100 cm. The largest

Table 6. Values of the Leaching Fraction at Steady State and Simulation Time When Steady State Was Reached
Days to Steady StateRAIJ ET AL.

LF Values at Steady State

Seepage Face Pressure Head

LFW

LFCl

LFEC

LFW

LFCl

LFEC

0 cm
230 cm2100 cm

29
1
16

192
169
53

229
225
198

0.34
0.34
0.34

0.34
0.34
0.33

0.33
0.34
0.35

WATER FLOW AND MULTICOMPONENT SOLUTE TRANSPORT

6569

Water Resources Research

10.1002/2016WR018930

h = 0 cmh = −30 cm

h = −100 cm
ESP %
30

Time = 50 days

28
Z
Y

Z
Y

Z
Y

24

Time = 100 days
Z
Y

26

Z
Y

Z
Y

22
20

Time = 150 days
Z
Y

18
Z
Y

Z
Y

16
14

Time = 200 days
Z
Y

12
Z
Y

Z
Y

10

Figure 7. Exchangeable sodium percentage (ESP) spatial distributions for simulations with seepage face values of 0, 230, and 2100 cm atfour times in a cross section 50 cm wide (Z) and 70 cm high (Y) at X 5 20.26 cm (see Figure 1).

discrepancies in ESP were consistently located below the root zone, close to the LBC. It is well known that
increased ESP may inﬂuence hydraulic conductivity [McNeal and Coleman, 1966; McNeal, 1968; Frenkel et al.,
1978; Shainberg et al., 1980; Crescimanno et al., 1995] depending on the soil type, water quality, and pH[Simůnek et al., 2012]. Varying ESP values within a lysimeter proﬁle could affect local hydraulic conductivities, leading to altered water ﬂow and solute transport patterns. If studies of reclamation of sodic soils or irrigation with varying water qualities are to be assessed using lysimeters, then understanding the adsorbed
and dissolved solute distributions throughout the lysimeter proﬁle becomes critical.
Different salinities and water contents due to different BCs could also affect root water uptake, depending on
the uptake and stress parameters related to particular crops. Soilless media used for food production in containers or constructed root zones using imported media in the soil itself introduce salts that may necessitateleaching before plant establishment or during the growing season. While water ﬂow in various growing media
has been studied and characterized [Raviv and Lieth, 2008; Caron et al., 2015b], there are still knowledge gaps
concerning solute transport, exchange processes and nutrient availability in detached, or otherwise imported,
growing systems designed to increase crop production [Heinen, 2001; Ityel et al., 2011; Caron et al., 2015a].
Salinity in speciﬁc areas of the root zone, produced by a combination of irrigation water quality, delivery
method, discharge rate, evaporative demand, and container or root zone medium geometry and size can limit
root growth and therefore plant development [Raviv and Lieth, 2008]. The calibrated model presented in thisstudy can be used in order to compare the solute distribution in containers or constructed root zones of different sizes and with different irrigation strategies, as well as their differences with the distributions expectedunder more traditional ﬁeld conditions. The approach of the current study could prove beneﬁcial to othercases necessitating three-dimensional consideration, for example, drip irrigation of particularly variable soils,
orchards irrigated with low-quality waters [Russo et al., 2013], as a method to evaluate solute transport asaffected by root zone proﬁle depth due to presence of a shallow water table or impermeable soil horizon.
In addition, the model can be used to study the differences between conditions produced in artiﬁcial rootzones including lysimeters or growing containers, and parallel ﬁeld conditions and to analyze design and management criteria of the constructed root zones, particularly with respect to the choice of geometry, the bottom
BC, and dripper location and discharge rate. Once the design criteria are chosen, the model can be further
applied in order to assess different methods of transient LF calculations for prediction of salinity management.

7. Conclusions
HYDRUS (2D/3D) coupled with UNSATCHEM was successfully applied to simulate and describe threedimensional transient water ﬂow and multicomponent solute transport in lysimeters. Monitored

RAIJ ET AL.

WATER FLOW AND MULTICOMPONENT SOLUTE TRANSPORT

6570

Water Resources Research

10.1002/2016WR018930

concentrations in drainage and soil pore water solutions from lysimeter experiments, in which variable
water qualities were used for irrigation, correlated well with simulated values, with R2 5 0.94 and a slope of
0.92 for all solutes measured in the drainage and soil solution samples (n 5 209). For 60 pairs of measured
and modeled EC values, R2 was 0.94 with a slope of 0.87. Models with speciﬁc chemistry, such as
UNSATCHEM coupled with HYDRUS (2D/3D), are less ﬂexible than general geochemical models [Simůnek
et al., 2014], but they are simpler to use and implement. In addition, it should be noted that long computational times are required for multicomponent multidimensional modeling, and it may be advantageous tostreamline the system by using simpler geometries and models of lower dimension (3-D/2-D/1-D).
While active and passive root solute uptake are included in the standard HYDRUS (2D/3D) [Simůnek et al.,2016], solute uptake is not considered when the multicomponent solute transport module is activated[Simůnek et al., 2012]. Considering speciﬁc solute uptake could help describe processes such as ion toxicity
and subsequent water uptake reductions. Solute uptake could also be important in quantifying moderatechanges in the soil solution composition, especially for irrigation with varying water qualities [Groenveld
et al., 2013; Silber et al., 2015]. Combining multicomponent root uptake models, such as the one presented
by Silberbush and Ben-Asher [2001], in soil-based models such as HYDRUS (2D/3D), could improve simula€der et al., 2012; Raij et al., 2013].
tions of cropping systems [SchroThe calibrated model was used to study the effects of lysimeter design on transient leaching fractions andon ESP proﬁle distribution. Different lysimeter outlet designs led to slightly different steady state leaching
fractions calculated using EC. Transient responses of the LFs between initial conditions and steady state
were a function of the LBC. Disparities in proﬁle ESP were found as a function of lysimeter design, demonstrating the potential of the calibrated model for studying how LBCs affect the proﬁle chemistry and its consequence on hydraulic conductivity.
As discussed by Simůnek et al. [2016], HYDRUS has been widely used for simulating soil water, solute andnutrient distribution for different drip irrigation with varying levels of salinity. However, the past studies
were either one-dimensional [Gonçalves et al., 2006; Ramos et al., 2011] or considered only a single solute or
general salinity as a tracer [Ajdary et al., 2007; Ramos et al., 2012; El-Nesr et al., 2013]. When using drip irrigation with water containing salts, ions are leached from the root zone to the surrounding areas, generating
favorable conditions for plant roots [Hanson et al., 2009; Phogat et al., 2012; Shan and Wang, 2012; Chenet al., 2014; Mguidiche et al., 2015]. This process can be modeled in a 2-D axisymmetrical domain as long as
the wetted areas of each dripper do not overlap [Kandelous et al., 2011]. The model presented in this study
can describe water ﬂow and solute movement in a complex three-dimensional root zone, where water contents and dissolved, adsorbed and precipitated species are highly variable. HYDRUS (2D/3D) coupled withUNSATCHEM can be used to analyze and optimize design of ﬁeld lysimeters for in situ monitoring of soluteleaching, for prevention of soil salinization, and for monitoring of salts and agrochemicals moving out of
the root zone and into deep soils and groundwater. In addition, the model can assist in designing lysimeters, commercial growing containers or constructed root zones for better management of irrigation with
low-quality waters, providing information concerning spatial distribution of general salinity, speciﬁc soluteions, and adsorbed and precipitated species.

Acknowledgments
This research was partially supportedby the Israel Ministry of Agriculture
and Rural Development (Eugene
Kandel Knowledge centers) as part of
the Root of the Matter - The root zone
knowledge center for leveraging
modern agriculture, I-CORE Program of
the Planning and Budgeting
Committee and The Israel Science
Foundation (grant 152/11), by theGoldinger Trust and the USDA ISE
Program (grant 2010-51160-21070).
The data of this paper are available
upon request (iael@post.bgu.ac.il).

RAIJ ET AL.

References
Abdou, H. M., and M. Flury (2004), Simulation of water ﬂow and solute transport in free-drainage lysimeters and ﬁeld soils with heterogeneous structures, Eur. J. Soil Sci., 55(2), 229–241, doi:10.1046/j.1365-2389.2004.00592.x.
Ajdary, K., D. K. Singh, A. K. Singh, and M. Khanna (2007), Modelling of nitrogen leaching from experimental onion ﬁeld under drip fertigation, Agric. Water Manage., 89, 15–28, doi:10.1016/j.agwat.2006.12.014.
Ali, M. H. (2010), Fundamentals of Irrigation and On-Farm Water Management, vol. 1, Springer, N. Y.
Allen, R. G., L. S. Pereira, D. Raes, and M. Smith (1998), Crop Evapotranspiration: Guidelines for Computing Crop Water Requirements, FAO Irrig.
Drain. Pap. 56, Food and Agric. Organ. of the U. N., Rome.
Ayers, R. S., and D. W. Westcot (1994), Water Quality for Agriculture, FAO Irrig. Drain. Pap. 29, Food and Agric. Organ. of the U. N., Rome.
Bear, J. (1972), Dynamics of Fluids in Porous Media, Am. Elsevier Publ. Co., N. Y.
Ben-Gal, A., and U. Shani (2002), A highly conductive drainage extension to control the lower boundary condition of lysimeters, Plant Soil,
239, 9–17.
Ben-Gal, A., H. Borochov-Neori, U. Yermiyahu, and U. Shani (2009), Is osmotic potential a more appropriate property than electrical conductivity for evaluating whole-plant response to salinity?, Environ. Exp. Bot., 65, 232–237, doi:10.1016/j.envexpbot.2008.09.006.
Bergstr€om, L. (1990), Use of lysimeters to estimate leaching of pesticides in agricultural soils, Environ. Pollut., 67(4), 325–347.
Bhantana, P., and N. Lazarovitch (2010), Evapotranspiration, crop coefﬁcient and growth of two young pomegranate (Punica granatum L.)
varieties under salt stress, Agric. Water Manage., 97(5), 715–722, doi:10.1016/j.agwat.2009.12.016.

WATER FLOW AND MULTICOMPONENT SOLUTE TRANSPORT

6571

Water Resources Research

10.1002/2016WR018930

Bresler, E., B. L. McNeal, and D. L. Carter (1982), Saline and Sodic Soils, edited by D. F. R. Bommer, et al. Springer, Berlin.
Bryla, D. R., T. J. Trout, F. Collins, and J. E. Ayars (2010), Weighing lysimeters for developing crop coefﬁcients and efﬁcient irrigation practices for vegetable crops, HortScience, 45(11), 1597–1604.
Caron, J., R. Heinse, and S. Charpentier (2015a), Organic materials used in agriculture, horticulture, reconstructed soils, and ﬁltering applications, Vadose Zone J., 14(6), 1–6, doi:10.2136/vzj2015.04.0057.
Caron, J., J. S. Price, and L. Rochefort (2015b), Physical properties of organic soil: Adapting mineral soil concepts to horticultural growing
media and histosol characterization, Vadose Zone J., 14(6), 1–14, doi:10.2136/vzj2014.10.0146.
Chen, L. J., Q. Feng, F. R. Li, and C. S. Li (2014), A bidirectional model for simulating soil water ﬂow and salt transport under mulched drip
irrigation with saline water, Agric. Water Manage., 146, 24–33, doi:10.1016/j.agwat.2014.07.021.
Chen, W., Z. Hou, L. Wu, Y. Liang, and C. Wei (2010), Evaluating salinity distribution in soil irrigated with saline water in arid regions ofnorthwest China, Agric. Water Manage., 97, 2001–2008, doi:10.1016/j.agwat.2010.03.008.
Corwin, D. L. (2002), Miscible solute transport: Suction cups, in Methods of Soil Analysis, Part 4: Physical Methods, edited by J. H. Dane andG. C. Topp, pp. 1261–1266, Soil Sci. Soc. of Am., Madison, Wis.
Corwin, D. L., J. D. Rhoades, and J. 
Simůnek (2007), Leaching requirement for soil salinity control: Steady-state versus transient models,
Agric. Water Manage., 90(3), 165–180, doi:10.1016/j.agwat.2007.02.007.
Crescimanno, G., M. Iovino, and G. Provenzano (1995), Inﬂuence of salinity and sodicity on soil structural and hydraulic characteristics, Soil
Sci. Soc. Am. J., 59, 1701–1708, doi:10.2136/sssaj1995.03615995005900060028x.
Dudley, L. M., A. Ben-Gal, and N. Lazarovitch (2008a), Drainage water reuse: Biological, physical, and technological considerations for system management, J. Environ. Qual., 37, S25–S35, doi:10.2134/jeq2007.0314.
Dudley, L. M., A. Ben-Gal, and U. Shani (2008b), Inﬂuence of plant, soil, and water on the leaching fraction, Vadose Zone J., 7, 420–425, doi:
10.2136/vzj2007.0103.
El-Nesr, M. N., A. A. Alazba, and J. Simůnek (2013), HYDRUS simulations of the effects of dual-drip subsurface irrigation and a physical barrier on water movement and solute transport in soils, Irrig. Sci., 32(2), 111–125, doi:10.1007/s00271-013-0417-x.
Evett, S. R., R. C. Schwartz, T. A. Howell, R. L. Baumhardt, and K. S. Copeland (2012), Can weighing lysimeter ET represent surrounding ﬁeld
ET well enough to test ﬂux station measurements of daily and sub-daily ET?, Adv. Water Resour., 50, 79–90, doi:10.1016/
j.advwatres.2012.07.023.
Flury, M., M. V. Yates, and W. A. Jury (1999), Numerical analysis of the effect of the lower boundary condition on solute transport in lysimeters, Soil Sci. Soc. Am. J., 63, 1493–1499.
Frenkel, H., J. O. Goertzen, and J. D. Rhoades (1978), Effects of clay type and content, exchangeable sodium percentage, and electrolyte
concentration on clay dispersion and soil hydraulic conductivity, Soil Sci. Soc. Am. J., 42, 32–39, doi:10.2136/
sssaj1978.03615995004200010008x.
Gasser, M. O., J. Caron, M. R. Laverdière, and R. Lagac
e (2002), Solute transport modeling under cultivated sandy soils and transient water
regime, J. Environ. Qual., 31, 1722–1730.
Gonçalves, M. C., J. Simůnek, T. B. Ramos, J. C. Martins, M. J. Neves, and F. P. Pires (2006), Multicomponent solute transport in soil lysimetersirrigated with waters of different quality, Water Resour. Res., 42, W08401, doi:10.1029/2005WR004802.
Groenveld, T., A. Ben-Gal, U. Yermiyahu, and N. Lazarovitch (2013), Weather determined relative sensitivity of plants to salinity: Quantiﬁcation and simulation, Vadose Zone J., 12(4), 9, doi:10.2136/vzj2012.0180.
Groh, J., J. Vanderborght, T. P€
utz, and H. Vereecken (2016), How to control the lysimeter bottom boundary to investigate the effect of climate change on soil processes?, Vadose Zone J., 15(7), 1–15, doi:10.2136/vzj2015.08.0113.
Hannes, M., U. Wollschl€ager, F. Schrader, W. Durner, S. Gebler, T. P€
utz, J. Fank, G. von Unold, and H. J. Vogel (2015), High-resolution
estimation of the water balance of high-precision lysimeters, Hydrol. Earth Syst. Sci. Discuss., 12, 569–608, doi:10.5194/hessd-12-5692015.
Hanson, B. R., J. 
Simůnek, and J. W. Hopmans (2006), Evaluation of urea-ammonium-nitrate fertigation with drip irrigation using numericalmodeling, Agric. Water Manage., 86, 102–113, doi:10.1016/j.agwat.2006.06.013.
Hanson, B. R., J. W. Hopmans, and J. Simůnek (2008), Leaching with subsurface drip irrigation under saline, shallow groundwater conditions, Vadose Zone J., 7(2), 810–818, doi:10.2136/vzj2007.0053.
Hanson, B. R., J. W. Hopmans, and D. E. May (2009), Drip irrigation provides the salinity control needed for proﬁtable irrigation of tomatoes
in the San Joaquin Valley, Calif. Agric., 63(3), 131–136.
Hartz, T. K., and G. J. Hochmuth (1996), Fertility management of drip-irrigated vegetables, Hort Technol., 6(3), 168–172.
Heinen, M. (2001), FUSSIM2: Brief description of the simulation model and application to fertigation scenarios, Agronomie, 21, 285–296,
doi:10.1051/agro:2001124.
Ityel, E., N. Lazarovitch, M. Silberbush, and A. Ben-Gal (2011), An artiﬁcial capillary barrier to improve root zone conditions for horticultural
crops: Physical effects on water content, Irrig. Sci., 29(2), 171–180, doi:10.1007/s00271-010-0227-3.
Jacobsen, O. H., F. J. Leij, and M. T. van Genuchten (1992), Lysimeter study of anion transport during steady ﬂow through layered coarsetextured soil proﬁles, Soil Sci., 154(3), 196–205.
Jacques, D., C. Smith, J. 
Simůnek, and D. Smiles (2012), Inverse optimization of hydraulic, solute transport, and cation exchange parameters
using HP1 and UCODE to simulate cation exchange, J. Contam. Hydrol., 142–143, 109–125, doi:10.1016/j.jconhyd.2012.03.008.
Kaledhonkar, M. J., and A. K. Keshari (2006), Modelling the effects of saline water use in agriculture, Irrig. Drain., 55, 177–190.
Kandelous, M. M., J. 
Simůnek, M. T. van Genuchten, and K. Malek (2011), Soil water content distributions between two emitters of a subsurface drip irrigation system, Soil Sci. Soc. Am. J., 75(2), 488–497, doi:10.2136/sssaj2010.0181.
Lazarovitch, N., A. Ben-Gal, and U. Shani (2006), An automated rotating lysimeter system for greenhouse evapotranspiration studies,
Vadose Zone J., 5, 801–804, doi:10.2136/vzj2005.0137.
Letey, J., and G. L. Feng (2007), Dynamic versus steady-state approaches to evaluate irrigation management of saline waters, Agric. Water
Manage., 91, 1–10, doi:10.1016/j.agwat.2007.02.014.
Letey, J., G. J. Hoffman, J. W. Hopmans, S. R. Grattan, D. Suarez, D. L. Corwin, J. D. Oster, L. Wu, and C. Amrhein (2011), Evaluation of soil
salinity leaching requirement guidelines, Agric. Water Manage., 98(4), 502–506, doi:10.1016/j.agwat.2010.08.009.
Marek, T., G. Piccinni, A. Schneider, T. Howell, M. Jett, and D. Dusek (2006), Weighing lysimeters for the determination of crop waterrequirements and crop coefﬁcients, Appl. Eng. Agric., 22(6), 851–856.
McNeal, B. L. (1968), Prediction of the effect of mixed-salt solutions on soil hydraulic conductivity, Soil Sci. Soc. Am. J., 32, 190–193, doi:
10.2136/sssaj1968.03615995003200020013x.
McNeal, B. L., and N. T. Coleman (1966), Effect of solution composition on soil hydraulic conductivity, Soil Sci. Soc. Am. J., 30, 308–312, doi:
10.2136/sssaj1966.03615995003000030007x.

RAIJ ET AL.

WATER FLOW AND MULTICOMPONENT SOLUTE TRANSPORT

6572

Water Resources Research

10.1002/2016WR018930

MDH Engineered Solutions Corp. (2003), Evaluation of Computer Models for Predicting the Fate and Transport of Salt in Soil and Groundwater,
Edmonton, Alberta, Canada.
Mguidiche, A., G. Provenzano, B. Douh, S. Khila, G. Rallo, and A. Boujelben (2015), Assesing HYDRUS-2D to simulate soil water content (SWC)and salt accumulation under an SDI system: Application to a potato crop in a semi-arid area of central Tunisia, Irrig. Drain., 64, 263–274.
Neuman, S. P. (1990), Universal scaling of hydraulic conductivities and dispersivities in geologic media, Water Resour. Res., 26(8), 1749–
1758, doi:10.1029/WR026i008p01749.
NOAA-ESRL (2015), Trends in atmospheric carbon dioxide, Earth Syst. Res. Lab., Global Monit. Div., Boulder, Colo. [Available at http://www.
esrl.noaa.gov/gmd/ccgg/trends/global.html#global_data.]
Oren, O., Y. Yechieli, J. K. B€
ohlke, and A. Dody (2004), Contamination of groundwater under cultivated ﬁelds in an arid environment, central
Arava Valley, Israel, J. Hydrol., 290, 312–328, doi:10.1016/j.jhydrol.2003.12.016.
Phogat, V., M. Mahadevan, M. Skewes, and J. W. Cox (2012), Modelling soil water and salt dynamics under pulsed and continuous surface
drip irrigation of almond and implications of system design, Irrig. Sci., 30, 315–333, doi:10.1007/s00271-011-0284-2.
Poeter, E. P., M. C. Hill, E. R. Banta, S. Mehl, and S. Christensen (2005), UCODE_2005 and six other computer codes for universal sensitivity
analysis, calibration and uncertainty evaluation, U.S. Geol. Surv. Tech. Methods, 6-A11, 283 pp.
Poeter, E. P., M. C. Hill, D. Lu, C. R. Tiedeman, and S. Mehl (2014), UCODE_2014, with new capabilities to deﬁne parameters unique to predictions, calculate weights using simulated values, estimate parameters with SVD, evaluate uncertainty with MCMC, and more, Rep.
GWMI 2014-02, Integ. Groundwater Model. Cent., Golden, Colo.
Raij, I., Lazarovitch, N., Ben-Gal, A., Yermiyhau, U., and Jacques, D. (2013), Accounting for solution composition in a plant roots active nutrient uptake model, in Proceedings of the 4th International Conference on HYDRUS Software Applications to Subsurface Flow and Contaminant Transport Problems, edited by J. Simůnek, M. Th. van Genuchten, and R. Kodesova, pp. 299–306, Dep. of Soil Sci. and Geol., Czech
Univ. of Life Sci., Prague, Czech Republic.
Ramos, T. B., J. Simůnek, M. C. Gonçalves, J. C. Martins, A. Prazeres, N. L. Castanheira, and L. S. Pereira (2011), Field evaluation of a multicomponent solute transport model in soils irrigated with saline waters, J. Hydrol., 407(1–4), 129–144, doi:10.1016/j.jhydrol.2011.07.016.
Ramos, T. B., J. 
Simůnek, M. C. Gonçalves, J. C. Martins, A. Prazeres, and L. S. Pereira (2012), Two-dimensional modeling of water and nitrogen fate from sweet sorghum irrigated with fresh and blended saline waters, Agric. Water Manage., 111, 87–104, doi:10.1016/
j.agwat.2012.05.007.
Rasouli, F., A. K. Pouya, and J. Simůnek (2012), Modeling the effects of saline water use in wheat-cultivated lands using the UNSATCHEM
model, Irrig. Sci., 31, 1009–1024, doi:10.1007/s00271-012-0383-8.
Raviv, M., and H. J. Lieth (2008), Soilless Culture: Theory and Practice, Elsevier B.V., Amsterdam, Netherlands.
Rawlins, S. L. (1973), Principles of managing high frequency irrigation, Soil Sci. Soc. Am. Proc., 37, 626–629.
Rhoades, J. D., A. Kandiah, and A. M. Mashali (1992), The Use of Saline Waters for Crop Production, FAO Irrig. Drain. Pap. 48, Food and Agric.
Organ. of the U. N., Rome.
Rhoades, J. D., F. Chanduvi, and S. Lesch (1999), Soil Salinity Assessment: Methods and Interpretation of Electrical Conductivity Measurements,
FAO Irrig. Drain. Pap. 57, Food and Agric. Organ. of the U. N., Rome.
Roberts, T., N. Lazarovitch, A. W. Warrick, and T. L. Thompson (2009), Modeling salt accumulation with subsurface drip irrigation using HYDRUS-2D, Soil Sci. Soc. Am. J., 73(1), 233–240, doi:10.2136/sssaj2008.0033.
Ruiz-Pe~nalver, L., J. A. Vera-Repullo, M. Jim
enez-Buendıa, I. Guzman, and J. M. Molina-Martınez (2015), Development of an innovative low
cost weighing lysimeter for potted plants: Application in lysimetric stations, Agric. Water Manage., 151, 103–113, doi:10.1016/
j.agwat.2014.09.020.
Russo, D., A. Laufer, R. H. Shapira, and D. Kurtzman (2013), Assessment of solute ﬂuxes beneath an orchard irrigated with treated sewage
water: A numerical study, Water Resour. Res., 49, 657–674, doi:10.1002/wrcr.20085.
Schoen, R., J. P. Gaudet, and D. E. Elrick (1999), Modelling of solute transport in a large undisturbed lysimeter, during steady-state waterﬂux, J. Hydrol., 215(1–4), 82–93, doi:10.1016/S0022-1694(98)00263-7.
Schr€oder, N., M. Javaux, J. Vanderborght, B. Steffen, and H. Vereecken (2012), Effect of root water and solute uptake on apparent soil dispersivity: A simulation study, Vadose Zone J., 11(3), 1–14, doi:10.2136/vzj2012.0009.
Shainberg, I., J. D. Rhoades, and R. J. Prather (1980), Effect of low electrolyte concentration on clay dispersion and hydraulic conductivity ofsodic soils, Soil Sci. Soc. Am. J., 45, 273–277, doi:10.2136/sssaj1981.03615995004500020012x.
Shan, Y., and Q. Wang (2012), Simulation of salinity distribution in the overlap zone with double-point-source drip irrigation using HYDRUS-3D, Aust. J. Crop Sci., 6(2), 238–247.
Shenker, M., D. Harush, J. Ben-Ari, and B. Chefetz (2011), Uptake of carbamazepine by cucumber plants—A case study related to irrigation
with reclaimed wastewater, Chemosphere, 82, 905–910, doi:10.1016/j.chemosphere.2010.10.052.
Silber, A., Y. Israeli, I. Elingold, M. Levi, I. Levkovitch, D. Russo, and S. Assouline (2015), Irrigation with desalinated water: A step toward
increasing water saving and crop yields, Water Resour. Res., 51(1), 450–464, doi:10.1002/2014WR016398.
Silberbush, M., and J. Ben-Asher (2001), Simulation study of nutrient uptake by plants from soilless cultures as affected by salinity buildup
and transpiration, Plant Soil, 233, 59–69.
Simůnek, J., and J. W. Hopmans (2009), Modeling compensated root water and nutrient uptake, Ecol. Modell., 220(4), 505–521, doi:10.1016/
j.ecolmodel.2008.11.004.
Simůnek, J., and M. T. van Genuchten (2006), Contaminant transport in the unsaturated zone theory and modeling, in The Handbook of
Groundwater Engineering, edited by J. W. Delleur, pp. 22-1–22-46, CRC Press, Boca Raton.
Simůnek, J., M. T. van Genuchten, and M. Sejna (2008), Development and applications of the HYDRUS and STANMOD software packages
and related codes, Vadose Zone J., 7(2), 587–600, doi:10.2136/vzj2007.0077.
Simůnek, J., M. 
Sejna, and M. Th. van Genuchten (2012), The UNSATCHEM module for HYDRUS (2D/3D) simulating two-dimensional movement of and reactions between major ions in soils, version 1.0, 54 pp., PC Progress, Prague, Czech Republic.
Simůnek, J., D. Jacques, T. B. Ramos, and B. Leterme (2014), The use of multicomponent solute transport models in environmental analyses,
Chapter 16, in Application of Soil Physics in Environmental Analyses: Measuring, Modelling and Data Integration, edited by W. Teixeira,et al., pp. 377–402, Springer, Heidelberg, Germany, doi:10.1007/978-3-319-06013-2.
Simůnek, J., M. T. van Genuchten, and M. Sejna (2016), Recent developments and applications of the HYDRUS computer software packages, Vadose Zone J., 15(7), 25, doi:10.2136/vzj2016.04.0033.
Skaggs, T. H., and F. J. Leij (2002), Solute transport: Theoretical background, in Methods of Soil Analysis, Part 4: Physical Methods, edited byJ. H. Dane and G. C. Topp, pp. 1353–1378, Soil Sci. Soc. of Am., Madison, Wis.
Skaggs, T. H., T. J. Trout, J. Simůnek, and P. J. Shouse (2004), Comparison of HYDRUS-2D simulations of drip irrigation with experimental
observations, J. Irrig. Drain. Eng., 130(4), 304–310, doi:10.1061/(ASCE)0733-9437(2004)130:4(304)CE.

RAIJ ET AL.

WATER FLOW AND MULTICOMPONENT SOLUTE TRANSPORT

6573

Water Resources Research

10.1002/2016WR018930

Skaggs, T. H., D. L. Suarez, S. Goldberg, and P. J. Shouse (2012), Replicated lysimeter measurements of tracer transport in clayey soils: Effectsof irrigation water salinity, Agric. Water Manage., 110, 84–93, doi:10.1016/j.agwat.2012.04.003.
Skaggs, T. H., D. L. Suarez, and D. L. Corwin (2014), Global sensitivity analysis for UNSATCHEM simulations of crop production with degraded waters, Vadose Zone J., 13(6), 1–13, doi:10.2136/vzj2013.09.0171.
Smiles, D. E., and C. J. Smith (2004), Absorption of artiﬁcial piggery efﬂuent by soil: A laboratory study, Aust. J. Soil Res., 42, 961–975, doi:
10.1071/SR04008.
Suarez, D. L., and J. 
Simůnek (1997), UNSATCHEM: Unsaturated water and solute transport model with equilibrium and kinetic chemistry,
Soil Sci. Soc. Am., 61, 1633–1646.
Tripler, E., U. Shani, A. Ben-Gal, and Y. Mualem (2012), Apparent steady state conditions in high resolution weighing-drainage lysimeters
containing date palms grown under different salinities, Agric. Water Manage., 107, 66–73, doi:10.1016/j.agwat.2012.01.010.
Turkeltaub, T., D. Kurtzman, E. E. Russak, and O. Dahan (2015), Impact of switching crop type on water and solute ﬂuxes in deep vadosezone, Water Resour. Res., 51, 9828–9842, doi:10.1016/0022-1694(68)90080-2.
Vanclooster, D., M. Mallants, J. Vanderborght, J. Diels, J. Van Orshoven, and J. Feyen (1995), Monitoring solute transport in a multi-layered
sandy lysimeter using time domain reﬂectometry, Soil Sci. Soc. Am. J., 59, 337–344.
Vanderborght, J., and H. Vereecken (2007), Review of dispersivities for transport modeling in soils, Vadose Zone J., 6, 29–52, doi:10.2136/
vzj2006.0096.
van Genuchten, M. T. (1980), A closed-form equation for predicting the hydraulic conductivity of unsaturated soils, Soil Sci. Soc. Am., 44(5),892–898.
van Genuchten, M. T. (1987), A numerical model for water and solute movement in and below the root zone, Res. Rep. 121, Riverside, California, USDA, ARS.
van Genuchten, M. T., F. J. Leij, and S. R. Yates (1991), The RETC code for quantifying the hydraulic functions of unsaturated soils, EPA/600/
2-91/065, pp. 1–85, U.S. Department of Agriculture, Calif., doi:10.1002/9781118616871.
Vrugt, J. A., M. T. van Wijk, J. W. Hopmans, and J. 
Simůnek (2001), One-, two-, and three-dimensional root water uptake functions for transient modeling, Water Resour. Res., 37(10), 2457–2470, doi:10.1029/2000WR000027.
Wang, Q., K. Cameron, G. Buchan, L. Zhao, E. Zhang, N. Smith, and S. Carrick (2012), Comparison of lysimeters and porous ceramic cups formeasuring nitrate leaching in different soil types, N. Z. J. Agric. Res., 55(4), 333–345, doi:10.1080/00288233.2012.706224.
Warrick, A. W., and N. Lazarovitch (2007), Inﬁltration from a strip source, Water Resour. Res., 43(3), 1–5, doi:10.1029/2006WR004975.
Weiherm€
uller, L., J. Siemens, M. Deurer, S. Knoblauch, H. Rupp, a G€
ottlein, and T. P€
utz (2007), In situ soil water extraction: A review, J. Environ. Qual., 36(6), 1735–1748, doi:10.2134/jeq2007.0218.
Winton, K., and J. B. Weber (1996), A review of ﬁeld lysimeter studies to describe the environmental fate of pesticides, Weed Technol., 10(1),202–209.

RAIJ ET AL.

WATER FLOW AND MULTICOMPONENT SOLUTE TRANSPORT

6574