Journal of Hydrology 313 (2005) 257–273
www.elsevier.com/locate/jhydrol

Single-porosity and dual-porosity modeling of water flow andsolute transport in subsurface-drained fields using effective
field-scale parameters
Nathan W. Hawsa,*, P. Suresh C. Raoa,b, Jirka Simunekc, Irene C. PoyeraaSchool of Civil Engineering, Purdue University, West Lafayette, IN 47907-1150, USAb
Department of Agronomy, Purdue University, West Lafayette, IN 47907-2054, USAc
Department of Environmental Sciences, University of California, Riverside, CA 92521, USA
Received 28 July 2004; revised 14 March 2005; accepted 16 March 2005

Abstract
Water and solute flux from the subsurface drains of macroporous agricultural fields are simulated using two-dimensional
single-porosity and dual-porosity models. Field-averaged (i.e. effective) parameters are calibrated from drainage outflow andvalidated for water flow and the transport of non-reactive solutes applied at discrete locations on the field. Both the singleporosity and the dual-porosity simulations capture the observed trends in the drainage hydrographs, with the dual-porosity
model performing slightly better than the single-porosity model. The values of the effective hydrologic parameters, however,were not fully characteristic of macroporous soils. The physical meaning of the effective parameters was further questioned asneither the single-porosity nor the dual-porosity models could simulate the rapid transport of solutes to the subsurface drain.
The discrepancies between the simulated and observed solute flux indicate that the actual spatial area contributing to drainage
outflow in the field experiment was much greater than the integrated area of the simulated domain over which the effective
parameters were calibrated. Supplementary simulations using parameters calibrated from solute flux data (outflow soluteconcentration multiplied by the outflow water flux) also fail to match both water and solute fluxes. The failure of the simulationsis attributable to factors such as non-unique parameters and problems with representing a three-dimensional heterogeneous
domain as a two-dimensional homogeneous system, including a misrepresentation of macropore flow paths. This study shows
the fallacy of interpreting a hydrograph fit as evidence of the physical meaning of model parameters.
q 2005 Elsevier B.V. All rights reserved.
Keywords: Dual-porosity; HYDRUS-2D; Effective parameters; Subsurface-drains; Model calibration

1. Introduction* Corresponding author. Address: Department of Geography andEnvironmental Engineering, Johns Hopkins University, Baltimore,
MD 21218-2686, USA. Tel.: C1 410 516 4255; fax: C1 410 516
8996.
E-mail address: nhaws@jhu.edu (N.W. Haws).

0022-1694/$ - see front matter q 2005 Elsevier B.V. All rights reserved.
doi:10.1016/j.jhydrol.2005.03.035

Agricultural landscapes consist of multi-scale
heterogeneities that influence the movement ofwater and chemicals. Usually, a detailed quantification of spatial variability is either impractical or258

N.W. Haws et al. / Journal of Hydrology 313 (2005) 257–273

impossible. Consequently, models for predicting flow
and transport through the vadose zone must somehow
simplify the representation of the transport domain,
yet still effectively reproduce the flow and transport
response at a control plane or outlet point. In a
deterministic model, this is often accomplished byassuming a representative elementary volume (REV),
which proposes that a medium can be represented byhomogenized (i.e. effective) parameters provided
there is a separation of scales such that thecharacteristic volume of spatial heterogeneities ismuch smaller than macroscopic size of the transportdomain (Royer et al., 2002).
The REV assumption is commonly invoked tocalibrate effective model parameters from catchment
outflow hydrographs for rainfall-runoff models. Often
the legitimacy of these models is validated solely on
the basis of hydrograph prediction (Beven, 1993;
Beran, 1999). Using the outflow response as an
exclusive measure of performance, however, is being
increasingly questioned as these models fail torepresent the actual drainage pathways and thephysical processes occurring within the flow domain(Binley et al., 1991; Jain et al., 1992; Blazkova et al.,2002). Consequently, models become more sophisticated to better account for the internal system
structure. The increased sophistication results inincreased parameterization requirements. Because it
is typically infeasible to make detailed distributed
measurements, the added parameters must be estimated based on the same hydrograph information,
which leads to a higher likelihood of non-uniqueparameter sets and an ill-posed inverse calibration
problem—even for simple, small-scale systems(Hopmans and Šimůnek, 1997; Durner et al., 1997;
Madsen et al., 2002; Doherty and Johnston, 2003).
While the REV approach to finding effective
parameters is simple in concept, verifying its
existence can be complicated. (Baveye and Sposito,
1984; Berkowitz et al., 1988; Neuman and Orr,
1993; Lugo et al., 1998; Tartakovsky and Neuman,
1998a,b,c; Indelman, 2002; Fernandez-Garcia et al.,
2002; Royer et al., 2002). Neuman and Orr (1993)
demonstrated mathematically that while effective
values may be found for infinite subsurface domains,
the flux for bounded systems becomes non-local, andthus, an average flow model generally does not exist.
For radially converging flow systems, however, they

showed that effective parameters could be founddepending on the correlation structure and physicaldimensions of the domain.
Temporal variability also creates difficulties infinding effective parameters. Tartakovsky and
Neuman (1998b) reported that for porous media, aneffective hydraulic conductivity will exist in the ‘strict
sense’ only when the mean head and the residual flux
are constant in space and time. Though this is notpossible under transient flow—except if the storativity
is zero or time approaches infinity—there are anumber of cases, where effective hydraulic conductivities can be approximated in real, Laplace and/or
Fourier spaces (Tartakovsky and Neuman, 1998b,c).
Some of the more complex flow and transportdomains are subsurface-drained agricultural fields of
the US Midwest region. Analogous to catchmenthydrographs, the subsurface drain response is anaggregation of the spatial and temporal variability
across that part of the field that contributes to drainage
outflow. Consequently, the drainage hydrograph may
be useful for estimating effective parameters if thetransport domain can be approximated as a radiallyconverging flow system, spatial correlation lengths ofvariability in the soil properties are small compared to
the field dimensions, and rainfall events are limited toa realistic range of magnitudes such that the field
volume contributing to the outflow does not appreciably change between events. Still, even when these
conditions are met, the non-linear processes in thevadose zone, coupled with flow and transport throughmacropore networks that are abundant in these soils,
may violate the ergodic hypothesis (Yeh, 1997) andthe representations of the driving processes.
As in catchment modeling, using the outlet
response to derive effective hydraulic properties for
subsurface-drained fields is not new (Hoffman and
Schwab, 1964; Skaggs, 1976), and automated inverse
procedures have been developed and tested for avariety of systems of different scales and for a numberof boundary conditions (Durner et al., 1997; Hopmans
and Šimůnek, 1997; Jacques et al., 1997; Šimůnek
et al., 1998). De Vos et al. (1997) used the inverse
procedure in the HYDRUS-2D flow and transportmodel to estimate effective field-scale parametersfrom measurements of drainage outflow and groundwater level fluctuations. The effective parametermodel could satisfactorily reproduce the overall

N.W. Haws et al. / Journal of Hydrology 313 (2005) 257–273

water balance and the observed water fluxes;
however, it failed to predict the rapid arrival at thedrain outlet of a bromide tracer that was surface
applied to the field. De Vos et al. (1997) attributed themodel’s failure in the solute simulations to the soil
macropores that created a preferential transportdomain that was not well described by a single
saturated hydraulic conductivity value. They further
surmised that a more complex, dual-domain model
would better simulate solute transport. More recently,
Gerke and Köhne (2004) modeled the one-dimensional leaching of bromide in a subsurface-drained
field using both the single-porosity HYDRUS-1Dmodel and a one-dimensional dual-porosity/dual
permeability (i.e. two mobile domains) model.
Again, they found that both models could satisfactorily explain water flow, yet only the dual-permeability model could reproduce the bromide
concentrations in the outflow.
With the supposition that the subsurface-drain
integrates field-scale heterogeneity, this study follows
the work of De Vos et al. (1997) and Gerke and Köhne(2004) by using a version of the HYDRUS-2D code
(Šimůnek et al., 1999) that was modified to include adual-porosity/single-permeability (i.e. mobile–immobile) flow and transport domain (Šimůnek
et al., 2003). Furthermore, the solute simulations are
conducted for tracers applied on discrete locations of
the field, so that the solute concentrations in thedrainage signal a two-dimensional field location andtransport pathway. Effective parameters are calibrated
for both single- and dual-porosity models based on
subsurface-drainage outflow data. The effective
parameter models are then evaluated for both waterand solute flux. The objectives of these modelingstudies is to evaluate effective hydraulic parametersthat are calibrated based on drainage data, to assess
the benefits of increasing model complexity (singlevs. dual-porosity), and to test the physical representativeness of parameters by their ability to simulate
solute flux and reproduce the contributing drainage
areas the flow domain.

2. Field experiments
Replicate multi-solute experiments were conducted on two subsurface-drained cornfield plots

259(labeled Plots 51 and 52) at the Purdue Agronomy
Water Quality Field Station (WQFS), West Lafayette,
Indiana. The soil at the field station is classified asDrummer silty clay loam (fine-silty, mixed mesic
Typic Endoquoll), has a primarily sub-angular blocky
structure in the upper 2 m profile, and an A horizon
that typically extends to at least 40 cm (USDA, 1998).
Although the soil is classified as poorly drained,
macropores that are formed by root channels, earthworm burrows, and other bio-processes can rapidly
convey water and solute from the soil surface to thedrain outlet.
Plots 51 and 52 each measure 48.5!60 m2(L!W), with crop rows running along the 48.5 m
length, parallel to the subsurface drain. Each plot isalso drained by two perimeter drains and a central
drain that bisects the perimeter drains, giving a 30 mhorizontal drain spacing. The subsurface drains areapproximately 1 m below the ground surface, and a
relatively impermeable glacial till lies at about 2 m
below ground. The central subsurface drain of eachplot routes effluent to collection huts equipped withtipping bucket flow meters and automatic water
sampling devices. Bucket tips were recorded withdata loggers (Campbell Scientific, Model CR510),
triggered via a magnetic sensor switch located on thebottom of the tipping bucket. Drainage effluent was
sampled with automatic samplers (ISCO, Model
3700) that were programmed to draw a 300 ml
water sample every 100 tips of the bucket meter(approximately 130 l depending on the specifictipping bucket calibration). The flow-weighted
sampling provided a greater sampling intensity duringstorm peaks and reduced sampling frequencies during
the hydrograph recessions. The 300 ml samples were
sub-sampled in 40 ml volatile organic vials, andstored at 4 8C until analyzed.
The experiments were conducted under naturalrainfall and atmospheric conditions. Rainfall andevaporation values, used as input boundary fluxes inthe HYDRUS-2D model, were obtained frommeasurement locations at or near the field site. Hourly
rainfall amounts were recorded with two rain gauges
located at the field site. Daily evaporation values were
measured from a Class A evaporation pan locatedabout 0.8 km south of the field plots.
A water budget was computed to verify that themeasured rainfall, evaporation, and drainage outflow

260

N.W. Haws et al. / Journal of Hydrology 313 (2005) 257–273

30
20

Evaporation
Plot 52 Outflow
Plot 51 Outflow

100Plot 52 Balance
Plot 51 Balance–10–20
Rainfall

2
/0
14
6/

2

02
4/
6/

/0

2
25

/0
15

5/

02
5/

5/
5/

4/

25

/0

/0
15
4/

5/
4/

2

2

–30
02

Cumulative Outward Flux (cm)accounted for all the water fluxes to and from theplots. For the water budget, a constant pan factor of
0.7 was used to adjust the potential evaporation values
measured in the Class A pan to actual fieldevaporation values. It should be noted that in reality,
pan factors are not constant, but can vary from about
0.35, for dry conditions, to about 0.85, for very
wet/humid conditions, and the factor of 0.7 is anaverage value suggested by Brouwer and Heibloem(1986). Another complication with computing thewater budget is that no water content data wereavailable so that there could be no direct accounting
for the change in water storage within the flow
domain. Consequently, the water budget computations began at the cessation of drainage so that theinitial soil condition could be assumed to be athydrostatic equilibrium, and only boundary fluxes(rainfall, evaporation, and drainage) were consideredin the water budget calculations. The cumulative
water budget is shown in Fig. 1. In this figure, water
fluxes entering the domain (rainfall) are taken as
negative values, and water fluxes leaving the domain(outflow and evaporation) are positive. Theoretically,a budget value of zero at the end of a drainage event(i.e. when the soil again reaches hydrostatic equilibrium) means a perfect balance between inflow andoutflow fluxes. In general, a good water balance isachieved before and after mid-May. The larger
negative errors during mid-May indicate a likely

Date (MM/DD/YY)

Fig. 1. Cumulative water fluxes for rainfall, evaporation, anddrainage for Plots 51 and 52. Also shown is the water balance for
Plots 51 and 52, computed as rainfallCevaporationCoutflow.
Water flux out of the domain is taken as positive.

underestimation of evaporation (or underestimation of
the pan factor) from about May 10–25. (Note, that any
underestimate cannot be attributed to an additional
demand from cropwater uptake since corn was notplanted until June.) Some of the imbalance in thebudget also may result from changes in water storage
that were not explicitly considered.
The field experiments began in April, 2002. Prior
to the solute application, potassium chloride andpotassium bromide were dissolved individually inwater, and stored separately. On April 6, 2002, theconcentrated solutions were surface applied toseparate 1 m wide ‘strips’ that ran parallel to thesubsurface drain. The chloride solution was applied to
a 1 m wide strip centered at a lateral offset of 0.5 m
from the drain. The bromide solution was applied to a
1 m wide strip centered 5 m from the drain. The strip
widths were chosen based on two criteria: (1) to bethin enough so that solute outlet response would
signal a distinct field origin, and (2) wide enough toallow a convenient and accurate solute application.
Although both Plots 51 and 52 have drain half-spacingof 15 m, the outer strip was centered only 5 m fromthe subsurface drain in order to negate any boundary
effects and prevent the chance of solute capture in theperimeter drains. Several spray passes were madealong the length of each strip to ensure uniform tracer
application. For each application, the spray nozzle
was kept close to the soil surface to minimize drift ofthe spray droplets. To prevent cross-contamination,different backpack sprayers were used for each tracer
mixture. Water flux and solute concentration in thedrain outflow were monitored through the end of
2002; however, only data obtain prior to the planting
date of June 3, 2002 were used in the modeling study.
Chloride and bromide concentrations in the drain
effluent samples were analyzed utilizing a Dionex
DX600 Ion Chromatography system, equipped withan ED50 electrochemical detector, an AG40 potassium hydroxide eluent generator, a GP50 pump, and an
AS50 autosampler. Analytes were separated using a
Dionex IonPac 4!50 mm guard column in tandemwith an IonPac AS17 4!250 mm analytical column
and a potassium hydroxide eluent gradient ramped
from 12 to 25 mM. Background anion suppression
was maintained to less than 2.0 mS (micro-Siemens)with an ASRS Ultra 4 mm suppressor. The method
detection limit was 0.25 mg/l.

N.W. Haws et al. / Journal of Hydrology 313 (2005) 257–273

3. Flow and transport modeling and parameter
estimation
Flow and transport simulations were conductedwith the HYDRUS-2D finite-element code (Šimůnek
et al., 1999), modified to include a dual-porosity/
single-permeability module (Šimůnek et al., 2003).
The finite-element grid was constructed to represent a
two-dimensional, vertical section of the WQFS plots(Fig. 2). The modeled field domain was 15!2 m2(width!depth), assuming symmetry about the vertical plane of the drain and a no-flux boundary at the
midplane between the central and perimeter drains.
For the water flow simulations, the field was
represented with 60!43 finite element cells with acell spacing that tapered toward the drainage outlet
located at zZ100 cm at the right end of the domain.
This mesh configuration was sufficiently refined togive mass balance errors of less than one percent forwater flux. The sides and bottom of the domain were
specified as no-flux boundaries, and the top of thedomain had an atmospheric boundary. The subsurface
drain was represented as a nodal sink (Šimůnek et al.,1999) with an effective drain diameter, de, of 1 cm anda reduction of hydraulic conductivities of the square in
the finite element mesh surrounding the drain using
the correction factor, Cd, of four (Mohammad andSkaggs, 1983; Fipps et al., 1986).
In the single-porosity model, the entire flow
domain conducts water according to Richards’
equation and transports solute according to theadvection–dispersion equation.
Water:vq
Z V½KðhÞVðh K zÞ
vt

(1a)

q Z KKðhÞVðh K zÞ

(1b)KðhÞ ZKs

qKqr
qs Kqr

1 (
2

261





qKqr
1K 1K
qs Kqr)2
 n  nK1
n
nK1

(1c)
qZ

qs K qr
C qr½1 C jahjn 1Kð1=nÞ
q Z qsfor h! 0(1d)

for hO 0

Solute:vðqCÞ Z VðDVC K qCÞ
vt

(2a)

q
q

(2b)

D Z D0 C l

where q[L3/L3] is the water content (subscripts s and rdenote saturated and residual, respectively), K [L/T] isthe hydraulic conductivity, h [L] is the pressure head,a [1/L] and n [K] are the van Genuchten (1980)
pressure-saturation parameters, and t [T] and z [L] aretemporal and spatial coordinates (z positive downward), respectively. Also, C [M/L3] is the aqueous
solute concentration, q [L/T] is the water flux, and D is
the hydrodynamic dispersion coefficient [L2/T],
which is defined by the dispersivity, l [L], and the
diffusion coefficient, D0 [L2/T].
The dual-porosity model divides the flow domain
into mobile (subscript m) and immobile (subscript im)regions with first-order advection and diffusion
between the two domains (Šimůnek et al., 2003).
Water:
vqm
Z V Kðhm ÞV hm K z K Gw
vt

(3a)

vqim
Z Gw
vt

(3b)

q Z KKðhm ÞVðhm K zÞ

(3c)

Gw Z uw ðhÞðhm K him Þ(3d)Solute:Fig. 2. The HYDRUS-2D finite element mesh used in the water flow
simulations.

v
ðq C Þ Z VðDm VCm K qm Cm Þ K Gs
vt m m

(4a)

qm
qm

(4b)

Dm Z Dm0 C lm

262

N.W. Haws et al. / Journal of Hydrology 313 (2005) 257–273

v
qim Cim Z Gs
vt

(4c)

Gs Z Gw C C us ðCm K Cim Þ(4d)where Gw[1/T] is the water transfer term, Gs [M/L3T]
is the solute transfer term and C*ZCm for GwO0 andCim for Gw!0. The solute transfer rate coefficient, us[1/T], is a constant that accounts for the diffusive
solute transport due to concentration gradients. The
water transfer rate coefficient, uw [1/LT], is defined asfollows (Gerke and van Genuchten, 1993; Šimůnek
et al., 2003):
uw ðhÞ Zb
Ka ðhÞgw
d2(5)

where b [K] is a geometric shape factor, d [L] is thecharacteristic length of the soil matrix, gw [K] is ascaling coefficient, and Ka is the effective hydraulicconductivity of the fracture–matrix interface, calculated using Eq. (1c) with immobile zone retention
parameters. Ka is commonly defined as an arithmetic
average involving both hm and him as follows (Gerkeand van Genuchten, 1996):
Ka ðhÞ Z 0:5½Km ðhm Þ C Kim ðhim Þ(6)

In HYDRUS-2D, uw is a lumped coefficient thatcombines all the right-hand side terms in Eq. (5). The
input value for uw is based on Ka at soil saturation,and the actual value of uw varies with pressure/saturation level according to Eq. (3d).
Five parameters are needed to describe the variably
saturated water flow in the single-porosity domain:
saturated water content (qs), residual water content(qr), the van Genuchten (1980) soil–water characteristic parameters (a, n), and the saturated hydraulicconductivity (Ks). Water flow in the dual-porosity
model requires the same parameters as in the singleporosity model in the mobile region plus five otherparameters for the immobile region: qsim, qrim, aim,
nim, and uw. The additional parameters of the dualporosity model, while complicating the inverse
procedure, were expected to better represent themacropore flow processes than the single-porosity
model.
To find and test the effective hydraulic parameters,
the field data were parsed into two sections. The firstsection, used for calibration, began on April 3 at

0:00 h for Plot 51 and April 6 at 12:00 h for Plot 52
and ended on April 23 at 19:00 h for both plots—a
total of 500 and 416 h for each plot, respectively. The
beginning dates and times were selected to correspondto the completion of a previous drainage event foreach of the plots to justify the initial condition of ahorizontal water table located at the drain invert andhydrostatic equilibrium within the domain. The rainfall–evaporation–outflow data used for parameter
validation began at the ending date and time of thecalibration sets and continued until just prior to theplanting date on June 3, 2002 (results are only shown
up to May 24, 2002 since no drainage events occurred
after this date). Ideally, a longer data set would havebeen chosen for model calibration; however, experimental data for water and solute were simultaneously
available only for the 2002 season. It is possible that
performing the calibration over a longer period would
lead to different parameter values; nevertheless, thetwo rainfall events in the calibration period and thedrying interlude between the events provide a goodrange of transient conditions for estimating the model
parameters.
Parameters were calibrated using the inverse
module within the HYDRUS-2D software package.
HYDRUS-2D employs the Marquart–Levenberg
optimization algorithm to aid in finding best-fit
parameters. The objective of this routine is tominimize the sum of the squared deviations (SSD)
between the observed and simulated drainage outflow
values. In addition to using the minimum SSD
objective, the parameter sets were also evaluatedwith the correlation coefficient
rx;y Z½COVðx; yÞ½sx ½sy (7)where COV is the covariance functions, s is thestandard deviation, and x and y are the observed andsimulated outflows, respectively. The additional useof the correlation coefficient was helpful as anindication of how well the model predictions matched
the timing and general shape of the observed waterand solute outflow fluxes while the SSD was used as ameasure of the distance between simulated andobserved outflow fluxes.
In the inverse procedure, all parameters were
calibrated except the residual water content forN.W. Haws et al. / Journal of Hydrology 313 (2005) 257–273

263

Table 1
Range of parameter values for (a) single-porosity and (b) dual-porosity models
Parameter
(A)
qr
qs
a
n
Ks
Dp
lL, lT(B)
qrm (mobile)qsm (mobile)am(mobile)
nm(mobile)Ks
qrim (immobile)qsim (immobile)aim (imsmobile)nim (immobile)
uw
us
Dmo
lmL, lmTab

Value or range

Definition

0.07–0.2 (cm3/cm3)a
0.3–0.6 (cm3/cm3)a0.003–0.4 (1/cm)
1.01–6 (K)1–2000 (cm/h)
1 (cm2/h)
15, 1.5 (cm)Residual water content
Saturated water content
Water retention parameter (van Genuchten, 1980)
Water retention parameter (van Genuchten, 1980)
Saturated hydraulic conductivity
Diffusion coefficientLongitudinal, tranverse dispersivity

00.01–0.15 (cm3/cm3)b0.003–0.4 (1/cm)
1.01–6 (K)1–2000 (cm/h)
0.07–0.2 (cm3/cm3)
0.25–0.45 (cm3/cm3)
0.003–0.4 (1/cm)
1.01–5 (K)0–0.1 (1/cm-h)0–0.1 (1/h)
1 (cm2/h)
15, 1.5 (cm)Residual water content in the mobile (macropore) domain
Saturated water content in the mobile (macropore) domain
Water retention parameter in mobile domain
Water retention parameter in mobile domain
Saturated hydraulic conductivity
Residual water content in the immobile (matrix) domain
Saturated water content in the immobile (matrix) domain
Water retention parameter in the immobile domain
Water retention parameter in the immobile domainFirst-order advective exchange termFirst-order diffusion exchange term
Diffusion coefficientLongitudinal, tranverse dispersivity

USDA (1998).
Logsdon (2002), Casey et al. (1997, 1998, 1999) and Jaynes et al. (1995).

the mobile domain, qrm and the values for lT, lL, andD0 (Table 1). The value of qrm was held constant at 0.
The values of lT and lL were kept at fairly high values
to account for heterogeneity in the large field domain
and to help maintain numerical stability. The value ofD0 was also high, though its effect on the overall
dispersion coefficient is likely to be dominated byhydrodynamic dispersion. The potential values of theother model parameters were only weakly constrained
to be within a realistic range based on literature values(Table 1). For parameters, where literature referenceswere unavailable or inconclusive, a relatively wide
range of values was permitted to account for the nonideal behavior of macroporous soils. In particular, uw
was permitted to fluctuate from relatively large
values, as would be expected for short transport
lengths and large macropore conductivities to smallvalues that might result from clay-organic coatings on
the interfaces between the soil matrix and themacropores (Gerke and van Genuchten, 1993). By
not rigidly constraining parameter ranges, non-uniqueparameter sets are possible, especially with the large
number of terms in each model. The possibility of

non-uniqueness is further addressed in the discussionof the modeling results. Following the calibration of
all water flow parameters, the value of us was
adjusted separately for chloride and bromide during
the solute flux simulations.
The numerical simulations for solute flux coupledwith water flow were analogous to the simulations
with water flux alone except that the finite-element
grid was refined to 86!53 finite elements (rather than
the 60!43 grid used in the water flux study) in orderto obtain acceptable mass-balance errors for thesolutes. The solute applications were simulated by
specifying an initial solute concentration to the meshelements approximately corresponding to the top
2 cm of the domain.
Each solute simulation was first run using theeffective parameters estimated from the water flow
simulations, and from these model runs, us was
adjusted based on the cumulative solute flux data(measured solute concentration multiplied by themeasured water flow). Following these simulations,effective parameters were recalibrated based solely on
the cumulative solute flux. The purpose of calibrating

264

N.W. Haws et al. / Journal of Hydrology 313 (2005) 257–273

Water

Plot 51
7.0E-03(A)
A(

Drainage Outflow/Total Cumulative Rainfall

6.0E-03

7.0E-0351 Obs
51-1W

5.0E-03

4.0E-03

4.0E-033.0E-033.0E-032.0E-032.0E-03

1.0E-03

1.0E-03

4/4
7.0E-03

4/14

4/24

5/4

C
(C)

0.0E+00
5/24
4/4
7.0E-03

5/14
51 Obs
51-2W

6.0E-035.0E-03

4.0E-03

4.0E-033.0E-033.0E-032.0E-032.0E-03

1.0E-03

1.0E-03

4/4
7.0E-03

4/14

4/24

5/4

(E)
E

52-1W

4/14

4/24

5/4

5/14

D
(D)

0.0E+00
5/24
4/4

5/14

52 Obs

5/24

52 Obs
52-2W

6.0E-035.0E-03

0.0E+00

BB)

6.0E-035.0E-03

0.0E+00

Plot 52

4/14

4/24

5/4

5/14

5/24

51 Obs51-2Br

6.0E-035.0E-03
4.0E-033.0E-032.0E-03
1.0E-03
0.0E+00

Relative Cumulative Drainage

4/4
0.40
0.35

4/14

4/24

5/4

5/14

5/24
0.40

F
(F)

(G)
G

0.35

52-Obs

0.30

0.30
51-2W

51-1W

0.25
0.20

0.25
0.20

51-2Br

0.15

0.15

0.10

0.10

0.05

0.05
51-Obs

0.00
4/4

4/14

4/24

5/4

5/14

52-1W

52-2W

0.00
5/24

4/4

4/14

4/24

5/4

5/14

5/24

Time (month/day)Fig. 3. Simulated (single-porosity [1W], dual-porosity [2W], and dual-porosity with parameters calibrated from bromide flux [2Br]) andobserved normalized water fluxes for Plots 51 and 52. Cumulative normalized water fluxes shown in (E) and (F).

N.W. Haws et al. / Journal of Hydrology 313 (2005) 257–273

based on solute flux was to further constrain the
inverse solution because the solute flux data comprised both the water flux data and the measured
solute concentrations. Due to the lengthy simulation
times, only the dual-porosity model for Plot 51 wasinvestigated for this part of the study. Plot 51 waschosen for these simulations because of the greater
confidence in the accuracy of the estimates of the
background chloride concentration (see Section 4
below). In this re-calibration of parameter values, theeffective parameter values calibrated in the water flow
data were used as the initial parameter estimates.
Initially, Ks, qsm, and uw were calibrated with theother parameters held constant; then, qsim and the van
Genuchten (1980) parameters in both regions were
calibrated before a final adjustment of uw. The
calibration period for chloride and bromide wasmaintained as defined for the water flux simulations;however, an additional constraint was imposed for thebromide calibration in that the initial breakthrough of
the simulation had to also correspond with theobserved initial bromide breakthrough. As will bedemonstrated, this additional constraint was incorporated due to the difficulty in reproducing the initialarrival of bromide in the drain outflow.

4. Results
4.1. Hydrograph simulations using effective
parameters computed from water flux
Both the single-porosity and dual-porosity simulations successfully capture the main attributes of theflow hydrographs for Plots 51 and 52 (Fig. 3(A)–(D)).
The correlation between simulated and observedhydrographs is generally better for Plot 52 than 51,
as reflected in the values of rxy that are above 0.9 forthe validation period for Plot 52 and only 0.75 and0.72 for the single- and dual-porosity model simulations, respectively, for Plot 51 (Table 2). One
apparent reason for the stronger correlation with the
Plot 52 hydrograph is the difficulty in simulating a
sharp, high peak flow followed by a short recession
limb. Since the peaks are generally smaller and therecession limbs longer in Plot 52, the simulations ofdrainage outflow in Plot 52 fit the observed data
better. The sharp peaks and short recessions are

265

Table 2
Correlation coefficients and sum of squared deviations for waterflow in single-porosity simulations (1W), dual-porosity simulations(2W) and dual-porosity simulation using parameters calibrated frombromide flux (2Br)
Simulation

51-1W
51-2W
52-1W
52-2W
51-2BrCorrelation coefficient(rxy)

Sum of squared deviations(SSD)

Calibration

Validation

Calibration

Validation

0.75
0.72
0.92
0.94
0.64

5.6!10
5.6!10K6
3.8!10K5
3.8!10K6
1.6!10K4

2.6!10K4
1.2!10K4
8.7!10K5
4.7!10K5
4.0!10K4

0.82
0.94
0.91
0.98
0.56

K5

indicative of flow through macropores that dominate
the drainage response. Macropores quickly convey
water to the drain outlet upon filling, and then rapidly
drain as the soil becomes even slightly unsaturated(e.g. Mohanty et al., 1996). That the simulations havedifficulty matching this response is illustrated in thelarge over-predictions of the cumulative discharge(Fig. 3(F) and (G)) and manifests a deficiency inrepresenting the field-scale flow processes. The overpredictions of cumulative discharge in the latter partof the simulations also may be due, in part, to adiscrepancy in the water budget during May. Yet, thisdiscrepancy in the water budget is reconciled by theend of May while the over-predictions in the
cumulative water fluxes continue to increase throughout the duration of the simulation, indicative of aconsistent failure in the simulations to match the rapid
drainage recessions of the measured outflows.
The single-porosity and the dual-porosity simulations have a comparable ability to predict the timing
and shape of the outflow hydrographs - manifested bytheir similar values of rxy (Table 2); however, thesuperiority of the dual-porosity model is evident incomparison to the single-porosity model because of
its ability to more closely match peak discharge
values and recession limbs. The qsm value of the dualporosity model (Table 3) permits an abrupt termination in the recession limbs, and water exchange
between the mobile and immobile regions allows
water to be stored within the flow domain so a waterbalance can be met without prolonged flow throughthe subsurface drain. The advantage of the dualporosity model in replicating both the hydrographbehavior and the cumulative drainage is reflected in266

N.W. Haws et al. / Journal of Hydrology 313 (2005) 257–273

Table 3
Calibrated parameters from HYDRUS-2D inverse modeling of water flow data
Mobile domain

51-1W
51-2W
52-1W
52-2W
51-2BrSand
SCLoamImmobile domain

qrm

qsm

am
(cmK1)

nm

Ks(cm/h)

qrim

qsim

aim
(cmK1)nim

uw (hK1)

us (hK1)

0.25
0
0.25
0
0
0.05
0.09

0.40
0.073
0.40
0.070
0.035
0.43
0.43

0.010
0.036
0.010
0.032
0.01
0.15
0.01

1.75
1.30
1.75
1.51
1.30
2.68
1.23

6.0
4.8
5.0
3.6
25
30
0.07

–
0.10
–
0.10
0.1
–

–
0.30
–
0.50
0.4
–
–

–
0.020
–
0.016
0.01
–
–

–
6.0
–
6.0
6.0
–
–

–
5.0!10K5–
1.6!10K5
1.75!10K4–––
2.5!10K4a–
1.5!10K3a
5.4!10K2a––

Typical values, estimated from the Rosetta database (Schaap, 1999), for sand and silty-clay loam (SCLoam) are given for comparison.
aThis value was used only for the chloride flux simulations. For all simulations of the bromide flux, usZ0.

the SSD values of the dual- porosity simulations thatare a full order of magnitude smaller than the singleporosity simulations for the calibration period andabout twice as small for the validation period(Table 2).
Despite the good reproduction of the mainhydrograph features, the values of the effective
parameters (Table 3) cast ambiguity on how wellthese parameters physically represent the flow
domain. The fitted values for the saturated hydraulicconductivity (Table 3) are in the range of, though in
general much lower than (by a factor of 2–3 times insome cases) steady-state infiltration rates measured
for macroporous soils at a nearby location (Haws
et al., 2004). The precipitation events had rainfall
intensity rates lower than the calibrated Ks value
resulting in a calibrated value beyond the direct
measurement range. Still, these smaller-than-expectedvalues for the hydraulic conductivity, while notunrealistic, do raise concern regarding the uniquenessand physical meaning of the parameters.
The values for van Genuchten (1980) parameters in
the mobile and immobile domains are also unexpected. In theory, the water content in macroporesshould rapidly drain to a residual value at small
tensions, resulting in a hydraulic conductivity function that almost immediately decreases from an
extremely large saturated value to 0 (Logsdon andJaynes, 1993; Mohanty et al, 1996, 1997). Thesemacropore characteristics are comparable to thedrainage behavior of coarse sand. In contrast, thewater release characteristics of the soil matrix regions
should better mimic silty-clay loam soil, having

a much slower release of water and retaining a
relatively high degree of saturation at large tension
values. The more non-linear water retention andhydraulic conductivity functions of sands are represented in the van Genuchten (1980) model withlarger values of a and n, such that these parameters
would be expected to be greater for the macropore
domain in comparison to the soil matrix. The
optimized a values are, as expected, larger for themobile domain than for the immobile domain; yet incontrast to intuition, the calibrated n parameter
exhibits the opposite trend (Table 3, Fig. 4). For thedual-porosity simulations, the n parameter is smaller
for the mobile (macropore) regions and larger for theimmobile (clayey) domain. A possible reason thatthese n parameters are able to reproduce the observed
water flow is that the van Genuchten (1980) modelcannot accurately describe the characteristic curve formacropores at pressure heads near saturation (De Vos
et al., 1997; Mohanty et al., 1997). This nearsaturation region, though a relatively small portionof the full range of pressure heads encountered in
the simulation (0 to K400 cm), would govern the flowof water during the major rainfall events. Due to thelocation of the drain at 100 cm below the soil surface
and the fluctuation of the groundwater level above thisdepth, the pressure heads encountered in the soil
profile are mostly between 0 and K100 cm. In thesimulations, only in the narrow soil surface layer(S30 cm depth) did the calculated pressure heads
drop below K100 cm. Therefore, the parameter
calibrations are dominated by the flow events in the
narrow pressure head range between 0 and K100 cm.

N.W. Haws et al. / Journal of Hydrology 313 (2005) 257–273
1.0

Relative SaturationAinternal flow and transport pathways within the
integrated hydrologic system.

SCLoam

0.8
52-2W (m)4.2. Solute simulations using effective parameters
computed from water flux

0.6
52-1W

0.4
52-1W (im)

0.2
Sand

0.0
0

–50

–100 –150 –200 –250 –300 –350 –400

1.0

B
Relative Saturation

267

SCLoam

0.8
0.6

52-1W
52-2W (m)0.4
0.2

52-1W (im)
Sand

0.0
0

–50

–100 –150 –200 –250 –300 –350 –400

Pressure Head, h (cm)Fig. 4. Water retention characteristics from the soil hydraulicparameters estimated from the inverse modeling for (a) Plot 51 and(b) Plot 52.

Within this range the calibrated parameters produce
water retention curves that would be expected formobile and immobile regions, i.e. with the effective
saturation being larger for the immobile domain ascompared to the mobile domain for the same pressure
heads (Fig. 4).
Another likely reason for the large n values of thesoil matrix region is the possibility of correlations
between the parameters. With nine parameters in
the dual-porosity calibration (or 10 if qrm is not setto 0), there is a high likelihood of correlation. Hence, asmall hydraulic conductivity might be compensatedby unrealistic values of the other parameters. The nparameter is often found to be correlated with othersoil hydraulic parameters, mainly Ks and qr (Šimůnek
and van Genuchten, 1996; Šimůnek et al., 1998).
Consequently, an impressive hydrograph fit could
be achieved by artificially relegating the storage andrelease of water rather than matching the actual

The simulations of the chloride fluxes further show
the superiority of the dual-porosity model over the
single-porosity model. The dual-porosity simulations
better match the general timing and shape of thechloride fluxes, particularly after the first event(Fig. 5(A)–(D)), resulting in reasonable rxy values of
0.67 for Plot 51 and 0.50 for Plot 52 (Table 4).
Because the dual-porosity model can transport waterand solute through a small fraction of the total domain
and also transfer water to, and store water in, the soil
matrix, it can produce higher transport velocities (i.e.
shorter peak arrival times) with smaller overall waterand solute mass flux than the single-porosity model.
Despite the reasonable chloride predictions withthe dual-porosity simulations, the solute flux simulations reveal many of the inadequacies of theeffective parameters that were not apparent in thewater outflow simulations. Even though the initialarrival of chloride is predicted well and early
cumulative fluxes are only slightly under-predicted,
it appreciably over-predicts the cumulative chloride
fluxes in both plots after the calibration period(Fig. 5(F) and (G)). Again, some of this overprediction may be attributable to the potential underbalance in the water budget during May. However, theinitial slight under-prediction of chloride fluxes
followed by large over-predictions indicates that thesimulated mean arrival time for chloride is later thanin the field.
Model deficiencies are more evident with thebromide simulations (Fig. 6). The simulated bromide
breakthrough curves for both the single- and dualporosity models significantly lag the observed bromide fluxes. The simulations with the single-porosity
model never show bromide breakthrough, and eventhe dual-porosity simulations do not predict bromide
breakthrough until the final drainage event. Virtuallyno correlation is evident between observed andsimulated results (Table 4). It should be noted thatthe experimental data for Plot 52 show a slightly
larger relative amount of bromide drained than ofchloride (compare Figs. 5(G) and 6(D)). This is likely

268

N.W. Haws et al. / Journal of Hydrology 313 (2005) 257–273

Chloride

Plot 51
2.0E-02

Plot 52

2.0E-02A

B

51 Obs52-1W

1.6E-0251-1W

1.5E-0252 Obs

1.2E-02
1.0E-028.0E-03

Chloride Mass Flux/Total Chloride Applied

5.0E-03

4.0E-03

0.0E+00
4/4
2.0E-02

4/14

4/24

C

5/4
Peak
=0.025

1.5E-020.0E-00
5/24
4/4
2.0E-02

5/14

51 Obs

4/14

4/24

5/4

5/24

52 ObsDD)
1.6E-02

51-2W

5/14

52_2W

1.2E-02
1.0E-028.0E-035.0E-03

4.0E-03

0.0E+00
4/4
2.0E-02(E)E
1.6E-02

4/14

4/24

Peak
=0.0305/4
Peak=0.095

1.2E-02

0.0E+00
5/24
4/4

5/14

4/14

4/24

5/4

5/14

5/24

51 Obs51-2BrPeak=0.045

8.0E-03

Relative Cumulative Chloride Flux

4.0E-03
0.0E+00
4/4

4/14

4/24

5/4

5/14

5/24

0.40
0.35

0.40

(F)
F

0.30

51-1W
(Total on
5/24=0.53)

51-2Br(Total on
5/24 =0.68)

0.25

52-1W
(Total on
5/24=0.56)

G

0.35
0.30

52-2W

0.25

0.20

0.20

51-2W

0.15

0.15

0.10

0.10

0.05

0.05

51-Obs0.00

52-Obs

0.00
4/4

4/14

4/24

5/4

5/14

5/24

4/4

4/14

4/24

5/4

5/14

5/24

Time (month/day)Fig. 5. Simulated (single-porosity [1W], dual-porosity [2W], and dual-porosity with parameters calibrated from bromide flux [2Br]) andobserved chloride mass flux for Plots 51 and 52. Cumulative chloride flux shown in shown in (E) and (F).

N.W. Haws et al. / Journal of Hydrology 313 (2005) 257–273

269

Table 4
Correlation coefficients and sum of squared deviations statistics for solute transport simulations using effective parameters calibrated from waterflow for single-porosity simulations (1W) dual-porosity simulations (2W) and dual-porosity simulation using parameters calibrated frombromide flux (2Br)
Simulation

Correlation coefficient (rxy)

51-1W
51-2W
52-1W
52-2W
52-2BrSum of squared deviations (SSD)

Chloride

BromideChloride

Bromide

0.27
0.67
0.38
0.50
0.70

0.054
0.069
K0.13
K0.12
K0.18

3.7!10
6.1!10K3
9.1!10K3
2.2!10K3
1.1!10K3

2.7!10K4
2.4!10K4
5.2!10K4
1.3!10K3
2.8!10K2

K3

due to an overestimation of the background chloride
concentration in Plot 52. While it was not expectedthat a model based on the assumption of homogeneityand effective properties would be able to reproduce alarger cumulative chloride flux than bromide, theinability to match the bromide breakthrough timeclearly evidences that the models do not accurately

capture the internal flow and transport processes eventhough they provided good matches with the observed
water flux.
These findings are comparable with those ofprevious studies. In the study of De Vos et al.
(1997), solute arrival times predicted using effective
parameters for a single-porosity model were
BromidePlot 51

Bromide Mass Flux/
Total Bromide Applied

5.0E-03

(A)
A

Peak
=0.032Peak=0.0284.0E-033.0E-0351 Obs51-1W
51-2W
51-2Br(B)
B
4.E-0351-1W(No arrival)

2.0E-0352 Obs
52-1W
52-2W

52-1W
(No arrival)

3.E-0351-2W

2.E-0351-2W

51-2Br

1.0E-03

1.E-030.0E+00

Relative Cumulative Bromide Flux

Plot 52

5.E-03

4/4

4/14

4/24

5/4

5/14

0.E+00
4/4
5/24

0.20

4/14

4/24

5/4

5/14

5/24

0.20

(C)
C

(D)
D
51-2Br(Total on
5/24 =0.67)

0.15

0.15

51-2W

0.10

0.10

51-Obs51-2W51-Obs

0.05

0.05
52-1W
(No arrival)

51-1W(No arrival)

0.00

0.00
4/4

4/14

4/24

5/4

5/14

5/24

4/4

4/14

4/24

5/4

5/14

5/24

Time (month/day)Fig. 6. Simulated (single-porosity [1W], dual-porosity [2W], and dual-porosity with parameters calibrated from bromide flux [2Br]) andobserved bromide mass flux for Plots 51 and 52. Cumulative chloride flux shown in (C) and (D).

270

N.W. Haws et al. / Journal of Hydrology 313 (2005) 257–273

considerably later than the observed results. The
authors surmised that significant improvement could
be made using a dual-domain model. Later, Gerke andKöhne (2004) demonstrated that a dual-domain model
can improve the prediction of solute arrival times.
Though, in the present study, the dual-porosity model
does perform better than the single-porosity model,
the overall failure of both models implies that ahomogeneous, two-dimensional representation cannot
adequately capture the complex flow and transport
pathways of a heterogeneous three-dimensional
domain. Specifically, the lag in the simulated bromide
breakthrough in this study indicates that the areas
contributing water to outflow drainage in the actual
field plots are much greater than those of thesimulated transport domains. This discrepancy precludes affixing a rigid physical basis to the effective
parameters since they were essentially calibrated
across different spatial support areas as opposed towhat contributed to the drainage outflow measured in
the subsurface drains. Hence, even though theeffective parameter models did successfully predictwater flow, it must be conceded that the parameters
themselves are simply empirical fits.
4.3. Effective parameters derived from solute data
The parameters from the inverse optimization
using the cumulative solute flux are included withthe parameters for calibration using only water flow
data in Table 3. In the calibration using the cumulative
chloride flux, no improvement could be made over theparameters calibrated from the water flux (51–2 W).
Therefore, Table 3 does not report any additionalparameters for the chloride flux calibration. In thecalibration using the cumulative bromide flux, thevalues of qsm, as, Ks, and uw were most dramaticallyaltered (Table 3). As expected, the value of qsm
decreased and the value of Ks increased to allow
increased transport velocities. The increase in uw may
be partly due to correlations with changes in the water
retention parameters.
The simulations using the calibrated parameters
based on cumulative bromide flux reproduce theinitial arrival of the applied bromide exceptionallywell. However, after the first drainage event, bromide
flux is dramatically over-predicted (Fig. 6(C)). In
addition, the simulated total water flow is over fourtimes greater than observed for the water flux-based
parameter simulations (Fig. 3(F)). Thus, the rapid
bromide arrival is reproduced only by increasing thewater flux and not by accurately simulating the field
preferential flow and dual-domain transport processes. Before separately adjusting the value of us
for chloride, the chloride flux simulations, using theeffective parameters calibrated from the bromide flux
data, dramatically over-predicted the observed flux
for all events. The simulated chloride flux is kept close
to the observed chloride flux for the first event only
after adjusting us for chloride to a very large value(0.054 hK1). However, after the calibration period,
where the simulated chloride flux is greater than theobserved flux (Fig. 5(E) and (F)), the large differencein the value of us for two similar solutes placed on thesame field is further indication that the parameterslack the ability to mechanistically explain thecomplicated multi-domain processes.

5. Discussion and conclusions
The failure of the single-porosity model, and eventhe dual-porosity model, to simultaneously predictboth water and solute flux at the subsurface drain
outlet using effective parameters can be credited toseveral factors other than the underlying conceptualmodels. A more rigorous and better-constrained
inverse procedure may have improved parameter
optimization. Specifically, the simultaneously fittingof nine parameters is prone to correlated parameters,
non-unique solution sets, and local minima in theobjective function (Hopmans and Šimůnek, 1997). In
addition, the final calibrated values for the hydraulicparameters were certainly influenced by the initial
parameter estimates. Thus, it should be stressed thatthe inability to calibrate realistic parameter values
does not implicate a failure in the conceptual model
per se.
Beyond the limitations of the calibration procedure, however, a fundamental issue contributing to
the deficiencies of the simulations is imposing
homogeneity on a domain, where flow and transport
is highly dependent on the geomorphology of intrinsic
preferential networks. Disregard of preferential pathways often leads to significant disagreement between
the simulated and observed solute arrival times

N.W. Haws et al. / Journal of Hydrology 313 (2005) 257–273(Grayson and Bloschl, 2000) because preferential
pathways likely ‘short-cut’ the flow lines predictedfrom a Richards-based model. This shortcutting of‘ideal’ pathways may also explain why chloride wasbetter simulated than bromide. With chloride placedalmost directly over the drain, its modeled and actual
streamlines would both be fairly unidirectional andthus have less pronounced differences than those forbromide, whose ‘ideal’ streamlines would be longer
and more curvilinear, and therefore, more susceptible
to short-cutting. Similarly, this indicates why a dualdomain model might satisfactorily simulate solutetransport in a one-dimensional flow domain, where
short-cutting is not manifested, but fail to predict
solute arrival in a two-dimensional domain, where
short-cutting can be important.
The over-prediction of water flow in the simulations using parameters calibrated from bromide flux
data further highlight the misrepresentation of the true
field processes. The shortcutting of the ideal streamlines results in quick water and solute arrival times
with relatively little drainage. Because even thehomogeneous dual-porosity model transports waterand solute through ideal streamlines, it can only matchthe observed solute fluxes by flushing a larger volume
of water. The observed rapid arrival times withrelatively little water drainage also suggest a smaller
saturated mobile water content, qsm, than used in theeffective parameter simulations. Although the qsm
values used in the simulations agreed with fieldmeasured values, recent studies report that thecontinuous macropore fraction contributing to theoutflow response may be much less than values
measured at the soil surface (Jorgensen et al., 2002;
Deurer et al., 2003).
A general conclusion of this study is that a model’s
success or failure to represent the flow and transport
processes internal to the transport domain should not
be judged solely by the response prediction at a single
outlet point. As demonstrated in this study, parameters that produced a seemingly good hydrographfit may not accurately predict solute fluxes (representative of internal transport pathways). Thus a goodhydrograph fit does not denote a physical meaning of
the model parameters. Extrapolating the significanceof a model’s success in predicting flow at an outlet
point should be exercised with caution.

271

Acknowledgements
Sincere appreciation is expressed to Dr M. T. van
Genuchten at the USDA Soil Salinity Laboratory,
Riverside, CA for his encouragement to use dualporosity models for field-scale applications. The
authors also thank Dr Sylvie Brouder and Jim Beaty
in the Department of Agronomy at Purdue University
for providing access to the Water Quality Field Station.

References
Baveye, P., Sposito, G., 1984. The operational significance of thecontinuous hypothesis in the theory of water movement throughsoils and aquifers. Water Resour. Res. 20 (5), 521–530.
Beran, M., 1999. Hydrograph prediction—how much skill? Hydrol.
Earth Syst. Sci. 3 (2), 305–337.
Berkowitz, B., Bear, J.C., Braester, C., 1988. Continuum models forcontaminant transport in fractured porous media. Water Resour.
Res. 24 (8), 1225–1236.
Beven, K., 1993. Prophecy, reality and uncertainty in distributedhydrological modeling. Adv. Water Resour. 16 (1), 41–551.
Binley, A.M., Beven, K., Calver, A., Watts, L.G., 1991. Changing
responses in hydrology—assessing the uncertainty in physically
based model predictions. Water Resour. Res. 27 (6), 1253–1261.
Blazkova, S., Beven, K., Tacheci, P., Kulasova, A., 2002. Testing
the distributed water table predictions of Topmodel (allowing
for uncertainty in model calibration): the death of Topmodel?
Water Resour. Res. 38 (11) (Art. no. 1257).
Brouwer, C., Heibloem, M., 1986. Irrigation water management,
training manual no. 3: irrigation water needs, Food andAgricultural Organization of the United Nations, Rome, Italy
1986.
Casey, F.X.M., Logsdon, S.D., Horton, R., Jaynes, D.B., 1997.
Immobile water content and mass exchange coefficient in a fieldsoil. Soil Sci. Soc. Am. J. 61, 1030–1036.
Casey, F.X.M., Logsdon, S.D., Horton, R., Jaynes, D.B., 1998.
Measurement of field hydraulic and solute transport parameters.
Soil Sci. Soc. Am. J. 62, 1172–1178.
Casey, F.X.M., Logsdon, S.D., Horton, R., Jaynes, D.B., 1999.
Comparing field methods that estimate mobile–immobile model
parameters. Soil Sci. Soc. Am. J. 63, 800–806.
Deurer, M., Green, S.R., Clothier, B.E., Bottcher, J.,
Duijnisveld, W.H.M., 2003. Drainage networks in soils. A concept
to describe bypass-flow pathways. J. Hydrol. 272, 148–162.
De Vos, J.A., Simunek, J., Raats, P.A.C., Feddes, R.A., 1997.
Identification of the hydraulic characteristics of a layered silt
loam. In: van Genuchten, Th, M., Leij, F.J., Wu, L. (Eds.),
Characterization and Measurement of Hydraulic Properties of
Unsaturated Porous Media. University of California, Riverside,CA, pp. 783–798.

272

N.W. Haws et al. / Journal of Hydrology 313 (2005) 257–273

Doherty, J., Johnston, J.M., 2003. Methodologies for calibration andpredictive analysis of a watershed model. J. Am. Water Resour.
Assoc. 39 (2), 251–265.
Durner, W., Schultz, B., Zurmühl, T., 1997. State-of-the-art in inverse
modeling of inflow/outflow experiments. In: van Genuchten,Th, M., Leij, F.J., Wu, L. (Eds.), Characterization andMeasurement of Hydraulic Properties of Unsaturated PorousMedia. University of California, Riverside, CA, pp. 661–682.
Fernandez-Garcia, D., Sanchez-Vila, X., Illangasekare, T.H., 2002.
Convergent-flow tracer tests in heterogeneous media: combined
experimental-numerical analysis for determination of equivalent transport parameters. J. Contam. Hydrol. 57, 129–145.
Fipps, G., Skaggs, W.R., Neiber, J.L., 1986. Drains as boundary
conditions in finite elements. Water Resour. Res. 22 (11),1613–1621.
Gerke, H.H., Köhne, J.M., 2004. Dual-permeability modeling ofpreferential bromide leaching from a tile-drained glacial tillagricultural field. J. Hydrol. 289, 239–257.
Gerke, H.H., van Genuchten, M.T., 1993. Evaluation of a first-orderwater transfer term for variably saturated dual-porosity flow
models. Water Resour. Res. 29 (4), 1225–1238.
Gerke, H.H., van Genuchten, M.T., 1996. Macroscopic representationof structural geometry for simulating water and solute movementin dual-porosity media. Adv. Water Resour. 19 (6), 343–357.
Grayson, R., Bloschl, G., 2000. Spatial processes, organisation andpatterns. Spatial patterns in catchment hydrology: observations
and modeling. In: Grayson, R., Bloschl, G. (Eds.), Spatial
Patterns in Catchment Hydrology: Observations and Modelling.
Cambridge University Press, Cambridge, UK, pp. 3–16.
Haws, N.W., Liu, B., Boast, C.W., Rao, P.S.C., Kladivko, E.J.,
Franzmeier, D.P., 2004. Spatial variability and measurement
scale of infiltration rate on an agricultural landscape. Soil Sci.
Soc. Am. J. 68, 1818–1826.
Hoffman, G.J., Schwab, G.O., 1964. Tile spacing prediction basedon drain outflow. Trans. ASAE 7, 444–447.
Hopmans, J.W., Šimůnek, J., 1997. Review of inverse estimation ofsoil hydraulic properties. In: van Genuchten, Th, M., Leij, F.J.,
Wu, L. (Eds.), Characterization and Measurement of Hydraulic
Properties of Unsaturated Porous Media. University of California, Riverside, CA, pp. 643–660.
Indelman, P., 2002. On mathematical models of average flow inheterogeneous formations. Transport Porous Media 48, 209–224.
Jacques, D., Feyen, J., Mallants, D., 1997. Determination of
hydraulic properties using tension infiltrometer data and inverse
optimization. In: van Genuchten, Th, M., Leij, F.J., Wu, L.
(Eds.), Characterization and Measurement of Hydraulic Properties of Unsaturated Porous Media. University of California,
Riverside, CA, pp. 749–760.
Jain, S.K., Storm, S.B., Bathurst, J.C., Refsgaard, J.C., Singh, R.D.,
1992. Application of the SHE to catchments in India. Part 2.
Field experiments and simulation studies with SHE or the kolar
subcatchment of the Naramada river. J. Hydrol. 140, 25–47.
Jaynes, D.B., Logsdon, S.D., Horton, R., 1995. Field method formeasuring mobile/immobile water content and solute transfer
rate coefficient. Soil Sci. Soc. Am. J. 59, 352–356.
Jorgensen, P.R., Hoffman, M., Kistrup, J.P., Bryde, C., Bossi, R.,
Villholth, K.G., 2002. Preferential flow and pesticide transport

in a clay-rich till: field, laboratory, and modeling analysis.
Water Resour. Res. 38 (11), 1246.
Logsdon, S.D., 2002. Determination of preferential flow model
parameters. Soil Sci. Soc. Am. J. 66, 1095–1103.
Logsdon, S.D., Jaynes, D.B., 1993. Methodology for determining
hydraulic conductivity with tension infiltrometers. Soil Sci. Soc.
Am. J. 57, 1426–1431.
Lugo, J.E., Rio, J.A.D., Taguena-Martinez, J., 1998. Influence of
nonlinear local properties on effective transport. Transport
Porous Media 31, 89–108.
Madsen, H., Wilson, G., Ammentorp, H.C., 2002. Comparison ofdifferent automated strategies for calibration of rainfall-runoff
models. J. Hydrol. 261 (1–4), 48–59.
Mohammad, F.S., Skaggs, W.R., 1983. Drain tube opening effects
on drain inflow. Irrig. Drain. Div., Am. Soc. Civ. Eng. 109 (4),393–404.
Mohanty, B.P., Horton, R., Ankeny, M.D., 1996. Infiltration andmacroporosity under a row crop agricultural field in a glacial till
soil. Soil Sci. 169 (4), 205–213.
Mohanty, B.P., Bowman, R.S., Hendrickx, J.H.M., vanGenuchten, M.T., 1997. New piecewise-continuous hydraulicfunctions for modeling preferential flow in an intermittentflood-irrigated field. Water Resour. Res. 33 (9), 2049–2063.
Neuman, S.P., Orr, S., 1993. Prediction of steady state flow innonuniform geologic media by conditional moments: exact
nonlocal formalism, effective conductivities, and weak approximation. Water Resour. Res. 29 (2), 341–364.
Royer, P., Auriault, J.-L., Lewandowska, J., Serres, C., 2002.
Continuum modelling of contaminant transport in fractured
porous media. Transport Porous Media 49, 333–359.
Schaap, M., 1999. Rosetta Help File. US Salinity Laboratory,
Riverside, CA.
Šimůnek, J., van Genuchten, M.T., 1996. Estimating unsaturated
soil hydraulic properties from tension disc infiltrometer data bynumerical inversion. Water Resour. Res. 32 (9), 2683–2696.
Šimůnek, J., Wendroth, O., van Genuchten, M., Th, A., 1998. A
parameter estimation analysis of the evaporation method for
determining soil hydraulic properties. Soil Sci. Soc. Am. J. 62(4), 894–905.
Šimůnek, J., Sejna, M., van Genuchten, M.T., 1999. The HYDRUS2D software package for simulating the two-dimensional movement of water, heat, and multiple solutes in variably saturated
media. Version 2.0, IGWMC-TPS-53, International Groundwater
Modeling Center, Colorado School of Mines, Golden, CO.
Šimůnek, J., Jarvis, N.J., van Genuchten, M.T., Gardenas, A., 2003.
Review and comparison of models for describing nonequilibrium and preferential flow and transport in the vadose zone.
J. Hydrol. 272, 14–35.
Skaggs, W.R., 1976. Determination of the hydraulic conductivity—drainable porosity from water table measurements. Trans.
ASAE 19, 73–84.
Tartakovsky, D.M., Neuman, S.P., 1998a. Transient flow inbounded randomly heterogeneous domains 1. Exact conditional
moment equations and recursive approximations. Water Resour.
Res. 34 (1), 1–12.
Tartakovsky, D.M., Neuman, S.P., 1998b. Transient flow inbounded randomly heterogeneous domains 2. Localization of

N.W. Haws et al. / Journal of Hydrology 313 (2005) 257–273
conditional mean equations and temporal nonlocality effects.
Water Resour. Res. 34 (1), 13–20.
Tartakovsky, D.M., Neuman, S.P., 1998c. Transient effective
hydraulic conductivities under slowly and rapidly varying
mean gradients in bounded three-dimensional media. Water
Resour. Res. 34 (1), 21–32.
United States Department of Agriculture (USDA), 1998. Soilsurvey of Tippecanoe county, Indiana. Purdue University
Agricultural Experiment Station, Indiana. Indiana Department of Natural Resources, State Soil Conservation Board

273

and the Division of Soil Conservation. West Lafayette,
Indiana.
van Genuchten, M.T., 1980. A closed-form equation for predicting
the hydraulic conductivity of unsaturated soils. Soil Sci. Soc.
Am. J. 44, 892–898.
Yeh, T.-C.J., 1997. A site characterization method for the vadosezone. In: van Genuchten, Th, M., Leij, F.J., Wu, L. (Eds.),
Characterization and Measurement of Hydraulic Properties of
Unsaturated Porous Media. University of California, Riverside,CA, pp. 1377–1380.