Lateral Water Diffusion in an Artificial Macroporous System:
Modeling and Experimental Evidence
P. Castiglione,* B. P. Mohanty, P. J. Shouse, J. Simunek, M. Th. van Genuchten, and A. Santini
ABSTRACTused to interpret results from solute transport experiments. For certain conditions (e.g., negligible matrix
domain flux, steady-state water flow and instantaneousadsorption), the transfer rate coefficients for solute (␣s)
can be determined from moment analysis of the concentration breakthrough curve (BTC) (e.g., Valocchi, 1990).
This methodology was tested by Hu and Brusseau (1995)
in a series of leaching experiments through artificial soilcolumns. The authors explored different aggregate geometries, for which theoretical expressions of the transfer coefficients are available, and with most of the parameters of the two-region Mobile–Immobile Model(MIM) (van Genuchten and Wierenga, 1976) determined independently of the data.
In more general cases, for example when water moves
through both the matrix and the macropore regions,experimental measurements of the mass exchange are
difficult to make since more complex models with alarge number of parameters are generally required. Matrix and macropore flows cannot be easily discriminatedduring water flow experiments, thus making the hydraulic characterization of the two domains very uncertain.
On the basis of an analogy between water and solutetransfer processes, Gerke and van Genuchten (1993b)
suggested that the transfer rate coefficient for water(␣w) can be derived from the corresponding coefficientfor solute (␣s), as determined using the MIM hypotheses. To date, this assumption has not been validated
by experimental evidence, mostly because of technicallimitations in monitoring the lateral exchange process.
We present results from infiltration experiments conducted on a soil column with an artificial macropore.
Detailed knowledge of the geometry of the system andsoil matrix properties, as well as the controlled laboratory conditions, allowed us to address some of the uncertainties mentioned above. For instance, we were able
to discriminate matrix flow from macropore flow, andto measure the depth-integrated water exchange flux
between the two domains. The experimental results
were described reasonably well with a dual-permeability
model, with most of the parameters determined independently of the data. A numerical analysis of the watertransfer process also allowed us to determine the transfer rate coefficient for water, which otherwise would bea fitting parameter.

In two-domain schematizations of macroporous soils or fractured
rock systems, lateral mass exchange between macropores and the soil
matrix is generally modeled as an apparent first-order process. With
respect to lateral diffusion, the system is thus characterized by a single
parameter, the transfer rate coefficient, which is difficult to estimate
a priori. We conducted water infiltration experiments in a laboratory
column with an artificial macropore. The novel design of the experimental setup allowed us to discriminate between matrix flow andmacropore flow, from which we could estimate the water exchange flux
between the two domains. Most of the parameters in a dual-permeability model could be determined independently of the experimental
data. In particular, a theoretical expression for the transfer rate coefficient was derived by assuming lateral water and solute diffusion tobe similar processes. Numerical analysis of the water exchange process
revealed that the transfer coefficient depended also on the macropore
conductivity. When this dependency was taken into account, the model
reproduced the experimental data reasonably well.

M

ost deterministic approaches for modeling preferential flow in macroporous soils or unsaturated
fractured rock rely on two- or multidomain continuum
assumptions (Smettem and Kirkby, 1990; Cushman,
1990; Othmer et al., 1991; Pruess, 1999; National Resource Council, 2001). Typically, heterogeneous media
are modeled as two overlapping porous systems, a highpermeability domain associated with the macropore or
fracture network and a low permeability domain associated with the matrix blocks. Exchange of mass (water
and/or solute) between the two systems will occur in
response to nonequilibrium conditions, thus reproducing lateral exchange phenomena often recognized instructured media (Beven and Germann, 1982). Interdomain transfer processes are diffusion based and assuch depend on the diffusivity of the matrix blocks andthe geometric configuration of the preferred flow paths.
When these processes are modeled by means of a firstorder approximation (FOA), the system can be characterized by one unique parameter, the mass transfer
rate coefficient. Despite extensive efforts in the past two
decades, how closely the actual mass exchange process
is captured by the first-order approximation and how
to determine the transfer rate coefficients remain
poorly understood.
First-order mass transfer models have been widely
P. Castiglione, B.P. Mohanty, P.J. Shouse, J. Simunek, M.Th. van
Genuchten, George E. Brown, Jr. Salinity Laboratory, 450 W. Big
Springs Road, USDA, ARS, Riverside, CA; and A. Santini, Department of Agricultural Engineering, University of Naples “Federico
II”, Italy. B.P. Mohanty, currently, Biological and Agricultural Engineering, Texas A&M University, College Station, TX. Received 14
June 2002. Original Research Paper. *Corresponding author (paoloc@
uidaho.edu).

THEORYThe dual-permeability model (DPM) used in thisstudy to simulate vertical water flow through a singlemacropore column was described in detail by Gerke
Abbreviations: BTC, breakthrough curve; DPM, dual-permeability
model; FOA, first-order approximation; MIM, mobile–immobile model;
PAM, polyacrylamide.

Published in Vadose Zone Journal 2:212–221 (2003).

212

213

www.vadosezonejournal.org

and van Genuchten (1993a). The governing equations
for one-dimensional flow are

冢
冣
⌫
⳵␪
⳵
⳵h
K
⫽
⫺K 冣⫹
冢
⳵t
⳵z
⳵z
1⫺w
⳵␪f
⌫
⳵h
⳵
Kf f ⫺ K f ⫺ w
⫽
⳵t
⳵z
⳵z
wf
m

m

m

[1]

w

[2]

m

f

in which the subscript “f” refers to the macropore domain, and “m” to the matrix domain; h (L) is the matricpotential, K (L T⫺1) the hydraulic conductivity, z (L) is
the spatial coordinate, t (T) is time, ⌫w (T⫺1) is thewater transfer term, representing the flux of water being
exchanged between macropores and the soil matrix per
unit volume of bulk soil, and wf is a dimensionlessweighting factor. The volumetric water content ␪f in Eq.
[1] is defined as the ratio of the volume of water of thefracture domain (Vw,f ) to the total volume of fractures(Vt,f ) ␪f ⫽ Vw,f /Vt,f . Analogously, ␪m ⫽ Vw,m /Vt,m . The
volumetric weighting factor (wf) is defined as the ratioof the fracture domain volume relative to the total soilvolume (Vt):wf ⫽ Vt,f /Vt[3]

While some macroscopic quantities, such as watercontent or flux, for a dual-permeability medium aredefined as weighted averages of the corresponding properties of the individual domains, this cannot be done
for the hydraulic conductivity, since in general no uniquepressure head gradient exists (Beven and Germann,
1982). When the system is in hydraulic equilibrium (hf ⫽hm), however, it is possible to define the hydraulic conductivity as (Peters and Klavetter, 1988):
K ⫽ wf Kf ⫹ (1 ⫺ wf)Km[4]

Equation [4] holds in particular for saturated conditions
and thus permits us to estimate the macropore saturated
conductivity (Kf,sat) from Ksat , Km,sat , and wf .
The interdomain water exchange flux is modeled bymeans of the first-order approximation (FOA):⌫w ⫽ ␣w(hf ⫺ hm)[5]in which hf and hm are averaged over some representative elementary volume, and ␣w (L⫺1 T⫺1) is the masstransfer rate coefficient for water. Assuming analogy
between water and solute lateral exchange processes,
Gerke and van Genuchten, (1993b) derived ␣w from theequivalent transfer rate coefficient for solute. The FOA
for solute exchange flux (⌫s) is expressed as:⌫s ⫽ ␣s(cf ⫺ cm)[6]

where ␣ (T⫺1) is the transfer rate coefficient for solute,sand cf and cm (M L⫺3) are the solute concentrations
for the fracture and the matrix domain, respectively.
Assuming applicability of MIM (van Genuchten and
Wierenga, 1976), ␣s is fully determined by geometry
and diffusivity of the soil aggregates, and is generally
approximated by the expression (van Genuchten and
Dalton, 1986):
␣s ⫽H
␪mDa,m␭2[7]

Fig. 1. Scheme of the soil column with a vertical cylindrical macropore.

where H is a dimensionless geometry-dependent coefficient and ␭ (L) is the characteristic length of the matrix
structure (e.g., the radius for spherical or cylindricalaggregates, or one-half of the fracture spacing for parallel rectangular matrix blocks); Da,m (L2 T⫺1) is the effective diffusion coefficient of solute in the soil matrix(notice that the subscript “m” represents here the soil
“matrix,” not the “mobile” domain as is commonly done
in literature). For simple geometries of the aggregates,
theoretical values for the geometrical parameters in Eq.
[7] can be derived (Valocchi, 1990; van Genuchten and
Dalton, 1986). In particular, several expressions are reported in literature for cylindrical macropores (vanGenuchten and Dalton, 1986; Gerke and van Genuchten, 1996; Young and Ball, 1995), which we found to
give very similar results. Below we adopt the expression
proposed by Young and Ball (1995) to be used in Eq. [7]:
␭⫽b⫺a

[8]

8共␰ 20 ⫺ 1兲共␰ 0 ⫺ 1兲2

H⫽

4␰ 40ln␰ 0 ⫺ 3␰ 40 ⫹ 4␰ 20 ⫺ 1, with ␰ 0 ⫽ b .
a[9]

In Eq. [8] and [9], a is the radius of the macropore,and b is the radius of the soil mantle surrounding themacropore (see Fig. 1).
Due to their common diffusive nature, lateral waterand solute transfers are often considered similar process(Dykhuizen, 1990; Gerke and van Genuchten, 1993a;
Zimmerman et al., 1990). Assuming that the water flux
divergence is zero within the soil matrix, Gerke andvan Genuchten (1993b) found that the geometry of theaggregates affects water and solute mass transfers inthe same way. They therefore proposed the following
expression for the transfer rate coefficient for water:

214

VADOSE ZONE J., VOL. 2, MAY 2003

␣w ⫽

H
Ka␥w
␭2

[10]

which is derived from Eq. [7] with the term Ka␥w replacing ␪mDa,m. In Eq. [10], ␥w is a dimensionless scaling
factor and Ka (L T⫺1) the effective hydraulic conductivity function of the macropore–matrix interface, evaluated in terms of both hm and hf as follows:
Ka ⫽ 0.5[Ka(hm) ⫹ Ka(hf)][11]

In the absence of organic matter or other deposits along
the macropore wall (which may alter the permeabilityof the aggregates at the interface), Ka is best represented
by the matrix conductivity function Km (Gerke and van
Genuchten, 1993a).
Gerke and van Genuchten (1993b) and Gerke andvan Genuchten (1996) tested the applicability of Eq. [10]for several aggregate geometries (spherical, rectangularslab, and cylindrical macropores) by comparing simulation results of the one-dimensional DPM with thosegenerated using equivalent two-dimensional geometrybased transport models. For simplicity, the effects ofgravity were neglected in their calculations. These numerical studies revealed that, in general, Eq. [10] performs well as long as the dimensions of the matrix blocks
are not too large. Their results also suggest that, provided Eq. [11] is used, an average value of 0.4 for thescaling factor ␥w should be reasonably applicable to allof the geometries investigated. This value, however, may
not be appropriate for describing water exchange whenthe vertical water flow in the matrix is considered, orfor initial and boundary conditions different from thoseused in the original studies of Gerke and van Genuchten(1993b). To take this scenario-related dependency into
account, we performed a numerical analysis of the waterexchange process in a single macropore system for thesame initial and boundary conditions reproduced in our
experiments. Our calculations allowed us to determine
an optimal value for the scaling factor ␥w, which subsequently will be used in the analysis of our experimental data.

Numerical Analysis of the Water Transfer Process
We generated series of synthetic data for lateral waterexchange by simulating water infiltration through a soilcolumn with a single macropore, as depicted in Fig. 1.
Since the process is axial-symmetric, a modified versionof the two-dimensional HYDRUS-2D code (Simunek
et al., 1999) could be used. The synthetic water transfer
data were subsequently compared with results of the onedimensional dual-permeability model to test the applicability of Eq. [8], and to determine a value of the scaling
factor ␥w that best describes the exchange process for the
initial and boundary conditions used in our experiments.
For both the two- and the one-dimensional simulations,and for both the matrix and macropore domain, theinitial pressure head profile varied linearly from ⫺155cm at the top to ⫺135 cm at the bottom. The pressure
head was equal to zero at the top of the column, while
seepage conditions were imposed at the bottom bound-ary. In our simulation we explicitly considered verticalwater flow in the soil matrix.
Two-Dimensional Simulations
The HYDRUS-2D code solves the axial-symmetricflow equations in the r–z plane using a mass lumped
Galerkin finite element scheme. The flow domain consisted of two connecting porous media: a single macropore (r ⱕ 0.5 mm) at the center and the surrounding
matrix region. A very fine finite element grid, as wellas small time steps (⌬tmin was 10⫺9 h), were used toincrease the accuracy of the calculations. We initiallyexperienced some numerical problems, even with a veryfine discretization near the macropore–matrix interface.
This was caused by the fact that in HYDRUS-2D the
independent variables and hydraulic properties are assigned to nodes, and smoothed linearly between neighboring nodes, as is traditionally done in finite element
schemes. This approach causes flow and transport properties to be smoothed across the interfaces between
different materials. Most of the numerical problems
were solved by modifying the code to allow assignmentof material properties to elements, rather than to individual nodes.
The van Genuchten–Mualem model (van Genuchten,1980) was used to describe the retention and conductivity functions of both the matrix and the macropore domain:h(Se) ⫽

冤

冥

1/n
1
1/m
1 ⫺ 共S ⫺
⫺ 1兲
e
␣

冤

[12]

冥

m 2

K(Se) ⫽ KaS le 1 ⫺ 共1 ⫺ S e1/m兲

[13]

where Se ⫽ (␪ ⫺ ␪r)/(␪s ⫺ ␪r) is the effective saturation, ␣(L⫺1), n and l are empirical parameters, ␪r is the residual
water content, ␪s is the saturated water content, and Ks(L T⫺1) is the saturated hydraulic conductivity. We used
the m ⫽ 1 ⫺ 1/n constraint so that each domain ischaracterized by a total of six parameters (hereafterreferred to as “VGM parameters”). For our calculations, the matrix domain was assigned the VGM parameters estimated (as will be shown later) through inverse
analysis of data from an infiltration experiment without
a macropore (␪s ⫽ 0.38; Ks ⫽ 0.08 cm h⫺1; l ⫽ 0.5; ␪r ⫽0.20, ␣ ⫽ 0.005, and n ⫽ 1.664). In contrast, three different hydraulic functions were considered for the macropore domain. The saturated conductivity (Kf,sat) waschosen equal to 2000, 20 000, and 100 000 cm h⫺1, while
the remaining parameters were in all cases fixed at ␪r ⫽0, ␪s ⫽ 0.5, ␣ ⫽ 0.1, n ⫽ 2, and l ⫽ 0.5. Two somewhat
conflicting requirements were considered when choosing the VGM parameters for the macropore domain:(i) having a small air-entry tension (high ␣ value), below
which the hydraulic conductivity drops sharply to negligible values (high n value), and (ii) facilitating a stablenumerical solution of the flow equations. Each simulation was repeated for three different soil mantle radii:
1.7, 2.7, and 12 cm.

215

www.vadosezonejournal.org

Fig. 3. Partitioning of the bottom of the column in six pie sectors.
The metal dividers are 15 cm high. The stainless-steel nets support
the soil when the ceramics are removed.

Fig. 2. Experimental setup.

One-Dimensional Simulations
Infiltration processes were simulated for the same
scenarios (geometry, hydraulic properties, initial andboundary conditions) considered in the previous section
using the one-dimensional DPM (as incorporated into
HYDRUS-1D; Simunek et al., 1998). We next performed
one-dimensional inverse analyses of the synthetic datato optimize ␥w, using cumulative mass transfer data inthe objective function. Since optimized values for thetransfer coefficient depend somewhat on the characteristic time scale of the synthetic data (Young and Ball,
1995), we homogenized the different series by extending
each set of data to well after equilibrium was reached.
MATERIALS AND METHODS
Experimental Setup
A sandy loam (Typic Haploxeralf) with a small clay fraction(6%, mostly kaolinite) was used in this study because of its
low shrink–swell capacity. The soil, collected from the Angeles
National Forest (California) at 20 to 40 cm depth, was airdried, sieved (2 mm), and packed into a large column in 3-cm
increments to a dry bulk density of 1.56 g cm⫺3. The column(Fig. 2) was made from an acrylic cylinder (75 cm deep; i.d. ⫽24 cm), resting on a specifically designed support for controlling the bottom boundary condition and collecting the drainage effluent (Fig. 3). The lower portion of the column waspartitioned into six pie-shaped chambers using vertical dividers (Fig. 3). The dividers were made of 15-cm-high stainlesssteel sheets, glued into grooves in the base plate and held
rigidly by means of a vertical support placed in the center ofthe base plate. The base plate, also made of acrylic, housed
six ceramic plates (bubbling pressure ⫽ 5 ⫻ 104 Pa), one ineach of the six chambers (Fig. 3). The ceramic plates weresealed into place with a system of rubber gaskets and removable cover panels. Drainage effluent from each chamber wascollected separately in a system of six acrylic tubes (Fig. 2) that
permitted automated measurement of the cumulative outflow
using pressure transducers. The pressure head at the bottomof the soil column was regulated by adjusting the air pressure

within the outflow tubes. A 3-cm layer of diatomaceous earthmaterial (having a negligible hydraulic resistance) was placedin each chamber between the soil and the ceramic plates.
Collecting drainage effluent from each chamber independently allowed us to detect horizontal heterogeneities in thewetting front. In particular, we were able to verify the stabilityof the matrix flux in the absence of macropores. Most importantly, when a single artificial macropore was made and connected to one outflow chamber, separate measurements ofthe matrix and macropore fluxes could be obtained by contrasting the effluent fluxes from the different chambers. Thismethodology is applicable assuming preferential flow is fully
confined within the artificial macropore. We prevented preferential flow along the column walls by inserting small glue rings(“wall flow dams”; see Fig. 3) at several depths on the internal
wall of the column. The rings acted like a series of dams forcing any wall flow to reenter the soil matrix. Soil water pressure
heads were measured using pencil-size tensiometers placed at
5-cm depth increments, while horizontal heterogeneities inthe soil water pressure head were monitored using two setsof six tensiometers placed around the circumference of thesoil column at depths of 50 and 75 cm below the soil surface.
We also installed six time domain reflectometry three-rod
probes (10 cm long) at 10-cm depth increments to measure
the matrix water content profile.

Artificial Macropores
Different techniques for creating artificial macropores aredescribed in the literature, depending on the objectives of theresearch, soil type, and column dimensions (Hu and Brusseau
1995; Li and Ghodrati, 1997; Czapar et al., 1992). In our study
we created a single macropore by inserting a hollow stainlesssteel tube (1-mm o.d.) from the surface into the soil matrix.
The macropore was located approximately at the center ofthe column and reached the diatomaceous earth layer at thebottom of one of the chambers.
To help stabilize the artificial macropore, we applied awater-soluble polymer (polyacrylamide, PAM) along the macropore walls. This was done by pouring a PAM solution (75 g
L⫺1) into the stainless-steel tube once it was fully inserted into
the soil. The tube was subsequently slowly removed from thesoil such that the PAM left a coating on the macropore walls.
To enhance diffusion into the surrounding matrix, the PAM
was applied when the column was dry. After application, thecolumn was left for 48 h before a new experiment was performed. While the technique proved effective, the macropore

216

VADOSE ZONE J., VOL. 2, MAY 2003

itself was not always perfectly stable, and new applicationsof the polymer were needed periodically between separate
flow experiments.

Experimental Procedure
Using a 0.33 g L⫺1 CaCl2 solution, several infiltration anddrainage experiments were performed before making the artificial macropore. These initial experiments were used to test
the functionality of the experimental setup, especially the active drainage system, and to verify the uniformity of waterflow through the soil column. A first-type boundary condition
was imposed at the soil surface during the infiltration experiments by means of a tension infiltrometer (Ankeny et al.,
1988; Perroux and White, 1988); this also permitted us tomeasure the infiltration rate. We used a disk (Soil Measurement Systems, Tucson, AZ) with the same diameter as thecolumn, so that flow was always one-dimensional. A thin sandlayer was placed on the soil surface to ensure perfect contact
between soil and the infiltrometer. The head loss through thesand was measured by means of a very small tensiometer
placed at the sand–soil interface.
One major challenge of the experimental setup was accurate
control of the bottom boundary condition. When a tension is
imposed through the ceramic plates, the presence of entrapped
air, as well as the existence of fluxes greater than the conductivity of the ceramic plates, may cause the effective bottom tension to diverge from the applied one. Such a situation could
lead to misinterpretations of the experimental results, particularly since the magnitude of the flux through the artificial
macropore is not known a priori. We therefore designed thebase of the column to allow us to remove or install different
ceramic plates as needed. For example, to facilitate drainage
from the column we applied tension at the bottom, whereas
the ceramics were removed during the infiltration experiments
to produce a seepage-type bottom boundary condition.
We present results for two infiltration experiments, one withand one without the artificial macropore, performed underidentical initial and boundary conditions. The effects of themacropore on the overall flow regime were deduced by contrasting the measured inflow and outflow rates. In both experiments, the column was initially drained by applying suction
at the bottom (h ⫽ ⫺250 cm). Pressure head profiles during
the drainage showed fairly uniform values in the upper threequarters of the column, and a steep decrease in the pressure
head in the lower portion. Once a nearly uniformly dry profile
was obtained in the upper three-quarters of the column, theceramic plates were removed and the column was left to equilibrate for about 1 wk. Having a relatively low initial watercontent distribution favored extensive lateral exchange of water during the macropore infiltration experiment. In both experiments, the top pressure head equaled zero, while a seepage
boundary condition was applied at the bottom by removing
the ceramic plates.

RESULTS AND DISCUSSION
Results from Numerical Simulations
Results from the two-dimensional simulations highlighted several factors affecting the water exchange between the macropore and matrix. In general, water flowin the two regions did not evolve independently (i.e.,based on their specific hydraulic properties), but werevery much affected by each other. For example, thewetting front in the macropore had a tendency to slow
down because of diffusion of water into the surrounding

Fig. 4. Two dimensional simulations of pressure head distribution
during infiltration processes for a 2.7-cm soil mantle. All figures
show snapshots at 0.05 h.

matrix. This lateral exchange process is evident in Fig. 4,
which shows snapshots (t ⫽ 0.05 h) of the infiltration
process simulated for three values of Kf,sat. Notice thatthe relative position of the two wetting fronts determined the longitudinal extent of the region affected bytransfer process (indicated with d in Fig. 4). Since thisregion is large for highly conductive macropores, thedepth-integrated water exchange rate is higher for larger
Kf,sat values. This result is also evident from Fig. 5, where
the cumulative transfer flux integrated over the entire
column is shown for the different scenarios. For eachradius of the soil mantle, the water exchange rate ishigher for higher macropore conductivity. Notice alsothat the total amount of exchanged water depends only
on the size of the soil mantle for the smallest radius(Fig. 5b), whereas for larger mantles it also depends on
Kf,sat (Fig. 5a). This is due to vertical flow of water entering the matrix region from the top. In general, the maximum value of the cumulative transfer rate is smallerwhen vertical flow in the matrix is relatively fast compared with that of the macropore (smaller Kf,sat), since
the pressure difference between the two domains (i.e.,the driving force for water exchange) will vanish sooner.
This effect is negligible for the smallest radii, for which
equilibrium is reached only in a relatively short time.
Comparing the results from two- and one-dimensional simulations, we found that when the scaling factor
␥w in Eq. [8] was set equal to 1, water transfer rates
predicted with the dual-permeability model were higher
than those calculated with the two-dimensional simulations, as shown in Fig. 6 for a 12-cm soil mantle andKf,sat ⫽ 100 000 cm h⫺1. These results are consistent withthose obtained by Gerke and van Genuchten (1993b).
When optimized values of the scaling factor were used,the DPM reproduced the synthetic data more closely,
not only in terms of the cumulative water exchange(Fig. 6), but also for the inflow and outflow rates. The
optimized scaling factors are listed in Table 1 for the217

www.vadosezonejournal.org

Fig. 6. Cumulative water exchange predicted by the one-dimensional
dual-permeability model (DPM) compared with the transfer calculated by HYDRUS 2-D.

the scaling factor ␥w is correlated to the ratio of the twodomain’s conductivities.

Infiltration Experiment without MacroporesFig. 5. Summary of the cumulative water transfer for nine two-dimensional simulations.

different scenarios considered. In particular, data for the
12-cm soil mantle radius (the radius of our experimental
column) allowed us to determine an appropriate value
for the transfer coefficient in subsequent analyses ofthe experimental results. The data are plotted in Fig. 7
against Kf,sat. The dependency of the scaling factor on themacropore saturated conductivity appears pronounced
for the lowest Kf,sat, whereas ␥w approaches a fairly constant value for high Kf,sat values. We note that Gerke andvan Genuchten (1993b), who explored different hydraulic properties for the matrix domain in their numericalanalysis, also found a slight variability of ␥w with thematrix hydraulic conductivity. Although our results cannot be directly compared with the findings of Gerkeand van Genuchten (1993b) (since different boundary
conditions were adopted in the two studies, and since
we explicitly took into account the vertical water movement in the matrix domain), it appears that, in general,

Water flow during infiltration without the macropore
appeared to occur under equilibrium conditions. Thiswas evident from visual inspection of the sharp wetting
front and from the analysis of the data, which were welldescribed using the Richards equation (see Fig. 8 and 9).
The cumulative inflow curve displayed a pronounced sorption phase at the beginning of the infiltration (Fig. 8).
During this stage, the infiltration rate was high, being
driven by both capillary and gravity forces. As the wetting front moves down, pressure head gradients eventually become negligible and the infiltration rate approaches
a constant value, roughly equal to the saturated hydraulic conductivity. Outflow started about 13 h later andrapidly became constant and equal to the inflow rate.
The total effluent was evenly distributed among the six
chambers (data not shown), and the column at this time
appeared thoroughly saturated (Fig. 9).
Data from the infiltration experiment were used tocharacterize the soil matrix hydraulic properties. We
used the HYDRUS-1D code (Simunek at al., 1998) for
inverse analysis (e.g., Hopmans et al., 2002) using inflow,
outflow, and tension data. Of the six VGM parameters,
we measured ␪s independently of the infiltration experiment (␪s ⫽ 0.38), determined Ks from the steady-state
Table 1. Optimized values of the scaling factor ␥w for different soil
mantle radii and macropore saturated hydraulic conductivities.
Mantle radius
cm
1.7
2.7
12

Kf,satOptimized ␥w

cm h⫺12 000
20 000100 000
2 000
20 000100 000
2 000
8 000
20 000100 000

0.04870.19910.60620.10080.24020.68740.1026
0.27920.40330.6348

218

VADOSE ZONE J., VOL. 2, MAY 2003

Fig. 7. Scaling factor dependency on the macropore saturated conductivity.

infiltration rate (Ks ⫽ 0.08 cm h⫺1), and fixed l to 0.5,
whereas the remaining three parameters were optimized(␪r ⫽ 0.20, ␣ ⫽ 0.005, and n ⫽ 1.664). Good correspondence was obtained between the fitted and measured
data (Fig. 8 and 9). The uniqueness of the solution was
empirically verified by repeating the inverse analysis
with different initial values of the unknown parameters.
We note, however, that the estimated hydraulic properties probably were somewhat biased by a lack of information of soil matrix hydraulic behavior in the dry range(because of the dimension of the column, it was notpossible to achieve very dry conditions). Nevertheless,the hydraulic characteristics in the wet range, which areof interest in this study, were successfully estimatedfrom the available data. In particular, the early stage ofthe infiltration process should define the soil matrix
diffusivity, which regulates interdomain exchange, andultimately determines the strength of macropore flow.

Infiltration Experiment with One
Artificial Macropore
From capillary theory, a top boundary pressure headgreater than ⫺3 cm was needed to initiate flow throughFig. 8. Measured and fitted inflow and outflow curves during infiltration in the homogeneous soil column (without macropore). Simulated results are obtained using Richards’ equation.

Fig. 9. Pressure head profiles during the infiltration experiments without the macropore. Simulated results are obtained using Richards’ equation.

the 1-mm-diameter macropore. The occurrence of macropore flow was immediately evident, since the infiltration rate was considerably higher than that observed
without the macropore (Fig. 10). Furthermore, after
about 45 min a roughly constant flux was observed in
the chamber where the macropore was located, while
no outflow was observed in the other chambers (datanot shown). This flux, 0.082 cm h⫺1, was therefore assumed to be all macropore outflow (Of). This experiment allowed us to discriminate between macropore
flow and matrix flow, and to measure the interdomain
exchange flux. We reproduced the same initial andboundary conditions as for the infiltration experiment
without a macropore. In particular, a zero pressure head
was maintained at the top. During infiltration, the fol-Fig. 10. Cumulative infiltration rates measured with and without
the macropore.

219

www.vadosezonejournal.org

lowing relationship holds:
Itot(t) ⫽ Im(t) ⫹ If(t)[14]where Im is the flux of water infiltrating into the soil
matrix, If is the flux entering the macropore, and Itot
represents the total infiltration rate, as measured withthe tension infiltrometer. An analogous relationship canbe written for the outlet flux:
Otot(t) ⫽ Om(t) ⫹ Of(t)[15]where Om and Of are the contributions of matrix andmacropore, respectively, to the total outflow (Otot). Because of lateral diffusion, in general If ⬎ Of, the difference between the two fluxes representing the total water
flux exchanged between macropore and matrix:⌫w(t) ⫽ If(t) ⫺ Of(t)[16]Note the different symbols for the local (⌫w) and thedepth-integrated exchange flux (⌫w). The flux entering
the macropore (If) could not be measured directly.
Nonetheless, assuming that Im in Eq. [14] equals theinfiltration rate measured during the experiment without the macropore, and hence that the presence of themacropore does not affect the matrix flux, If could be
estimated from the difference between the inflow rates
observed during the two experiments (Fig. 10). This
assumption seems very reasonable at the early stage ofthe process, when lateral diffusion is still confined toonly a relatively small portion of the column. As theinfiltration process proceeds, however, the assumption
will become unrealistic, especially if lateral equilibrationis relatively fast. The validity of the above assumption
was tested by suddenly interrupting the macropore flow
process (before steady-state conditions were reached)and verifying that the cumulative infiltration curve proceeded parallel to the corresponding nonmacropore
curve (see Fig. 10). Macropore flow was interrupted by
reducing the top pressure head to below the macropore
air-entry-value (⫺3 cm). Several steady-state experiments without macropores had previously revealed thatsuch a small pressure drop did not affect the matrix
flow rate. Following the interruption, bubbling in the
Mariotte and the reservoir towers of the tension infiltrometer stopped for several minutes until the new pressure
head was established at the top boundary. During thisperiod the instrument could not record the matrix water
flow, which explains the short horizontal segment in the
cumulative infiltration data shown in Fig. 10. The total
exchange flux during the macropore infiltration process
can therefore be measured as⌫w ⫽ (Itot ⫺ Im) ⫺ Of

[17]

The cumulative lateral flux measured in this way isshown in Fig. 11, and will be discussed in more detailbelow. The experiment was repeated three times, withthe three replications giving nearly identical results (notshown here). We emphasize here that steady-state conditions could never be reached completely because ofslow, yet persistent clogging of the macropore during infiltration.

Fig. 11. Cumulative infiltration curves measured in the soil column
with a macropore, and calculated using the dual permeability
model, for ␥w optimized and ␥w ⫽ 1.

Comparisons of Model Results with
Experimental Data
Application of the DPM involves two sets of VGM
parameters, one for each domain, plus the mass transfer
rate coefficient (␣w) and the volumetric weighting factor(wf). In our case, wf is known from the geometry of thesystem (wf ⫽ 1.8 ⫻ 10⫺5). The characteristics of thematrix domain were determined from an analysis of theinfiltration data without the macropore. Estimation of
the macropore parameters was much more uncertain,even for the controlled conditions of our experiment.
A single macropore remains either completely full ornearly empty, except possibly for some film flow, sothat its behavior during variably saturated flow cannot
be easily investigated. Still, this behavior suggests thatthe macropore conductivity function, Kf(hf), should becharacterized by a sharp drop to negligible values aftera very small air-entry suction is exceeded. This representation of a cylindrical macropore as an equivalent porous medium reduced the range of possible values forthe macropore domain hydraulic parameters.
Since it was not possible to maintain macropore flow
until steady-state conditions were reached, Eq. [4] could
not be used to estimate the saturated hydraulic conductivity. We optimized Kf,sat by fitting this parameter tothe inflow data using the parameter estimation featureof HYDRUS-1D. Cumulative infiltration rates were
found to be sensitive to Kf,sat, but not affected by theother hydraulic parameters. Our measured data wereclosely reproduced when using Kf,sat ⫽ 25 000 cm h⫺1(Fig. 11). The contribution of the macropore to thecomposite soil saturated conductivity (Ksat) can now be
estimated from Eq. [4]: for Kf,sat ⫽ 25 000 cm h⫺1 andwf ⫽ 1.8 ⫻ 10⫺5, Ksat was found to be equal to 0.45 cm
h⫺1, as compared with a value of only 0.08 cm h⫺1 asobserved in the absence of the macropore. The scaling
factor corresponding to this value was obtained from
Fig. 7 (␥w ⫽ 0.441). Three of the six VGM parameters
for the macropore domain were fixed (␪r ⫽ 0, ␪s ⫽ 0.5,l ⫽ 0.5), while the remaining parameters were determined by fitting lateral exchange and outflow data, on

220

VADOSE ZONE J., VOL. 2, MAY 2003

Fig. 12. Observed vs. predicted cumulative water transfer between
the macropore and matrix. In the dual permeability model, both␥w optimized and ␥w ⫽ 1 were used.

which they appear to have the largest impact. Leastsquares fitting gave ␣ ⫽ 0.1 cm⫺1; n ⫽ 2.
To highlight the model sensitivity to the scaling factor,
we compared the experimental data with the model
results obtained with both the optimized value of the ␥w
parameter, and for ␥w ⫽ 1. In both cases, the cumulative
water exchange predicted with the DPM did not matchthe experimental data perfectly (Fig. 12). However, thecalculated water transfer rate (given by the slope of the
cumulative curve) reproduced the observations reasonably well when using the optimized scaling factor, while
it persistently overestimated the experimental resultswith ␥w ⫽ 1. Even with the optimized ␥w, an offset ofabout 80 cm3 was still observed between the experimental and simulated cumulative curves. Assuming a porosity of 0.45 for the fine sand placed on top of the column,
this amount corresponds to water adsorbed in about
0.4 cm of sand after infiltration started. Therefore, thedisagreement between measured and calculated data in
Fig. 12 may be attributed to a slight difference in the
thickness of the sand layer between macropore and nonmacropore experiments.
Finally, we show in Fig. 13 the impact of the scaling
factor on cumulative outflow. The relative poor comparison with the experimental results indicates several limitations of the DPM model. For example, we observed
a constant effluent flux until the preferential flow wasinterrupted, whereas the simulation predicted a progressively increasing outflow rate. This behavior resulted
from having a constant flux at the inlet and progressively
less exchange of water between the macropore and thematrix as equilibrium is approached. An excellent match
was obtained between the experimental and simulated(with optimized ␥w) time intervals between the beginning of the infiltration process and the appearance ofthe first outflow. For ␥w ⫽ 1, the DPM failed to predict
the breakthrough time.

CONCLUSIONS
Experiments in artificial macroporous systems provide insight into preferential flow phenomena thatwould be difficult to achieve otherwise. In this study

Fig. 13. Observed vs. predicted cumulative outflow. In the dual permeability model, both ␥w optimized and ␥w ⫽ 1 were used.

we discriminate macropore flow and matrix flow. In
addition, the water exchange flux integrated over thecolumn was measured. The experimental data were described reasonably well with the dual-permeability
model, with most of the parameters determined independently of the observations. For the investigated geometry, the first-order approximation closely described
the water transfer data as long as an appropriate value
for the transfer rate coefficient (or equivalently for thecorrecting scaling factor ␥w) was used.
Similarly to the previous study by Gerke and van
Genuchten (1993b), we derived the transfer rate coefficient for water from the equivalent coefficient for solutetransfer through numerical simulations of the lateral
diffusion process. Our analysis showed that the scaling
factor ␥w, in general, is dependent on the hydraulic conductivity of the macropore domain.
ACKNOWLEDGMENTS
The research was conducted at the George E. Brown Jr.
Salinity Laboratory, Riverside, CA. P. Castiglione would liketo acknowledge the financial support from the University of
Naples, Italy, during his stay at Riverside as part of his Ph.D.
research. The authors also acknowledge the technical supportof J. Fargerlund and J.A. Jobes. This study was partially supported by NSF grant EAR-0106956.

REFERENCES
Ankeny, M.D., T.C. Kaspar, and A.A. Hussen. 1988. Designed for an
automated tension infiltrometer. Soil Sci. Soc. Am. J. 52:893–895.
Beven, K.J., and P. Germann. 1982. Macropores and water flow insoils. Water Resour. Res. 18:1311–1325.
Cushman, J.H. 1990. An introduction to hierarchical porous media.
p. 1–6. In J.H. Cushman (ed.) Dynamics of fluids in hierarchicalporous media. Academic Press, London.
Czapar, G.F., R. Horton, and R.S. Fawcett. 1992. Herbicide and tracer
movement in soil columns containing an artificial macropore. J.
Environ. Qual. 21:110–115.
Dykhuizen, R.C. 1990. A new coupling term for dual-porosity models.
Water Resour. Res. 26:351–356.
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 1storder water transfer term for variably saturated dual-porosity flow
models. Water Resour. Res. 29:1225–1238.

www.vadosezonejournal.org

Gerke, H.H., and M.Th. van Genuchten. 1996. Macroscopic representation of structural geometry for simulating water and solute movement in dual-porosity media. Adv. Water Resour. 19:343–357.
Hopmans, J.W., J. Simunek, N. Romano, and W. Durner. 2002. Inverse
modeling of transient water flow. p. 963–1008. In J.H. Dane andG.C. Topp (ed.) Methods of soil analysis. Part 4. SSSA Book Ser. 5.
SSSA, Madison, WI.
Hu, Q.H., and M.L. Brusseau. 1995. Effect of solute size on transport
in structured porous media. Water Resour. Res. 31:1637–1646.
Li, Y.M., and M. Ghodrati. 1997. Preferential transport of solute
through soil columns containing constructed macropores. Soil Sci.
Soc. Am. J. 61:1308–1317.
National Research Council. 2001. Conceptual models of flow andtransport in the fractured vadose zone. National Academy Press,
Washington, DC.
Othmer, H., B. Diekkruger, and M. Kutilek. 1991. Bimodal porosity
and unsaturated hydraulic conductivity. Soil Sci. 152:139–150.
Perroux, K.M., and I. White. 1988. Designs for disc permeameters.
Soil Sci. Soc. Am. J. 52:1205–1214.
Peters, R.R., and E.A. Klavetter. 1988. A continuum model for water
movement in an unsaturated fractured rock mass. Water Resour.
Res. 24:416–430.
Pruess, K. 1999. A mechanistic model for water seepage through thickunsaturated zones in fractured rocks of low matrix permeability.
Water Resour. Res. 34:1039–1051.
Simunek, J., M. Sejna, and M.Th. van Genuchten. 1998. The HYDRUS1D software package for simulating the one-dimensional movement

221

of water, 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.
Simunek, J., M. Sejna, and M.Th. van Genuchten. 1999. The HYDRUS2D software package for simulating two-dimensional movementof water, heat, and multiple solutes in variably saturated media.
Version 2.0. IGWMC-TPS-53. International Ground Water Modeling Center, Colorado School of Mines, Golden, CO.
Smettem, K.R.J., and C. Kirkby. 1990. Measuring the hydraulic properties of a stable aggregated soil. J. Hydrol. (Amsterdam) 117:1–13.
Valocchi, A.J. 1990. Use of temporal moment analysis to study reactive solute transport in aggregated porous media. Geoderma 46:
233–247.
van Genuchten, M.Th., and P.J. Wierenga. 1976. Mass transfer studies
in sorbing porous media. I. Analytical solutions. Soil Sci. Soc. Am.
J. 40:473–480.
van Genuchten, M.Th. 1980. A closed-form equation for predicting
the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am.
J. 44:892–898.
van Genuchten, M.Th., and F.N. Dalton. 1986. Models for simulating
salt movement in aggregated field soils. Geoderma 38:165–183.
Young, D.F., and W.P. Ball. 1995. Effects of column conditions on thefirst-order rate modeling of non-equilibrium solute breakthrough.
Water Resour. Res. 31:2181–2192.
Zimmerman, R.W., G.S. Bodvarsson, and E.M. Kwicklis. 1990. Absorption of water into porous blocks of various shapes and sizes.
Water Resour. Res. 26:2797–2806.