ELSEVIER

Advances in Water Resources, Vol. 19, No. 6, pp. 343-357, 1996
Copyright 0 1996 Published by Elsevier Science Ltd
Printed in Great Britain. All rights reserved
PII: S0309-1708(96)00012-7
0309-1708/96/$15.00 + 0.00

Macroscopic representation of structural geometry
for simulating water and solute movement indual-porosity media
Horst H. Gerke
Department of Soil Landscape Research, Center for Agricultural Landscape and Land Use Research, Eberswalder StraJe 84.
D-15374 Miincheberg, Germany

Martinus Th. van Genuchten
U.S. Salinity Laboratory, United States Department of Agriculture, Agriculture Research Service, 450 West Big Springs Road,
Riverside, California 92507-4617, USA(Received 6 October 1995; accepted 16 April 1996)

The structure of macroporous or aggregated soils and fractured rocks is generally
so complex that it is impractical to measure the geometry at the microscale (i.e., thesize and the shape of soil aggregates or rock matrix blocks, and the myriad offissures or fractures), and use such data in geometry-dependent macroscale flow
and transport models. This paper analyzes a first-order type dual-porosity model
which contains a geometry-dependent coefficient, p, in the mass transfer term tomacroscopically represent the size and shape of soil or rock matrix blocks. As areference, one- and two-dimensional geometry-based diffusion models were used tosimulate mass transport into and out of porous blocks of defined shapes. Estimatesfor ,0 were obtained analytically for four different matrix block geometries. Values
for ,0 were also calculated by directly matching analytical solutions of the diffusion
models for a number of selected matrix block geometries to results obtained withthe first-order model assuming standard boundary conditions. Direct matching
improved previous results for cylindrical macropore geometries, especially whenrelatively small ratios between the outer soil mantle and the radius of the innercylinder were used. Results of our analysis show that /3 is closely related to the ratioof the effective surface area available for mass transfer, and the soil matrix volume
normalized by the effective characteristic length of the matrix system. Using valuesof /3 obtained by direct matching, an empirical function is derived to estimate
macroscopic geometry coefficients from medium properties which in principle aremeasurable. The method permits independent estimates of p, thus allowing thedual-porosity approach eventually to be applied to media with complex and mixed
types of structural geometry. Copyright 0 1996 Published by Elsevier Science Ltd
Key words: preferential flow, variably-saturated structured media, dual-porosity
model, mass transfer coefficient, geometry coefficient.

conditions in the pressure head or solute concentration
perpendicular to the flow direction. Preferential flow
phenomena have been extensively studied8”7*49 because of
their implications on the hydrologic and contaminanttransport properties of soil/rock systems. Dual- (ordouble-) porosity approaches’15’ have been applied tovarious problems of water flow in fissured reservoirs,
solute transport in aggregated soils, and transient flow
and transport in variably-saturated structured media (for

INTRODUCTION
Structured soils and fractured rock formations often
contain macropores or fractures through which water andsolutes can move at considerably larger velocities than in
the porous matrix. The resulting microscopic or porescale heterogeneities in flow velocities may eventually
manifest themselves at the macroscopic scale in terms ofpreferential flow involving transient nonequilibrium
343

344

H. H. Gerke. M. Th. van Genuchten

an extensive review, see Gerke and van Genuchten15).
Conceptually, a dual-porosity medium is generally
assumed to involve two overlapping but interacting
pore domains having different hydraulic and transport
properties. Water and solute movement in dual-porosity
models (e.g., Dykhuizen’j and Jarvis et a1.24) is simulated
separately for the fracture and matrix pore systems, while
coupling terms are used to describe the transfer of waterand solutes between the two pore systems.
A commonality of dual-porosity approaches is thatthey require information about the geometry of the soil
or rock system being simulated, i.e., the size and shapeof matrix blocks or soil aggregates, including thegeometry of the interface area through which water andsolutes are exchanged between the two pore systems.
The interface area may be associated with the surface
structure of either the fracture pore system (e.g., planarcracks or cylindrical macropores), or of the matrix poresystem itself (e.g., rock matrix blocks or soil aggregates).
Naturally structured media often consist of a mixture of
matrix blocks having different shapes and sizes, thereby
leading to potentially very complex interfacial geometries as has been demonstrated recently by means ofcomputed tomography.20925’50
Mass transfer in dual-porosity models is frequently
approximated using a first-order differential equation.‘0*‘695’ More sophisticated formulations that consider a nalytical solutions of flow and transport into or
out of porous matrix blocks of known size and shape
are generally limited to relatively simple situations, suchas for well-defined macropore systems or idealized
matrix block geometries,34 steady-state water fl~w,~~
and for single pulses of water or solutes moving along
a fracture.14y52l53 Microstructure-type mode1s22s23135@
describing transport at both scales of the porousm e d i u m , i.e.,, macroscopic transport in the fracture
pore s yst e m1 and local-scale transport inside the aggregates or matrix blocks, mostly assume a single equivalent matrix block geometry (e.g., a sphere) for which
local-scale solutions of diffusion-type flow and transportmodels are available. Equivalent geometries may beobtained by deriving a volume-weighted average sphereradiu:-,32 using surface-to-volume ratios and distanceprobability functions for matrix diffusion obtained
from fractured rocks or aggregated soi1s,29’33 invoking
scaling methods involving shape factors,41,44 using
scaling methods based on ratios between the volume
and surface area of matrix blocks,54 or by means of‘block -geometry functions’ associated with planar, cylindrical, or spherical block shapes,2 among other methods.
Although conceptually of considerable interest, theapproaches above are still difficult to apply to fieldproblems since the microscopic model parameters are
not easily measured. A related issue, mostly stillunresolved, is the question of how to obtain adequate
microstructural information required to extend suchmodels from simple or artificial media like glass beads tofield problems.’ Consequently, the first-order masstransfer coefficient is at present still treated mostly as
an empirical parameter which must be calibrated byfitting to observed field data, or to solute breakthroughcurves obtained in displacement experiments.45 One of
the challenges of dual-porosity models is to identify
those local-scale processes that have the most impact at
the macroscale and to incorporate those processes in amacroscopic model in such a way that the model can beapplied to realistic field problems. With respect to a
structured dual-porosity medium, this means that thegeometric effects somehow must be represented by alimited few lumped parameters related to the bulk
volume of the medium. The aim of this paper is to derive
expressions that appropriately describe the structuralgeometry of a dual-porosity medium using parametersthat can be related to measurable properties.

DUAL-POROSITY MODEL AND FIRST-ORDERMASS TRANSFER
The dual-porosity model used in this study to simulate
one-dimensional vertical water flow and solute transport
in a variably-saturated structured medium was describedin detail by Gerke and van Genuchten.” The model
assumes that the properties of the bulk porous medium
can be characterized by two sets of local-scale properties: one set associated with the fracture pore system(subscript f), and the other with the matrix pore system(subscript m), as follows
E = It+& + (1 - Wf)&(la)

0 - M’& + (1 - wr)f&(lb)

4 = VIf + (1 - Qqn,(lc)where E is the porosity (L3Lp3), 6’ is the water content(L3L-3), q is the fluid flux density (LT-‘), and wf is therelative volumetric proportion of the fracture poresystem (0 < wf < 1). Assuming applicability of Darcy’s
law,45 water flow in the fracture and matrix pore regions
are described by a coupled pair of Richard’s equations:

@I
I‘,
+_-

1 - Wf(2b)where h is the pressure head (L), C is the specific watercapacity df?/dh (L-t), K is the hydraulic conductivity(LT-I), z is depth taken to be positive downward (L), t
is time (T), and Fw is the water transfer term (T-l) takenasrw = ‘Yw (hf - h,)(3)

in which cr, is a first-order mass transfer coefficientSimulation of water and solute movement in dual-porosity media(L-‘T-l) for water flow. The following expression for(Y, was derived by Gerke and van Genuchten”’ using ascaling procedure:Pow = ;;z %v&(~>(4)

where p is a dimensionless geometry-dependent coefficient, a is the characteristic length of the matrix
structure (L) (e.g., the radius of a spherical or solid
cylindrical aggregate, or half the fracture spacing in thecase of parallel rectangular voids), T,,, is a dimensionlessscaling coefficient for which an average value of 0.4 was
obtained empirically assuming parallel rectangular slabs(Fig. la), and K, is the effective hydraulic conductivity
function of the fracture-matrix interface evaluated interms of both h, and hf as follows
K, = @5K.&) + K,(&Jl(5)
Analogous to eqn (2), the transport of solutes in avariably-saturated dual-porosity medium is describedusing two coupled convection-dispersion equations:(6a)where D is the dispersion coefficient (L*T-‘), R is adimensionless retardation factor accounting for linear
equilibrium so tion, and Fs is the solute mass transfer
term (MLV3T- ‘p) formulated here in a slightly different
way as done previously in Gerke and van Genuchten.”rs = %(l - Wfh%l(Cf - CnJ +

1

r&f
r c
wm

rw L 0r,<o(7)

in which Q, is the first-order solute mass transfer
coefficient (T-l) of the form
o.J = $Da(8)

where Da is an effective diffusion coefficient (L*T-‘)
which represents the diffusion properties of the fracturematrix interface as well as other parameters and which,however, has not yet been evaluated analogously to
yWKa in eqn (4) for water transfer.
For steady-state water flow with qm = 0 (no flow inthe matrix pore system), the above variably-saturated
dual-porosity model reduces to the two-region mobileimmobile solute transport model of van Genuchten and
Wierenga.43

4IJG A=;t a:(cf-c,)(9b)where the water contents, br and r9,, as well as theretardation factors, RF and R& are expressed in terms of345

the bulk soil volume, i.e.,6, = ( 1‘19f = WBfi

-w,)&

(lOa; b)

The first-order mass transfer coefficient, (Y: (T-l), in (9b)characterizing diffusive exchange of solutes between
the mobile and immobile liquid phases is of the general
form

P%D:,
a: = --g-(11)

where 0: is the effective diffusion coefficient for solute
mass transfer into the immobile pore water region or thematrix pore system. We refer to Gerke and van
Genuchten” and van Genuchten4* for additional
discussions of the dual-porosity models presented here,
including possible ways of how the different parameters,other than p, may be measured or estimated.
GEOMETRY REPRESENTATION OF THE MASS
TRANSFER COEFFICIENT
Analytical evaluation

In order to relate the mass transfer coefficient, ai, to the
diffusion properties of soil aggregates or rock matrix
blocks, the first-order formulation given by eqns (9)and (11) may be compared with more comprehensive
geometry-based microstructure-type models for diffusion
into porous blocks having spherical, cylindrical, or
rectangular slab-type geometries. A variety of techniques
have been used to derive expressions for CX: in eqn (1 l),
including Laplace transform comparisons,31~39~44 moment
analysis,30l40 and the use of algebraic expansions.4 The
Laplace transform method may also be used to derive
expressions for the first-order water transfer coefficient,ow, in e n (4) by linearizing the horizontal flow
equation. 18 Such analyses lead to first-order mass transfer
coefficients in eqns (4), (8), and (11) which are of thegeneral form

cu,*

P

(12)
6, a*’
with P = D& in the two-region mobile-immobile typetransport model (9), D, in the dual-porosity solutetransport model (6), and -yWKa in the dual-porosity
water flow model (2). Thus, we hypothesize that valuesof @ can be derived that are identical for these threemodels, and hence a function only of the geometry.
Mass transfer coefficients obtained analytically using
Laplace transform comparisons differ only in the valueof the geometry-dependent coefficient, p. Most of thestudies mentioned above derived values for ,0 of 3 forrectangular slabs, 8 for solid cylinders, and 15 forspheres, although other studies obtained slightly different values.“y3’ For a hollow cylindrical dual-porosity
medium involving macropores surrounded by cylindrical soil or rock matrix mantles, van Genuchten and
Q,3(Y,~---_-

346

H. H. Gerke. M. Th. van Genuchten

Dalton4 derived the following expression for 0:2(Co - u2(13)

where co = (a + b)/b, in which b is the radius of thecylindrical macropore and a is the thickness of the soil
matrix mantle (Fig. lb).

For hollow cylindrical geometries the flux across thefracture-matrix interface (the surface of the innercylinder) is very similar to the soil water flux into aroot located in the center of a cylindrical soil segment.

The latter process has been widely studied.28 Analogous
to the water transfer term, rw, root water uptake isgenerally included as a sink term in the Richard’s
equation (2), and often formulated as an effective firstorder approach. Such an approach was actually used by
Moldrup et aL2’ who derived a soil resistance term byintegrating the flow equation for a cylindrical coordinate system between the inner and outer radius of thesoil mantle and assuming dO/dt to be constant (>O).

This approach leads to another equivalent expression
for the geometry-dependent coefficient, p. Using thenotation of this paper, this alternative expression is

b

a

/

/

/

(b)
a

b

(4
Fig. 1. Schematic illustration of porous blocks having (a) rectangular slab-type, (b) hollow cylindrical, (c) solid cylindrical, and(d) spherical geometries. Widths of 2a and 26 are associated with the matrix and fracture pore systems, respectively.

Simulation of water and solute movement in dual-porosity media

given by

d 0.5
P = lnK0) -c _
1

[

]

-1

;

b>1

(14)

Analytically derived mass transfer coefficientsare limited to basically four geometries for which
solutions of the diffusion equations are available.
Additional limitations are due to simplifying assumptions made during mathematical derivation of thecoefficients, e.g., neglecting higher-order terms in theseries expansions.

the soil mantle surrounding the macropore). Solving
eqn (15) for a uniform initial solute concentration,ca(r, 0), of zero, a constant concentration of unity at the
inner, ca(b, t) = 1, and a zero concentration gradient (noflow) at the outer boundary (at r = a + b), and substituting that solution into eqn (16) leads to9c,(Th) =’ - ““-2[; _ n=,1

x exp(-k&-b)
Direct matching
Alternatively, estimates for the mass transfer coefficient,crf, may also be obtained by directly matching solutionsof the first-order approach, i.e., of eqn (9b), with calculated average solute concentrations of soil aggregates orrock matrix blocks of known geometry. For this purpose
one could use analytical solutions of the equations
governing diffusion into spherical, rectangular (slabtype), or solid cylindrical aggregates, or those pertaining
to diffusion from hollow cylindrical macropores into thesurrounding soil matrix. The direct matching technique
was previously used by van Genuchten41 to compare
the diffusion properties of a soil aggregate with given
geometry and size with those having different shapes andsizes in a nonflowing system. His method is employed
here to obtain geometry-dependent coefficients forvarious matrix block sizes for use in the first-ordertransfer terms. The direct matching technique is demonstrated below for the case of a hollow cylindricalmacropore geometry.
Radial diffusion from the macropore through thesurface of the inner cylinder (Fig. lb) into thesurrounding mantle (assuming an infinite longitudinallength) is described bybLr<(a+b)(15)

where c,(r, t) is the local solute concentration inside thecylindrical soil mantle, r is the local radial coordinate,
and the subscript m indicates that Rk and 0:
represent matrix pore system properties which areconsistent with the transport model in eqn (9) at themacroscopic scale. We assume here that the soluteconcentration of the matrix pore system, c,, defines theaverage solute concentration, c,, of the liquid phase in
the hollow cylindrical (subscript h) porous matrix asfollows2

atb

E,(f) = (a+b)2 -62 bsrca(r’t)dr(16)

in which a is the characteristic length (the radial width of

347

(17)

where Jo, J1, Y,, and Y1 are Bessel functions and p,, arepositive roots ofJo(&) YI (&lCo) - J1 MO) Yo(&) = 0(18)
and where C~ is expressed in terms of dimensionless time,Th, as follows

J&f‘% = R$a;and where the subscript h on T and a is used to indicate ahollow cylindrical macropore.
Analogously, solutions for the average concentrations, i;,, as a function of dimensionless time, T, areavailable for spherical, solid cylindrical, and rectangular
geometries.9”’ Plots of the resulting curves of &(T)versus r for different aggregate geometries exhibit
relatively similar, but not identical, sigmoidal shapes.41
In the original paper, van Genuchten41 introduced a‘shape factor’, f, for conversion of one geometry into
an equivalent other geometry by matching the F,( )curves of two geometries through modification of thecharacteristic length a in the dimensionless time (e.g.,eqn (19)) pertaining to one of those geometries. The
shape factor was evaluated at an average value of T, =0.5, thus ignoring any observed dependency off on &.
A similar approach can be used for estimating shapefactors in the first-order transfer term.4’ Integratingeqn (9b) for c,(O) = 0 and cf( t) = 1 gives
T

I?~ E C,(To) = 1 - eXp(-To)

(20)

where
cY*t
7-O = $--&

(21)is the dimensionless time in the first-order model(subscript 0). The shape factor fh,o for converting ahollow cylindrical into the equivalent first-order rate
model was defined by van Genuchten4’ as

fh$ = ($

“2/4=o.1(22)in which the dimensionless times Th and 7. were both
evaluated at C~ = c,,, = 0.5. Substituting eqns (19) and348

H. H. Gerke. M. Th. van Genuchten(21) into eqn (22) shows that the shape factor,f, is related

to the geometry coefficient, /3, in the first-order transfer
term of the dual-porosity model by means of p = l/fz.
By using the direct matching method, we obtained
values for 0 of 3.5 for rectangular slabs, 11 for solid
cylinders, 22.7 for spheres (see also van Genuchten4’),
and the following empirical regression function(2 = 0.996) for a hollow cylindrical geometry:’ = 10.19lnil6Q]‘;

1 < <s < 100(23)

Figure 2 compares the geometry coefficients forhollow cylindrical matrix geometries as obtained bymeans of analytical derivation (eqns (13) and (14)) withthose obtained using eqn (23) which is based on the
direct matching approach. Equation (13) closely agrees
with the direct matching results in the limit of relatively
thick cylindrical soil mantles (large b), whereas eqn (14)approximates the direct matching results more closely inthe limit when <s -t 1. The results in Fig. 2 are aconsequence of simplifying assumptions invoked in thederivations of eqns (13) and (14), i.e., by having systems
with relatively few macropores per unit area (& > lo),or relatively dense root systems (<a < 5). The effect ofusing different values of fl in eqn (4) on soil water flow

simulations will be shown later by directly comparing
the dual-porosity model with equivalent two-dimensional models for different matrix geometries. We notethat the scaling procedure described by Gerke and van
Genuchten16 was employed to evaluate the hydraulicconductivity of the fracture-matrix interface, K,(h), in
eqn (4) using the analytically-derived geometry coefficient for parallel rectangular slabs u = 3).
COMPARISONS USING EQUIVALENTTWO-DIMENSIONAL MODELS
To study the geometry coefficients in the mass transfer

term (4) as evaluated by different methods, considertransient unsaturated water flow in both the fracture
and matrix pore systems. For convenience, and to
facilitate a direct comparison, vertical movement in thematrix pore system has been neglected. Results obtainedwith the dual-porosity model using the first-ordertransfer term will be compared with those generatedusing the SWMS_2D finite element code.” The flow
region was designed to imitate rectangular, hollow
cylindrical, and solid cylindrical geometries of thematrix blocks (Fig. 1). Hydraulic properties of relatively

j ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,(,, I ,,,,,,,”

0

10

20

30

40

50

60

70

80

“,““,

90

100

Outer-to-Inner-Radius Ratio, c0
Fig. 2. The geometry coefficient, @, as a function of the ratio, CO, between the outer and inner radius, assuming hollow cylindricalgeometry. Plotted are values and curves obtained by (a) analytical derivation using eqn (13), (b) direct matching using eqn (23), and(c) analytical derivation using eqn (14).

Simulation of water and solute movement in dual-porosity media

349

Table 1. Model parameters and initial conditions used in the twodimensional comparisons
n

(Y

4

4

(1/cm)

MatrixFracture

0.005
0.100

1.5

2.0

0.11
0.0

0.5.
0.5

coarse- and fine-textured soil materials (Table 1) wereassigned to the fracture and matrix pore systems,
respectively. Water flow for the parallel rectangular
fracture (slab) geometry was simulated by using a twodimensional cross-section (x-z plane) where x is thehorizontal and z is the vertical spatial direction. For
cylindrical geometries we assumed radial symmetry, andsimulated flow in a quasi three-dimensional region (r-zplane) using the flow equation for a cylindricalcoordinate system. In order to compare the results ofthe two-dimensional calculations with those obtainedwith the dual-porosity model (2), nodal water contents
at given simulation times were averaged over thehorizontal or radial distance and transformed into
equivalent pressure head values of the fracture andmatrix pore systems using the inverse of the water
retention function (e.g., van Genuchten et a1.47):
h = i (s;‘/” _

4 (2)
(cm)

1.05
2000

- 1000
-1000

4 (2)

0.277
0.005

for solid cylinders (subscript c); for rectangular prisms(subscript r) we assumed a square base area of (2a)2.
Flow systems containing a spherical geometry wereapproximated by assuming that the spherical matrix
blocks or soil aggregates were surrounded by thin
spherical mantles imitating the fracture pore system.
In case the spheres are completely surrounded by thefractures, the equivalent characteristic length of thefracture system, b:, isb: = a,[(1 - wf)-lj3 - l]

(28)
Assuming that the spherical matrix blocks are in contact
with each other (Fig. Id), the volume, V,,, of the top andbottom sphere-caps (subscript SC) of the fracture mantle
V, = : bi (3a, + 2b,)(29)

has to be distributed in the fracture pore system
resulting in a value of b, which is slightly larger thanbl as follows

1)“”’ - wfb3 _ 3a ‘a5 - wf 62 _ 3a2 -bb,+az
=0
s
s
s

where Q and n are empirical parameters, m = 1 - l/n,s, = (0 - &>l(& - 0,) is effective saturation, and 0, and8, are the residual and saturated water contents,
respectively. The procedure used for comparing modeling results is described in detail by Gerke and van
Genuchten.16
For the idealized geometries (Fig. 1) used here, thecharacteristic length, b (e.g., the half-width in the caseof a rectangular slab, or the radius in the case of acylinder or sphere) of the fracture pore system is related
to the relative proportion of the fracture system, wf, andthe characteristic length, a, of the matrix system by

Wf

s

Wf

(30)

Equation (30) can either be solved exactly, using aclosed-form trigonometric solution for cubic equations,or solved approximately in an iterative manner beginning with 6: obtained from eqn (28) as the initial
estimate. For example, assuming a = 1 cm and wr = 0.05,
one obtains bf = 0.017245 cm and b, M 0.0174cm. The
values of b for the other geometries (Table 2) were used todesign equivalent two-dimensional flow regions andassociated numerical finite element discretizations.
To illustrate physical nonequilibrium during transient
flow in an unsaturated dual-porosity-type structured
porous medium, the simulation example used previously’5~‘6 was selected. We assumed the applicationof water at a constant rate of 50 cm/day to a 40-cm-deepdual-porosity medium having an initial uniform pressure
head of h,i = h,,i = -1000 cm. Water was allowed toinfiltrate exclusively into the fracture pore system(wr = 0.05), such that a local infiltration rate ofqf,o = l000 cm/day was imposed at the upper and afree-drainage condition at the lower boundary of the

Bt = al&f

for rectangular slabs (subscript l), by

for hollow cylinders (subscript h), and bybcir = a,;,[(1 - wf)-li2 - l]KS
(cm/day)(27)

Table 2. Model parameters in first-order transfer term and two-dimensional representation

Rectangular slab
Solid cylinder
Hollow cylinderP
analytic
3
8
2.42

P
matching
3.5
11.0
1.5

a

(cm)
1
1

1

xv

(A)
0.02563
0.02598
0.28800

0.4
0.4
0.4

KW

(cm/day)
0.01
0.01
0.01

0.1
0.1
0.1

Simulation of water and solute movement in dual-porosity media
fracture pore system. Parameters of the soil water
retention function ((Y, n, O,, O,), the hydraulic conductivity at saturation, K,, the initial pressure head, hi, andwater content, ei, and the first-order mass transfer
coefficients can be found in Tables 1 and 2, respectively.
Figure 3 shows calculated two-dimensional pressure
head distributions at t = 0.04 days for (a) a parallel
rectangular slab, (b) a hollow cylindrical geometry, and(c) a solid cylindrical geometry assuming relatively low(K,,a = 0.01 cm/day) water transfer rates between thefracture and matrix pore systems. Note that thehorizontal axis is enlarged 20 times. The simulation
examples in Fig. 3 illustrate significant preferential flow
and pressure head nonequilibrium between the fracture
and matrix pore domains for all three geometries. Afterabout 1 hour (t = 0.04 days) of constant infiltration,water moved downward to a depth of 23cm in therectangular pore system (Fig. 3a), to z = 35 cm in thehollow cylindrical pore system (Fig. 3b), and toz = 15cm in the fracture pore system surrounding thesolid cylindrical soil matrix (Fig. 3c). Since all model
parameters, including the relative volume, wf, of thefracture pore system and the initial and boundary
conditions were identical, differences in calculatedpressure distributions were caused only by the geometrydependent rates of water transfer from the fracture into351

the matrix pore system. Two-dimensional model simulations were also carried out for the case of relatively highwater transfer rates assuming K,+ = 0.1 cm/day (results
not shown here).
Horizontally averaged two-dimensional simulation
results are compared with results obtained from thedual-porosity model for relatively small water transfer
rates (I& = 0.0 1 cm/day) assuming matrix structures
involving rectangular slabs (Fig. 4), hollow cylinders(Fig. 5), and solid cylinders (Fig. 6). As was discussed
previously,‘6 the first-order term generally underestimates the transfer rate at the infiltration front due tofailure of the linear first-order term to capture the highly
nonlinear transient flow process near the front (see also
Zimmerman et al.“). However, except in areas close to
the infiltration front, the proposed linear term yielded
relatively good approximations for the transfer of water(parts b of Figs 4-7), so that the calculated pressure
head (parts a) and water content profiles (parts c) stillclosely matched the reference results. The match
generally improved when allowing higher water transfer
rates by using KS,, = 0.1 cm/day (cfX, Figs 5 and 7 for ahollow cylindrical geometry and p = 1.5). The results in
Figs 4-6 reveal similar dependencies of the transfer rate
on p among the different geometries, i.e., the larger thegeometry-dependent mass transfer rates along the

0

10

(a)
Pressure Head, h (cm)

Water Transfer Rate, I’, (l/day)

Volumetric Water Content, 0

Fig. 5. Pressure head, h (a), water transfer rate, Pw (b), and total volumetric water content, 0 (c), versus depth, z, at t = 0.04 days as

calculated with the one-dimensional dual-porosity (Dual) and two-dimensional (2-D) models, assuming a hollow cylindricalgeometry and relatively small water transfer rates (K,+ = 0.01 cm/day); results are for geometry coefficients derived analytically(p = 2.4) and by direct matching (p = 1.5).

352

H. H. Gerke, M. Th. van Genuchten

sequence: ‘hollow cylinders < rectangular slabs < solid
cylinders’, the closer the match with the reference
results. Thus, the closest correspondence was obtained
for the solid cylindrical geometry (Fig. 6) assuming
p=g
In t ontrast to the above, the analysis becomes more
complicated when, for a particular geometry, different
methods for evaluating the geometry coefficient, p, arecomp ared. Because of the approximate nature of thefirst-order approach, and the nonlinear flow regimes
involved, the comparisons will differ for each particular
set o f hydraulic parameters and at each particularsimul ation time. We again present results at t = 0.04 foran infiltration event, and selected a situation with fairly
typical results. In case of a hollow cylindrical geometry(Figs 5 and 7), the geometry coefficient obtained bydirect matching, i.e., p = 1.5 calculated using (23),substantially improved the correspondence with the twodimen sional reference results, as compared to theanaly tically derived value (p = 2.42) using eqn (13).
These results confirm that eqn (13) is limited to dualporosity systems with relatively large values of Q., On
the other hand, results obtained with different values offl (as estimated analytically or using direct matching)were quite similar for rectangular slabs (Fig. 4), whereas
for solid cylinders the results using the analytically derived

value turned out to be superior as compared to thedirect matching method (Fig. 6). The results in Figs 4-6
only in part reflect the different sensitivities of the masstransfer rate on 0 for different matrix geometries. The
results also support the assumption that a single
constant value of 0.4 for the scaling coefficient, T,,,, ineqn (4) in combination with an arithmetically averaged
K, (eqn (5)) can be reasonably well applied to all threegeometries.
EXTENSION TO MEDIA WITH COMPLEXSTRUCTURAL GEOMETRIESTo be able to extend the approach to porous media withmore complex structural geometries, we assume thatsomehow a macroscopic geometry-dependent coefficient
may be defined which effectively captures the overall
diffusion properties of the structured system. So far, /3 in
eqn (4) was quantitatively related to the shape of a welldefined matrix structure. Alternatively, /3 could be
related also to the ratio between the surface area of
the interface between the fracture and matrix poresystem and the volume of the matrix blocks. Both ofthese quantities are, at least in principle, measurable.
The surface-to-volume ratio has been used in various

(c)

(b)

wJltlI)IlIIJIIIl~lIII’

-100

-10

Pressure Head, h (cm)

-1 0

5

10

15

Water Transfer Rate, Tw, (l/day)

0.30

0.35

0.40

0.45

0.50

Volumetric Water Content, 8Fig. 6. Pressure head, h (a), water transfer rate, I’, (b), and total volumetric water content, fI (c), versus depth, z, at t = 0.04 days as
calculated with the one-dimensional dual-porosity (Dual) and two-dimensional (2-D) models, assuming a solid cylindrical geometry
and relatively small water transfer rates (& = 0.01 cm/day); results are for geometry coefficients derived analytically (J = 8) and bydirect matching @ = 11).

Simulation of water and solute movement in dual-porosity media
approaches to obtain geometric information of soil orrock structures?9z3’” mostly from microscopic points ofview.
The ratio between the surface area, A, (L2), of thematrix-fracture interface (e.g., the soil aggregate surface)
and the volume, Vm (L3), in the case of a matrix blockformed by parallel rectangular slabs (plates), is given by

where ah = bh(CO - 1) is the characteristic length of thematrix structure.
For application to media with arbitrary geometries
and various characteristic diffusion lengths, we introduce a dimensionless surface-area-to-volume ratio, c(-),
normalized per unit length of the matrix structure forrectangular and solid cylindrical geometry
*

A
1
m,l_
v, -z

(31)

which assumes that the contributions of the sides (topand bottom ends) can be ignored. For solid cylinders
and rectangular prisms we obtain similarly(32)

while for spheres the surface-area-to-volume ratio is
Am,s
3
v,=Q,

(33)

For hollow cylinders, assuming that water and solute
mass transfer is restricted to the surface area of the innercylinder, one obtains
Am h
2
2
A=
vm

353

bh(&f - 1) = ah(cO f 1)’

co > 1

(34)

c=$a*

(35)

m

and for hollow cylindrical structures<=&Tyco > 1

where Ah is the effective surface area and a* is theeffective length of the matrix pore system, both being
lumped macroscopic scale parameters which integrate
over all geometries of a defined volume of the structured
medium.
By assuming lumped parameters, we may plot thegeometry coefficient, p, versus the normalized surfacearea-to-matrix-volume ratio, r;, obtained from thedirectly-matched shape factors for various sizes andshapes for conversion into the equivalent first-order rate
model using eqn (22). Results are shown in Fig. 8 where,
for comparison, the analytically-derived coefficients for

15-

(a)
20-sm. 1 . . ,..“‘I . . ,““‘. . 1
-1000
-100
-10
-10

Pressure Head, h (cm)

I
5

(b)
” ’ ‘1””
10

Water Transfer Rate, I’, (l/day)

Volumetric Water Content, 8Fig. 7. Pressure head, h (a), water transfer rate, lTW (b), and total volumetric water content, 8 (c), versus depth, z, at t = 0.04 days as
calculated with the one-dimensional dual-porosity (Dual) and two-dimensional (2-D) models, assuming a hollow cylindricalgeometry and relatively large water transfer rates (KS+ = 0.1 cm/day); results are for geometry coefficients derived analytically(@ = 2.4) and by direct matching (p = 1.5).

H. H. Gerke, M. Th. van Genuchten

354

which is simpler but also somewhat less accurate thaneqn (37).
The proposed functions in eqns (37) or (38) relate thegeometry-dependent coefficient, p, of the macroscopic
scale mass transfer terms in eqn (12) to macroscale
properties of an arbitrarily structured medium. The
functional interpretation of structure used here focuses
only on the size and shape of the fracture-matrix
interface area through which diffusive mass transfer
actually takes place. An important lumped parameter inthis context is the effective surface area, A;, which
integrates all ‘internal’ structural surfaces over thevolume of the porous medium; its value also contains
contributions from imperfectly developed (sub-) structures such as small fissures in spherical or cylindricalmatrix blocks. Values of the normalized surface-area-tovolume ratio smaller than unity represent structured
soils or fractured rocks in which the ‘internal’ surface
area for diffusive mass transfer is primarily concave(e.g., hollow cylinders). For porous media having valuesof 5 > 1, the fracture-matrix interface areas are predominantly convex (e.g., spheres). The effective characteristic length parameter, a*, may be interpreted to be
the smallest average distance for diffusion between theplane sheets, solid cylinders, and spheres have been
included. Figure 8 shows that the difference in the valueof p between structures that have the same reduced
surface-to-volume ratio, are much smaller than differences in 0 with respect to the method of evaluation(e.g., ,0 = 22.7 for a sphere and 20.7 for a cube). This
fact permits one to describe the geometry coefficient
over a wide range of c-values by the fitted empiricalexpression[0191n{16(~-1)}]-z
3.6~“’Ul + u2c +0.0198<<<1

(37)

lLc12
2<s110u3c2

in which u, = 11.4275, u2 = -7.4438, and u3 = 3.5473
and where 2 = 1. For 1 < < < 2, the appropriate
geometry coefficients may be obtained by using theinterp olation function suggested in eqn (37). Alternatively. the dependency may be described also by[

0.0198<<<[0491n{16(f-1))1-2
0.65-20.09 + <1

_
_’

1

(38)

I<r<lO

_

,
_‘_

.
_‘_

(

..,..
:

‘.

,

,.
:

.
‘.

.
_’

,.

‘.

1

.’

7 Equation

.

37

1

Normalized Surface-to-Volume Ratio, sFig. 8. Plot of the geometry coefficient 0 as a function of the normalized effective surface-to-volume ratio, c = a*Ak/ If,,,, as given byeqn (37), as well as p-values obtained by analytical evaluation and direct matching.

Simulation of water and solute movement in dual-porosity media

fracture-matrix interface and the center of the matrix
blocks perpendicular to the macroscopic scale flow
direction. The results in Fig. 8 suggest that the effect of pon flow and transport calculations will be most apparent
when c is close to 1.0. The effect of ,LI declines for c < 1
where /3 approaches a relatively constant value of 0.5(i.e., the asymptotic behavior of fl for c < 0.1 isdemonstrated in Fig. 2) and total mass transfer rates
decrease because of a reduced effective surface area fortransfer, Ah. For c > 1, the effect of p may also decline
if, in the case of relatively large mass transfer rates, thesituation of a homogeneous (unstructured) medium isapproached which will not exhibit any preferential flow.
We emphasize that the coefficient /I as estimated witheqn (37) or eqn (38) accounts for the shape dependencyof the mass transfer process, and indirectly still is afunction of the characteristic length, a*. The macroscaleparameters a*, b*, and wr, however, are no longer related
to each other in well-defined ways as for idealized
geometries using eqns (31)-(34); they have to bemeasured independently.
POSSIBLE PARAMETERIZATION OF
STRUCTURAL GEOMETRY
Most approaches involving microstructure-type mobileimmobile models assume prior knowledge of the type ofgeometry, or rely on detailed micromorphologicalmeasurements.29V32V33,54 In addition, an equivalent diameter concept has been proposed for tubular macroporesin saturated media.12 However, we believe that realistic
simulations of preferential flow require the use of
macroscopic scale flow and transport parameters forboth pore systems.
The macroscopic approach presented in this paper
permits structural effects to be included without assuming any type of geometry a priori, and hence follows amore functional view as suggested previously by Boumaet al7 The proposed macroscopic scale representation of
the structural geometry leading to diffusive masstransfer relies on one of the most basic properties of astructured medium, i.e., the fracture-matrix interface
area. The extent, size and shape of the interface, as wellas its physical, chemical, and biological properties (e.g.,
those of fracture coatings) are major factors determining
transport processes in aggregated soils or fractured
rocks. Of relevance here are, as discussed earlier,t6 not
only the interfacial permeability, K,, but also theeffective diffusion coefficient, D,, the chemical composition and sorption sites, biological activities along
fractures, as well as other properties occurring at ornear such interfaces.
Unfortunately, the effective fracture-matrix interface
area is not easily determined experimentally. Undercertain conditions the area and shape of the interface
may actually depend on the water content and flow rate,355

while also partial wetting of aggregate surfaces orswelling and shrinking of the matrix may occur.
Recently, methods have been suggested to calculate
the surface areas of structured soils using dyes that stain
macropore area&‘* or by using an image analyzer.6
Measurements of air permeability have also beensuggested36 to determine the effects of pore geometry
on transport processes in macroporous soils. Birgersson
et al.’ used two different tracers at the same time toestimate the flow-wetted surface of fractured rocks into
which the dyes can diffuse and sorb. In another study,
McKay et a1.26 used nonreactive solute and colloid-sizedbacteriophage tracers to reveal large differences in transport rates between the two tracers because of exclusion of
the large tracers from diffusing into the matrix.
The second macroscopic property, the effective
characteristic length of the matrix structure, a*, may
be determined using diffusion experiments. As a reasonable initial estimate for a*, one-half of the volumesurface mean diameter d (as defined by Gregg andSing”) may be used:&,A 1 c(d)a* =2=2C(n&)(39)

where ni are the number of particles (cubes) of edge
length li of a specific volume. Often, the values of thelengths li in eqn (39) have to be estimated, for example,
by using the mean projected diameter,” which is thediameter of a circle having the same area as the
projected image of the particle when viewed in a
direction normal to the plane of greatest stability. One
could apply a similar approach to aggregate mixtures so
as to obtain a weighted average of the geometry andlength coefficient. Mass transfer rates could also be
obtained by multiple calculations using observed shapeand size distributions of the matrix pore system, andapplying a volumetric weighing procedure.
The above modeling concept offers also the possibility
to combine mass transfer with preferential flow insituations involving partial wetting of matrix surfaces ofcracked clays or fractured rocks, for instance, byassuming
A:, = A,$,

OIS<l

(40)

where s(t) is a dimensionless wetting parameter which
may be a function of time or certain soil properties.2’
The effects of swelling and shrinking in soils in thisapproach could be described directly by means of theparameter wf.

CONCLUSIONS
Preferential flow of water in a dual-porosity system ofvarious matrix geometries was simulated using a firstorder type water transfer term. Direct comparisons with

356

H. H. Gerke, M.. Th. van Genuchten

reference two-dimensional model simulations demonstrate (1) inherent limitations in the approach because of
simplification of highly nonlinear transfer processes into
an approximate first-order rate model, and (2) sensitivityof the mass transfer coefficient in the dual-porosity
model to the exact geometry of the matrix pore system(e.g., rectangular slabs, and hollow or solid cylinders).
The comparisons revealed similar tendencies for all
three geometries when using an appropriate value of thecoefficient, & describing the specific geometry of thematrix pore system. By scaling the diffusion properties
and relating the geometry coefficient to the ratio
between effective surface area at the fracture-matrix
interface and matrix volume, the approach may be
extended to media with arbitrary and complex internal
geometries.
The methodology in this paper gives qualitatively agood approximation of more realistic two-dimensional
representations of the mass transfer process. An
important advantage of the scaling approach is thatthe method can be readily incorporated in multidimensional numerical models (e.g., Vogel et ~1.~‘).
Also, the model may be used as a tool to represent or
better understand preferential flow and transportphenomena in a mechanistic (process-based) fashion
with a reasonable numerical effort. The approachenhances the utility of dual-porosity models in that
otherwise empirical fitting parameters now can be
expressed in terms of measurable soil/rock properties.
Also, the approach is shown to be not limited to porous
media containing a single type of aggregate shape andsize (geometry) or uniformly-sized aggregates, nor is ana priori qualitative description of the geometry of theporous matrix required. The approach is also notlimited to experimentally determined aggregate size
distributions, which often are difficult to obtain for
specific soils. Nevertheless, the study of interfacial
macroscale properties will require new or improved
experimental methods for determining effective surface
areas, i.e., those of soil aggregates or fracture walls.
The approach may help to further increase theapplicability of the dual-porosity approach to analyzetransport in media with complex and widely-differing
structural geometries.
REFERENCES
Barenblatt, G. I., Zheltov, I. P. & Kochina, I. N., Basicconcepts in the theory of seepage of homogeneous liquids
in fissured rocks. J. Appl. Math. Mech. Engl. Transl., 24(1960) 1286-1303.
Barker, J. A., Block-geometry functions characterizing
transport in densely fissured media. J. Hydrol., 77 (1985)
263-279.

Birgersson, L., Moreno, L., Neretnieks, I., Widen, H. &Agren, T., A Tracer Migration Experiment in a small
fracture zone in granite. Water Resour. Res., 29 (1993)
3867-3878.

4. Bolt, G. H., Movement of solutes in soil: Principles ofadsorption/exchange chromatography. In Soil Chemistry,
B. Physico-Chemical Models, Dev. Soil Sci. 5B, ed. G. H.

Bolt, 1979, pp. 295-348.
5. Booltink, H. W. G. & Bouma, J., Physical and morphological characterization of bypass flow in a well-structured
clay soil. Soil Sci. Soc. Am. J., 55 (1991) 1249-1254.
6. Booltink, H. W. G., Hatano, R. & Bouma, J., Measurement and simulation of bypass flow in a structured clay
soil: A physico-morphological approach., J. Hydrol., 148(1993) 149-168.
7. Bouma, J., Belmans, C. F. M. & Dekker, L. W., Water
infiltration and redistribution in a silt loam subsoil withvertical worm channels. Soil Sci. Soc. Am. J., 46 (1992)
917-921.
8. Brusseau. M. L. & Rao. P. S. C.. Modeling solutetransport in structured soiIs: A review. Geoderma, 46(1990) 169-192.
9. Carslaw, H. S. & Jaeger, J. C., Conduction of Heat in Solids,
2nd edition. Oxford University Press, London, 1959.
10. Coats, K. H. & Smith, B. D., Dead end pore volume anddispersion in porous media. Soc. Petrol. Eng. J., 4 (1964)
73-84.

11. Crank, J., The Mathematics of Diffusion. Oxford Uuiversity
Press, London, 1956, 347 pp.
12. Dunn, G. H. & Phillips, R. E., Equivalent diameter ofsimulated macropore systems during saturated flow. Soil

Sci. Soc. Am. J., 55 (1991) 1244-1248.
13. Dykhuizen, R. C., Transport of solutes through unsaturated fractured media. Water Res., 21 (1987) 1531-1539.
14. Dykhuizen, R. C., A new coupling term for doubleporosity models. Water Resour. Res., 26 (1990) 351-356.
15. Gerke, II. H. & van Genuchten, M. T., A dual-porosity
model for simulating the preferential movement of waterand solutes in structured porous media. Water Resour.
Res., 29 (1993), 305-319.
16. Gerke, H. H. & van Genuchten, M. T., Evaluation ofa first-order water transfer term for variably-saturated
dual-porosity models. Water Resour. Res., 29 (1993) 12251238.
17. Gish, T. J. & Shirmohammadi, A.(eds.), Preferential flow.
Proceedings of the national symposium, 16- 17 December
1991, Chicago, Illinois, American Society of AgriculturalEngineers, St. Joseph, Michigan, 1991, 408 pp.
18. Ghodrati, M. & Jury, W. A., A field study using dyes tocharacterize preferential flow of water. Soil Sci. Soc.
Am. J., 54 (1990) 1558-1563.

19. Gregg, S. J. & Sing, K. S. W., Adsorption, Surface Area andPorosity, 2nd edition. Academic Press, London, 1982.

20. Heijs, A. W. J., de Lange, J., Schoute, J. F. T. & Bouma, J.,

Computed tomography as a tool for non-destructive
analysis of flow patterns in macroporous clay soils.

Geoderma, 64 (1995) 183-196.
21. Hoogmoed, W. B. & Bouma, J., A simulation model for
predicting infiltration into cracked clay soil. Soil Sci. Soc.
Am. J., 44 (1980) 458-461.
22. Hornung, U., Homogenization of miscible displacementin unsaturated aggregated soils. In Workshop on Composite Media and Homogenization Theory, ICTP, Tricste,
1990, Progress in Nonlinear Differential Equations andApplications, eds., DalMaso, G. & Dell’Antonio, G. F.
BirkhHuser, Boston, 1991, pp. 143-153.
23. Hornung, U. & Showalter, R. E., Diffusion models for
fractured media. J. Math. Anal. Appl., 147 (1990) 69--80.
24. Jarvis, N. J., Jansson, P.-E., Dik, P. E. & Messing, I.,
Modelling water and solute transport in macroporous soil.
I. Model description and sensitivity analysis. J. Soil Sci.,
42 (1991) 59-70.

Simulation of water and solute movement in dual-porosity media
25. Joschko, M., Graff, O., Muller, P. C., Kotzke, K.,

Lindner, P., Pretschner, D. P. & Larink, O., A nondestructive method for the morphological assessment ofearthworm burrow systems in three dimensions by X-ray
computed tomography. Biol. Fertil. Soils, 11(199 1) 88-92.
26. McKay, L. D., Gilham, R. W. & Cherry, J. A., Field
experiment in a fractured clay till. Water Resour. Res., 29357

42. van Genuchten, M. T., New issues and challenges in soil
physics research, Transactions, 15th World Congress ofSoil Science, Acapulco, Mexico, July 10-16, 1994, Volume
1, 5-27, Int. Soc. of Soil Sci. and Mexican Soc. of Soil Sci.,
1994.
43. van Genuchten, M. T. & Wierenga, P. J., Mass transfer
studies in sorbing porous media. I. Analytical Solutions.

(1993) 3879-3890.
27. Moldrup, P., Rolston, D. E., Hansen, J. A. & Yamaguchi, T.,

Soil Sci. Soc. Am. J., 40 (1976) 473-480.
44. van Genuchten, M. T. & Dalton, F. N., Models for

A simple, mechanistic model for soil resistance to plantwater uptake. Soil Sci., 153 (1992) 87-93.
28. Molz, F. J., Models of water transport in the soil-plantsystem: A review. Water Resour. Res., 17 (1981) 1245-1260.
29. Neretnieks, I. & Rasmuson, A., An approach to modellingradionuclide migration in a medium with strongly varying
velocity and block sizes along the flow path. Water Resour.

simulating salt movement in aggregated field soils.

Res., 20 (1984) 1823-1836.
30. Parker, J. C. & Valocchi, A. J., Constraints on the validityof equilibrium and first-order kinetic transport models in
structured soils. Water Resour. Res., 22 (1986) 399-407.
31. Raats, P. A. C., Transport in structured porous media.
Proc. Euromech. (Delft), 143 (198 1) 221-226.

32. Rao, P. S. C., Jessup, R. E. & Addison, T. M.,

Experimental and theoretical aspects of solute diffusionin spherical and non-spherical aggregates. Soil Sci., 133(1982) 342-349.
33. Rappoldt, C., The application of diffusion models to an
aggregated soil. Soil Sci., 150 (1990) 645-661.
34. Rasmuson, A. & Neretnieks, I., Exact solution for
diffusion in particles and longitudinal dispersion inpacked beds. J. Am. Inst. Chem. Eng., 26 (1980) 686-690.
35. Rasmuson, A., Gimmi, T. & Fliihler, H., Modelingreactive gas uptake, transport, and transformation inaggregated soils. Soil Sci. Soc. Am. J., 54 (1990) 1206-

1213.
36. Roseberg, R. J. & McCoy, E. L., Measurement of soil
macropore air permeability. Soil Sci. Soc. Am. J., 54 (1990)
969-974.
37. Simunek, J., Vogel, T. N. & van Genuchten, M. T., The

SWMS_2D code for simulating water flow and solutetransport in two-dimensional variably-saturated media,
Research Report No. 126, U. S. Salinity Laboratory,
USDA, ARS, Riverside, CA, 1992.
38. Skopp, J., Gardner, W. R. & Tyler, E. J., Miscible
displacement in structured soils: Two region model withsmall interaction. Soil Sci. Soc. Am. J., 45 (1981) 837-842.
39. Sudicky, E. A., The Laplace Transform Galerkin Technique for efficient time-continuous solution of solutetransport in double-porosity media. Geoderma, 46 (1990)
209-232.

40. Valocchi, A. J., Use of temporal moment analysis to study
reactive solute transport in aggregated porous media.
Geoderma, 46 (1990) 233-247.

41. van Genuchten, M. T., A general approach for modelingsolute transport in structured soils. Memoires Znt. Assoc.
Hydrogeol., 17 (1985) 513-526.

Geoderma, 38 (1986) 165-183.

45. van Genuchten, M. T. & Gerke, H. H., Dual-porosity
models for simulating solute transport in structured media,
Proceedings, Colloquium ‘Porous or Fractured Unsaturated Media: Transports and Behaviour’, October 5-9,
1992, Monte Verita, Ascona, Switzerland, 1993, pp. 182205.
46. van Genuchten, M. T., Tang, D. H. & Guennelon, R.,
Some exact and approximate solutions for solute transport
through soils containing large cylindrical macropores.
Water. Resour. Res., 20 (1984) 335-346.

47. van Genuchten, M. T., Leij, F. J. & Yates, S. R., The
RETC-code for quantifying the hydraulic functions of
unsaturated soils, Report EPA/600/2-91/065, Washington,
DC, 1991.
48. Vogel, T. N., Zhang, R., Gerke, H. H. & van Genuchten,
M. T., Modeling two-dimensional water flow and solutetransport in heterogeneous soil systems. In Hydrologic
Investigation, Evaluation, and Ground Water Modeling, eds.
Y. Eckstein & A. Zaporozec. Proceedings of Industrial andAgricultural Impacts on the Hydrologic Environment, The
Second USA/CIS Joint Conference on Environmental
Hydrology and Hydrogeology, 16-21 May 1993,
Washington, D. C., 1993, pp. 279-302.
49. Wang, J. S. Y., Flow and transport in fractured rocks.
Reviews of Geophysics, Supplement (1991) 254-262.
50. Warner, G. S., Nieber, J. L., Moore, I. D. & Geise, R. A.,
Characterizing macropores in soil by computed
tomography. Soil Sci. Soc. Am. J., 53 (1989) 653-660.
51. Warren, J. E. & Root, P. J., The behaviour of naturally
fractured reservoirs. Soc. Petrol. Eng. J., 3 (1963) 245-255.
52. Zimmerman, R. W. & Bodvarsson, G. S., Integral method
solution for diffusion into a spherical block. J. Hydrol.,111 (1989) 213-224.
53. Zimmerman, R. W. & Bodvarsson, G. S., An approximate solution for one-dimensional absorption in unsaturated porous media. Water Resour. Res., 25 (1989) 1422-

1428.
54. Zimmerman, R. W., Bodvarsson, G. S. & Kwicklis, E. M.,

Adsorption of water into porous blocks of various shapes
and sizes. Water Resour. Res., 26 (1990) 2797-2806.
55. Zimmerman, R. W., Chen, G., Hadgu, T. & Bodvarsson,
G. S., A numerical dual-porosity model with semianalytical treatment of fracture/matrix flow. Water Resour. Res.,
29 (1993) 2127-2137.