Journal of Hydrology 290 (2004) 63–79
www.elsevier.com/locate/jhydrol

Two-dimensional simulation of water flow and solute transport
below furrows: model calibration and validation
Fariborz Abbasia, Jan Feyena,*, M.Th. van Genuchtenb,1aInstitute for Land and Water Management, Katholieke Universiteit Leuven, Vital Decosterstraat 102, 3000 Leuven, BelgiumbGeorge E. Brown Jr., Salinity Laboratory, USDA-ARS, 450 W. Big Springs Road, Riverside, CA 92507-4617, USA
Received 7 October 2002; revised 7 November 2003; accepted 21 November 2003

Abstract
In this study a two-dimensional numerical flow/transport model (HYDRUS-2D) was calibrated and experimentally validated
using data from long furrow irrigation experiments. The model was calibrated using data from an experiment carried out
assuming free-draining (FD) outlet conditions, and subsequently validated against data from three experiments assuming
blocked-end conditions. The data were analyzed using the Richards’ equation for variably saturated flow and either the
traditional convection– dispersion equation (CDE) or the physical nonequilibrium mobile – immobile (MIM) model for solutetransport. Optimization was accomplished by means of Levenberg – Marquardt optimization in combination with the
HYDRUS-2D model. Simultaneous and two-step optimization approaches were used to estimate the soil hydraulic and solutetransport parameters near the FD furrow inlet and outlet sites. First, the saturated hydraulic conductivity ðKs Þ and CDE or MIM
solute transport parameters were estimated simultaneously. We also used sequential (two-step) estimation in which we firstestimated the soil hydraulic parameters followed by estimation of the solute transport parameters. In the two-step method, thesaturated soil water content ðus Þ; the n parameter in van Genuchten’s soil water retention model, and Ks values were estimatedduring the first step, and the CDE or MIM solute transport parameters during the second step.
Estimated soil hydraulic and solute transport parameters were found to vary substantially between the inlet and outlet sites.
Estimated CDE and MIM transport parameters were very similar for both optimization approaches. The two-step method
significantly improved predictions of the soil water content during model calibration, while the solute concentration predictions
were nearly the same for both approaches, with both not providing a good description of the observed concentrations. Solutedata were also analyzed using horizontal averages to somewhat lessen the effects of spatial variability. Horizontally averaged
concentration distributions showed better agreement with the predictions. Soil water contents for the three blocked-end
experiments during model validation were well predicted. The two-step method produced slightly better agreement withobserved data. However, both optimization approaches produced relatively poor agreement between measured and predicted
solute concentrations and deep percolation rates.
q 2003 Elsevier Ltd. All rights reserved.
Keywords: Water flow; Solute transport; Furrow irrigation; Parameter estimation; Simultaneous and sequential parameter estimation* Corresponding author. Tel.: þ32-16-329756; fax: þ 32-16329760.
E-mail addresses: jan.feyen@agr.kuleuven.ac.be (J. Feyen),
frabbasi@hotmail.com (F. Abbasi), rvang@ussl.ars.usda.gov (M.T.
van Genuchten).
1Tel.: þ1-909-369-4847.
0022-1694/$ - see front matter q 2003 Elsevier Ltd. All rights reserved.
doi:10.1016/j.jhydrol.2003.11.028

1. Introduction
Numerical models are increasingly being used forpredicting or analyzing water flow and chemical

64

F. Abbasi et al. / Journal of Hydrology 290 (2004) 63–79

transport in soils and groundwater. The governing
flow and transport equations are now reasonably wellestablished, while a large number of approximate
analytical or more comprehensive numerical solutions
have been published during especially the last three
decades. Still, a number of numerical and conceptual
difficulties remain to be solved, especially for largescale transient, multi-dimensional field applications.
For example, hydrologists and engineers lack precise
methodologies for identifying values of the parameters in the governing flow and transport equations.
Identification and accurate simulation of field-scale
processes subject to natural boundary conditions
remains a complex task since many input parameters
are difficult to measure at the desired scale. Recent
studies (Jacques et al., 2002) indicate that a combination of numerical solutions of the governing flow
and transport equations with inverse optimization
algorithms and detailed measurement of different
variables is a promising approach for process andparameter identification. Simultaneous estimation of
the soil hydraulic and solute transport parametersfrom measured flow and transport data appears to be
the most beneficial since the method uses cross-overeffects between state variables and parameters (Sunand Yeh, 1990), and hence takes advantage of allavailable information. Simultaneous estimation of soil
hydraulic and solute transport parameters often also
yields smaller estimation errors than sequential
estimation (Mishra and Parker, 1989; Šimůnek et al.,2002). Carrera and Neuman (1986) for groundwaterstudies similarly found that using steady-state andtransient data jointly may lead to less uncertainty inmany (but not all) cases, as opposed to using either
type of data independently.
Notwithstanding the fact that subsurface waterflow and solute transport processes have been actively
studied for several decades, relatively few studies
exists in which inverse methods have been used toestimate solute transport parameters using data frompractical field experiments, especially when
conducted during variably saturated flow. Main
reasons for this are a number of practical difficulties
often encountered during data collection (e.g. extreme
variability of the subsurface, the time-consuming
nature of many field experiments, available resources,
and required labor), as well as numericaland conceptual difficulties inherent to accurately

simulating field-scale solute transport. For example,
one of the problems of simultaneous estimation offlow and transport parameters stems from the fact thatsolutions may become unstable for certain combinations of the estimated parameters (Carrera, 1987).
The main objective of this study was to estimate
soil hydraulic and solute transport properties from apractical two-dimensional (2D) transient furrow
irrigation experiment using simultaneous and sequential parameter optimization approaches. Optimizations were carried out using the Levenberg –Marquardt optimization in combination with the
HYDRUS-2D variably saturated flow/transport software package of Šimůnek et al. (1999). The optimized
parameters were subsequently used to predict waterflow and solute transport in three other furrow
irrigation experiments conducted at the same scalein the same field.

2. Furrow irrigation experiments
Details about the experiments can be found in
Abbasi et al. (2003a,d). Only those parts of theexperiments related to this study are briefly summarized here. Four large-scale furrow irrigation experiments (Fig. 1), one under free-draining (FD) andthree under blocked-end conditions, were conductedat the Maricopa Agricultural Center, Phoenix, AZ, in
2001 on a bare sandy loam soil. The experiments werecarried out on 115-m long furrows, spaced 1 m apart.
Each experiment included three furrows, one monitored non-wheel furrow in the middle, and two guard
wheel furrows, one at each side. The FD experiment
was run for two successive irrigation events 10 days
apart. The first irrigation lasted 275 min and thesecond 140 min. The average inflow rates were 1.07
and 1.03 l s21 for the first and second irrigations,
respectively. Bromide in the form of CaBr2 wasinjected at a constant rate (6.3 g Br l21) during theentire first irrigation. The second irrigation was
carried out with unamended water (without bromide)
to investigate the distribution of moisture and Br in
the soil profile below and adjacent to the furrow.
Two sets of neutron probe access tubes were
installed at x ¼ 5 and 110 m along the monitored
furrow. Hereafter we refer to these locations as theinlet and outlet sites, respectively. Each set included

F. Abbasi et al. / Journal of Hydrology 290 (2004) 63–79

65

Fig. 1. Plane view of the furrow irrigation field experiments at Maricopa Agricultural Center, Phoenix, AZ (axes are not to scale).

five 2.2-m long neutron probe access tubes, installedin two rows at different locations in a cross-section
perpendicular to the furrow, spaced 50 cm apart toavoid mutual interference of the readings (Fig. 2). A
site-calibrated neutron probe was used for repeatedmeasurement of the soil water content at depths of 20,
40, 60, 80, 100, 140 and 180 cm. In addition to initial
readings before each irrigation event, water contents
were recorded 6 and 12 h after each irrigation and thendaily. Measurements were taken up to 15 days after
the second irrigation. Water contents of the surface
layer (0 – 30 cm) were also measured with Time
Domain Reflectometry at the same times as neutronprobe readings were taken.
Three experiments were conducted on blocked-end
furrows using the following three solute application

options (referred to as the 100%, FH and SH
experiments, respectively):bromide applied during the entire irrigation(100%),
bromide applied during only the first half of theirrigation (FH),bromide applied during the second half of theirrigation (SH).
All three experiments involved one 140-min long
irrigation. The average inflow rates for the three runs
were 1.29, 1.32, and 1.28 l s21, respectively. Averageapplied bromide concentrations for the 100%, FH, andSH experiments were 2.36, 2.79, and 5.35 g l21,
respectively. We tried to apply the same amount ofFig. 2. Position of neutron probe access tubes at different locations in the furrow cross-section. Numbers relate to access tubes installed in two
different rows; the first row includes tubes 2 and 4 along the sides and the second row includes tubes 1, 3, and 5.

66

F. Abbasi et al. / Journal of Hydrology 290 (2004) 63–79

solute (not the same concentration) during eachexperiment. However, because of several unforeseen
practical problems, this objective was not fully met.
Below we summarize the remaining measurements
for both the FD and blocked-end experiments.
Soil samples for analyzing bromide concentrations
and gravimetric soil water contents were collected
five days after each irrigation at three different
locations: at 10, 60 and 100 m from the inlet for theblocked-end experiments, and at 5, 60 and 110 m for
the FD experiment (the inlet, middle, and outlet sites,
respectively). The samples were taken from one sideof the monitored furrows at three locations (top, side,and bottom of the furrows; e.g. at locations 1, 2 and 3
in Fig. 2) in a cross-section perpendicular to thefurrow axis at similar depths as used for the neutronprobe measurements. Soil samples from the FD
experiment were also taken five days after the second
irrigation. Soil water extractions (1:1 weight:volume
ratios) were taken and analyzed for bromide with a
Lachat QuikChem flow injection analyzer using
standard colorimetric procedures. Overland water
samples (for analysis of bromide concentrations)and water flow depths in the furrows were also
taken at the inlet, middle, and outlet sites every few
minutes as soon as water reached a particular station.
Geometries of the furrows were measured with aprofilometer before and after each irrigation at the
inlet, middle, and outlet sites. These measured
geometries served as the upper boundaries for our
numerical calculations using HYDRUS-2D. We also
took soil bulk densities at several locations using a
Madera soil sampler up to a depth of 180 cm. In
addition, 38 undisturbed soil samples (6-cm long, 5.4cm diameter soil cores) were collected at selected
locations and depths (up to 100 cm) for laboratory
analyses of the soil water retention curve. The first FD
irrigation was carried out on March 23, 2001, and thesecond irrigation 10 days later, while the 100% andFH experiments were carried out on the same day(April 3, 2001) and the SH experiment one day later.
An 18-mm 3-h long rainfall event occurred on April 7.

the physical nonequilibrium mobile – immobile(MIM) model for solute transport (van Genuchten
and Wagenet, 1989; Clothier et al., 1998). Detailed
descriptions of the governing equations for water flow
and solute transport and the invoked parameter
optimization procedure are given by Abbasi et al.
(2003c).
2.2. Initial and boundary conditions
Measured bromide concentrations and soil watercontents before the experiments were used as initialconditions within the flow domain. Time-spacedependent flow depths (surface ponding, hðx; tÞ in
Fig. 3) and overland bromide concentrations were
specified as the upper boundary condition in the
furrow during irrigation, while averages of measured
pan evaporation rates from the nearest weather station(about 150 m away from the experimental field) andestimated reference evapotranspiration rates obtainedusing the Penman – Monteith method (as reported by
Abbasi et al., 2003a) were used as atmosphericboundary conditions during redistribution phase. As
indicated earlier, in each experiment flow depths at
the inlet, middle, and outlet sites were frequently

2.1. Governing equations
The data were analyzed in terms of the Richards’
equation for variably saturated flow and either the
traditional convection– dispersion equation (CDE) or

Fig. 3. Water and solute boundary conditions used for the
HYDRUS-2D calculations.

F. Abbasi et al. / Journal of Hydrology 290 (2004) 63–79

67

measured. A Cauchy condition was used for the upperboundary condition for solute transport, while freedrainage conditions for both water and solute wereapplied to the lower boundary of the domain (Fig. 3).
No-flux boundary conditions were applied to bothsides of the flow domain. Measured furrow crosssections were used to define the upper section of theflow domain for each experimental site. A relatively
fine grid (, 0.5 cm) was used near the soil surface,and a much coarser grid (, 4 cm) at and near thebottom of the domain.

(model validation). In this study we assumed
equivalent homogeneous soil profiles because of
insufficient field data to distinguish between the
different soil layers. A previous, more detailed study
on short furrows (3-m long) at the same field showedthat considering soil heterogeneity did not
significantly improve the model predictions (Abbasi
et al., 2003b).

2.3. Parameter estimation

3.1. Estimation of soil water retention properties

Soil hydraulic and CDE and MIM solute transportparameters were inversely estimated using theLevenberg – Marquardt optimization (Marquardt,
1963; Šimůnek and Hopmans, 2002) in combination
with the HYDRUS-2D numerical code (Šimunek
et al., 1999) using two optimization approaches. First,
we estimated the saturated hydraulic conductivity ðKs Þand the CDE or MIM solute transport parameterssimultaneously. Of the soil hydraulic parameters only
Ks was included in the optimization since we did nothave any direct measurements of Ks : Also, it isgenerally not recommended to estimate too many
parameters simultaneously because of possible correlation among the parameters and resulting uniquenessproblems (Hopmans et al., 2002). To further improvethe model predictions, a two-step optimization
approach was also carried out in which we firstestimated the saturated soil water content ðus Þ; theparameter n in van Genuchten’s soil hydraulicproperty model, and Ks as the most sensitive soil
hydraulic parameters (Abbasi et al., 2003b), followedby estimation of transport parameters. Soil watercontent readings and solute concentrations were usedin the objective function during the optimization
process.
As indicated earlier, two sets of neutron probe
access tubes at the inlet and outlet sites of themonitored furrow were available during the FD
experiments for detailed water content measurements.
The FD experiment hence produced relatively more
information for the purpose of inversely estimating
the soil hydraulic and solute transport parameters(model calibration). FD parameter estimates were
subsequently used to simulate the other experiments

The soil water retention curve, uðhÞ; was describedusing the closed-form equation of van Genuchten(1980), and the unsaturated soil hydraulic conductivity
function, KðhÞ; using the statistical pore-size distribution model of Mualem (1976). Laboratory-measured
soil water retention curves from 38 soil samples werescaled using linear scaling relationships proposed by
Vogel et al. (1991) to reduce and quantify spatial
variability in the soil hydraulic properties. In this study
we fixed the saturated ðups ¼ 0:411 cm3 cm23 Þ andresidual ðupr ¼ 0:106 cm3 cm23 Þ soil water contents at
the averages of the fitted local us and ur values, while lin Mualem’s equation was assumed to be 0.5 (Mualem,1976). Ks values were inversely estimated from themeasured soil water contents since no hydraulicconductivity measurements were made.
Measured (un-scaled), scaled, and reference soil
water retention curves are shown in Fig. 4. The effectof scaling is clearly noticeable on the retention curves,
and reflects the considerable spatial variability in thesoil hydraulic properties at the field site. Scaled a andn values in van Genuchten’s (1980) equation were0.039 cm21 and 1.55, respectively. The soil hydraulicparameters thus obtained (including us and ur ) weresimilar (except possibly a) to those estimated usingneural network-based pedotransfer functions (theROSETTA code) derived by Schaap et al. (2001).
ROSETTA predicts the van Genuchten parameters and
Ks in a hierarchical manner from soil textural
fractions, bulk density, and one or two water retention
points. Values for ups ; upr ; a; and n estimated with
3
23
ROSETTA (assuming l ¼ 0:5) were 0.407 cm cm ,323
21
0.066 cm cm , 0.055 cm , and 1.50, respectively.

3. Results and discussion

68

F. Abbasi et al. / Journal of Hydrology 290 (2004) 63–79

Fig. 4. Measured (un-scaled), scaled, and reference soil water retention curves.

4. Model calibration
Simultaneously optimized Ks values and CDE andMIM solute transport parameters for the inlet andoutlet sites of the FD experiment are listed in Table 1.
Estimated Ks values were similar for both sites along
the furrow, and comparable to those estimated with
the ROSETTA code (Schaap et al., 2001) as well withthose obtained earlier for short furrows at the same
field (Abbasi et al., 2003b,c). The Ks value estimatedwith ROSETTA was 0.097 cm min21.
Notice that the optimized CDE and MIM solutetransport parameters (Table 1) are very similar at agiven location, but differ substantially between theinlet and outlet sites. The considerable differences
between the longitudinal ðlL Þ and transverse ðlT Þdispersivities of the inlet and outlet sites reflect theconsiderable spatial variability in the soil hydraulicand transport characteristics in spite of scaling of thewater retention data. This discrepancy may be partly
due also to slightly different imposed flow conditions(flow depths and water and solute application times) at
the two sites. The lL and lT values are comparable to

the range of values obtained for the short furrows(Abbasi et al., 2003b,c). The lT value for the inlet site
was somewhat larger than those reported earlier forthe short furrows. We believe that this is caused bymore lateral water flow and solute spreading at the
inlet site, which had a higher water level in the furrow
and longer water and solute application times ascompared to the outlet site. Also contributing to thehigher lT value of the inlet site likely was its
relatively low initial water content and the larger
scale of the FD experiment. Fig. 7a and b show thatlateral solute spreading was very pronounced at the
inlet site. Fried (1975) and Gelhar et al. (1992), among
others, showed that the dispersivity increases with thescale of observation, while studies by White et al.
(1986) and Steenhuis and Muck (1988) and others
showed that the initial soil water content can also havean important effect. However, the latter is in contrast
to a study by Flury et al. (1994) who showed that of 14
different field sites the initial soil water content had
relatively little impact at several of the sites and noeffect at all at other sites. Numerical studies also
showed a dependency of lateral solute spreading on

Table 1
Summary of the simultaneously optimized saturated hydraulic conductivity, Ks ; and CDE and MIM solute transport parameters for the FD inlet
and outlet sites
Site

Ks (cm min21)

uim (cm3 cm23)

v (min21)lL (cm)lT (cm)SSQ ( –)

R2 (–)

Inlet
OutletAveragea

0.0835 (0.0847)
0.0953 (0.0945)
0.0892 (0.0895)

0.030
0.045
0.0367

0.028
0.020
0.024

22.2 (24.2)
9.1 (9.2)
14.3 (14.9)

4.4 (4.2)
0.01 (0.07)
0.21 (0.54)

5.975 (5.797)
1.413 (1.328)–

0.829 (0.828)
0.849 (0.847)–

aOptimized values for the CDE runs are given in parentheses.
Geometric average of estimated parameters for the inlet and outlet sites used for model validation.

F. Abbasi et al. / Journal of Hydrology 290 (2004) 63–79
Table 2a
Summary of optimized soil hydraulic parameters during the firststep of the sequential (two-step) parameter estimation method
Site

n
(–)

us
(cm3 cm23)

Ks(cm min21)

SSQ
(–)

R2
(–)

Inlet
OutletAveragea

1.53
1.31
1.41

0.301
0.387
0.341

0.04240.1066
0.0679

0.317
0.795
–

0.733
0.315
–

a
Geometric average of estimated parameters for the inlet andoutlet sites used for model validation.

the initial soil water content between and below the
furrows (Abbasi, 2003).
Similar to the small-scale experiments of Abbasi
et al. (2003c), the FD experiment revealed relatively
low immobile water contents ðuim Þ; thus indicatingthat advection and dispersion were the main transport
processes in this study. Estimated first-order exchange
coefficients ðvÞ were generally somewhat lower thanthose obtained for the small-scale experiments, but
comparable with those reported in the literature for
one-dimensional (1D) water flow and solute transport
studies (Nkedi-Kizza et al., 1983; Seyfried and Rao,
1987; Maraqa et al., 1999; among others).
Estimated soil hydraulic and solute transportparameters using the sequential (two-step) optimization method are given in Table 2. Notice that the sumof the squared residuals (SSQ) are much lower thanthose obtained with the simultaneous approach.
However, the R2 values for the two-step method are
generally lower than those for the simultaneous
approach. The soil hydraulic parameters (Table 2a)
are consistent with those obtained earlier for the shortfurrows (Abbasi et al., 2003b). The value of n is also
comparable with laboratory-measured values, as wellas with those estimated with ROSETTA , while us values
are somewhat smaller than those estimated withROSETTA and measured in the laboratory on small

69

soil samples. Estimated Ks values differed considerably between the inlet and outlet sites, and werelower/higher than those obtained with simultaneous
optimization method (Table 1). However, they arestill within range of estimated values obtained for the
short furrows.
Estimated solute transport parameters (Table 2b)
for the CDE and MIM models were again similar;
however, they varied substantially between the twosites. Since they produced very similar results, theCDE model is preferred because it has fewer
parameters, is less vulnerable to parameter uniquenessproblems, and requires less computational time. The
estimated uim and v values for the inlet site were
substantially larger than those for the outlet site andthose estimated using the simultaneous optimization
method (Table 1). The lL and lT parameter values
differed also considerable for the two sites. Overall,estimated lL values using the simultaneous and twostep approaches for both inlet and outlet sites of theFD experiment were consistent with those obtained
for the short furrow experiments assuming homogeneous soil profiles (Abbasi et al., 2003c), butsomewhat lower than those estimated for the same
experiments assuming layered profiles (Abbasi et al.,2003b). Discrepancies between the estimated parameters for the long FD experiment and the shortfurrows may be due to the inherent spatial variability
in the soil hydraulic and solute transport properties,
and partly due to the imposed flow conditions (flowdepths and water and solute application times).
Measured soil water contents and solute concentration values for both the inlet and outlet sites arecompared in Fig. 5 with the predicted values. Twostep optimization produced better agreement withthe observed water contents, whereas solute concentrations obtained with both optimization methods
were more or less the same, with both of them giving

Table 2b
Summary of optimized solute transport parameters during the second step of the sequential (two-step) parameter estimation methodSiteuim (cm3 cm23)

v (min21)lL (cm)lT (cm)SSQ ( –)

R2 (–)

Inlet
OutletAveragea

0.10
0.036
0.06

0.34750.03140.104520.05 (19.85)
1.74 (1.84)
5.91 (6.04)

4.34 (4.34)
0.04 (0.03)
0.422 (0.36)

0.308 (0.308)
0.159 (0.159)–

0.703 (0.702)
0.850 (0.849)–

aOptimized values for the CDE runs are given in parentheses.
Geometric average of estimated parameters for the inlet and outlet sites used for model validation.

70

F. Abbasi et al. / Journal of Hydrology 290 (2004) 63–79

relatively poor descriptions of the data. Notice that theR2 values in Fig. 5 represent regression of measured
versus predicted results for the inlet and outlet sites
whereas the R2 values given in Tables 1 and 2 hold forspecific data types used in the individual inlet andoutlet optimizations.
Measured and calibrated (using the simultaneousand two-step optimization approaches) soil watercontents and solute concentrations at the inlet andoutlet sites of the FD experiment are presented in Figs.
6 and 7. Results are given by means of 1D curves toprovide a better visual comparison between themeasured and calculated distributions. Results are
given at two different times (12 h and five days after
the first irrigation, being representatives of relatively
wet and dry conditions, for soil water contents andfive days after each irrigation for concentrations)and for three different locations in the furrowcross-section (bottom, side and top of the furrow) up
to a depth of 100 cm below the soil surface. Results
are plotted versus depth (instead of versus lateral
distance) since considerable more data were available

versus depth. The CDE and MIM model simulations
were always very similar (Tables 1 and 2). For thisreason, for both model calibration and validation, weare showing here only results obtained with the CDEtransport model.
The simultaneous optimization approach overestimated soil water contents of the different inlet
sampling locations (Fig. 6a and b). We believe thathysteresis, which often occurs in field-scale studies,and air entrapment, may have played a major role inthis study since the calculated water contents immediately after each irrigation quickly reached saturation(us ¼ 0:411 cm3 cm23 as measured on small soil
samples in the laboratory) in almost the entire soil
profile. By comparison, the maximum measured watercontent was only about 0.34 cm3 cm23 (about 83% ofus ), which occurred in the surface layer below
the furrow bottom (results not further shown here).
These results agree with findings by Klute (1986),
among others, who showed that field-measured us
values are often much lower than the porosity becauseof entrapped or dissolved air. Using laboratory

Fig. 5. Measured versus predicted (using simultaneous and two-step parameter optimizations) soil water contents and solute concentrations (forboth the CDE and MIM models) for the inlet and outlet sites of the free-draining experiment.

F. Abbasi et al. / Journal of Hydrology 290 (2004) 63–79

71

Fig. 6. Measured and calibrated (using simultaneous and two-step parameter optimizations) soil water contents for the inlet and outlet sites ofthe free-draining experiment (measured: symbols, simultaneous: solid lines, two-step: dashed lines).

estimated soil water retention parameters to describe
soil water status at the field scale is another possible
reason for the discrepancy between the measured
and simulated soil water contents using the simultaneous parameter estimation approach. Estimation

and application of soil hydraulic properties should bemade in the same scale. Agreement between measured
and calculated water contents was much better for theoutlet site than the inlet site (Fig. 6c and d). Hysteresis
likely was far less important at the outlet site because

72

F. Abbasi et al. / Journal of Hydrology 290 (2004) 63–79

Fig. 7. Measured and calibrated (using simultaneous and two-step parameter optimizations, and the CDE solute transport model) bromide
concentrations for the inlet and outlet sites of the free-draining experiment at different times (measured: symbols, simultaneous: solid lines, twostep: dashed lines).

of shorter water and solute application (opportunity)
times, which caused predicted water contents never toreach us : For instance, opportunity times for the inlet
and outlet sites during the first irrigation were 300 and
63 min, respectively. Our data were not sufficient totake into account the hysteresis.

Simulations obtained using the two-step method
were satisfactory compared to the simultaneous
optimization approach. This was expected since more
hydraulic parameters were optimized with the two-stepmethod (in particular us ). Agreement between
measured and simulated soil water contents during

F. Abbasi et al. / Journal of Hydrology 290 (2004) 63–79

the second irrigation was more or less the same as forthe first irrigation (results not further shown here).
Overall, agreement between measured and calculated solute concentrations (Fig. 7) was less satisfactory, likely because of spatial and temporal variability
in the soil hydraulic and solute transport properties.
Temporal variability was caused by physical deterioration of the furrow surface during the two irrigation
events. Both optimization methods produced poor
results at several locations (Fig. 7a and b; bottom of
the furrow, and Fig. 7d; bottom and side of the
furrow), but better agreement with the data in othercases (Fig. 7c; all locations).
As noted earlier, flow and solute transportparameters can display substantial spatial and temporal variations that are difficult to characterize since
perfect knowledge about the variations is generally
impossible to acquire. Many methods, often forsaturated media, have been employed to approximate
spatial distributions, such as zonation, pilot points,

73

splines, and finite elements (Yeh, 1986; Carrera,1987). In a recent transient study, Forrer et al. (1999)
processed their 2D solute distributions by means ofhorizontal and vertical averaging to reduce spatial
variability and produce 1D concentration profiles forcomparisons with model predictions. They found thathorizontally averaged solute distributions were muchsmoother than vertically averaged distributions.
To further lessen the effects of spatial variability(in addition to scaling of the soil water retention data),
solute concentration data were subjected to horizontal
and vertical averaging. Horizontal averages of themeasured and predicted bromide concentrations forthe different FD sampling sites are shown in Fig. 8.
Following Forrer et al. (1999), we used simple
arithmetic averaging. Vertically averaged results are
not presented because of limited data at each site (onlythree points were available). Horizontal averaging
significantly improved agreement between measured
and model predictions for the first irrigation of the FD

Fig. 8. Horizontally averages of measured and calibrated (using simultaneous and two-step parameter optimizations, and the CDE solutetransport model) bromide concentrations for the different sites of the free-draining experiment (measured: symbols, simultaneous: solid lines,
two-step: dashed lines).

74

F. Abbasi et al. / Journal of Hydrology 290 (2004) 63–79

experiment (Fig. 8). However, the expected improvement was not achieved for the second irrigation (fivedays after the second irrigation). Note that results for
the FD middle sites, predicted by geometric means of
the obtained values for the inlet and outlet sites given
in Tables 1 and 2, were also included in Fig. 8.

5. Model validation
Optimized soil hydraulic and CDE solute transportparameters obtained with the simultaneous and twostep optimization methods from the FD experimentwere used to experimentally test (validate) the
HYDRUS-2D model for the other three experiments.
For this purpose, we used geometric averages of the
optimized parameters from the inlet and outlet sites of
the FD experiment (Tables 1 and 2) to predict soil
water contents, solute concentrations, and deep
percolation rates of water and solutes from the
100%, FH, and SH experiments. Horizontally averaged solute concentrations were presented because oflarge spatial variability between the different sites in
the furrow cross-section.
Measured and simulated soil water contents of thedifferent sites (inlet, middle, and outlet) of the 100%

Fig. 9. Measured and predicted (using simultaneous and two-step parameter optimizations) soil water contents for the different sites of the 100%experiment, five days after irrigation (measured: symbols, simultaneous: solid lines, two-step: dashed lines).

F. Abbasi et al. / Journal of Hydrology 290 (2004) 63–79

experiment are compared in Fig. 9. Despite considerable variability and uncertainty in the soil hydraulicparameters, typical of field-scale studies, agreement
between measured and predicted soil water contents
was relatively good. Similar to model calibration, thetwo-step optimization approach produced better
predictions with the measured values. Discrepancies
in some cases (Fig. 9b) were due to the some of
the factors discussed under model calibration(e.g. hysteresis and variability in the soil hydraulic75

properties). Results for the FH and SH experiments
were more or less like 100% experiment (not further
shown here).
Neither the simultaneous nor the two-step optimization methods described the solute concentrations
well. Concentrations were overestimated at some ofthe sites (Fig. 10a, middle and outlet sites). We
believe that this was mostly due to having insufficient
data for the optimization procedure to accurately
estimate the transport parameters, and because of

Fig. 10. Horizontally averages of measured and predicted (using simultaneous and two-step parameter optimizations, and the CDE solutetransport model) bromide concentrations for the different sites of the 100%, FH, and SH experiments, five days after irrigation (measured:symbols, simultaneous: solid lines, two-step: dashed lines). Note that water did not reach the end in the SH experiment.

76

F. Abbasi et al. / Journal of Hydrology 290 (2004) 63–79

Table 3
Measured and predicted (using the simultaneous and two-stepparameter optimization methods and the CDE solute transportmodel) deep percolation rates of water and solutes at the different
sites of the various experiments
Experiment

Measured

Predicted(simultaneous)Predicted(two-step)

Water
(l/m)

Solute
(g/m)

Water
(l/m)

Solute
(g/m)

Water
(l/m)

Solute
(g/m)

FD
Inlet
MiddleOutlet

82.7
50.2
24.1

139.7
95.9
8.4

129.2
126.0
5.9

382.5
363.6
0.0

65.1
39.4
4.6

101.0
30.9
0.0

100%
Inlet
MiddleOutlet

57.0
71.6
46.6

13.0
3.2
1.1

91.4
49.7
22.2

78.7
31.2
8.7

68.9
39.5
19.6

31.7
10.1
2.1

FH
Inlet
MiddleOutlet

68.6
66.9
45.9

16.1
11.2
0.7

28.9
18.3
1.6

11.8
4.2
0.0

10.7
4.8
0.0

1.1
1.1
0.0

SH
Inlet
Middle

40.3
23.6

23.0
5.8

34.4
2.9

16.0
0.5

14.1
0.0

0.4
0.0

averaging of the optimized transport parameters
obtained from the inlet and outlet sites (which showed
relatively large differences in dispersivities; see
Tables 1 and 2). Also, point measurements often leadto poor representations of water flow and transport

behavior at field scale (Tsang et al., 1996). Point
measurements additionally can cause considerable
underestimation of chemical fluxes as shown byRitsema and Dekker (1996), and generally involve
uncertainties induced by human, instrumental, andhydrogeological errors (Xiang et al., 1992). Bothoptimization methods predicted the concentrations
near the soil surface relatively well, except for the
inlet site of the 100% experiment where they were
underestimated. Notice that the irrigation water in the
furrow never reached the end of the SH experimentbecause of high infiltration and roughness properties.
Measured and predicted deep percolation rates forwater (WDP) and solute (SDP) five days after
irrigation are presented in Table 3. WDP and SDP
are defined as the amounts of water and solute,respectively, percolating below a depth of 100 cm. In
order to calculate WDP and SDP, soil profiles below
the depth of 100 cm down to 180 cm were divided
into a network of rectangular elements, each being
approximately 25 cm by 30 cm. Volumes of water andmasses of bromide were determined for eachrectangle, and then summed to obtain the total
amounts of water/bromide. We used for this purpose
measured gravimetric soil water contents, soil bulk
densities, and bromide concentrations at different
depths.
Deep percolation rates of water and solutes
were overestimated/underestimated for the various
experiments (Table 3). The simultaneous approach

Fig. 11. Measured and predicted (using simultaneous and two-step parameter optimizations, and the ROSETTA code) soil water contents at thebottom of the inlet, middle, and outlet sites of the free-draining experiment, five days after the first irrigation.

F. Abbasi et al. / Journal of Hydrology 290 (2004) 63–79

overestimated the rates for the inlet and middle sitesof the FD experiment, while the two-step method
predicted the values somewhat better for those sites.
However, both methods underestimated the rates forthe outlet site. Both optimization approaches underestimated deep percolation rates of water and solutes
for the different sites of the various blocked-end
experiments, except those estimated with the simultaneous method for the inlet site of the 100%experiment. The two-step predictions were slightly
better for the inlet site of the 100% experiment,
whereas the simultaneous predictions were somewhat
better for other sites. The rates for the FH and SH
experiments were also predicted somewhat better withthe simultaneous approach.
Finally, we compared predicted soil water content
distributions obtained using the simultaneous andtwo-step optimization methods, against predictions
using pedotransfer functions obtained with the
ROSETTA code of Schaap et al. (2001). Essentiallythe same water content distributions were obtainedwith ROSETTA and the two optimization approaches,
exemplified by the plots in Fig. 11. This shows thatROSETTA can be a relatively fast, easy, and costeffective method for estimating soil hydraulic properties, as compared to using direct measurement
methods.

6. Conclusions
In this study we calibrated and experimentallytested (validated) the HYDRUS-2D numerical code
against four real furrow irrigation experiments. The
soil hydraulic and solute transport parameters were
estimated using simultaneous and sequential optimization of data from two sites near the inlet and outlet
sites of a FD furrow experiment assuming homogeneous soil profiles. The optimized parameters were
subsequently used for model validation to produce
soil water contents and solute concentrations in threeother blocked-end furrow experiments conducted in
the same field.
The optimized parameters were very similar to
those reported earlier for short furrows at thesame experimental site. Minor differences wererelated to the presence of soil heterogeneity acrossthe field and versus depth, different observation77

scales, and different imposed flow conditions (flowdepths and water and solute applications times).
Optimization resulted in similar parameters for the
CDE and MIM transport models. In spite of scaling of
the laboratory-measured soil water retention data,
estimated parameters along the monitored furrow,particularly the longitudinal and transverse dispersivities, were found to be vary spatially because ofinherently different soil properties and differently
imposed flow and solute transport boundary conditions (e.g. varying flow depths, and irrigation andsolute application times). Agreement between
measured and predicted soil water contents was
satisfactory with the sequential optimization approach
whereas the simultaneous method overestimated thesoil water contents. Observed solute concentrations
along the furrow cross-sections displayed more spatial
variability than the soil water contents and were notpredicted as well as the soil water contents. Spatial
averaging across the furrow cross-sections somewhat
improved agreement between observed and calculated
solute concentrations. Both optimization approaches
predicted the solute concentrations near the soil
surface reasonably well. Deep percolation rates ofwater and solutes were overestimated/underestimatedwith the simultaneous method, and underestimatedwith the sequential approach. Two-step optimization
produced better agreement between predicted andmeasured soil water contents, but required morecomputational time as compared to simultaneous
optimization.
This study provides further evidence that theunsaturated soil hydraulic properties should be estimated and used in model predictions at the same scale.
Estimated soil hydraulic parameters from small soilcores in the laboratory may not adequately describe
field-scale soil water status in view of considerable
variability, uncertainty, and a number of experimental
complications (especially air entrapment and macropore flow). Also, considering temporal variability inthe soil hydraulic and transport parameters could well
improve the model predictions further.

Acknowledgements
This material is based on work supported in part by
SAHRA (Sustainability of semi-Arid Hydrology

78

F. Abbasi et al. / Journal of Hydrology 290 (2004) 63–79

and Riparian Areas) under the STC Program of the
National Science Foundation, Agreement No. EAR9876800. The authors would like to thank the threeanonymous referees for their useful comments andsuggestions.

References
Abbasi, F., 2003. Field-scale analyses of water flow and solutetransport in furrows. PhD Thesis, No. 555, K.U. Leuven, Fac.
Landbouwkundige en Toegepaste Biologische Wetenschappen,
Leuven, Belgium.
Abbasi, F., Feyen, J., Roth, R.L., Sheedy, M., van Genuchten,
M.Th., 2003a. Water flow and solute transport in furrowirrigated fields. Irrig. Sci. 22 (2), 57 –65.
Abbasi, F., Jacques, D., Šimůnek, J., Feyen, J., van Genuchten,
M.Th., 2003b. Inverse estimation of the soil hydraulic andsolute transport parameters from transient field experiments:
heterogeneous soil. Trans. ASAE 46 (4), 1097–1111.
Abbasi, F., Šimůnek, J., Feyen, J., van Genuchten, M.Th., Shouse,
P.J., 2003c. Simultaneous inverse estimation of the soil hydraulicand solute transport parameters from transient field experiments:
homogeneous soil. Trans. ASAE 46 (4), 1085–1095.
Abbasi, F., Šimůnek, J., van Genuchten, M.Th., Feyen, J., Adamsen,
F.J., Hunsaker, D.J., Strelkoff, T.S., Shouse, P., 2003d. Overland water flow and solute transport: model development andfield data analysis. J. Irrig. Drain. Engng 129 (2), 71–81.
Carrera, J., 1987. State of the art of the inverse problem applied to theflow and solute transport equation, in analytical and numericalgroundwater flow and quality modelling. In: Custodio, E., et al.
(Eds.), NATO-ARW Series C: Mathematical and PhysicalSciences, vol. 224. Reidel, Norwell, Mass, pp. 549 –583.
Carrera, J., Neuman, S.P., 1986. Estimation of aquifer parameters under transient and steady-state conditions: application to synthetic and field data. Water Resour. Res. 22 (2),
228–242.
Clothier, B.E., Vogeler, I., Green, S.R., Scotter, D.R., 1998.
Transport in unsaturated soil: aggregates, macropores, andexchange. In: Selim, H.M., Ma, L. (Eds.), Physical Nonequilibrium in Soils: Modeling and Application, Ann Arbor
Press, Chelsea, MI, pp. 273– 295.
Flury, M., Fluhler, H., Jury, W.A., Leuenberger, J., 1954.
Susceptibility of soils to preferential flow of water: a fieldstudy. Water Resour. Res. 30, 1945–1954.
Forrer, I., Kasteel, R., Flury, M., Fluhler, H., 1999. Longitudinal andlateral dispersion in an unsaturated field soil. Water Resour. Res.
35, 3049–3060.
Fried, J.J., 1975. Groundwater Pollution. Elsevier, New York.
Gelhar, L.W., Welty, C., Rehfeldt, K.R., 1974. A critical review ofdata on field-scale dispersion in aquifers. Water Resour. Res. 28(7), 1955–1955.
van Genuchten, M.Th., 1980. A closed-form equation for predicting
the hydraulic conductivity of unsaturated soils. Soil Sci. Soc.
Am. J. 44, 892–898.

van Genuchten, M.Th., Wagenet, R.J., 1989. Two site/two region
models for pesticide transport and degradation: theoreticaldevelopment and analytical solutions. Soil Sci. Soc. Am. J. 53(5), 1303–1310.
Hopmans, J.W., Šimůnek, J., Romano, N., Durner, W., 2002.
Inverse methods. In: Dane, J.H., Topp, G.C. (Eds.), Methodsof Soil Analysis. Part 4. Physical Methods, SSSA Book Series
5, Soil Science Society of America, Madison, WI, pp.
963 –1008.
Jacques, D., Šimůnek, J., Timmerman, A., Feyen, J., 2002.
Calibration of Richards’ and convection–dispersion equations
to field-scale water flow and solute transport under rainfall
conditions. J. Hydrol. 259, 15 –31.
Klute, A., 1986. Water retention: laboratory methods. In: Klute, A.,(Ed.), Methods of Soil Analysis. Part 1. Physical andMineralogical Methods, 2nd ed., Agronomy, 9(1). American
Society of Agronomy, Madison, WI, pp. 635 –662.
Maraqa, M.A., Wallace, R.B., Voice, T.C., 1999. Effects ofresidence time and degree of water saturation on sorption
nonequilibrium parameters. J. Contam. Hydrol. 36, 53 –72.
Marquardt, D.W., 1963. An algorithm for least squares
estimation of non-linear parameters. J. Ind. Appl. Math. 11,
431 –441.
Mishra, S., Parker, J.C., 1989. Parameter estimation for coupledunsaturated flow and transport. Water Resour. Res. 25 (3),385 –396.
Mualem, Y., 1976. A new model for predicting the hydraulicconductivity of unsaturated porous media. Water Resour. Res.
12 (3), 513–522.
Nkedi-Kizza, P., Biggar, J.W., van Genuchten, M.Th., Wierenga,
P.J., Selim, H.M., Davidson, J.M., Nielsen, D.R., 1983.
Modeling tritium and chloride36 transport through an aggregated oxisol. Water Resour. Res. 19, 691 –700.
Ritsema, C.J., Dekker, L.W., 1996. Influence of sampling strategy
on detecting preferential flow paths in water-repellent sand.
J. Hydrol. 177, 33–45.
Schaap, M.G., Leij, F.J., van Genuchten, M.Th., 2001. ROSETTA :a computer program for estimating soil hydraulic parameterswith hierarchical pedotransfer functions. J. Hydrol. 251,
163 –176.
Seyfried, M.S., Rao, P.S.C., 1987. Solute transport in undisturbed
columns of an aggregated tropical soil: preferential flow effects.
Soil Sci. Soc. Am. J. 51, 1434–1444.
Šimůnek, J., Hopmans, J.W., 2002. Parameter estimation andnonlinear fitting. In: Dane, J.H., Topp, G.C. (Eds.), Methodsof Soil Analysis. Part 4. Physical Methods, SSSA Book
Series 5, Soil Science Society of America, Madison, WI, pp.
139 –157.
Šimůnek, J., Sejna, M.Th., van Genuchten, M., 1999. The
HYDRUS-2D Software Package for Simulating the Twodimensional Movement of Water, Heat, and Multiple Solutes
in Variably Saturated Media, Version 2.0, IGWMC-TPS-70,
International Ground Water Modeling Center, Colorado Schoolof Mines, Golden Colorado.
Šimůnek, J., van Genuchten, M.Th., Jacques, D., Hopmans, J.W.,
Inoue, M., Flury, M., 2002. Solute transport during variablysaturated flow: inverse methods. In: Dane, J.H., Topp, G.C.

F. Abbasi et al. / Journal of Hydrology 290 (2004) 63–79(Eds.), Methods of Soil Analysis. Part 4. Physical Methods,
SSSA Book Series 5, Soil Science Society of America,
Madison, WI, pp. 1435– 1449.
Steenhuis, T.S., Muck, R.E., 1988. Preferred movement of
nonadsorbed chemicals on wet, shallow, sloping soils.
J. Environ. Qlty 17 (3), 376–384.
Sun, N.Z., Yeh, W.W.G., 1990. Coupled inverse problems
in groundwater modeling. 1. Sensitivity analysis andparameter identification. Water. Resour. Res. 26 (10),2507–2525.
Tsang, Y.W., Tsang, C.F., Hale, F.V., Dverstorp, B., 1996. Tracer
transport in a stochastic continuum model of fractured media.
Water Resour. Res. 32 (10), 3077–3092.

79

Vogel, T., Cislerova, M., Hopmans, J.W., 1991. Porous media withlinearly variable hydraulic properties. Water Resour. Res. 27,
2735–2741.
White, R.E., Dyson, J.S., Gerstl, Z., Yaron, B., 1986. Leaching of
herbicides through undisturbed cores of a structured clay soil.
Soil Sci. Soc. Am. J. 50, 277 –283.
Xiang, Y., Sykes, F., Thomson, N.R., 1992. A composite L1parameter estimator for model fitting in groundwater flow andsolute transport simulation. Water Resour. Res. 29 (6),1661–1673.
Yeh, W.W.-G., 1986. Review of parameter identification procedures in groundwater hydrology: the inverse problem. Water
Resour. Res. 22 (2), 95 –108.