SIMULTANEOUS INVERSE ESTIMATION OF SOIL HYDRAULICAND SOLUTE TRANSPORT PARAMETERS FROM TRANSIENTFIELD EXPERIMENTS: HOMOGENEOUS SOIL
F. Abbasi, J. Simunek, J. Feyen, M. Th. van Genuchten, P. J. Shouse
ABSTRACT. Inverse estimation of unsaturated soil hydraulic and solute transport properties has thus far been limited mostly
to analyses of one–dimensional experiments in the laboratory, often assuming steady–state conditions. This is partly becauseof the high cost and difficulties in accurately measuring and collecting adequate field–scale data sets, and partly becauseof difficulties in describing spatial and temporal variabilities in the soil hydraulic properties. In this study, we estimated soil
hydraulic and solute transport parameters from several two–dimensional furrow irrigation experiments under transientconditions. Three blocked–end furrow irrigation experiments were carried out, each of the same duration but with different
amounts of infiltrating water and solutes resulting from water depths of 6, 10, and 14 cm in the furrows. Two more experiments
were carried out with the same amounts of applied water and solute, and hence for different durations, on furrows with water
depths of 6 and 10 cm. The saturated hydraulic conductivity (Ks ) and solute transport parameters in the physical equilibrium
convection–dispersion (CDE) and physical nonequilibrium mobile immobile (MIM) transport models were inversely
estimated using the Levenberg–Marquardt optimization algorithm in combination with the HYDRUS–2D numerical code. Soilwater content readings, cumulative infiltration data, and solute concentrations were used in the objective function during
the optimization process. Estimated Ks values ranged from 0.0389 to 0.0996 cm min–1, with a coefficient of variation of 48%.
Estimated immobile water contents (qim ) were more or less constant at a relatively low average value of 0.025 cm3 cm–3,
whereas the first–order exchange coefficient (w) varied between 0.10 and 19.52 min–1. The longitudinal dispersivity (DL )ranged from 2.6 to 32.8 cm, and the transverse dispersivity (DT ) ranged from 0.03 to 2.20 cm. DL showed some dependency
on water level and irrigation/solute application time in the furrows, but no obvious effect was found on Ks and other transportparameters, most likely because of spatial variability in the soil hydraulic properties. Agreement between measured andpredicted infiltration rates was satisfactory, whereas soil water contents were somewhat overestimated, and soluteconcentrations were underestimated. Differences between predicted solute distributions obtained with the CDE and MIM
transport models were relatively small. This and the value of optimized parameters indicate that observed data were
sufficiently well described using the simpler CDE model, and that immobile water did not play a major role in the transport
process.
Keywords. Furrow irrigation, Homogeneous soil, Inverse solution, Solute transport, Water flow.

While considerable progress has been madeduring the past several decades in describing
and modeling water flow and solute transport
processes under controlled, often steady–state, conditions at the laboratory scale, detailed analyses offield–scale experiments remain limited, mostly because oflabor and cost requirements, but also because of inherent
complications posed by field–scale heterogeneity. For
example, the unsaturated soil hydraulic properties canchange by several orders of magnitude, even over relatively

Article was submitted for review in July 2002; approved for publicationby the Soil & Water of ASAE in May 2003.
The authors are Fariborz Abbasi, Soil Physicist, and Jan Feyen,
Professor, Institute for Land and Water Management, Katholieke
Universiteit Leuven, Belgium; and Jirka Simunek, Soil Physicist,
Martinus Th. van Genuchten, Research Leader, and Peter J. Shouse, Soil
Physicist, USDA–ARS, George E. Brown Jr. Salinity Laboratory,
Riverside, California. Corresponding author: Jan Feyen, Institute forLand and Water Management, Vital Decosterstraat 102, 3000–Leuven,
Belgium; phone: +32–16–329756; fax: +32–16–329760; e–mail: jan.
feyen@agr.kuleuven.ac.be.

short distances (Biggar and Nielsen, 1976; Bresler et al.,1984), thus making it difficult to assign values that are
applicable at the field scale.
A wide range of models exists for simulating water flow
and solute transport in the vadose zone. Addiscott andWagenet (1985) gave a review of one–dimensional (1–D)
transport models. They compared and classified the models
as being deterministic or stochastic, mechanistic or functional, numerical or analytical, and research or management
oriented. They also discussed the degree of complexity,
flexibility, transferability, and usefulness of these models forfield conditions. One–dimensional water flow and solutetransport models (such as those developed by Jarvis et al.,
1991; Sulekha and Duijnisveld, 1998; Shao et al., 1998;
Simunek et al., 1998; among others) are useful for many
applications. While requiring less expertise in programming,
fewer input data, and less computer time and memory than
two–dimensional (2–D) models, 1–D approaches generally
do not accurately describe geometry and thus the dynamics
of water and solutes during irrigation using furrows, drippers,
perforated tubes, and porous lines, or some other method. Insuch cases, flow and transport are two– or even three–dimen-

Transactions of the ASAE
Vol. 46(4): 1085–1095

2003 American Society of Agricultural Engineers ISSN 0001–2351

1085

sional (3–D), thus requiring a multi–dimensional simulation
model.
A number of 2–D (e.g., Neuman et al., 1975; Simunek etal., 1999a) and 3–D (Wu and Chieng, 1995a, 1995b; Russo
et al., 1998) models currently exist, with the governing
equations generally being solved using a variety of finite
difference, finite element, or other numerical methods (e.g.,
Segol, 1994; Forsyth et al., 1995; van Genuchten and
Sudicky, 1999, Kees and Miller, 2001). Multi–dimensionalmodels are more suitable for studying the effects of soil
spatial variability on water flow and solute transport at thefield scale. Unfortunately, field applications of such models
remain limited since they generally require large amounts of
input data and field measurements, and consequently can
become quite costly to implement. For example, Wu andChieng (1995a, 1995b) developed a three–dimensional
multi–component solute transport model. The finite element
model was used successfully to simulate chemical reactions
such as acid–base reaction, complexation, ion exchange,
precipitation and dissolution, as well as heat transfer andvapor movement. Loretta and Wu (1999) later used thismodel to study the transport of four heavy metals through aclay barrier. Russo et al. (1998) similarly presented a 3–D
model for simulating water flow and solute transport within
the unsaturated zone. They used a fully implicit Euler scheme
to solve the water flow equation, while solute transport wasapproximated by an operator–splitting approach. Application to bromacil transport showed that the patterns ofmoisture contents, pressure heads, pressure gradients, hydraulic conductivities, and root–water uptake were considerably affected by spatial heterogeneity in the soil hydraulicproperties.
Of special interest to this study is the HYDRUS–2D codedeveloped by Simunek et al. (1999a). HYDRUS–2D is a
Microsoft Windows–based modeling environment for analyzing water flow and solute transport in 2–D variably
saturated porous media, as well as for inverse estimation ofsoil hydraulic and/or solute transport properties. The governing equations are numerically solved using Galerkin–typelinear finite elements, while the Levenberg–Marquardt
algorithm (Marquardt, 1963) is used for parameter optimization. The model has been successfully used in manylaboratory and field studies (e.g., Gribb, 1996; Simunek andvan Genuchten, 1996, 1997; Kodesova et al., 1998, 1999) tooptimize soil hydraulic properties. For example, Simunek
and van Genuchten (1996) showed that measured infiltration
data using a tension disc infiltrometer at constant tension did
not provide enough information to inversely estimate soil
hydraulic properties, whereas using infiltration data at
multiple tensions yielded good estimates of the soil hydraulicproperties (Simunek and van Genuchten, 1997). The model
was also successfully used to inversely estimate soil
hydraulic properties from cone permeameter experiments(Gribb, 1996; Kodesova et al., 1998, 1999; Simunek et al.,1999b). In other studies, Wang et al. (1997) used theCHAIN–2D model (Simunek and van Genuchten, 1994), an
early version of the HYDRUS–2D, to investigate the effects
of different irrigation methods and spatial variability in thesaturated hydraulic conductivity on subsurface solute transport. Mohanty et al. (1997, 1998) similarly used CHAIN–2D
to study water flow and nitrate transport in a soil profile with
subsurface tile drains. The model predicted nitrate concentrations in the soil profile reasonably well, but underesti-1086

mated nitrate flux concentrations in the subsurface tile drain,
and overestimated those in piezometers.
Despite considerable efforts in studying solute transport at
the field scale during transient flow, particularly during thepast 20 years, the implications of using different irrigation
methods on subsurface solute transport have received
relatively little attention. The limited studies available in theliterature suggest that transport processes can be significantly
affected by irrigation regime and application method (Bowman and Rice, 1986; Jaynes et al., 1988; Wang et al., 1997).
In an extensive furrow irrigation study, Abbasi et al. (2003)
recently showed that water flow depth and irrigation/solute
application time both play a major role in transporting anddistributing water and solutes below the furrows. The main
purpose of this study was to quantify soil hydraulic and solutetransport properties from various 2–D field experiments on
blocked–end furrows under transient conditions. The parameter optimization features of HYDRUS–2D were used toinversely estimate Ks and the transport parameters in theclassical convection–dispersion (CDE) and physical nonequilibrium mobile–immobile (MIM) solute transport models assuming the presence of a homogeneous soil profile.
While the CDE model holds for relatively uniform flow
conditions, the MIM model also considers nonequilibrium
and/or preferential flow associated with the presence ofrelatively immobile liquid regions in the soil. The MIM
model reduces to the CDE model when nonequilibrium is notpresent.

MATERIALS AND METHODSFIELD EXPERIMENTS
Five field experiments on short blocked–end furrows wereconducted at the Maricopa Agricultural Center (MAC) in
Phoenix, Arizona, on a non–vegetative Casa Grande sandy
loam soil (fine–loamy, mixed, hyperthermic Typic Natrargids). Experiments were carried out in plots made up of threefurrows approximately 20 cm deep and spaced 1.0 m apart.
The monitored center furrow of each plot was a non–wheel
furrow, with a wheel track furrow on each side of themonitored furrow. The blocked portion of the furrows was
3 m in length. Details about the experiments can be found in
Abbasi et al. (2003). Below we briefly discuss only thoseparts of the experiments that are directly related to this study.
Two series of experiments were carried out. We first
performed same–duration (SD) experiments involving threeflow depths (6, 10, and 14 cm) used in two successive
irrigations, 10 days apart. The first irrigation took place in
two phases. During the first phase, water was applied for
60 min (the same for all three experiments) to wet the soil
profile. Water was then removed from the furrows andmeasured using an electric scale. During the second phase,water amended with CaBr2 was applied for 30 min again thesame for the three indicated plots. Water in the furrow wasagain removed and the amount determined. The second
irrigation utilized the same depths of unamended water as
employed during the first irrigation, and lasted 90 min. As forthe first irrigation, water stored in the furrow was pumped out
and measured at the end of the experiment.
The second series of experiments involved similar
amounts of applied water and solutes (SWS). The amount ofwater applied was the same as that infiltrated for the 14 cm

TRANSACTIONS OF THE ASAE

depth treatment of the SD experiments. Water levels of 6 and10 cm were used in this case. The first irrigation was carriedout in two steps, similarly to the SD scenarios, with unamended water being applied first followed by bromide–
amended water, except that the irrigation times were adjusted
for each of the water levels such that predefined amounts ofwater, obtained from the previous experiment with 14 cm
flow depth, infiltrated. The second irrigation used unamended water applied with the same water levels as for the
first irrigation, but with times adjusted so that again the sameamount of water infiltrated during the second irrigation as forthe 14 cm depth treatment in the first set of experiments. As
stated above, the 14 cm flow depth experiment plays as a transient experiment between the SD and SWS scenarios.
A set of five neutron probe access tubes (each 3.3 m indepth) was installed to measure soil water contents at different locations perpendicular to the axis of the monitored furrows. The neutron tubes were installed in two rows 50 cm
apart in order to avoid mutual effects on the readings (see Abbasi et al., 2003, for details). A site–calibrated neutron probe
was used to measure soil water contents. The readings weretaken at depths of 20, 40, 60, 80, 100, 140, 180, 220, and 260
cm before each irrigation to provide initial conditions, immediately after each irrigation, then hourly up to 6 hours after
each irrigation, and subsequently each 3 hours up to 24 hours.
Measurements were later taken 3 to 4 times per day up to 3
days after irrigation, and then 1 to 2 times per day up to thenext irrigation. Since neutron probe readings are generally
imprecise near the soil surface, water contents of the surface

layer (0 to 30 cm) were measured using a site–calibrated time
domain reflectometry (TDR) probe. The TDR and neutronprobe readings were taken at the same times as indicatedabove.
Soil samples for analyzing bromide concentrations weretaken manually, at depths and locations corresponding to theneutron probe access tube measurements, four times during
the experiments, i.e., prior to the experiments as initial values, 5 days after the first irrigation, and 6 and 20 days after
the second irrigation. The samples were air dried and crushed
to pass a 2 mm sieve. Soil extractions (1:1 weight:volume)
were made and analyzed for bromide with a Lachat QuikChem flow injection analyzer using standard colorimetric
procedures. In addition, 38 undisturbed soil samples (6 cm
long, 5.4 cm diameter soil cores) were collected randomly at
different locations and depths (up to 100 cm) for laboratory
analyses of the soil water retention curve and the saturated
hydraulic conductivity.
The same bromide concentration (10 g L–1) was used in all
experiments. Flow depths were kept constant during the irrigation by adjusting the water levels to the desired heights
as determined from staff gauges placed at the bottom of the
furrows. Geometries of the experimental furrows were determined before each irrigation at two locations along the furrows in order to calculate the volume of water needed to fill
the furrow section and to infer geometry parameters required
for numerical calculations using HYDRUS–2D. Measured
and fitted geometries of the furrows before each irrigation are
shown in figure 1.

Furrow Width (cm)
0

10

20

30

40

50

60

70

80

90

100

70

80

90

100

Furrow Depth (cm)

0
2.5
5
7.5
10
12.5
15
17.5
20
22.5
25

(a)

Furrow Width (cm)
0

10

20

30

40

50

60

0
Furrow Depth (cm)

2.5
5
7.5
10
12.5
15
17.5
20

(b)

Figure 1. Geometries of the experimental furrows and the fitted profiles. Symbols are the measured data and solid lines are the best fits: (a) before thefirst irrigation, and (b) before the second irrigation.

Vol. 46(4): 1085–1095

1087

6

4

2

0.3
0.2
0.1
0

0

20
Jan

30
Jan

9
Feb

19
Feb

1
Mar

11
Mar

21
Mar

31
Mar

Experimental Period
Figure 2. Reference pan evapotranspiration rates and rates estimated using the Penman–Monteith method (dashed line = Penman–Monteith, andcontinuous line = pan).

All plots were covered with plastic sheets during possible
rainfall events to make sure that irrigation was the only
source of water during the experiments. Reference evapotranspiration rates from the nearest weather station (approximately 150 m from the experimental field) and estimatedevapotranspiration rates using the Penman–Monteith method(Allen et al., 1994) are given in figure 2. The SD and SWS
experiments started on 30 January and 26 February 2001,
respectively, and each lasted 30 days.

MODEL DESCRIPTIONWATER FLOW
Two–dimensional isothermal Darcian flow of water in avariably saturated rigid porous medium is given by thefollowing form of the Richards equation (Richards, 1931):∂θ∂∂h[ K ( K ijA=+ K izA )] − Γ∂t ∂xi
∂xi

(1)

where
q
= volumetric soil water content (L3 L–3)h= pressure head (L)G= sink term (e.g., root water uptake; T–1)
= spatial coordinates (L; i = 1, 2)xit= time (T)Kij A and Kiz A = components of the dimensionlessanisotropy tensor KAK= unsaturated hydraulic conductivity
function (L T–1) given by:
K(h) = Ks Kr (h)

(2)

where
Ks = saturated hydraulic conductivity (L T–1)
Kr = relative hydraulic conductivity.
In this study, KA was assumed to be isotropic. The sinkterm (G) in equation 1 was set to zero since the soil in thisstudy was left bare (no root water uptake).
SOIL HYDRAULIC PROPERTIESThe soil water retention curve, q(h), was described using
the closed–form equation of van Genuchten (1980), and theunsaturated soil hydraulic conductivity function, K(h), wasdescribed using the capillary model of Mualem (1976) as
follows:

1088

0.4(cm 3 cm–3 )Soil Water Content

Reference Evapotranspiration(mm/day)

0.5

8

0

1

2
3
Log Pressure Head (cm)4

5

Figure 3. Measured and fitted soil water retention curves (symbols = measured, dashed line = RETC and Rosetta optimization, and continuousline = RETC optimization).

θ( h ) = θ r +θs − θr(3)(1 + αh n ) m

K(h) = Ks Se l[1 – (1 – Se 1/m)m]2

(4)

where
Se =

θ − θr ,θs − θr1
m = 1− ,
n

n>1

(5)

and
qr
= residual soil water contents (L3 L–3)
= saturated soil water contents (L3 L–3)
qs
Se
= relative saturation (dimensionless)m, n, l, a = empirical parameters (m, n, and l dimension–
less; a in L–1).
Soil water retention data (fig. 3) were determined on small
undisturbed soil samples in the laboratory using Tempe cells(for pressures below 1 bar) and pressure plates (for pressuresover 1 bar). The van Genuchten (VG) retention parameters qr,
qs , a, and n were estimated from the measured q(h) data using
the RETC non–linear optimization program of van Genuchten et al. (1991). The optimized parameters, sum of squared
residuals (SSQ), and regression between the observed andestimated values (R2) are listed in table 1. The optimized
value for qr (0.10 cm3 cm–3) was found to be slightly higher
than the water contents (approximately 0.090 cm3 cm–3) of
the soil surface layers (0 to 0.20 m) measured at the end of theexperiments. For this reason, we also estimated the qr usingneural network–based pedotransfer functions derived by
Schaap et al. (2001) and incorporated into their Rosetta code.
Rosetta predicts VG retention parameters and Ks in ahierarchical manner from soil textural class information, soil
textural fractions, the bulk density, and one or two water
retention points as input (Schaap et al., 2001).
Table 1. van Genuchten soil hydraulic properties obtainedwith the RETC and Rosetta optimization codes.
θr
θs
α
n
SSQ[a]
R2
Method(cm3 cm–3) (cm3 cm–3) (cm–1) (–)
(–)
(–)
RETC
0.10
Rosetta
0.065
Rosetta+RETC 0.065[b]

0.433
0.407
0.407 [b]0.0758 1.478
0.0545 1.503
0.0689 1.37

0.273
––
0.3

0.92
––
0.91

[a] SSQ = sum of squared residuals.
[b] Estimated with Rosetta and used as fixed value for the RETC optimiza-tion.

TRANSACTIONS OF THE ASAEVG parameters estimated with Rosetta (assuming l = 0.5
in eq. 4) are given in table 1. Rosetta predicted qr, qs , and avalues that were somewhat lower than those obtained with
RETC, while both approaches produced almost the samevalues for n. The Rosetta values for both qr and qs agreed wellwith those reported by Carsel and Parrish (1988). As a result,
we used a combination of Rosetta and RETC to finalize thesoil hydraulic parameters: qr and qs were chosen from
Rosetta and used as fixed values in RETC to estimate a andn from the measured retention data (the third option in table
1 and fig. 3).
The saturated conductivity (Ks ) was measured on undisturbed soil samples (6 cm long and 5.4 cm diameter) using theconstant head method. The average Ks from 38 samples was
725 cm day–1, with a coefficient of variation of 73%. The
average Ks was much larger than those of most agriculturalsoils, even coarse–textured soils (e.g., Carsel and Parrish,
1988; Schaap et al., 1998). The large Ks values could havebeen caused by possible disturbance of the samples duringcollection and transportation (from the field site in Arizonato the Salinity Laboratory in Riverside, California), measurement errors that are sometimes as large as one order ofmagnitude, and size effects or altered boundary conditions of
the samples. Direct measurements of Ks are generally
difficult for undisturbed samples taken in the field. In
addition, Ks is usually very sensitive to the saturated watercontent, since small errors in the water content nearsaturation can lead to large errors in Ks . Mallants et al. (1997)
compared saturated hydraulic conductivities of three different undisturbed soils and reported average value of 334, 127,
and 11.8 cm day–1 for small (5 cm long, 5 cm diameter),intermediate (20 cm long, 20 cm diameter), and large(100 cm long, 30 cm diameter) soil columns, respectively,
taken from the same field. Because of the above difficulties,
we decided to estimate Ks from the field–measured soil watercontent and infiltration data using the parameter estimation
features of HYDRUS–2D.
SOLUTE TRANSPORT
The physical non–equilibrium convection–dispersion
model, also known as the mobile–immobile (MIM) or
two–region model, for transport of non–reactive solute invariably saturated porous media is given by (van Genuchten
and Wagenet, 1989; Clothier et al., 1998):∂(θmCm )∂C∂Cm ∂(qi Cm )
∂
(θ D
)−
+ θim im =∂t∂t
∂xi m ij ∂x j
∂xi
θim

∂Cim
= ω(Cm − Cim )
∂t

(6)

(7)

where
qm = mobile water content (L3 L–3)
qim = immobile water content (L3 L–3)
Cm = resident solute concentrations in the mobile region(ML–3)
Cim = resident solute concentrations in the immobile
region (ML–3)
w = first–order exchange coefficient (T–1)
qi = volumetric flux (LT–1)
Dij = dispersion coefficient tensor (L2 T–1) defined asfollows:

Vol. 46(4): 1085–1095

θmDij = DT q δij + (DL − DT )qi q j+ θmDo τoδij
q

(8)

where
Do = ionic or molecular diffusion coefficient in freewater (L2 T–1)
to = tortuosity factor (dimensionless)dij = Kronecker delta function (dimensionless)
DL = longitudinal dispersivity (L)DT = transverse dispersivity (L).
Please note that equation 6 simplifies to the standard
convection–dispersion equation (CDE) when no immobile
water is present (qim = 0).
INITIAL AND BOUNDARY CONDITIONS
Measured bromide concentrations and soil water contents,
the latter taken with a neutron probe, before the experiments
were used as initial conditions within the flow domain. A
constant pressure head (surface ponding, ho in fig. 4) was
specified as the upper boundary condition in the furrow
during irrigation, while an atmospheric boundary condition
was used after irrigation (during redistribution) as:K(KijA∂h+ KizA )ni ≤ E
∂x j

(9)

hA ≤ h ≤ hs

(10)

where
E
= potential evaporation rate
= components of the outward unit vectornihA and hs = minimum and maximum pressure heads,
respectively, allowed under the prevailingsoil conditions.
The value for hA is determined from the equilibrium
conditions between soil water and atmospheric water vapor,whereas hs is usually set to zero. In this study, hA was fixedat –15 bar, being the equilibrium value with a soil surface
water content of 0.09 cm3 cm–3 as measured at the end of theexperiments. The average of the reference and estimatedevapotranspiration rates was used as the atmospheric boundary condition.
A Cauchy condition was used for the upper boundary
condition for solute transport, while free–drainage conditions for both water and solute were applied to the lower
boundary of the domain (fig. 4). No–flux boundary conditions were applied to both sides of the flow domain.
Measured furrow cross–sections as shown in figure 1, not thebest fits, were used to define the upper section of the transportdomain for each experimental plot. A finer grid (~0.5 cm)
was used near the soil surface, and a much coarser grid(~4 cm) was specified at the bottom of the domain.
INVERSE OPTIMIZATION
The inverse problem was based on numerical solutions of
the Richards equation (eq. 1) and the CDE or MIM
convection–dispersion model (eqs. 6 and 7) using theLevenberg–Marquardt optimization procedure (Marquardt,
1963). In the inverse procedure, the unknown parameters are
optimized by minimizing the objective function F(q,b)
defined as (Simunek et al., 1999b):

1089

0
ho (water level)
Soil Profile Depth (cm)Atmosphericboundary
condition–100

No–flux boundaries
Free–draining
boundary0

Furrow width (cm)100

Figure 4. Applied water and solute boundary conditions for the specified
soil profiles.

mconcentrations, and the many transport parameters to be
estimated, no guarantee exists that a unique set of Ks andtransport parameters was obtained. In this study, we used both
the CDE and MIM approaches to see whether or not the
two–region transport model would improve the predictions.
Ks and solute transport parameters were estimated simultaneously since this approach takes advantage of crossover
effects between the water flow and solute transport statevariables and parameters (Sun and Yeh, 1990). This method
uses simultaneously all available information and hence
should yield smaller estimation errors than sequential
estimation (Mishra and Parker, 1989; Simunek et al., 2002).
The current version of HYDRUS–2D is capable of simultaneously estimating up to 15 different soil hydraulic, root
water uptake, and/or transport parameters. However, werecommended to rarely, if ever, optimize such a large numberof parameters. The number of parameters that can be
simultaneously optimized depends on the availability ofmeasured data used in the objective function and on theuncertainty in the estimated parameters.

nΦ(q, b) = ∑v j ∑ wij[q*j (x, z, ti ) − q j (x, z, ti , b)]2
j =1 i =1

where
n

(11)

RESULTS AND DISCUSSION= number of observations for the jth measure–ment set (e.g., water contents, infiltrations,concentrations, ...)
qi *(x,z,ti ) = specific measurements at time ti , location x,
and depth z
qi (x,z,ti ,b) = corresponding model predictions obtainedwith the vector of optimized parameters b =(DL , DT, w, qim , ...)
vj and wij = weights associated with a particularmeasurement set or point, respectively.
Weighting coefficients wij were assumed to be equal to 1,
while the vj values were determined according to Clausnitzer
and Hopmans (1995).
Measured soil water contents, infiltration data, andbromide concentrations were used in the inverse optimization to estimate simultaneously Ks and two CDE (i.e., DL andDT ) or four MIM transport parameters (i.e., DL ,DT, w, andqim ). Other parameters (notably the VG soil hydraulicparameters) during the optimization were fixed at theirvalues given in table 1. Each inverse simulation was run atleast three times with different initial estimates to increase
the probability of finding the global minimum of theobjective function. Still, we understand that with the limited
amount of measured data available, particularly bromide

Final optimized values of the saturated hydraulic conductivity and CDE and MIM transport parameters for thedifferent experiments are summarized in table 2. Notice thatthe optimized parameters and sum of squared residuals (SSQ)
for the CDE and MIM transport models are very similar forthe different plots (except for DL in plot 1).
The estimated values for Ks were more or less the same forall experiments, except for plot 5, whose Ks was about two
times larger than for the other experiments. These results are
consistent with field observations indicating that during thesecond irrigation the same amounts of water infiltrated inalmost the same time in plots 3 and 5, which had water levelsof 14 and 10 cm, respectively, whereas a shorter application
time was observed for plot 5 notwithstanding its lower waterlevel of 10 cm. The estimated Ks values are close to thosereported by Carsel and Parrish (1988) and Schaap et al.
(1998) for sandy loam soils. The optimization results
indicated very little correlation between Ks and the otherparameters. Correlation coefficients were generally less thanabout 0.05 for the different plots.
Contrary to recent studies (e.g., Jacques, 2000; Ventrella
et al., 2000), which suggest that the presence of immobile
water is important for modeling field–scale solute transport,
our analysis produced relatively low values for immobile
water content (qim ) (table 2), thus indicating that immobile

Table 2. Summary of the optimized saturated hydraulic conductivity (Ks ) and CDE and MIM transportparameters for the different experiments. Optimized values for the CDE runs are given in parentheses.
Ks
θim
ω
DL
DT
SSQ
Experiment(cm min–1)(cm3 cm–3)(min–1)
(cm)
(cm)
(–)

R2
(–)

SD[a]

Plot 1 (6 cm)
Plot 2 (10 cm)
Plot 3 (14 cm)

0.0394 (0.0389)
0.0396 (0.0392)
0.0497 (0.0497)

0.029
0.010
0.057

19.52
0.10
4.11

32.8 (18.3)
16.8 (17.0)
13.5 (12.5)

0.07 (0.03)
0.81 (0.40)
0.05 (0.04)

1.69 (1.69)
3.48 (4.22)
6.40 (6.32)

0.99 (0.99)
0.99 (0.99)
0.98 (0.98)

SWS[b]

Plot 4 (6 cm)
Plot 5 (10 cm)

0.0456 (0.0436)
0.0957 (0.0996)

0.010
0.020

7.87
11.20

7.9 (8.3)
3.7 (2.6)

0.04 (0.07)
2.10 (2.20)

8.73 (8.93)
3.23 (3.50)

0.99 (0.99)
0.99 (0.99)

43.9 (47.5)

76.8

86.4

74.7 (55.1)

145 (170)––––

CV (%)[a] Same duration (SD) experiments.
[b] Same amount of applied water and solute (SWS) experiments.

1090

TRANSACTIONS OF THE ASAE

water did not play a major role in our study, and hence thatadvection and dispersion were the main transport processes.
This may be due in part to the coarse texture (sandy loam) of
the soil studied. Values for qim for the various plots varied
between 0.010 and 0.057 cm3 cm–3. Unlike qim , thefirst–order exchange coefficient (w) values (table 2) were
found to vary highly between the different experiments(CV = 86.4%), with the estimated values being about one totwo orders of magnitude higher than those reported in theliterature for both laboratory (e.g., Jacobsen et al., 1992;
Mallants et al., 1996) and field studies (Jaynes et al., 1995;
Jacques et al., 2002). The high values reflect rapid soluteexchange between the mobile and immobile regions. Thissuggests that even if some immobile water is present, soluteexchange between the two regions would be so rapid that thetransport process would become macroscopically indistinguishable from a CDE–type process in which advection anddispersion dominate.
The highly transient nature of the 2–D infiltration
processes may well have contributed to the relatively
standard convective–dispersive nature of the transport processes in our study. Other studies in which mobile–immobiletransport processes were investigated and/or MIM parameters were derived usually involved 1–D studies, mostly
during flow conditions at or close to steady–state (Nkedi–
Kizza et al., 1983; Seyfried and Rao, 1987; Jaynes et al.,
1995; Mallants et al., 1996; among others). Additional 2–D
water flow and transport studies under transient infiltration
conditions may be needed to confirm our results; we areunaware of any other studies of this type in the literature. We
note that our optimizations never converged when we usedvery low initial guesses for w (less than 0.001 min–1 for someexperiments, and less than 0.1 min–1 for others). Also, spatial
variability in the soil hydraulic properties made it difficult todraw conclusions about the effects of the water level on qim
and w. Additional measurements, particularly of soluteconcentrations at more frequent time levels, perhaps could
have produced more precise optimized parameters.
The optimized DL values ranged from 2.6 to 32.8 cm,
while DT values varied between 0.03 and 2.20 cm (table 2).
DL values were somewhat larger than those typicallyobtained from 1–D laboratory transport experiments during
steady–state flow. However, they were fairly similar tovalues previously reported by Bowman and Rice (1986) andJaynes et al. (1988) for field experiments at the same field site
involving intermittent and continuous flood irrigation regimes, respectively. For instance, DL values (as estimatedwith the one–dimensional CDE) obtained by Jaynes et al.
(1988) varied between 13.8 to 22.8 cm. Estimated DL and DT
values were comparable with those reported earlier by Forrer
et al. (1999) for a 2–D field–scale study under steady–stateflow conditions. The relatively large optimized values for DL
and DT observed for several plots were likely due to inherent
2–D behavior of the solute transport process and more lateralflow spreading in the furrows. They also reflect the natural
heterogeneity of our field site, spatial variability in soil watercontents, and perhaps our somewhat lower water contents ascompared to the experimental conditions of Forrer et al.
(1999). The relatively low DT values for some of our plots
could also have been compensated somewhat by the large DL
values (correlation coefficients between DL and DT were
found to be about 0.25). Overall, however, the relatively large
confidence intervals for the optimized DT values indicate thatVol. 46(4): 1085–1095

our experiments did not provide enough information toreliably estimate this parameter.
Overall, the obtained DL values for the drier plots(SD experiments) were higher than those for the SWS plots
with the higher soil water contents. This is consistent withtheoretical analyses by Russo (1993, 1998), who showed thatthe dispersivity should increase with decreasing soil watercontents and increasing spatial variability in the soil watercontent. DL also showed some dependency on the water levelin the furrow and the water/solute application time. Higher
values for DL were obtained for the SD experiments. It isdifficult to discern if this difference is caused by differencesin the water level and/or application time, or due tovariability in the soil properties.
The relatively large estimated DT value for plot 5 may
have been caused by increased lateral flow in this plot. As
demonstrated by Abbasi et al. (2003), observed concentrations, particularly in the SWS plots, showed considerable
lateral solute spreading. In addition, DT values showed larger
CVs as compared to the other optimized parameters (table 2).
Previous sensitivity analyses indicate that numerical models
are more sensitive to changes in DT (Forrer et al., 1999) thanin DL (Forrer et al., 1999; Ventrella et al, 2000). Furthermore,both theoretical and experimental investigations have shownthat field–scale dispersivities are several orders of magnitudes larger than laboratory–scale values for the same porous
medium (Gelhar et al., 1992) and increase with the scale ofobservation (Fried, 1975; Gelhar et al., 1992). Observed
longitudinal dispersivities for 59 different soils and aquifers
ranged from 10–2 to 104 m for scales ranging from 10–1 to 105m (Gelhar et al., 1992).
Measured cumulative infiltration rates, soil water contents, and solute concentrations from all five experimental
plots are compared in figure 5 with predicted values ascalculated with the optimized parameters. Measured infiltration data matched the predicted values reasonably well
(R2 was about 0.95), with agreement between measured andpredicted water contents and concentrations being somewhat
less. Water contents were mostly overestimated, and concentrations were underestimated. Overestimation of watercontents could be due to possible overestimation of thelaboratory–measured qs values). Including qs and other soil
hydraulic parameters in the optimization could well have
improved the model predictions. However, as indicatedearlier, we felt it is important to limit as much as possible thenumber of unknown parameters in view of the limited
amount of measured data and the relatively high sensitivityof the parameter estimation technique to the number of
unknowns (Simunek et al., 1999b). Only Ks was included in
the optimization procedure since we did not have direct fieldmeasurements of this parameter. The MIM and CDEtransport models showed more or less the same agreement
between observed and predicted values (fig. 5). Note that theR2 values in figure 5 represent regressions of measured versus
predicted results for all of the data pertaining to the five
optimized plots, whereas the R2 values in table 2 hold forspecific data types (infiltration, soil water content, andconcentration) used in the individual optimizations.
Selected comparisons between measured and model–predicted soil water contents and solute concentrations arepresented in figures 6 and 7, respectively. Plot 1 was selected
for these comparisons because of relatively large differences
between the optimized CDE and MIM transport parameters

1091

Infiltration(L/m)

Concentration(g/L)

Soil water content(cm 3 cm –3 )

200

120
80
40

40

0.30

80 120 160 200

0.20

200

0
40

1
1

3

4

5

6

6
1:1 Line
R 2 = 0.430

0.30

0.20
0.10
0.10

80 120 160 200

2Observed

40

0

2

0

5

1:1 Line
R 2 = 0.544

4
Predicted

80

Predicted

1:1 Line
R 2 = 0.909

120

3

0

0.40

0.40

160
Predicted

0.30

1:1 Line
R 2 = 0.521

4Observed

Observed

MIM model

5

0.20
0.10
0.10

0
0

1:1 Line
R 2 = 0.425
Predicted

1:1 Line
R 2 = 0.897

Predicted

Predicted

CDE model

6

0.40

160

3
2
1
0

0.20

0.30

0.40

0

1

Observed

Observed

2

3

4

5

6

Observed

Figure 5. Measured versus predicted (using the CDE and MIM transport models) cumulative infiltration rates, soil water contents, and solute concentrations for all experimental plots.

for this plot (table 2). Relatively good agreement was foundbetween observed and calibrated soil water contents in theredistribution phase (figs. 6b and 6c), whereas water contents
were overestimated during infiltration (fig. 6a) and also in thedrier conditions (fig. 6d). Results for the second irrigation
event (not shown) were more or less the same as those for the
first irrigation shown in figure 6.
Agreement between measured and predicted solute concentration distributions was somewhat poor near the soil

Soil Profile Depth (cm)Observed(a) Immediately

surface, where the model underestimated concentrations, but
it was much better in deeper layers (fig. 7). In spite ofunderestimating peak concentrations, the model estimatedthe observed solute concentration fronts at different times
very well. We believe that one–dimensional soil watercontent and solute concentration breakthrough curves could
have provided a better and more cogent comparison between
model predictions and field measurements. Also, notwith–standing the large discrepancy between the optimized MIM(b) 6 hours(c) 1 day(d) 5 days

0

0

0

0

– 25

–25

–25

–25

– 50

–50

–50

–50

– 75

–75

–75

–75

– 100
0

25

50

75

–100
100 0

25

50

75

–100
100 0

25

50

75

–100
100 0

25

50

75

100

25

50

75

100

Soil Profile Depth (cm)Calibrated

Furrow Width (cm)
0

0

0

0

–25

–25

–25

–25

–50

–50

–50

–50

–75

–75

–75

–75

–100
0

25

50

75

–100
100 0

25

50

75

–100
100 0

25

50

75

–100
100 0

Furrow Width (cm)Figure 6. Measured and calibrated soil water contents in plot 1 at different times after the first irrigation event.

1092

TRANSACTIONS OF THE ASAE

Soil Profile Depth (cm)Observed(a) 5 days after1st irrigation(c) 20 days after 2nd irrigation(b) 6 days after 2nd irrigation

0

0

0

–25

–25

–25

–50

–50

–50

–75

–75

–75

–100
0

25

50

75

–100
100 0

25

50

75

100

–100
0

25

50

75

100

25

50

75

100

Soil Profile Depth (cm)Calibrated

Furrow Width (cm)
0

0

0

–25

–25

–25

–50

–50

–50

–75

–75

–75

–1 00
0

25

50

75

–100
100 0

25

50

75

–100
100
0

Furrow Width (cm)Figure 7. Measured and calibrated solute concentrations in plot 1 at different soil sampling times (dashed lines = CDE, and continuous lines = MIM).

and CDE dispersivity values for plot 1 (table 2), no obvious
differences were found between predictions with the CDE
and MIM transport models (fig. 7). This again indicates
relatively low sensitivity of solute transport models to thelongitudinal dispersivity, as was previously shown also by
Ventrella et al. (2000) and Forrer et al. (1999) for one– andtwo–dimensional field–scale studies, respectively.

CONCLUSIONSThe HYDRUS–2D code in combination with the Levenberg–Marquardt optimization algorithm was used to inversely estimate the saturated hydraulic conductivity and CDE andMIM transport parameters from several field–scale waterflow and solute transport experiments carried out undertransient conditions. The optimized parameters, particularly
the saturated hydraulic conductivity and the longitudinal andtransverse dispersivities, corresponded well with values
reported in the literature for field–scale studies. Results
indicate that immobile water did not play an important role,
perhaps because of the relatively coarse texture of our fieldsite. Immobile water contents were more or less the same forthe different experiments, giving an average value of only
0.025 cm3 cm–3. The first–order mass transfer coefficient was
higher than those reported for most 1–D analyses on soilcolumns (Seyfried and Rao, 1987; Mallants et al., 1996) andfield studies (Jaynes et al., 1995; Jacques, 2000). This reflects
rapid exchange of solute with any immobile water that ispresent, and hence indicates that the effects of immobile
water were relatively minimal in our study. These results are
consistent with the fact that only minor differences were
found between the CDE and MIM solute concentration
predictions.
Estimated values for the longitudinal dispersivity werelarger than those reported earlier for laboratory soil columns,

Vol. 46(4): 1085–1095

but similar to those measured for field studies (Bowman and
Rice, 1986; Jaynes et al., 1988; Forrer et al., 1999). We
believe that our somewhat higher values are related to having
slightly lower water contents and more irregular flow
patterns than previous studies and also uncertainty inestimating several parameters, simultaneously. Longitudinal
dispersivities showed some dependency on the water level inthe furrows and the water/solute application time, but noclear effects were found of the water level on the transverse
dispersivity, the immobile water content, and the masstransfer coefficient. Agreement between model predictions
and measured infiltration data was generally satisfactory,
while it was relatively poor for soil water contents and soluteconcentrations. Soil water contents were somewhat overestimated during the infiltration and later stages of redistribution
when the soil dried out, while solute concentrations were
underestimated near the soil surface.
Simultaneous estimation of soil hydraulic and transportparameters has the advantage that this can be done in a single
step. The method, however, may sometimes lead to inaccurate transport parameters since the objective function isdefined in terms of both hydraulic and transport data. Large
hydraulic residuals could then compensate (or overshadow)
small transport residuals, and vice versa. The final outcomeof the optimization process then depends very much on what
weights are given to the particular data sets.
The present simulations were carried out assuming ahomogeneous soil profile. Accounting for soil layering andimplementing a scaling method to describe spatial variability
in the soil hydraulic and transport properties and consequently estimating the invoked parameters with the classicalsequential estimation of hydraulic parameters followed bysolute transport parameters will further improve the predictions. One concern we had during analysis of our data was therequired computational time for the 2–D water flow andsolute transport calculations, especially during inverse1093

analyses of the heterogeneous case. Depending on the initial
estimates of the unknown parameters, optimization required
up to several days on a Pentium III PC. In addition, accurate
evaluation of soil hydraulic and solute transport properties
from field–scale data, which are subject to considerable
spatial and temporal variability, requires extensive datacollection. For these reasons, it is important to measure as
many parameters as possible independently, especially when
limited or highly variable field data are available.
ACKNOWLEDGEMENTS
The authors are grateful to Jack Jobes of the George E.
Brown Jr. Salinity Laboratory, Riverside, California, for his
help in analyzing soil samples for soil water retention curves
and saturated hydraulic conductivity. This study was partially supported by SAHRA (Sustainability of Semi–Arid
Hydrology and Riparian Areas) under the STC Program of
the National Science Foundation, Agreement No.
EAR–9876800.

REFERENCES

Abbasi, F., F. J. Adamsen, D. J. Hunsaker, J. Feyen, P. Shouse, and
M. Th. van Genuchten. 2003. Effects of flow depth on waterflow and solute transport in furrow irrigation: Field data
analysis. J. Irrig. and Drainage Eng. 129(4): 237–246.
Addiscott, T. M., and R. J. Wagenet. 1985. Concepts of soluteleaching in soils: A review of modeling approaches. European J.
Soil Science 36(3): 411–424.
Allen, R. G., M. Smith, A. Perrier, and L. S. Pereira. 1994. An
update for the definition of reference evapotranspiration. ICID
Bulletin 43(2): 1–34.
Biggar, J. W., and D. R. Nielsen. 1976. Spatial variability of theleaching characteristics of a field soil. Water Resources Research
12(1): 78–84.
Bowman, R. S., and R. C. Rice. 1986. Transport of conservative
tracers in the field under intermittent flood irrigation. Water
Resources Research 22(11): 1531–1536.
Bresler, E., G. Dagan, R. J. Wagenet, and A. Laufer. 1984.
Statistical analysis of salinity and texture effects on spatial
variability of soil hydraulic conductivity. SSSA J. 48(1): 16–25.
Carsel, R. F., and R. S. Parrish. 1988. Developing joint probability
distributions of soil water retention characteristics. Water
Resources Research 24(5): 755–769.
Clausnitzer, V., and J. W. Hopmans. 1995. Non–linear parameter
estimation: LM–OPT. General–purpose optimization code basedon the Levenberg–Marquardt algorithm. Land, Air, and Water
Resources Paper No. 100032. Davis, Cal.: University ofCalifornia.
Clothier, B. E., I. Vogeler, S. R. Green, and D. R. Scotter. 1998.
Transport in unsaturated soil: Aggregates, macropores, andexchange. In Physical Nonequilibrium in Soils: Modeling andApplication, 273–295. H. M. Selim and L. Ma, eds. Chelsea,
Mich.: Ann Arbor Press.
Forrer, I., R. Kasteel, M. Flury, and H. Fluhler. 1999. Longitudinaland lateral dispersion in an unsaturated field soil. Water
Resources Research 35(10): 3049–3060.
Forsyth, P. A., Y. S. Wu, and K. Pruess. 1995. Robust numericalmethods for saturated–unsaturated flow with dry initialconditions. Advances in Water Resources 18(1): 25–38.
Fried, J. J. 1975. Groundwater Pollution. New York, N.Y.: Elsevier.
Gelhar, L. W., C. Welty, and K. R. Rehfeldt. 1992. A critical review
of data on field–scale dispersion in aquifers. Water Resources
Research 28(7): 1955–1974.
Gribb, M. M. 1996. Parameter estimation for determining hydraulicproperties of a fine sand from transient flow measurements.
Water Resources Research 32(7): 1965–1974.

1094

Jacobsen, O. H., F. J. Leij, and M. Th. van Genuchten. 1992.
Parameter determination for chloride and tritium transport inundisturbed lysimeters during steady flow. Nordic Hydrology
23(2): 89–104.
Jacques, D. 2000. Analysis of water flow and solute transport at thefield scale. PhD diss. No. 454. Leuven, Belgium: Katholieke
Universiteit Leuven, Fac. Landbouwkundige en Toegepaste
Biologische Wetenschappen.
Jacques, D., J. Simunek, A. Timmerman, and J. Feyen. 2002.
Calibration of Richards’ and convection–dispersion equations tofield–scale water flow and solute transport under rainfall
conditions. J. Hydrology 259: 15–31.
Jarvis, N. J., P. E. Jonsson, P. E. Dik, and I. Messing. 1991.
Modeling water and solute transport in macroporous soil: I.
Model description and sensitivity analysis. European J. Soil
Science 42(1): 59–70.
Jaynes, D. B., R. S. Bowman, and R. C. Rice. 1988. Transport of aconservative tracer in the field under continuous flood irrigation.
SSSA J. 52(3): 618–624.
Jaynes, D. B., S. D. Logsdon, and R. Horton. 1995. Field method
for measuring mobile/immobile water content and solute transfer
rate coefficient. SSSA J. 59(2): 352–356.
Kees, C. E., and C. T. Miller. 2001. Higher order integration
methods for two–phase flow. Advances in Water Resources25(2): 159–177.
Kodesova, R., M. M. Gribb, and J. Simunek. 1998. Estimating soil
hydraulic properties from transient cone permeameter data. Soil
Science 163(6): 436–453.
Kodesova, R., S. E. Ordway, M. M. Gribb, and J. Simunek. 1999.
Estimating of soil hydraulic properties with cone permeameter:
Field studies. Soil Science 164(8): 527–541.
Loretta, Y., and G. Wu. 1999. Numerical simulation of four heavy
metals in kaolinite clay. J. Environ. Eng. 125(4): 314–324.
Mallants, D., M. Vanclooster, and J. Feyen. 1996. Transect study on
solute transport in a macroporous soil. Hydrological Processes
10(1): 55–70.
Mallants, D., B. P. Mohanty, A. Vervoort, and J. Feyen. 1997.
Spatial analysis of saturated hydraulic conductivity in a soil withmacropores. Soil Technology 10(2): 115–131.
Marquardt, D. W. 1963. An algorithm for least squares estimation of
non–linear parameters. J. Ind. Appl. Math. 11: 431–441.
Mishra, S., and J. C. Parker. 1989. Parameter estimation for coupledunsaturated flow and transport. Water Resources Research 25(3):385–396.
Mohanty, B. P., R. S. Bowman, M. H. Hendrickx, and M. Th. van
Genuchten. 1997. New piecewise–continuous hydraulicfunctions for modeling preferential flow in anintermittent–flood–irrigated field. Water Resources Research
33(9): 2049–2063.
Mohanty, B. P., R. S. Bowman, M. H. Hendrickx, J.Simunek, and
M. Th. van Genuchten. 1998. Preferential transport of nitrate toa tile drain in an intermittent–flood–irrigated field: Modeldevelopment and experimental evaluation. Water Resources
Research 34(5): 1061–1067.
Mualem, Y. 1976. A new model for predicting the hydraulicconductivity of unsaturated porous media. Water Resources
Research 12(3): 513–522.
Neuman, S. P., R. A. Feddes, and E. Bresler. 1975. Finite element
analysis of two–dimensional flow in soils considering water
uptake by roots: I. Theory. SSSA Proc. 39(2): 224–230.
Nkedi–Kizza, P., J. W. Biggar, M. Th. van Genuchten, P. J.
Wierenga, H. M. Selim, J. M. Davidson, and D. R. Nielsen.
1983. Modeling tritium and chloride36 transport through anaggregated oxisol. Water Resources Research 19(3): 691–700.
Richards, L. A. 1931. Capillary conduction of liquids throughporous medium. Physics 1: 318–333.
Russo, D. 1993. Stochastic modeling of macrodispersion for solutetransport in a heterogeneous unsaturated porous formation.
Water Resources Research 29(2): 383–397.

TRANSACTIONS OF THE ASAE

Russo, D. 1998. Stochastic analysis of flow and transport inunsaturated heterogeneous porous formations: Effects ofvariability in water saturation. Water Resources Research 34(4):569–581.
Russo, D., J. Zaidel, and A. Laufer. 1998. Numerical analysis offlow and transport in a three–dimensional partially saturated
heterogeneous soil. Water Resources Research 34(6):1451–1468.
Schaap, M. G., F. J. Leij, and M. Th. van Genuchten. 1998. Neuralnetwork analysis for hierarchical prediction of soil hydraulicproperties. SSSA J. 62(4): 847–855.
Schaap, M. G., F. J. Leij, and M. Th. van Genuchten. 2001.
ROSETTA: A computer program for estimating soil hydraulicparameters with hierarchical pedotransfer functions. J.
Hydrology 251: 163–176.
Segol, G. 1994. Classic Groundwater Simulations: Proving and
Improving Numerical Models. Englewood Cliffs, N.J.:
Prentice–Hall.
Seyfried, M. S., and P. S. C. Rao. 1987. Solute transport inundisturbed columns of an aggregated tropical soil: Preferential
flow effects. SSSA J. 51(6): 1434–1444.
Shao, M., R. Horton, and R. K. Miller. 1998. An approximate
solution to the convection–dispersion equation of solutetransport in soil. Soil Science 163(5): 339–345.
Simunek J., and M. Th. van Genuchten. 1994. The CHAIN–2D
code for simulating the two–dimensional movement of water,
heat, and multiple solutes in variably saturated porous media.
Research Report 136. Riverside, Cal.: USDA–ARS, George E.
Brown Jr. Salinity Laboratory.
Simunek, J., and M. Th. van Genuchten. 1996. Estimating
unsaturated soil hydraulic properties from tension disc
infiltrometer data by numerical inversion. Water Resources
Research 32(9): 2683–2696.
Simunek, J., and M. Th. van Genuchten. 1997. Estimating
unsaturated soil hydraulic properties from multiple tension disc
infiltrometer data. Soil Science 162(6): 383–398.
Simunek, J., M. Sejna, and M. Th. van Genuchten. 1998. The
HYDRUS–1D software package for simulating theone–dimensional movement of water, heat, and multiple solutes
in variably saturated media, Version 2.0. IGWMC–TPS–70.
Golden Colo.: Colorado School of Mines, International Ground
Water Modeling Center.
Simunek, J., M. Sejna, and M. Th. van Genuchten. 1999a. The
HYDRUS–2D software package for simulating thetwo–dimensional movement of water, heat, and multiple solutes
in variably saturated media, Version 2.0, IGWMC–TPS–70.
Golden Colo.: Colorado School of Mines, International Ground
Water Modeling Center.

Vol. 46(4): 1085–1095

Simunek, J., R. Kodesova, M. M. Gribb, and M. Th. van
Genuchten. 1999b. Estimating hysteresis in the soil water
retention function from cone permeameter experiments. Water
Resources Research 35(5): 1329–1345.
Simunek, J., D. Jacques, J. W. Hopmans, M. Inoue, M. Flury, and
M. Th. van Genuchten. 2002. Solute transport during
variably–saturated flow – Inverse methods. Chapter 6.6 in
Methods of Soil Analysis: Part 1. Physical Methods,
1435–1449. 3rd ed. J. H. Dane and G. C. Topp, eds. Madison,
Wisc.: SSSA.
Sulekha, D. K., and W. H. M. Duijnisveld. 1998. Two–phasereactive solute transport simulation in unsaturated zone. J. Irrig.
Drainage Eng. 124(3): 140–147.
Sun, N. Z., and W. W. G. Yeh. 1990. Coupled inverse problems ingroundwater modeling: 1. Sensitivity analysis and parameter
identification. Water Resources Research 26(10): 2507–2525.
van Genuchten, M. Th. 1980. A closed–form equation forpredicting the hydraulic conductivity of unsaturated soils. SSSA
J. 44(5): 892–898.
van Genuchten, M. Th., and E. A. Sudicky. 1999. Recent advances
in vadose zone flow and transport modeling. In Vadose Zone
Hydrology: Cutting Across Disciplines, 155–193. M. B.
Parlange and J. W. Hopmans, eds. New York, N.Y.: Oxford
University Press.
van Genuchten, M. Th., and R. J. Wagenet. 1989.
Two?site/two?region models for pesticide transport anddegradation: Theoretical development and analytical solutions.
SSSA J. 53(5): 1303?1310.
van Genuchten, M. Th., F. J. Leij, and S. R. Yates. 1991. The RETC
code for quantifying the hydraulic functions of unsaturated soils.
Report No. EPA/600/2–91/065. Ada, Okla.: U.S. Environmental
Protection Agency, R. S. Kerr Environmental ResearchLaboratory.
Ventrella, D., B. P. Mohanty, J. Simunek, N. Losavio, and M. Th.
van Genuchten. 2000. Water and chloride transport in a
fine–textured soil: Field experiments and modeling. Soil Science
165(8): 624–631.
Wang, D., S. R. Yates, J. Simunek, and M. Th. van Genuchten.
1997. Solute transport in simulated conductivity fields under
different irrigations. J. Irrig. Drainage Eng. 123(5): 336–343.
Wu, G., and S. T. Chieng. 1995a. Modeling multicomponent
reaction chemical transport in non–isothermal
unsaturated/saturated soils: Part 1. Mathematical model
development. Trans. ASAE 38(3): 817–826.
Wu, G., and S. T. Chieng. 1995b. Modeling multicomponent
reaction chemical transport in non–isothermal
unsaturated/saturated soils: Part 2. Numerical simulation. Trans.
ASAE 38(3): 827–838.

1095

1096

TRANSACTIONS OF THE ASAE