Published online December 16, 2005

Inverse Dual-Permeability Modeling of Preferential Water Flow in a Soil Column and
Implications for Field-Scale Solute Transport
J. Maximilian Köhne,* Binayak P. Mohanty, and Jirka Šimůnek
flow and solute transport. While the soil hydraulic functions for DPM can be independently estimated fromwater retention and hydraulic conductivity measurements with double-hump curve shapes (e.g., Mohanty
et al., 1997; Köhne et al., 2002), this approach cannot
provide the parameters of water transfer between matrix and macropore domains that are also required. Inverse simulation is an alternative way to estimate DPM
parameters that are otherwise difficult to obtain. However, questions about the uniqueness of such optimized
parameters still exist (Šimůnek et al., 2001). Several
investigations on preferential water flow mechanisms
in soils with PFPs were reported (e.g., Booltink andBouma, 1991; Tofteng et al., 2002; Logsdon, 1995), but
only few preferential water flow experiments have been
analyzed in detail for evaluation of the region-specifichydraulic parameters related to the matrix or the PFPs
(e.g., Šimůnek et al., 2001; Castiglione et al., 2003; Köhne
and Mohanty, 2005). While most studies employed a mixture of measurements, estimates, and manual calibrationof DPM parameters to simulate preferential solute transport (e.g., Jarvis et al., 1991; Jabro et al., 1994; Andreu
et al., 1994; Ahuja et al., 1995; Larsson and Jarvis, 1999a,
1999b; Malone et al., 2001; Gerke and Köhne, 2004), inverse parameter optimization schemes have been employed as well (e.g., Schwartz et al., 2000; Abbaspour
et al., 2001; Kätterer et al., 2001).
Traditionally, the parameters of (nonpreferential)
flow and transport models are estimated sequentially.
First the soil hydraulic parameters are estimated using
hydraulic observations, followed by the estimation oftransport parameters using solute concentration data
(Jacques et al., 2002; Abbasi et al., 2003). However,comparisons between this sequential fitting procedureand simultaneous estimation of hydraulic and transportparameters from combined data of water flow and solute concentrations rarely showed the sequential scheme
to be superior (e.g., Abbasi et al., 2003), while otherstudies found the simultaneous scheme to be more efficient (Medina and Carrera, 1996; Mishra and Parker,
1989; Inoue et al., 2000). For example, Šimůnek et al.
(2002) used experimental data of Inoue et al. (2000) andshowed that although the final parameter estimates forboth sequential and simultaneous optimizations werevery similar, the confidence intervals for the simultaneous optimization of soil hydraulic and solute transportparameters, as well as the final value of the objectivefunction, were smaller. It remains unclear whether thesequential parameter estimation scheme can be employed for the inverse simulation of preferential water

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

ABSTRACTThe question of whether or not soil hydraulic parameters of dualpermeability models (DPM) can be properly identified by inverse
analysis of preferential water flow data has not been resolved to date.
We applied a DPM based on two coupled Richards’ equations tocompare the performance of inverse and forward simulations of laboratory preferential flow data. Infiltration and drainage experiments
were conducted using a repacked loam soil column (80 cm long, 24-cm
diameter) containing a cylindrical sand region (2.4-cm diameter) as thepreferential flow path (PFP) along its central axis. The forward DPM
water flow simulations relied on hydraulic parameters for the matrix
and the PFP as determined by means of separate infiltration anddrainage experiments on loam and sand columns, respectively. One
inverse DPM approach relied on standard (lumped) observations ofinfiltration and outflow, while the other included outflow through thematrix and the PFP. Both inverse approaches provided accurate
matches of bulk infiltration and outflow, but the outflow out of thematrix and the PFP could only be described when fitting the DPM toregion-specific outflow data. The practical implication of this finding
for predicting solute transport in natural soils was evaluated. An observed tile-drainage hydrograph was used for inverse hydraulic DPM
parameter estimation, followed by fitting the solute transfer coefficientand the dispersivity for simulating Br2 tracer concentrations. Thissequential fitting procedure was successful for hydrograph simulation
but unsuccessful for Br2 breakthrough simulation. Simultaneous hydraulic and transport parameter estimation considerably improved
the approximation of Br2 concentrations. This study shows that ahydrograph alone is not sufficient for inverse identification of soil
hydraulic DPM parameters. Simultaneously employing hydrographand solute breakthrough data may facilitate identification of hydraulicand transport DPM parameters to characterize preferential flow andsolute transport.

D

UAL-PERMEABILITY MODELS are increasingly used
for analysis of preferential solute transport instructured soils that contain PFPs, such as fractures
or macropores, both on the laboratory column scale(e.g., Gwo et al., 1995, 1996; Allaire et al., 2002a, 2002b;
Castiglione et al., 2003; Greco, 2002; Kätterer et al., 2001)
and on the plot or field scale (e.g., Jarvis et al., 1991;
Andreu et al., 1994; Larsson and Jarvis, 1999a, 1999b;
Kohler et al., 2001). A problem receiving considerable
attention recently is the determination of the many DPM
parameters required for simulating preferential water

J.M. Köhne and B.P. Mohanty, Dep. of Biological & AgriculturalEngineering, Texas A&M Univ., Scoates Hall, College Station, TX
77843-2117; J. Šimůnek, Univ. of California, Dep. of Environmental
Sciences, 900 University Av., A135 Bourns Hall, Riverside, CA 92521;
J.M. Köhne, currently, Institute for Land Use, Univ. Rostock, Justusvon-Liebig Weg 6, 18059 Rostock, Germany. Received 14 Jan. 2005.
*Corresponding author (max.koehne@uni-rostock.de).
Published in Vadose Zone Journal 5:59–76 (2006).
Original Researchdoi:10.2136/vzj2005.0008ª Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA

Abbreviations: DPM, dual-permeability model; IC, initial conditions;
LBC, lower boundary conditions; PFP, preferential flow path; SWCC,
soil water characteristic curve; TDR, time domain reflectometry.

59

60

VADOSE ZONE J., VOL. 5, FEBRUARY 2006

Table 1. Particle size distribution, organic carbon content (OC), and bulk density (r) of the soil materials used in the homogeneous soilcolumns (i.e., loam, medium-sand, and coarse-sand columns) and in the columns with preferential flow paths (PFP).
Particle size fractions
Material (pore domain)†

Coarse sand
2–0.5 mm

Medium sand
0.5–0.2 mm

Fine sand
0.2–0.05 mm

0.0
87.8
16.9

0.6
11.3
55.1

52.5
0.9
27.4

Silt 0.05–0.002
mm

Clay,0.002 mm

OC

r†

35.8
0
0. 6

11.1
0
0

0.24
n.a.
n.a.

g cm23
1.66 (1.59)
1.65 (1.61)
1.64 (1.60)

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

%
Loam (matrix)
Coarse sand (PFP)
Medium sand (PFP)† Name or number in parentheses is for the column with PFP.

preferential solute transport prediction. Observations ofdrainage flux and preferential Br2 tracer breakthroughin a tile-drained field between 1991 and 1994 (Lennartz
et al., 1999) were simulated, and we compared the performances of sequential and simultaneous inverse estimationof soil hydraulic and solute transport DPM parameters.

flow and solute transport in the vadose zone, since waterflow observations may not always contain sufficient
information for proper inverse identification of hydraulic DPM parameters (Šimůnek et al., 2001). A recent
study showed that despite an excellent fit of the tiledrainage hydrograph, a two-dimensional dual-porosity(mobile–immobile) model may fail to describe observed
tracer breakthrough (Haws et al., 2005).
These findings from previous studies suggest that (i)hydrographs may not carry enough information to distinguish between local water flow paths with distinctly
different water flow velocities and (ii) while domainspecific flow velocities may not be critical to describewater flow, they may become significant for simulating
preferential solute transport. Although reasonable, theabove conclusions are only partly supported by experimental data, since these data were obtained rather
indirectly at the field scale, where many additional confounding unknowns are present, such as soil heterogeneity, off-site contributions, catchment boundaries, andinitial and boundary conditions. Hence, a closer examination of the underlying causes of the above findings
seems warranted.
The primary objective of this study was to investigate if
DPM hydraulic parameters for matrix and PFP regions
can be identified based on standard soil hydraulic observations, such as infiltration, outflow, pressure heads,
and water contents of a well-controlled column experiment, or if additional domain-specific outflow throughmatrix and PFP must be included. The secondary objective was to evaluate the practical implications of the
DPM hydraulic parameter identification for field-scale

MATERIALS AND METHODS
Physical Soil Properties
The particle size distributions of the sandy loam soil (henceforth referred to as loam) and coarse and medium sand arelisted in Table 1. The loam, sampled from arable land near
College Station, TX, was air-dried and ground. Additionally,
30 smaller soil columns (7.5 cm in diameter and height) were
packed with loam and the medium and coarse sand (10 columns per material) to measure the soil water characteristic
curves (SWCC) between saturation and 2500 cm pressure
head according to Klute (1986). Based on the 10 sets of data
for each material, the average water retention functions were
determined by fitting soil hydraulic functions (van Genuchten,1980) using the RETC model (van Genuchten et al., 1991).
Furthermore, the saturated hydraulic conductivity, Ks, wasmeasured on the same soil columns (Klute and Dirksen, 1986).
The results for Ks and the SWCC are provided in Table 2.

Flow Experiments with the Homogeneous Soil Columns
Laboratory soil columns 16 cm in diameter and 40 cm highwere packed with loam, coarse sand, and medium sand, respectively. These air-dry soil materials were packed using apiston compactor at 3-cm increments to obtain relatively homogeneous bulk densities between 1.6 and 1.7 g cm23. The
bulk densities theoretically correspond to porosities ranging

Table 2. Soil hydraulic van Genuchten’s (1980) parameters and the range of 95% confidence limits obtained by inverse solution of the
Richards’ equation (Šimůnek et al., 1998) based on experimental infiltration and drainage data for homogeneous loam and medium-sandand coarse-sand columns. The parameters were used in the forward dual-permeability simulations. Corresponding parameters are also
given for the soil water characteristic curve (SWCC) measured at static equilibrium and for Ks at steady-state flow (SSF).
ur

Exp.

us
3

cm cm

23

a

n

21

cm
Loam

Ks
cm d

R
21

d

Iinfiltration
Drainage
SWCC, SSF

0.208 6 0.007
0.195 6 0.016
0†0.340†
0.340†
0.391 6 0.006

0.0051 6 0.0001
0.0090 6 0.0000
0.0035 6 0.000
Medium sand2.67 6 0.223
1.27 6 0.07
2.40 6 0.32

1.7 6 0.1
7.4 6 3.7
11.2 6 2.6

n.a.¶
n.a.
n.a.

Infiltration
SWCC, SSF

0†
0.0

0.300†
0.370 6 0.012

0.0296 6 0.0015
0.0234 6 0.0038
Coarse sand1.93 6 0.16
2.04 6 0.14

609 6 34
850 6 274

n.a.
n.a.

Drainage
SWCC, SSF

0†0†0.340 6 0.0
0.370†

0.2500‡ 6 0.0140
0.1561 6 0.0421

1.76 6 0.027
2.24 6 0.31

10917 6 1476
15542 6 7091

0.0335 6 0.0013
n.a.

† Fixed (not fitted).
‡ At upper constraint.
¶ Not available.

61

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

between 0.35 and 0.4, if a particle density of 2.65 g cm23 isassumed. After slowly saturating the columns with deionizedwater from the bottom to the top to avoid air inclusions, outflow and infiltration experiments were conducted to determine
the soil hydraulic parameters. A tension infiltrometer with a
16-cm-diameter infiltration disc (i.e., matching the soil column)
was used for infiltration at zero pressure head. At the column
bottom the soil was in contact with a nylon membrane supported by a densely perforated PVC plate. The lower (seepage)
boundary was subject to adjustable suctions controlled by avacuum system with an automated outflow monitoring system.
Pressure heads were recorded across the soil column with tensiometers at 10-, 20-, and 30-cm depths, while water contents
at the 10- and 30-cm depths were measured using time domainreflectometry (TDR) probes of 8-cm length (the same probes
were used in the preferential flow experiments described below).
Experiments were conducted for the following initial andlower boundary conditions (IC and LBC): (i) infiltration into
the loam column (IC: 2140 cm, LBC: 0 cm), (ii) three-step outflow out of the initially saturated loam column (LBC: 0, 260,
and 2120 cm), (iii) infiltration into the medium-sand column(IC: 240 cm (top) to 220 cm (bottom), LBC: 215 and 0 cm),(iv) drainage from the initially saturated coarse-sand column(LBC: 0, 230, and 260 cm). Infiltration into coarse sand wasmeasured as well; however, unstable (fingering) flow was observed, and the results were excluded from further analysis.

Preferential Flow Experiments
Our soil column setup for studying preferential flow was
described in detail in a recent study (Köhne and Mohanty,2005). Only those setup features are mentioned here that are
relevant to the present study, or that were modified. The setup
consisted of a large loam soil column (24.4-cm diameter, 80-cm
length, henceforth referred to as matrix) with a cylindricalsand-filled PFP (2.4-cm diameter) at the center. A 5-cm-highcylindrical flow separator at the bottom was used to prevent
water exchange between the matrix and the PFP due toboundary effects. The PFP region was implemented by vertically attaching a metal tube of 2.4-cm diameter to the flow
divider, pouring coarse sand into the metal tube, and slightly
compacting the sand by knocking on the tube’s side-walls. The
loam was packed similarly as for the homogeneous columns,
leading to a comparable bulk density (Table 1). After both thematrix and PFP domains were established, the metal tube wascarefully pulled out of the column. The matrix and PFP domains were thus in direct contact with each other, but without
any lateral soil mixing or intrusions between loam and sand,as was verified in horizontal cross sections taken after completion of the experiments. The soil column was slowly saturated
from the bottom with deionized water before the first (drainage) experiment.
For infiltration under constant head boundary condition, atension infiltrometer with infiltration disc of matching diameter (24 cm) was employed. Column outflow through the bottom boundary (consisting of a perforated plate with nylon
membrane) was monitored separately for the matrix and the
PFP. Pressure heads within the PFP were measured withcustom mini-tensiometers (porous ceramic cup: 0.2-cm diam.,
1.5-cm length) at depths of 5, 20, 35, 50, and 65 cm. At the same
depths, another set of tensiometers with a 0.6-cm-diameter
ceramic cup was installed in the soil matrix at distances of 6
and 2 cm from the macropore wall, respectively.
Corresponding to the mini-tensiometers, 8-cm-long TDR
probes were used to measure the matrix water contents at thesame depths of 5, 20, 35, 50, and 65 cm. Conventional TDR
probes are too large to measure water contents in the PFP.

Nissen et al. (1998) presented coiled TDR probes with smaller
dimensions that allow measurements along short distances. We
found these to be more suitable to measure water contents in
the PFP, even if they are inherently less accurate than TDR
probes. Four TDR coil probes (0.3-cm diam., 2-cm length)were built according to the design of Nissen et al. (1998) andused to monitor the water content at depths of 20 and 35 cm
inside the sand-filled PFP and within the loam matrix at a
distance of 2 cm from the PFP. The TDR and TDR coil probes
were calibrated separately to measure water contents in thematrix loam and PFP sand.
Three preferential flow experiments were selected forthis study:
Infiltration under zero pressure head into the loam column
with a medium-sand PFP. Based on tensiometer readings, theinitial hydrostatic pressure head profile increased linearly from2100 cm at the top to 220 cm at the bottom of the soil column.
A constant pressure head of 220 cm was imposed at the bottom boundary.
Infiltration at zero applied pressure head into the loam column with a coarse-sand PFP for somewhat drier initial conditions. The hydrostatic pressure head profile varied between
2150 cm at the top and 270 cm at the bottom. The bottom outlets of the PFP and the matrix were open to the atmosphere.
Drainage of the initially fully saturated loam column witha coarse-sand PFP. A pressure head of 210 cm was applied at
the lower boundary until 138 min, when it was decreased to
245 cm.

Field Bromide Displacement Tests
For demonstrating the application and suitability of ourparameter estimation schemes, we used a field Br2 tracer
data set from Lennartz et al. (1999). The experimental fieldsite Bokhorst (northern Germany) had a monitored area of5010 m2. Tile lines were installed at a depth of 1 m and at a
spacing of 11 to 14.5 m. The loamy soil at the site was classified
as a poorly drained Dystric Gleysol. The 0.5-ha experimental
area was identified as heterogeneous in terms of particle size
distribution. While the clay content decreased with soil depth
from 18% to about 10% at the 1-m depth, bulk density increased from 1.46 to 1.9 g cm23. The high clay content of thesoil surface layer induced cracks and fissures during dry
periods. The conventional tilled experimental plot was under a3-yr crop rotation with winter rape (Brassica napus L. var.
napus) (1991/1992), winter wheat (Triticum aestivum L.) (1992/
1993), winter barley (Hordeum vulgare L.), (1993/1994), andwinter rape (1994/1995). In three Br2 tracer tests conducted
between 1991 and 1994, preferential flow was found to be amajor factor controlling solute transport. Bromide tracer wasapplied once every year in late fall as a KBr solution, and wasmonitored in the drainage discharge during the following 5 mo.
Applied mass was comparable among the three tracer tests(ranging between 16.75 and 19.3 kg). Readers are referred to
Lennartz et al. (1999) for further details about the site, soil,
and the tracer experiments.

Numerical Model Simulations
Flow in Homogeneous Soil Columns
The water flow experiments were described by fitting thenumerical solution of the Richards’ equation as implemented
in HYDRUS-1D (Šimůnek et al., 1998) to the observations ofinflow; outflow; pressure heads at the 10-, 20-, and 30-cm
depths; and water contents at the 10- and 30-cm depths. The
resulting van Genuchten (1980) soil hydraulic parameters rep-

62

VADOSE ZONE J., VOL. 5, FEBRUARY 2006

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

resented the matrix (loam) and two PFP materials (medium
sand and coarse sand). The residual water content, ur (L3 L23),
was set to zero for the medium and coarse sand. This assumption did not adversely affect the description of water flowin the relatively wet water content range as operated in thewater flow experiments. To facilitate the inverse parameter
identification for loam and medium sand, the saturated watercontent, us (L3 L23), was fixed at the measured value.
Preferential Flow in a Soil Column
Dual-permeability models describe preferential flow using
one equation for flow in the matrix and one for flow in thefracture pore system or PFPs, with the two flow regions
coupled by a water transfer term. The DPM of Gerke andvan Genuchten (1993a) as implemented into the modified
HYDRUS-1D (Šimůnek et al., 2003) uses the following two
mixed-type Richards equations to describe transient waterflow in fractures (Eq. [1a]) and the matrix (Eq. [1b]):

]uf
]
]hf
Gw
5
2Kf 2
Kf
]t
]z
wf
]z

[1a]

]um
]
]hm
Gw
5
2 Km 1
Km
]t
]z
12wf
]z

[1b]

1

2

1

2

where the subscript “f” defines a property of the fracture poresystem or PFP, the subscript “m” represents the matrix, u is thewater content (L3 L23), h is the pressure head (L), K denotes
the hydraulic conductivity function (L T21), wf is the fraction of
total soil occupied by the fracture system (0 , wf , 1), and
Gw is the water transfer term (T21)Gw 5b
Ka gw (hf 2 hm )a2[2]

where b is a dimensionless geometry coefficient, a is thecharacteristic length of the matrix structure (L), Ka is thehydraulic conductivity of the fracture–matrix interface (L T21),
and gw is a dimensionless scaling coefficient for which a valueof 0.4 was suggested (Gerke and van Genuchten, 1993b).
Because of the absence of depositional coatings at the PFP–
matrix interface in our repacked soil columns, we assumed that
Ka in the transfer term can be described as

Ka 5 0:5[Km (hf ) 1 Km (hm )][3]

The geometry factor b was calculated assuming hollow
cylindrical PFP geometry according to Gerke and van
Genuchten (1996) as1
b5;[0:19ln(161)]2

1 , 1 , 100[5]in which b is the cylindrical PFP radius. The volume fraction of
the PFP, wf, was calculated as

wf 5

b2
1
2 5
12
(a 1 b)Field-Scale Preferential Flow and Bromide Transport
For simplicity, the entire 0.5-ha tile-drained field was represented by a one-dimensional vertical single-layered model
domain, since assuming several layers or accounting for twodimensional heterogeneity was impractical for inverse DPM
simulation because of the increasing number of parameters.
Soil hydraulic and solute transport DPM parameters of the
one-dimensional homogeneous domain were assumed to represent effective soil properties lumped over heterogeneities insoil texture, structure, and effects of tile-induced flow-field.
Preferential water flow was described with Eq. [1] through[3]. However, since the parameters a, Ka, b, and gw in Eq. [2]are difficult to define on the field scale, they were lumped
into a single water transfer coefficient, vw (T21 L21), whichyields the following simplified expression for the water transfer term:

Gw 5 vw (hf 2 hm )[6]

The numerical grid consisted of 201 elements with a vertical
nodal spacing of 0.4 cm to represent the 80-cm-long column.
The values for mass-transfer related DPM parameters were

[7]

Convection–dispersion equations simulate transport in boththe fracture and matrix regions, again using the DPM of Gerkeand van Genuchten (1993a) as implemented into the modified
HYDRUS-1D (Šimůnek et al., 2003), as follows:]uf cf]]cf]qf cfGs
5
2
uf D f
2
]z
]t
]z
]z
wf

[8a]

]um cm]]cm]qm cm
Gs
5
1
um Dm
2
]t
]z
]z
12wf
]z

[8b]

1

1

2

2

where for the fracture and matrix regions, c is the solute concentration (L3 M21), q is the volumetric flux density (L T21),
and D is the dispersion coefficient (L2 T21) described with thefollowing equation:

uf Df 5 uf D0

u7/3
f
1 lf |qf |
u2s,f

[9a]

um Dm 5 um D0

u7/3m
1 lm |qm |
u2s,m

[9b]

[4]

where the parameter z was calculated as

a1b
15
b

derived from the experimental PFP–matrix geometry (a 511-cm matrix mantle radius, b 5 1.2-cm PFP radius, z 5 10.167,
b 5 1.0685, and wf 5 0.009675).
The initial condition was defined by the hydrostatic pressure
head profile interpolated between initial tensiometer readings
at various depths. The lower boundary in the model was represented by a seepage face boundary condition, with no flow
when the boundary was unsaturated, and a zero (atmospheric)pressure head once saturation was reached. For those experiments where suction was applied to the lower boundary, theseepage face boundary condition was modified by using thenegative value of applied suction as the threshold value between zero flux and constant pressure head.

where once again for the fracture and matrix regions, D0
represents the molecular diffusion (L2 T21), l signifies thedispersivity (L), and us denotes the saturated water content.
The rate of solute mass transfer between the fracture andmatrix regions, Gs, is given as the sum of diffusive and convective fluxes:

Gs 5 vs (1 2 wf )(cf 2 cm ) 1 Gw c*[10]where vs is a first-order diffusive solute mass transfer coefficient (T21), and c* is equal to cf for Gw . 0 and cm for Gw , 0.

63

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

The initial condition for water flow was assumed as a hydrostatic profile between zero pressure head at drainage depth(100 cm) and 2100 cm at the top. Measured rainfall rates lessthe actual evapotranspiration were applied as the upper boundary condition. A seepage face boundary condition was prescribed for water flow at the lower boundary. The initial
Br2 concentration was assumed to be zero in 1991 and 1993. In
1994, resulting from the prior tracer tests, an initial concentration of 4 mg L21 was assumed throughout the soil. The top
boundary condition for solute transport was represented by theprescribed Br2 concentration of the applied KBr solution,
while a zero Br2 concentration gradient was prescribed at thelower boundary.

Inverse Parameter Identification

plate had a higher perforation density below the PFP thanbelow the matrix.
Three different approaches were compared for estimating
the parameters of the DPM. First, the inverse-lumped approach relied on observations of pressure heads in the matrix
at 5-, 35-, and 65-cm depths, water contents in the matrix at
35-cm depth, and cumulative infiltration (across top boundary)or drainage (across bottom boundary) fluxes. Second, theinverse-local approach additionally utilized observations ofregion-specific (matrix and PFP) outflow, that is, data thatwere not available in earlier column studies or field experiments. Finally, in the forward approach we studied if independently determined soil hydraulic parameters for loam,
medium sand, and coarse sand could be successfully used forforward DPM simulation of the preferential flow data in therepacked soil columns with various matrix–PFP compositions.

Preferential Flow in a Soil Column
For describing preferential water flow in the soil column
using Eq. [1] to [6], the DPM requires 13 parameters (i.e., um,r,
um,s, am, nm, Km,s, uf,r, uf,s, af, nf, Kf,s, a, wf, b), of which threeparameters related to water transfer between domains are
known (a, wf, b). To additionally reduce the number of unknowns, the residual water content for the fracture poresystem, uf,r, was fixed at zero. This limited the total number of
unknown water flow parameters to nine. The DPM model
parameters were estimated by minimizing an objective function, F, (Šimůnek et al., 1998) representing the deviations between measured and calculated variables as follows:

O O w [O (z,t ) 2 E (z,t ,b)]m

F(b) 5

n

i, j

j

i

j

i

2

[11]

j51 i51

where m is the number of different sets of measurements,
n represents the number of observations in a particular measurement set, Oj(z,ti) represents observations at time ti for thejth measurement set at location (z), Ej(x, ti, b) are the corresponding estimated space–time variables for the vector b ofnine optimized van Genuchten (1980) parameters (um,r, um,s,am, nm, Km,s, uf,s, af, nf, Kf,s) and wi,j are weighting factorsthat here were manually calibrated in repeated inverse simulation runs until similar contributions of all observation sets
to the weighted least squares (Eq. [11]) of residuals between
observations and model estimates (Eq. [11]) were achieved.
This ad hoc procedure gave more evenly distributed deviations between individual data sets and model results
than using other weighting options for wi,j, such as the ratioof means or standard deviations of data sets. However,improved automated weighting procedures should be developed in future. The Levenberg-Marquardt algorithm was applied to minimize the objective function (Eq. [11]) (Šimůnek
et al., 1998).
For the case of drainage of the soil column initially completely saturated up to the top, a modification of the above
procedure became necessary to properly describe the gradual,
noninstantaneous decrease of positive initial pressure heads(see results). At the bottom of the matrix and PFP domains, athin layer corresponding to the bottom perforated-plate thickness (L 5 0.4-cm height) with a constant hydraulic resistancefactor of R 5 L/Ks was assumed. The factor R accounted forany potential hydraulic resistance of the bottom plate withnylon membrane during the early phase of drainage, characterized by positive pressure heads and large hydraulic gradients across the lower boundary. For technical reasons, R
could not be measured below matrix and PFP, and thus for thedrainage case, R for matrix and PFP had to be inversely estimated. Allowing for different magnitudes of R at the bottomof matrix and PFP was justified by the fact that the perforated-

Field-Scale Preferential Flow and Bromide Transport
Two inverse DPM approaches were compared to simulate
the field experiments of Lennartz et al. (1999). In the sequential DPM approach, eight hydraulic parameters (i.e., am,
nm, Km,s, uf,s, af, nf, Kf,s, vw) were estimated based on observations of tile-drainage outflow (both in terms of rate andcumulative values), followed by the estimation of one to threesolute transport parameters for different scenarios (vs; or vs,lf; or vs, lf, lm) based on the observed Br2 concentrations
in tile-drainage effluent. In the simultaneous DPM approach,all eight hydraulic and one to three solute parameters were
estimated at the same time using both tile-drainage and Br2
concentrations. In both approaches, the parameters were
estimated by minimizing the objective function (Eq. [11]) as
described above for simulations of soil column tests. The following parameters were fixed: the volumetric portion of thefracture pore system (wf 5 0.05), as estimated from dyestained soil fractions (Köhne, 1999), the residual water content
of the fracture pore system (uf,r 5 0), residual and saturated
matrix water contents (um,r 5 0.078, um,s 5 0.43) of a loam aspredicted by the neural network prediction option in
HYDRUS-1D (Šimůnek et al., 1998), and the matrix
dispersivity at an average magnitude (lm 5 0.5 cm) found inBr2 displacement tests conducted with Bokhorst soil columns(Meyer-Windel, 1998).

Goodness-of-Fit
The goodness-of-fit of the three parameter modelingschemes for preferential flow (inverse-local, inverse-lumped,and forward) was evaluated using the modified index ofagreement (MIA) as given in Legates and McCabe (1999).
For each measurement set j, the index MIAj was calculated
as follows:

O[|E (z,t ,b)2O (z,t )|]MIA 5 1:02
[12]
O[|E (z,t ,b)2O|]1O[O (z,t )2O]
n

j

i

j

i

i51

j

n

n

j

i51

i

j

i

i51

The MIA (here: MIAj) varies between 0.0 to 1.0, withhigher values indicating better agreement between the modeland observations, similar to the interpretation of the coefficient of determination, R2. Compared with R2, MIA is a morereliable statistical measure since it is more sensitive to differences in the observed and model simulated means and variances and is less sensitive to extreme values (Legates andMcCabe, 1999). There is no statistical basis to decide which
MIA value is a threshold separating “inaccurate” from

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

64

VADOSE ZONE J., VOL. 5, FEBRUARY 2006

“accurate” simulations. Here we assumed all column-scalewater flow simulations with MIA . 0.75 to be “accurate.”
However, this value would be too rigorous for assessing fieldscale preferential solute transport simulations, for which othercriteria are usually applied, such as the percentage of thesimulation results within a certain deviation from observations(e.g., Larsson and Jarvis, 1999a, 1999b). Hence, a common
MIA threshold value to separate accurate from inaccurate hydraulic and transport simulations at field scale is not applicable. However, visual comparison of the agreement between
model result and data may provide additional information(Green and Stephenson, 1986). Therefore, for the field-scale
simulations, MIA values were compared only by relative ranking, and additionally two subjective binary (yes/no) goodnessof-fit measures were employed based on visual comparison ofobserved and simulated Br2 breakthrough curves: (i) match offirst preferential Br2 peak and (ii) match of breakthroughcurve pattern (e.g., peaked vs. smooth, skewed vs. symmetric).

RESULTS AND DISCUSSION
Homogeneous Soil Columns
The experimental water flow observations and corresponding inverse simulation results for the Richards’
equation using inverse-HYDRUS-1D are documented
briefly (Fig. 1 and 2). The corresponding estimated van
Genuchten’s (1980) soil hydraulic parameters that were
subsequently used as region-specific parameters in the
DPM forward approach are given in Table 2. Figure 1
shows results of infiltration and drainage experiments
conducted in the loam column, such as cumulative inflow and outflow (Fig. 1a and 1b), pressure heads (Fig. 1cand 1d) and water contents (Fig. 1e and 1f). Figure 2
presents experimental and inverse simulation results for
infiltration into the medium-sand column and for drain-

age of the initially saturated coarse-sand column at three
applied bottom suctions of 0, 230, and 260 cm in sequence. For technical reasons, water content measurements could not be retrieved during the drainage
experiment in the coarse sand except for the saturated
initial water content (Fig. 2f). All simulations were
accurate according to the MIA values exceeding 0.75.

Infiltration into the Loam Column with
Medium-Sand Preferential Flow Path
The case of infiltration into the loam column containing a PFP of medium sand was selected to represent an
example of relatively small hydraulic nonequilibrium or
low degree of preferential flow. Figure 3 shows experimental and simulated results for the infiltration into theloam column with medium-sand PFP. The observed total
infiltration of 1067 cm3 of water (Fig. 3a) exceeded thetotal outflow of 665 cm3 (Fig. 3b) by 400 cm3; that is,there was an increase in water storage of about 40%.
The cumulative region-specific outflow via the matrixof 380 cm3 (Fig. 3c) and via the medium-sand PFP of
230 cm3 (Fig. 3d) was equivalent to a 60-fold local water
flux out of the PFP as compared with the matrix, considering the 100:1 volume ratio of matrix to PFP. The
corresponding larger Ks ratio of PFP to matrix materialsof 350 (609/1.7, see Table 2) suggests that outflow outof the PFP was considerably reduced by lateral watertransfer from the PFP into the matrix.
Both the inverse-local (Table 3) and inverse-lumped(Table 4) DPM approaches accurately matched the bulk
soil related observations of cumulative infiltration, qin(Fig. 3a), and cumulative outflow, qout (Fig. 3b), as revealed by the large MIA values (0.95–0.97, Table 5). The

Fig. 1. Experimental data and inverse HYDRUS-1D simulation results for the loam column, infiltration (left) and drainage (right): (a, b) cumulative
infiltration and outflow, (c, d) pressure heads at the 10-, 20-, and 30-cm depths, and (e, f ) water contents at the 10- and 30-cm depths.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

65

Fig. 2. Experimental data and inverse HYDRUS-1D simulation results for the infiltration into the medium-sand (left) and coarse-sand(right) columns: (a, b) cumulative infiltration and outflow, (c, d) pressure heads at 10-, 20-, and 30-cm depths, (e, f ) water contents at the 10- and30-cm depths.

forward DPM simulation also accurately (MIA 5 0.92)
predicted the measured outflow (Fig. 3b) without any
calibration, while the MIA for infiltration (0.72) was
only just below the assumed threshold value of 0.75
for an accurate simulation. The region-specific outflow
observations for the matrix (Fig. 3c) and PFP (Fig. 3d)regions were accurately described by inverse-local(MIA » 0.95) and forward DPM (MIA . 0.75) simula-Fig. 3. Infiltration into the loam column with a preferential flow path
(PFP) made of medium sand. Comparison of experimental
cumulative flow data and dual-permeability model simulations
using independently estimated parameters (forward), or inversely
identified parameters based on data of pressure heads and watercontents in the matrix, and domain-specific (inverse-local) or bulksoil related (inverse-lumped) outflow. (a) Infiltration, (b) outflow,(c) outflow from the matrix, and (d) outflow from the PFP.

tions, while the inverse-lumped approach overestimatedoutflow out of the matrix, qm,out, and compensated thisoverestimation by a corresponding underestimation of
outflow out of the PFP, qf,out (see Fig. 3c and 3d), yielding
MIA values for qm,out and qf,out far below 0.75 (Table 5).
Figure 4 compares experimental and simulation
results of pressure heads in the matrix and in the PFP
during the first 1.5 h at the 5- (Fig. 4a), 35- (Fig. 4b), and65-cm (Fig. 4b) depths. The pressure heads in the matrix
and the PFP were almost similar for a particular depth,
revealing the low degree of hydraulic nonequilibrium.
Both the inverse-local and the forward DPM gave a
succession of steeper pressure head increases withdepth, which suggested a somewhat sharper infiltration
front traveling through the column than would beassumed for more gradually increasing pressure headdata. Furthermore, the forward simulation gave pressure heads in the PFP that increased ahead of those in
the matrix. In the inverse-local simulation, this was thecase at 5 cm below the top of the column, while at the
65-cm depth the trend was reverse (i.e., matrix pressure
heads increased earlier than PFP). The reason for thisdiscrepancy is the difference between region-specifichydraulic conductivity functions: for inverse-local, thecrossover-point where the PFP hydraulic conductivity
falls below that of the matrix was at h 5 219 cm, while
for forward it was at h 5 2120 cm (Table 2). Consequently, in the inverse-local simulation the matrix had
a higher hydraulic conductivity in part of the unsaturated pressure head range observed in the experiment(240 cm , h , 0).
Although the inverse-local and the forward DPM
simulations were in qualitative agreement with pressure
heads (Fig. 4), the MIA values below 0.75 indicated thatpressure heads were not accurately described by any
approach. One reason for the deviations between model

66

VADOSE ZONE J., VOL. 5, FEBRUARY 2006

Table 3. Soil hydraulic van Genuchten’s (1980) parameters and their 95% confidence intervals obtained in the inverse dual-permeability
model simulation using region-specific water flow data (inverse-local DPM approach).
urFlow region

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

cm3 cm23n

usa

cm3 cm23

cm21
Infiltration (medium-sand PFP)

Ks

R

cm d21d

PFP
Matrix0†0.184 6 0.031

0.350‡ 6 0.
0.313 6 0.000

0.0969 6 0.0223
0.0059 6 0.0007
Infiltration (coarse-sand PFP)

2.85 6 1.01
2.84 6 0.39

550.5 6 22.3
3.9 6 0.2

n.a.¶
n.a.¶

PFP
Matrix0†0.000 6 0.044

0.380‡ 6 0.277
0.308 6 0.002

0.0554 6 0.1471
0.0026 6 0.0004
Drainage (coarse-sand PFP)

3.34 6 6.19
1.94 6 0.354343 6 71
3.1 6 0.4

n.a.¶
n.a.¶

PFP
Matrix0†0.230‡ 6 0.159

0.368 6 0.124
0.392 6 0.001

0.0873 6 0.0382
0.0053 6 0.0038

3.00§ 6 1.79
1.25§ 6 0.1438621 6 56567
16.6 6 9.7

0.02 6 0.10
7.15 6 31.43† Fixed (not optimized).
‡ At upper constraint.
§ At lower constraint.
¶ Not available.

sharp contrast with the inverse-local approach clearlymatching the matrix water contents (Fig. 5), such thatthe inadequacy of the MIA for this type of data scattering around an almost constant level became evident.
The inverse-lumped DPM simulation results (notshown) generally showed slightly stronger deviations
from the pressure head and water content data than the
inverse-local and the forward DPM approaches. The
results were omitted in Fig. 4 and 5 since the inverselumped approach already failed to describe the regionspecific outflow.

simulations and data may be the inherent inaccuracy of
the first-order water exchange term (Köhne et al., 2004).
Additionally, region-specific pressure head measurements may not have been as accurate as flow measurements. However, we note that matrix pressure heads(and matrix water contents shown below) in this study
are regarded as auxiliary information only, the primary
information being the infiltration and outflow, since
outflow was also the only hydraulic information available and used in the field-scale simulations.
Water contents in the matrix and in the PFP at the
35-cm depth are shown in Fig. 5. The matrix water contents stayed close to saturation without showing any
clear reaction to infiltration. The water contents in the
PFP increased even more gradually than pressure heads
at the 35-cm depth (see Fig. 4 and 5). However, it cannot
be entirely ruled out that part of the TDR coil probe was
not placed in the PFP, which would have resulted in ameasured water content averaged over the PFP andadjacent matrix. The forward simulation resulted in
the best prediction of water contents in the PFP. The
inverse-local approach produced a relatively sharp PFP
water content rise that was not characteristic of the observed PFP water contents (Fig. 5). The assumption ofzero residual water content of the medium-sand PFP,
made to reduce the number of unknowns, may not havebeen realistic. The small MIA values (Table 5) were in

Infiltration into the Loam Column with
Coarse-Sand Preferential Flow Path
The infiltration into the loam column containing a PFP
filled with coarse sand constitutes an example of relatively strong hydraulic nonequilibrium or pronounced
preferential flow. Figure 6 shows the observed and simulated cumulative flow at the top and bottom of thecolumn. Infiltration (Fig. 6a) and outflow out of bulk soil
(Fig. 6b) as well as outflow out of the PFP (Fig. 6d) weremuch faster, while outflow out of the matrix (Fig. 6c) was
less than for the column with medium-sand PFP (Fig. 3).
In this way, as expected, fast preferential flow in thecoarse-sand PFP was relatively less impeded by transfer
into the matrix than flow in the medium-sand PFP.

Table 4. Soil hydraulic van Genuchten’s (1980) parameters and their 95% confidence intervals obtained in the inverse dual-permeability
model simulations based on bulk soil water flow data (inverse-lumped DPM approach).
urFlow region

us
3

cm cm

23

a

n

21

cm
Infiltration (medium-sand PFP)

Ks
cm d

R
21

d

PFP
Matrix0†0.231 6 0.053

0.350‡ 6 0.120
0.314 6 0.001

0.1327 6 0.0777
0.0062 6 0.0028
Infiltration (coarse-sand PFP)

2.76 6 1.76
1.77 6 0.22

138.5 6 107.9
8.4 6 1.1

n.a.¶
n.a.¶

PFP
Matrix0†0.009 6 0.243

0.339 6 0.202
0.308 6 0.002

0.0171 6 0.0104
0.0030 6 0.0017
Drainage (coarse-sand PFP)

2.00 6 0.80
1.66 6 0.443380 6 447
11.6 6 4.6

n.a.¶
n.a.¶

PFP
Matrix0†0.230†

0.380‡ 6 0.075
0.390†

0.2000‡ 6 0.1806
0.0044 6 0.0119

3.00‡ 6 2.88
1.25†24290 6 43891
15.6 6 3.7

0.0055 6 0.006
2.37 6 1.77

† Fixed (not optimized).
‡ At upper constraint.
§ At lower constraint.
¶ Not available.

www.vadosezonejournal.org

67

Table 5. Modified index of agreement, MIA (Eq. [12]), for differences between various model results and experimental
observations.†
MIA
Model approach

qin

qout

qf,outqm,outhm

um

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

Infiltration (medium-sand PFP)
Inverse-local
Inverse-lumped
Forward

0.95‡
0.97‡
0.96‡
0.95‡
0.95‡
0.95‡
0.44
0.36
0.72
0.92‡
0.78‡
0.93‡
Infiltration (coarse-sand PFP)

0.38
0.36
0.38

0.22
0.18
0.06

Inverse-local
Inverse-lumped

0.97‡
0.96‡
0.96‡
0.07
0.97‡
0.96‡
0.80‡
0.01
Drainage (coarse-sand PFP)

0.78‡
0.79‡

0.81‡
0.80‡

Inverse-local
Inverse-lumped
Forward

n.a.§
n.a.
n.a.

0.89‡
0.88‡
0.40

0.48
0.44
0.04

Fig. 5. Same as Fig. 3, but showing water contents in the matrix and
preferential flow path (PFP) at the 35-cm depth.

† qin, cumulative infiltration; qout, qf,out, qm,out, cumulative outflow out of bulk
soil, preferential flow path, matrix; um, matrix water contents at 35 cm; hm,
matrix pressure heads at 5, 35, and 65 cm.
‡ MIA . 0.75 (i.e., ‘‘accurate’’ model simulation).
§ Not available.

example, an R2 value of 0.74 corresponds to the abovementioned MIA value of 0.07.
We did not include results of the forward DPM simulation because independently obtained parameters wereaffected by finger flow phenomenon in the homogeneous coarse-sand column, and we inferred that theseparameters would not be representative of the PFP.
The observed and simulated pressure heads in matrix
and coarse-sand PFP are shown in Fig. 7 for the depthsof 5 (Fig. 7a) and 35 cm (Fig. 7b). The pressure heads
showed more pronounced differences between the(coarse-sand) PFP and matrix than was observed between
the matrix and a PFP made of medium sand (Fig. 4). This

0.92‡
0.97‡
0.36

0.87‡
0.52
0.55

0.91‡
0.47
0.33

As evident from MIA values (Table 5) and visual assessment (Fig. 6), the inverse-local DPM simulations
once again gave an accurate description of both global(bulk soil related) and local (PFP) flows, and theinverse-lumped DPM approach provided a similar
match of global flow data, while PFP flow data werematched less. For the cumulative matrix outflow exhibiting a very flat increase, the low MIA value of 0.07
suggested failure even for the inverse-local simulation,
which according to visual evaluation was only slightly
overestimating the outflow (Fig. 6). Again, as alreadyfound above for the water content data, a conclusion
about goodness-of-fit of one subset of observations
should not be drawn from a statistical measure like the
MIA alone. We note also that R2 values (not shown)
were consistently higher than MIA values, which showsthat the MIA is a more rigorous statistical measure. For

Fig. 4. Same as Fig. 3, but showing pressure heads in matrix and preferential flow path (PFP) at depths of (a) 5, (b) 35, and (c) 65 cm.

Fig. 6. Infiltration into the loam column with the preferential flow
path (PFP) made of coarse sand. Comparison of experimental
cumulative flow data and dual-permeability model simulations
using independently estimated parameters (forward), or inversely
identified parameters based on data of pressure heads and watercontents in the matrix, and domain-specific (inverse-local) or bulksoil related (inverse-lumped) outflow. (a) Infiltration, (b) outflow,(c) outflow from the matrix, and (d) outflow from the PFP.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

68

VADOSE ZONE J., VOL. 5, FEBRUARY 2006

Fig. 7. Same as Fig. 6, but showing pressure heads at depths of (a) 5
and (b) 65 cm.

finding of stronger pressure head nonequilibrium can be
explained with the larger contrast between hydraulicproperties of loam and coarse-sand PFP than between
hydraulic properties of loam and medium-sand PFP.
Although within the pressure head range shown in Fig. 7
the inverse-local approach apparently better describedmatrix pressure head observations than the inverselumped approach, both gave accurate representations
according to similar MIAs of about 0.78 (Table 5).
Moreover, both inverse-local and inverse-lumped simulations slightly overestimated water contents in the matrix
at the 5-cm depth and underestimated them at the 35-cm
depth (Fig. 8).

Drainage of the Loam Column with Coarse-Sand
Preferential Flow Path
Finally, the drainage of the initially entirely saturated
loam column containing a PFP made of coarse sand was
conducted to see if findings for the infiltration cases
could be generalized to drainage flow conditions as well.
Findings are presented for a two-step drainage outflow
experiment. Figure 9 reveals that cumulative water outflow of the bulk soil (Fig. 9a), soil matrix (Fig. 9b), and
PFP (Fig. 9c) increased faster initially and showed slow
response to the step change of suction at the lower
boundary from 10 to 45 cm at 138 min. The inverse-local

Fig. 8. Same as Fig. 6, but showing matrix water contents at depths of
5 and 35 cm.

Fig. 9. Drainage from the loam column with the preferential flow path
(PFP) made of coarse sand. Comparison of experimental cumulative outflow data and dual-permeability model simulations using
independently estimated parameters (forward), or inversely identified parameters based on data of pressure heads and watercontents in the matrix and domain-specific (inverse-local) or bulksoil related (inverse-lumped) outflow. (a) Outflow, (b) outflow fromthe matrix, and (c) outflow from the PFP.

DPM approach again provided an accurate match ofboth lumped outflow and local outflow components,
yielding MIA values .0.75 (Table 5). The inverselumped DPM matched lumped outflow even better
according to the MIA (0.97 vs. 0.92), but once again
failed to describe local outflow components of the matrix and PFP (MIA , 0.75, Table 5). The forward DPM
approach underestimated outflow out of the matrix(Fig. 9b) and grossly underestimated outflow out of thecoarse-sand PFP (Fig. 9c). As a result of this behavior,outflow out of bulk soil was also underestimated by theforward DPM approach (Fig. 9a), and low MIA values
were obtained (Table 5).

Fig. 10. Same as Fig. 9, but showing pressure heads in the matrix andin the preferential flow path (PFP) at depths of (a) 5 and (b) 65 cm.

69

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

Fig. 11. Same as Fig. 9, but showing water contents in the matrix and in
the preferential flow path (PFP) at the 35-cm depth.

The pressure heads in the matrix and in the PFP at the
5-cm depth (Fig. 10a) decreased starting from slightly
negative values. Pressure heads at the 65-cm depth
started from initially positive values, decreased within
6 min (PFP) and 20 min (matrix) to negative values, andshowed a reaction to the suction step increase at 138 min(Fig. 10b). Differences between pressure heads in PFP
and matrix regions were less than for the infiltration case(compare Fig. 7 and 10), which reveals that the degree ofpressure head nonequilibrium during preferential flow
not only depends on the difference of hydraulic properties between PFP and matrix, but also on the flow direction (infiltration vs. drainage). The inverse-local andinverse-lumped simulations described pressure heads
accurately (MIA above 0.75; see Fig. 10) when accounting for a hydraulic resistance at the bottom. Pressure
heads simulated with the forward DPM approach did
not match with the observations (Fig. 10, Table 5).
Water contents at the 35-cm depth showed hardly any
reaction with time for the matrix, but a strong decreasein the water content of the PFP (Fig. 11) occurred during
the onset of drainage. According to visual evaluation,
the inverse models gave a satisfactory representationof the PFP water contents and excellent description of
matrix water contents not reflected in the MIA (Table 5),
while the forward simulation results deviated relatively
more from the experimental observations (Fig. 11).

Dual-Permeability Model Parameters
and Their Uncertainty
Model parameters and their 95% confidence intervals
for the forward, inverse-local, and inverse-lumped approaches are listed in Tables 2, 3, and 4. Results forsteady-state Ks and static retention curve measurements
are given for comparison (Table 2). The differences between parameter magnitudes for the various approaches
reflect the observed deviations among simulation results, which will be discussed for the Ks values as an
example. For the infiltration experiments, Ks for loam
was lower for the inverse-local approach, with values of
3.9 and 3.1 cm d21 for columns with coarse- and mediumsand PFPs, respectively (Table 3), than corresponding
loam Ks values for the inverse-lumped approach of 8.4

and 11.6 cm d21, respectively (Table 4). The large matrix
Ks values obtained with the inverse-lumped approach
resulted in the overestimation of the observed cumulative matrix outflow (Fig. 3 and 6). Accordingly, thevalues obtained in the inverse-local approach comparedbetter, with Ks of 1.7 cm d21 obtained for infiltration into
the homogeneous loam column. Assuming Ks equal to
1.7 cm d21 in the forward DPM approach still produced
accurate prediction of matrix outflow (Fig. 3). On theother hand, once again for the infiltration experiments,
the inverse-local approach gave higher Ks estimates for
mediums and coarse-sand PFPs of 550 and 4343 cm d21,
respectively (Table 3), than the inverse-lumped approach (138 and 3380 cm d21, Table 4), which for the
inverse-lumped approach resulted in the underestimation of the observed cumulative PFP outflow (Fig. 3 and
6). Consistent with these observations, for the mediumsand PFP the Ks of 550 cm d21 of the inverse-local approach was closer to Ks obtained for infiltration in themedium-sand column of 609 cm d21, than Ks of 138 cm
d21 for the inverse-lumped approach. Moreover, Ks
confidence limits obtained for infiltrations in matrix
and PFP were consistently smaller for the inverse-local
than inverse-lumped approach (Tables 3, 4). Furthermore, the correlation, r, between Ks,f and Ks,m values in
PFP and matrix was small for the inverse-local approach(r 5 20.33), but an almost perfect negative correlation(r 5 20.99) was obtained for the inverse-lumped approach (Table 6). This highly negative correlation of theinverse-lumped approach means that a higher Ks,f estimate is compensated by a lower Ks,m value and viceversa. In other words, Ks,f and Ks,m cannot be inversely
identified at the same time from bulk soil water flow
observations. In total, four “large” correlations (r .
0.75) were identified for the inverse-lumped approach,but only two values of r . 0.75 for the inverse-local
approach, between a and n of both matrix and PFP
Table 6. Correlation, r, between parameters for inverse dualpermeability model simulations (inverse-lumped and inverselocal approaches) of infiltration into loam with coarse-sand
PFP.
um,r

um,s

am

nm

Ks,m

uf,s

af

nf

Ks,f

Inverse-lumped
1
um,r
0.01
1
um,s
0.99‡ 0.00 1
am
nm0.28 20.04 0.39
0.01 0.10
Ks,m 0.13
uf,s
0.12 20.02 0.08
af 20.19 20.04 20.29
nf
20.05
0.00 20.19
Ks,f 20.16 20.01 20.13

1
20.57† 1
0.26 20.37
1
20.42
0.01 20.26
20.75‡ 0.23 20.20
0.57 20.99‡ 0.35
Inverse-local

1
um,r
1
um,s 20.06
0.08 1
am 20.41
nm 20.64† 0.03 0.81‡ 1
0.05 20.11 20.55† 1
Ks,m 0.18
uf,s
0.47 20.12 20.74† 20.45 20.30
0.02 0.29
0.25
0.04
af 20.14
20.09
0.00 0.19
0.19
0.02
nf
Ks,f 20.14
0.01 0.17
0.28 20.33† |r| . 0.5.
‡ |r| . 0.75.

10.84‡ 120.01 20.27

1

1
20.64† 1
20.54† 0.99‡ 1
0.07 20.12 20.21

1

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

70

VADOSE ZONE J., VOL. 5, FEBRUARY 2006

(Table 6). However, in the case of empirical curve-fit
parameters without much physical meaning, such as aand n, correlations may be acceptable if only the shapeof the estimated hydraulic functions is of interest, which
may be similar for different combinations of such correlated parameters.
All the above comparisons together suggest that theinverse-local approach not only accurately described
preferential water flow, but that the resulting inverse Ks
estimates (and other parameters except a and n) werealso generally well identified and physically meaningful,while all this was not the case for the inverse-lumped
approach. This is a key finding of this study, the practical
implications of which are shown below for DPM prediction of field-scale preferential solute transport basedon tile-drainage hydrographs.
For the drainage experiments and both inverse DPM
approaches, Ks,f and Ks,m and us,m were larger than for
infiltration (Tables 3 and 4). This discrepancy can be
explained when considering the two wetting fronts in thematrix during preferential infiltration, that is, one fromthe top and one from the side by water transfer from the

PFP. Hence, beneath the vertical wetting front in thematrix, there were coexisting wet (adjacent to PFP) anddry (far from PFP) regions, which resulted in lower Ks
and us estimates, as compared with values for waterdrainage coming from the total matrix cross section. Additionally, lateral water transfer during infiltration may
have entrapped air below the matrix wetting front, also
resulting in lower Ks and us estimates for the matrix. In
summary, the above findings are a manifestation of enhanced soil hydraulic hysteresis during preferential flow.
Indications of hysteresis during nonequilibrium flow
have been found before (e.g., Šimůnek et al., 2001;
Köhne and Mohanty, 2005). The relations between Ks
values found for the various approaches were not asclear for drainage as was found for infiltration, which we
attribute to the difficulties in modeling the observed
gradual positive pressure head dissipation, requiring theestimation of lower boundary resistivity factors formatrix and PFP. Moreover, in the drainage case, thehydraulic resistance factor greatly increased the uncertainty associated with the Ks values of matrix and PFP, asapparent from large confidence intervals (Tables 3 and 4).

Fig. 12. Soil hydraulic van Genuchten functions fitted to the soil water retention curve measured at equilibrium (static) and obtained by fitting thenumerical solution of Richards’ equation to water flow measurements performed on homogeneous soil materials (infiltration and drainage), andobtained by fitting the dual-permeability model to flow measurements conducted with heterogeneous soil systems consisting of loam (matrix) andmedium or coarse sand (preferential flow path, PFP): medium sand (a) water content and (b) hydraulic conductivity, coarse sand (c) watercontent and (d) hydraulic conductivity, and loam (matrix) (e) water content and (f) hydraulic conductivity.

71

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

Fig. 13. Rainfall, drain discharge flux, and cumulative drain outflow as observed at the Bokhorst field site for three runoff seasons between 1991 and1995 (data: Lennartz et al., 1999), and inverse dual-permeability model (DPM) simulation results for the sequential and simultaneous approaches.

For comparison, the magnitudes of Ks measured atsteady state for loam and medium and coarse sand werehigher than measured under transient flow conditions inhomogeneous columns, but the corresponding 95% Ks
confidence limits were still overlapping. The us values
obtained in static measurements were generally between those obtained for infiltration and drainage
experiments (Tables 2–4).
Those parameter values that were strongly correlated
are no longer amenable to interpretation and comparison, while the shapes of hydraulic functions, which may
be similar for different values of correlated parameters,
can still be compared across model approaches. Figure 12
compares water retention and hydraulic conductivity
functions obtained for homogeneous columns and steadystate or static conditions with those obtained for the matrix and PFP. For medium sand, inverse-local and-lumped
produced similar PFP retention functions (Fig. 12a) andhydraulic conductivity functions (Fig. 12c), but thesefunctions were higher close to saturation and lower in thedrier range than those obtained for infiltration into themedium-sand column. For coarse sand, the inverselumped water retention (Fig. 12b) and hydraulic conductivity (Fig. 12d) functions were vastly different, which
was unrealistic considering that retention functions obtained by static and infiltration measurements were quite
similar (Fig. 12b). For the loam, retention functions of all
approaches were almost parallel and revealed that estimated different ur values (and maybe n and a values)
obtained in different approaches are not necessarily
meaningful and probably cannot be used outside theexperimental pressure head range near saturation. The
different (convex or concave) shapes of the loam hydrau-

lic conductivity function for infiltration and drainage
were better reproduced by the inverse-local approach,disregarding differences in the absolute values (Fig. 12f).

Simulation of Field Tracer Tests
Results for simulated drain discharge rates and cumulative drain discharge using the sequential and simultaneous DPM simulation schemes are compared inFig. 13 with corresponding experimental observations
published by Lennartz et al. (1999). For completeness,
Table 7. Statistical goodness-of-fit: modified index of agreement,
MIA (Eq. [12]), for differences between model results andobservations (Lennartz et al., 1999) of drain discharge flux
(qout), cumulative drain discharge (qcumout), and Br2 concentrations in tile drainage water [c(Br2)]; and success (Yes) or
failure (No) of Br2 transport simulation based on visual comparison of observed and simulated Br2 breakthrough curves(Fig. 14 and 15); match of first preferential Br2 peak (Firstpeak); and match of general breakthrough curve pattern withregard to shape and secondary peaks (Pattern).
MIA†
Exp.

qout

1991/92
1993/94
1994/95

0.83
0.64
0.70

1991/92
1993/94
1994/95

0.71 (0.76)
0.61 (0.66)
0.65 (0.52)

qcumoutc(Br2)First peak

Patternc(Br2)c(Br2)No (No)No (Yes)
No (Yes)No (No)No (Yes)No (No‡)Yes (Yes)
Yes (Yes)
Yes (Yes)

No (Yes)
No (Yes)
No (Yes)Sequential
0.97
0.41 (0.42)
0.96
0.29 (0.75)
0.98
0.31 (0.61)
Simultaneous
0.97 (0.97)
0.96 (0.96)
0.91 (0.95)

0.61 (0.63)
0.40 (0.76)
0.46 (0.62)† Values in parentheses: MIA for simulations with fitted dispersivities.
‡ Flat matrix peak’’ of Br2 between 100 and 150 d was not approximated.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

72

VADOSE ZONE J., VOL. 5, FEBRUARY 2006

Fig. 14. Bromide tracer concentrations as observed in drain outflow at
the Bokhorst field site for three runoff seasons between 1992 and1995 (data Lennartz et al., 1999) and inverse dual-permeability
model (DPM) simulation results for the sequential and simultaneous approaches. The only parameter fitted in the solute transport
equation was the first-order solute transfer coefficient.

the natural daily rainfall used as the upper boundary
condition is provided as well (Fig. 13). The drain discharge rates closely reflect the rainfall pattern and werereproduced by both inverse DPM schemes, the simultaneous scheme yielding somewhat higher and narrower
discharge peaks than the sequential scheme. The resulting MIA values were only slightly larger for the sequential approach (Table 7). Cumulative drain discharge isreproduced equally well with both schemes, according to
MIA values around 0.97, except for 1994/1995 where thesequential approach yielded an MIA of 0.98, as compared with 0.91 for the simultaneous approach.
The observed and simulated Br2 concentrations areshown in Fig. 14 for the case with dispersivities in matrix, lm, and PFPs, lf, fixed at 0.5 cm as derived from
Br2 displacement experiments using soil columns(Meyer-Windel, 1998). The estimated soil hydraulicand transport DPM parameters are listed in Table 8. In
all 3 yr, there was a recurrent pattern of rapid Br2 concentration increase to the main preferential Br2 peakdirectly following the first rainfall event after solute
application, and secondary smaller Br2 concentration
peaks in response to later rainfall events. Only toward
the end of the 1994/1995 season, was there a slight increase in Br2 concentration considered the matrix Br2
peak (Fig. 14). Strikingly, only the simultaneous DPM
approach consistently reproduced the early preferential
peak for all three tracer tests, whereas the sequential
scheme overestimated the time lag between the rainfall
event and first peak, underestimated peak heights, andproduced no peak at all in the 1994/1995 runoff season(Fig. 14). Accordingly, MIA values for Br2 concentrations were consistently higher for the simultaneous
DPM approach (Table 7). However, neither approach
reproduced the breakthrough pattern of secondary
peaks or the matrix Br2 peak in the last season. Therefore, another set of simulations was run where lf wasadditionally inversely estimated (Fig. 15, Table 9). The
simultaneous DPM approach with fitted lf not only
simulated the initial peak but also the breakthroughpattern (Fig. 15), while water flow simulations wereabout as successful as the simultaneous DPM simulations with constant dispersivities (Table 7). The flat matrix peak in 1994/1995 could also be simulated, although
this required estimating lm (Fig. 15). The sequential
scheme improved, but still did not approximate thematrix peak in 1994/1995, even when both lf and lm
were fitted, and it did not capture the main peak in the

Table 8. Soil hydraulic van Genuchten’s (1980) and transport parameters as obtained by inverse dual-permeability model (DPM)
simulations of three field experiments between 1991 and 1995 (data of Lennartz et al., 1999). The sequential approach relied on tiledrainage fluxes for estimating hydraulic parameters, followed by estimating the solute transfer coefficient using Br2 effluent
concentrations. The simultaneous approach utilized drainage fluxes and Br2 concentrations together for identification of hydraulic andtransport parameters.
Pore region

Exp.

ur

us
cm3 cm23

a

n

cm21
SequentialKs

vw 3 1023

vs 3 1023

cm d21

cm21 d21

d21

4.48

2.2

4.61

9.2

0.72

150.0

0.10

166.6

0.04

55.6

0.02

18.9

Fracture pore system
MatrixFracture pore system
MatrixFracture pore system
Matrix1991/92
1991/92
1993/94
1993/94
1994/95
1994/95

0†0.082†
0†0.082†
0†0.082†

0.856
0.49†
0.875
0.49†
0.765
0.49†

0.018
0.011
0.018
0.011
0.028
0.030
Simultaneous

1.83
1.90
1.84
1.90
1.40
1.30

100.3
4.27
100.2
4.23
46.3
1.0§

Fracture pore system
MatrixFracture pore system
MatrixFracture pore system
Matrix1991/92
1991/92
1993/94
1993/94
1994/95
1994/95

0†0.082†
0†0.082†
0†0.082†

0.966
0.49†
0.220
0.49†
0.454
0.49†

0.023
0.047
0.008
0.001
0.030
0.029

1.70
2.50‡
2.34
2.50‡
1.60
1.30

66.2
50.0‡
47.9
9.01
40.0
1.82

† Fixed (not optimized).
‡ At upper constraint.
§ At lower constraint.

73

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

Fig. 15. Same as Fig. 14, but in the solute transport equation, both thefirst-order transfer coefficient and the fracture dispersivity (1994/
1995: additionally the matrix dispersivity) were fitted.

1991/1992 season, as reflected by lower MIA values forthe sequential than for the simultaneous approach. In
summary, visual and statistical criteria indicated preferential solute breakthrough was consistently better described using the simultaneous DPM approach.
Model parameters are discussed for the improved
simulations with fitted dispersivities (Table 9), whereparameters were assumed to better reflect the fieldsite. A caveat is that one should not assign too much
physical meaning to these parameters because of several
factors, such as parameter uniqueness (briefly discussed
below) and the effective nature of one-dimensionalmodel parameters.
Soil hydraulic van Genuchten (1980) parameters forboth the sequential and simultaneous approaches were
not constant between years, but varied within a factor of
2 (Table 9). No single parameter set could be identified

that gave an equally good representation of all years(not shown), which points to seasonal variability of soil
properties. One reason could be the variability of Ap
horizon structure in different years, since tracer application was done at different times (1 to 2 mo) after tillage.
With increasing time after tillage, some soil structure
may have developed PFPs, and different crops developed different root systems in each year. Moreover,initial soil moisture was known for 1991/1992 only, andwas assumed to be similar in the following years. However, Kätterer et al. (2001) found that a single DPM parameter set resulted in successful simulation of various
tracer elution curves except for dry initial soil moisture
conditions, which was attributed to water repellency developing under dry conditions.
The estimated Ks, a, and n values were generally
larger for the PFPs than for the matrix; hence, both inverse approaches generally produced plausible parameter results. For the 1991/1992 data-set, the simultaneous
approach estimated small values for the first-order coefficients for water transfer, vw, and solute transfer, vs.
In the sequential approach, vw estimated from drain discharge was large, and vs subsequently estimated fromBr2 concentrations approached zero (Table 9). Overall(advective and diffusive) Br2 transfer was still overestimated in the sequential approach because of the large
vw, which was the main reason for the failure in describing the preferential Br2 peak (Fig. 15).
Fitted lf values increased above the initial values by atleast one order of magnitude (Table 9), which suggests
heterogeneity of hydraulic PFP properties at the field
scale. The matrix dispersivity was not sensitive (notshown) except for the 1994/1995 breakthrough curvewith flat matrix peak, where a 10-fold increase of lm
fitted in the simultaneous DPM approach (lm 5 5.8 cm,
Table 9) was consistent with the large heterogeneityof soil texture observed at the Bokhorst field site
(Lennartz et al., 1999). In the sequential approach, lm
was at the upper constraint (100 cm), but still the matrix
Br2 peak was not approximated. The above findings
indicate that the failure to identify soil hydraulic func-

Table 9. Similar results as in Table 8, except that the dispersivity l was fitted for the fracture pore system (1991/92, 1993/1994 data-sets), or
for matrix and fracture pore systems (1994/95).
Pore region

Exp.

ur

us
cm3 cm23

a

n

Ks

cm21
Sequential

vw 3 1023l

vs 3 1023

cm21 d21

cm

d21

14.2
0.5†
100‡
0.5†
37.0
100‡

0.001§

Fracture pore system
MatrixFracture pore system
MatrixFracture pore system
Matrix1991/92
1991/92
1993/94
1993/94
1994/95
1994/95

0†0.082†
0†0.082†
0†0.082†

0.856
0.49†
0.875
0.49†
0.765
0.49†

0.018
1.83
0.011
1.90
0.018
1.84
0.011
1.90
0.028
1.40
0.030
1.30
Simultaneous

100.3
4.27
100.2
4.23
46.3
1.0†

4.48

Fracture pore system
MatrixFracture pore system
MatrixFracture pore system
Matrix1991/92
1991/92
1993/94
1993/94
1994/95
1994/95

0†0.082†
0†0.082†
0†0.082†

1.0§
0.49†
0.830
0.49†
0.529
0.49†

0.023
0.022
0.018
0.011
0.038
0.029

41.1
23.8
114.3
4.23
43.5
5.62

0.19

† Fixed (not optimized).
‡ At upper constraint.
§ At lower constraint.

1.50
2.50‡
2.12
1.77
2.82
2.30

4.61
0.72

4.81
2.51

20.4
0.5†
100‡
0.5†
100‡
5.78

6.1
23.9

1.0
17.6
7.2

74

VADOSE ZONE J., VOL. 5, FEBRUARY 2006

Table 10. Correlation, r, between parameters for inverse dual-permeability model (DPM) simulation of tile-drain outflow and Br2 effluent
concentrations as observed in 1994/95 (Lennartz et al., 1999). The sequential and simultaneous approaches are defined in the captionof Table 7.
am

nm

Ks,m

uf,s

nf

Ks,f

1
20.96‡
0.71†

1
20.68†

af

vw

lm

vs

lf

1
20.03
0.01

1
0.61†

1

1
20.56†
20.30

1
0.02

1

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

Sequential
am
nm
Ks,m
uf,s
af
nf
Ks,f
vw
lm
vs
lf

1
20.58†
0.75‡
20.65†
0.50†
0.28
20.16
0.65†

1
0.02
0.45
20.42
20.34
0.42
20.79

1
20.31
0.33
20.13
0.28
0.06

1
20.78‡
20.76‡
0.72
20.62

1
0.47
20.49
0.28

1

Simultaneous
am
nm
Ks,m
uf,s
af
nf
Ks,f
vw
lm
vs
lf

1
20.59†
0.93‡
20.14
0.66†
20.51†
0.20
0.35
20.10
0.08
0.07

1
20.38
0.33
20.64†
0.48
20.11
20.49
20.32
0.65†
20.29

1
0.01
0.57†
20.51†
0.18
0.19
20.29
0.29
0.11

1
0.20
20.29
0.66†
20.26
20.42
0.59†
0.33

1
20.83‡
0.57†
0.01
20.09
20.05
0.15

1
20.26
0.34
0.35
20.13
20.44

1
0.25
20.08
0.11
20.05

1
0.38
20.53†
0.02† |r| . 0.5.
‡ |r| . 0.75.

tions relevant to preferential flow in the sequential approach could not be compensated by estimating soluteparameters lm, lf, or vs.
Even for the simultaneous approach, not all 11 DPM
model parameters can be expected to be uniquely identified. However, Table 10 shows for 1994/1995 that
correlations between any two parameters of the simultaneous approach were mostly small (r # 0.5 for 41 of
55, or 75%) or moderate (r # 0.75 for 53 of 55, or 96%),
with only two (4%) highly correlated parameters. For
the sequential approach, a larger number of water flow
parameters were highly correlated (r . 0.75 for 5 of 28,
or 18%).
These above findings demonstrate two practical implications of experimental and simulation results obtainedat the column scale. First, the DPM hydraulic parameters that are identified based on bulk soil water flow data
may lead to failure in preferential solute transport simulation (sequential approach). Second, when using bothwater flow and solute concentration information, DPM
hydraulic parameters related to preferential flow canbe inversely estimated to produce improved preferential
solute transport simulations which are in agreementwith observations (simultaneous approach).
Our results may also add some explanations for previous findings, such as of poor results for two-dimensional two-domain (mobile–immobile) model simulationof preferential solute transport based on hydraulic parameters that gave excellent description of observed
tile-drain discharge (Haws et al., 2005) or for the findingthat almost identical simulations of tile-drainage water
discharge fluxes were obtained using single- and dualpermeability models, but only the DPM could reproduce
observed preferential solute breakthrough (Gerke andKöhne, 2004).

CONCLUSIONS
Several controlled infiltration and drainage experiments were conducted using a soil column setup withdesigned preferential flow paths. Hydraulic nonequilibrium increased with the contrast of hydraulic properties
between the matrix and PFP and was more pronounced
for infiltration than drainage conditions. Using independently determined soil hydraulic parameters, a dual permeability model provided mostly accurate predictions ofbulk soil related and region-specific flow data according
to the statistical modified index of agreement. An accurate match of the experimentally measured local flow
data could only be achieved with the DPM when utilizing region-specific water flow data for inverse parameter identification, as done in the inverse-local DPM
approach. Data of infiltration and preferential outflow
from bulk soil only, even when supplemented withmeasured pressure heads and water contents in the soil
matrix, did not contain enough information for proper
identification of region-specific parameters. This wasdemonstrated with the inverse-lumped approach, which
reproduced bulk soil flow data accurately, but notregion-specific flow. A practical consequence of thisfinding is that for natural soil where region-specific data
are not available, hydraulic DPM parameters characterizing preferential flow cannot be properly identified
from water flow observations alone. This was further
demonstrated using sequential inverse DPM simulations
that reproduced field-scale observations of tile-draindischarge but could not consistently describe preferential bromide breakthrough as observed in three tracer
tests. However, the successful field-scale preferential
solute transport simulation for concurrent estimation of
hydraulic and solute transport parameters (simultaneous

www.vadosezonejournal.org

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

DPM approach) suggests that effluent solute concentrations contain indirect information about preferential
water flow that can be used to identify hydraulic DPM
parameters. The principal implications of this study are
not limited by the basic one-dimensional DPM approach
adopted here, but should hold even when applying
further refined model representations of real soils.
ACKNOWLEDGMENTS
We appreciate the support by the National Science Foundation (grant EAR 0296158).

REFERENCES
Abbasi, F., J. Šimůnek, J. Feyen, M.Th. van Genuchten, and P.J.
Shouse. 2003. Simultaneous inverse estimation of soil hydraulic andsolute transport parameters from transient field experiments: Homogeneous soil. Trans. ASAE 46:1085–1095.
Abbaspour, K.C., A. Kohler, J. Šimůnek, M. Fritsch, and R. Schulin.
2001. Application of a two-dimensional model to simulate flow andtransport in a macroporous agricultural soil with tile drains. Eur. J.
Soil Sci. 52:433–447.
Ahuja, L.R., K.E. Johnson, and G.C. Heathman. 1995. Macropore
transport of a surface-applied bromide tracer: Model evaluation
and refinement. Soil Sci. Soc. Am. J. 59:1234–1241.
Allaire, S.E., S.C. Gupta, J. Nieber, and J.F. Moncrief. 2002a. Role ofmacropore continuity and tortuosity on solute transport in soils: 1.
Effects of initial and boundary conditions. J. Contam. Hydrol. 58:
299–321.
Allaire, S.E., S.C. Gupta, J. Nieber, and J.F. Moncrief. 2002b. Roleof macropore continuity and tortuosity on solute transport in soils:
2. Interactions with model assumptions for macropore description.
J. Contam. Hydrol. 58:283–298.
Andreu, L., F. Moreno, N.J. Jarvis, and G. Vachaud. 1994. Applicationof the model MACRO to water movement and salt leaching indrained and irrigated marsh soils. Agric. Water Manage. 25:71–88.
Booltink, H.W.G., and J. Bouma. 1991. Physical and morphologicalcharacterization of bypass flow in a well-structured clay soil. Soil
Sci. Soc. Am. J. 55:1249–1254.
Castiglione, P., B.P. Mohanty, P.J. Shouse, J. Šimůnek, M.Th. van
Genuchten, and A. Santini. 2003. Lateral water diffusion in an artificial maroporous system: Modeling and experimental evidence.
Available at www.vadosezonejournal.org. Vadose Zone J. 2:212–221.
Gerke, H.H., and M.Th. van Genuchten. 1993a. A dual-porosity model
for simulating the preferential movement of water and solutes instructured porous media. Water Resour. Res. 29:305–319.
Gerke, H.H., and M.Th. van Genuchten. 1993b. Evaluation of a firstorder water transfer term for variably saturated dual-porosity flow
models. Water Resour. Res. 29:1225–1238.
Gerke, H.H., and M.Th. van Genuchten. 1996. Macroscopic representation of structural geometry for simulating water and solute masstransfer in dual-porosity media. Adv. Water Resour. 19:343–357.
Gerke, H.H., and J.M. Köhne. 2004. Dual-permeability modeling ofpreferential bromide leaching from a tile drained glacial tillagricultural field. J. Hydrol. (Amsterdam) 289:239–257.
Greco, R. 2002. Preferential flow in macroporous swelling soil withinternal catchment: Model development and applications. J.
Hydrol. (Amsterdam) 269:150–168.
Green, I.R.A., and D. Stephenson. 1986. Criteria for comparison of
single event models. Hydrol. Sci. 31:395–411.
Gwo, J.P., P.M. Jardine, G.V. Wilson, and G.T. Yeh. 1995. A multiplepore-region concept to modeling mass transfer in subsurface media.
J. Hydrol. (Amsterdam) 164:217–237.
Gwo, J.P., P.M. Jardine, G.V. Wilson, and G.T. Yeh. 1996. Using a
multiregion model to study the effects of advective and diffusive
mass transfer on local physical nonequilibrium and solute mobility
in a structured soil. Water Resour. Res. 32:561–570.
Haws, N.W., P.S.C. Rao, J. Šimůnek, and I.C. Poyer. 2005. Singleporosity and dual-porosity modeling of water flow and solutetransport in subsurface-drained field using effective field-scaleparameters. J. Hydrol. (Amsterdam) (in press).

75

Inoue, M., J. Šimůnek, S. Shiozawa, and J.W. Hopmans. 2000.
Estimation of soil hydraulic and solute transport parameters fromtransient infiltration experiments. Adv. Water Resour. 23:677–688.
Jabro, J.D., J.M. Jemison, Jr., R.H. Fox, and D.D. Fritton. 1994. Predicting bromide leaching under field conditions using SLIM and
MACRO. Soil Sci. 157(4):215–223.
Jacques, D., J. Šimůnek, A. Timmerman, and J. Feyen. 2002.
Calibration of Richards’ and convection-dispersion equations tofield-scale water flow and solute transport under rainfall conditions.
J. Hydrol. (Amsterdam) 259:15–31.
Jarvis, N.J., P.-E. Jonasson, P.E. Dik, and I. Messing. 1991. Modelling
water and solute transport in macroporous soil. I. Model description and sensitivity analysis. J. Soil Sci. 42:59–70.
Kätterer, T., B. Schmied, K.C. Abbaspour, and R. Schulin. 2001.
Single- and dual-porosity modeling of multiple tracer transport
through soil columns: Effects of initial moisture and mode ofapplication. Eur. J. Soil Sci. 52:25–36.
Klute, A. 1986. Water retention: Laboratory methods. p. 635–660.
In A. Klute (ed.) Methods of soil analysis. Part 1. Agron. Monogr. 9.
ASA and SSSA, Madison, WI.
Klute, A., and C. Dirksen. 1986. Hydraulic conductivity and diffusivity:
Laboratory methods. p. 687–734. In A. Klute (ed.) Methods of soil
analysis. Part 1. Agron. Monogr. 9. ASA and SSSA, Madison, WI
Kohler, A., K.C. Abbaspour, M. Fritsch, M.T. van Genuchten, and R.
Schulin. 2001. Simulating unsaturated flow and transport in amacroporous soil to tile drains subject to an entrance head: Modeldevelopment and preliminary evaluation. J. Hydrol. (Amsterdam)
254:67–81.
Köhne, J.M. 1999. Analyse präferentiellen Wasserflusses und Stofftransports in strukturierten Böden mit Hilfe eines Dual-Porositätsmodells. Ph.D. dissertation (in German, with English summary).
In H. Roweck and B. Lennartz (ed.) Schriftenreihe Inst. Water
Management and Landscape Ecology 30. Christian-AlbrechtsUniversity, Kiel, Germany.
Köhne, J.M., S. Köhne, and H.H. Gerke. 2002. Estimating the hydraulic functions of dual-permeability models from bulk soil data,
Water Resour. Res. 38(7):1121. doi:10.1029/2001WR000492.
Köhne, J.M., and B.P. Mohanty. 2005. Water Flow Processes in a Soil
Column with a Cylindrical Macropore: Experiment and Hierarchical Modeling. Water Resour. Res. 41:W03010. doi:10.1029/
2004WR003303.
Köhne, J.M., B.P. Mohanty, J. Šimůnek, and H.H. Gerke. 2004. Evaluation of a second-order water transfer term for variably saturated
dual-permeability flow models. Water Resour. Res. 40:W07409.
doi:10.1029/2004WR003285.
Larsson, M.H., and N.J. Jarvis. 1999a. A dual-porosity model to quantify macropore flow effects on nitrate leaching. J. Environ. Qual. 28:
1298–1307.
Larsson, M.H., and N.J. Jarvis. 1999b. Evaluation of a dual-porosity
model to predict field-scale solute transport in a macroporous soil.
J. Hydrol. (Amsterdam) 215:153–171.
Legates, D.R., and G.J. McCabe, Jr. 1999. Evaluating the use of“goodness-of-fit” measures in hydrologic and hydroclimatic model
validation. Water Resour. Res. 35:233–241.
Lennartz, B., J. Michaelsen, W. Wichtmann, and P. Widmoser. 1999.
Time variance analysis of preferential solute movement at a tiledrained field site. Soil Sci. Soc. Am. J. 63:39–47.
Logsdon, S.D. 1995. Flow mechanisms through continuous and buried
macropores. Soil Sci. 160:237–242.
Malone, R.W., M.J. Shipitalo, L. Ma, L.R. Ahuja, and K.W. Rojas.
2001. Macropore component assessment of the root water quality
model (RZWQM) using no-till soil blocks. Trans. ASAE 44:843–852.
Medina, A., and J. Carrera. 1996. Coupled estimation of flow andsolute transport parameters. Water Resour. Res. 32:3063–3076.
Mishra, S., and J.C. Parker. 1989. Parameter estimation for coupledunsaturated flow and transport. Water Resour. Res. 25:385–396.
Meyer-Windel, S. 1998. Transport processes of bromide and herbicides
under steady-state and transient flow conditions. Ph.D. diss. (InGerman.) Schriftenreihe Institut für Wasserwirtschaft und Landschaftsökologie 28. Christian-Albrechts-University, Kiel, Germany.
Mohanty, B.P., R.S. Bowman, J.M.H. Hendrickx, and M.Th. van
Genuchten. 1997. New piecewise-continuous hydraulic functions
for modeling preferential flow in an intermittent-flood-irrigatedfield. Water Resour. Res. 33:2049–2063.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

76

VADOSE ZONE J., VOL. 5, FEBRUARY 2006

Nissen, H.H., P. Moldrup, and K. Henriksen. 1998. High-resolution
time-domain reflectometry coil probe for measuring soil watercontent. Soil Sci. Soc. Am. J. 62:1203–1211.
Schwartz, R.C., A.S.R. Juo, and K.J. McInnes. 2000. Estimating parameters for a dual-porosity model to describe non-equilibrium, reactive
transport in a fine-textured soil. J. Hydrol. (Amsterdam) 229:149–167.
Šimůnek, J., D. Jacques, J.W. Hopmans, M. Inoue, M. Flury, and M.Th.
van Genuchten. 2002. Solute transport during variably saturated
flow. Inverse methods. p. 1435–1449. In J.H. Dane and G.C. Topp(ed.) Methods of soil analysis. Part 4. SSSA Book Ser. 5. SSSA,
Madison, WI.
Šimůnek, J., N.J. Jarvis, M.T. van Genuchten, and A. Gärdenäs. 2003.
Review and comparison of models for describing nonequilibrium
and preferential flow and transport in the vadose zone. J. Hydrol.
(Amsterdam) 272:14–35.
Šimůnek, J., M. Sejna, and M.Th. van Genuchten. 1998. The

HYDRUS-1D software package for simulating the movement ofwater, heat, and multiple solutes in variably-saturated media.
Version 2.0. IGWMC-TPS-70. International Ground Water Modeling Center, Colorado School of Mines, Golden, CO.
Šimůnek, J., O. Wendroth, N. Wypler, and M.Th. van Genuchten. 2001.
Non-equilibrium water flow characterized by means of upward
infiltration experiments. Eur. J. Soil Sci. 52:13–24.
Tofteng, Ch., S. Hansen, and H.E. Jensen. 2002. Film and pulse flow in
artificial macropores. Nord. Hydrol. 334(4):263–274.
van Genuchten, M.Th. 1980. A closed-form equation for predicting thehydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:
892–898.
van Genuchten, M.Th., F.J. Leij, and S.R. Yates. 1991. The RETC codefor quantifying the hydraulic functions of unsaturated soils. EPA/
600/2-91/065. USEPA, Office of Research and Development, Washington, DC.