WATER RESOURCES RESEARCH, VOL. 40, W07409, doi:10.1029/2004WR003285, 2004

Numerical evaluation of a second-order water transfer term forvariably saturated dual-permeability models
J. Maximilian Köhne and Binayak P. Mohanty
Department of Biological and Agricultural Engineering, Texas A&M University, College Station, Texas, USA

Jirka Simunek
Department of Environmental Sciences, University of California, Riverside, California, USA

Horst H. Gerke
Center for Agricultural Landscape and Land Use Research, Institute of Soil Landscape Research, Müncheberg, Germany
Received 18 April 2004; accepted 26 May 2004; published 28 July 2004.

[1] Dual-permeability models contain a lumped mass transfer term that couples equations

for water flow in the soil matrix and fracture systems. Linear first-order transfer terms
cannot accurately calculate water transfer between the domains during early times ofpressure head nonequilibrium. In this study, a second-order equation for water transfer into
spherical rock matrix blocks [Zimmerman et al., 1993, 1996] was adapted and evaluatedcomprehensively for water transfer into and out of variably saturated soil matrix blocks ofdifferent hydraulic properties, geometries, and sizes, for different initial and boundary
conditions. Numerical solutions of the second-order term were compared with respective
results obtained with a first-order term and a one-dimensional horizontal flow equation.
Accurate results were obtained after implementing two modifications in the second-orderterm. First, the hydraulic conductivity was calculated as a weighted arithmetic averageof conductivities that used pressure heads in matrix and fracture. For a time-variable
pressure head boundary condition, a fixed weighting factor of 17 could be appliedirrespective of texture, initial condition, and matrix block size up to 10 cm. Second, if
direction of water transfer changed (to or from matrix), the initial pressure head was resetto the equilibrium pressure head at the time of transfer reversal. The modified secondorder term was implemented into a dual-permeability model, which closely approximatedreference results obtained with a two-dimensional flow model. For rectangular slab-type
or comparable geometry of soil matrix, the modified second-order term considerably
improves the accuracy of dual-permeability models to simulate variably saturated
INDEX TERMS: 1875 Hydrology: Unsaturated zone; 1866 Hydrology: Soil
preferential flow in soil.
moisture; 1829 Hydrology: Groundwater hydrology; KEYWORDS: dual-permeability model, preferential flow,second-order transfer term
Citation: Köhne, J. M., B. P. Mohanty, J. Simunek, and H. H. Gerke (2004), Numerical evaluation of a second-order water transfer
term for variably saturated dual-permeability models, Water Resour. Res., 40, W07409, doi:10.1029/2004WR003285.

1. Introduction[2] For simulating water exchange between regions ofcontrasting hydraulic properties in a fractured porous
medium, many dual-permeability models rely on a linear
first-order, ‘‘Warren-Root’’-type [Warren and Root, 1963]
water mass transfer term that in its most general form can be
written as

 dhm


¼ aKm hf  hm ;Gw ¼ 1  wf Cm
dt

ð1Þ

where Gw [T1] is the horizontal water flux (volumetric flux
per unit volume), wf [ ] is the volume fraction of the fracture
pore system, a [L2] is a first-order water transfer
Copyright 2004 by the American Geophysical Union.
0043-1397/04/2004WR003285$09.00

coefficient, Km [L T1] is the hydraulic conductivity of thematrix pore region, hf [L] is the pressure head of the fracture
pore region, hm [L] is the average pressure head of thematrix region, and Cm [L1] equal to dqm/dhm is the specificwater capacity of the matrix, where qm is the water content
in the matrix. A mathematical derivation of equation (1) was
given, for instance, by Zimmerman et al. [1993, 1996].
[3] A first-order transfer term such as equation (1) suffers
two drawbacks compared with the full description of flowusing, e.g., Richards’ equation. First, it is theoretically valid
only for quasi-steady water transfer [Zimmerman et al., 1993]
and therefore underestimates highly transient transfer flux forearly times [e.g., Zimmerman et al., 1993; Gerke and van
Genuchten, 1993a]. Second, in a lumped formulation such as
equation (1), Km needs to be evaluated as some blockaveraged, effective hydraulic conductivity function of thepressure head in the matrix and fractures that best represents

W074091 of 14

KÖHNE ET AL.: EVALUATION OF A SECOND-ORDER WATER TRANSFER TERM

W07409the conductivity of the actual pressure head profile at all
times. Gerke and van Genuchten [1993a] evaluated the
conductivity in equation (1) as an arithmetic average functionof pressure heads in both matrix and fractures and scaled thefirst-order transfer term to match the exact solution at 50%of the cumulative horizontal infiltration. Moreover, they
assumed an effective hydraulic conductivity at the fracturematrix interface, Ka [L T1], affected by fracture coatings[Gerke and van Genuchten, 1993a], to obtain following
modified first-order transfer term:

 dhm


¼ aGvG Ka hf  hm
Gw ¼ 1  wf Cm
dt
aGvG ¼ gw

b
a2

  

Ka ¼ 0:5 Ka hf þ Ka ðhm Þ ;

ð2aÞ
ð2bÞ
ð2cÞ

where the transfer coefficient aGvG includes the geometrydependent coefficient b [ ], a dimensionless scaling factor
gw , and the diffusion length a [L], representing the radius
for spheres (or cylinders) or the half width for cubes (orrectangular plates). On average, a scaling factor gw of 0.4
was deemed to give best results. However, best fit scaling
factors differed by up to 5 times for different initial pressure
head values between 30 and 3000 cm and varied by
20% or less between medium and fine textured soils[Gerke and van Genuchten, 1993a, Tables 2 and 3]. Evenfor the best fit scaling factors, the relative cumulative
horizontal infiltration was initially underestimated and later
overestimated, revealing the inherent limitations of the firstorder approach. A transfer term similar to equation (2) was
subsequently used by Gwo et al. [1995] and Ray et al.
[1997], among others.
[4] An accurate coupling term for a saturated dualpermeability system was suggested by Dykhuizen [1990].
The approach consisted of two different equations applicable to early- and late-time lateral imbibition. A nonlinear
ordinary differential equation accurate for all times waspresented by Zimmerman et al. [1993]. Their second-orderwater transfer term, derived by differentiating an analyticalsolution adapted from Vermeulen [1953], accurately approximates the exact solution for the pressure response of aspherical matrix block to a step function increase in thepressure head at its outer boundary [Zimmerman et al.,1993]. The term was subsequently modified for variably
saturated conditions [Zimmerman et al., 1996] and in thenotation of this paper is given by

2 
2
dhm bKm hm  hi  hf  hi


;
¼ 2
Gw ¼ 1  wf Cm
dt2ahm  hið3Þ

where t is time [T] and hi [L] is the initial pressure head, i.e.,a value that is assumed to be the same in the fractures andmatrix being initially in equilibrium. The second-order term(3) closely approximated reference solutions of horizontal
infiltration into saturated porous spheres and cubes[Zimmerman et al., 1993]. For unsaturated horizontal
infiltration into a 200-cm-thick rectangular rock slab, results
obtained with the dual-permeability transport of unsaturated
groundwater and heat (TOUGH) model with the secondorder mass transfer term were also in excellent agreement

W07409with the numerical reference simulation. However, forstructured soils, matrix aggregate sizes are much smaller
than 200 cm while matrix hydraulic conductivities are
generally much larger as compared with rock systems. It
was also shown that diffusion for fundamentally different
geometries like spheres and cylindrical macropores cannot
be described with a single, geometry-based model [vanGenuchten and Dalton, 1986; Young and Ball, 1997].
Hence it is not clear how the second-order term applies to an
aggregated soil and to the geometry of a cylindricalmacropore system.
[5] In this paper we will evaluate equation (3) with threemodifications for flow in macroporous soil systems. The firstmodification is that for variably saturated water flow, Km will
be evaluated as a block-averaged, effective function henceforth denoted K m . No hydraulic resistance at the fracturematrix interface will be assumed for K m as opposed to Ka inequation (2). Second, to allow for variable directions oftransfer between fracture and matrix with nonmonotonic
variation of hf, equation (3) is modified in a way that ensureshm ), as wasthat dhm /dt always has the same sign as (hf  shown and explained in detail by Zimmerman et al. [1993,
Appendix B] for the saturated system. Third, the initialpressure head in equation (3) is considered to be a dynamic
variable, hdi . Every time the transfer direction changes, hdi isreset to the pressure head value at the crossover point of
matrix and fracture pressure heads, i.e., the equilibriumpressure head at the time of transfer reversal. With thesemodifications, equation (3) changes to
 




 dhm bK m hf  hm hm  hdi  þ hf  hdi 


:
¼ 2
Gw ¼ 1  wf Cmhm  hd 
dt
2a
i
ð4Þ

[6] The primary objective of this study is to evaluate andenhance the accuracy of a second-order transfer term to
improve predictions of dual-permeability models for variably saturated preferential flow in structured soils. The
second-order mass transfer term (4) is evaluated in two
steps. In the first step, water transfer predictions obtainedwith the first- and second-order transfer terms are comparedwith a reference numerical solution of a horizontal flow
equation. Specific tasks with regard to the evaluation andimprovement of the second-order term are (1) to analyze theeffects of hydraulic properties, initial and boundary conditions, matrix size, and geometry, (2) to find the best single
estimation scheme for K m, and (3) to improve predictions of
the second-order term for variable transfer directions. In thesecond step, the performance of one-dimensional dualpermeability models with first- and second-order transfer
terms is assessed by comparing their results with numericalreference solutions of the two-dimensional Richards’ equation assuming different flow scenarios.

2. Comparison of the First- and Second-Order
Mass Transfer Terms With Horizontal
Flow Equation
2.1. Procedure[7] In numerical dual-permeability models, water transfer
between the matrix and fractures is assumed to be perpen-

2 of 14

W07409

KÖHNE ET AL.: EVALUATION OF A SECOND-ORDER WATER TRANSFER TERM

W07409[ 8 ] For numerical solution, the first-order transfer
term (2) and the second-order transfer term (4) are expressedin the following discretized forms (7) and (8) without theterm (1  wf ) so that they can be directly compared withequation (5).

Figure 1. Schematic representation of an elementary unitof a dual-permeability system consisting of parallel
rectangular matrix blocks of half-width a separated byfractures of half width b. White lines represent fictitious
pressure heads.

dicular to the direction of flow and has to be calculated foreach discretization element. Since preferential flow in soils
typically tends to be vertical, water transfer is essentiallyhorizontal. Therefore, as a starting point, water transfer
calculated with the first- and second-order mass transfer
terms will be compared with reference results obtained withthe following equation describing the horizontal water flow
in variably saturated rectangular matrix slabs:@qm ðx; hm Þ@
@hm ðx; t Þ¼
Km ðx; hm Þ
0  x  a;
@t
@x
@x

ð5Þ

where x [L] is the local horizontal coordinate between thefracture boundary at x = 0 and the center of a slab at x = a.
Figure 1 shows a schematic representation of a corresponding dual-permeability medium. For the purpose of comparing transfer terms and the horizontal flow equation, thevertical flow direction indicated by the vertical axis in
Figure 1 is neglected. A discretized form of equation (5)was numerically solved using HYDRUS-1D [Simunek etal., 1998] with a very small element size of 0.05 cm. On thebasis of results for qm in equation (5), effective values of thematrix pressure heads h were calculated from the average
saturation of the matrix block Se , using the inverse of van
Genuchten’s [1980] relation:
h n
i1n
h ¼ a1 S e1n  1
Za
Se ¼

qm  qm;rdx;
qm;s  qm;r
Dhm 0:4bK m ¼ 2
hf  hm
Dt
a Cm

ð7Þ

 



Dhm
bK m hf  hm hm  hdi  þ hf  hdi 


¼ 2
:
hm  hd Dt
2a C mi

ð8Þ

The nonlinear functions C m and K m were evaluated using
the van Genuchten [1980] equations of the soil hydraulicproperties. In equation (7), K m is evaluated according toequation (2c) with Ka replaced by Km. The first- andsecond-order mass transfer equations (7) and (8), respectively, were solved for hm iteratively at each time step using
Matlab (Mathworks, Inc.). Multiplying the solutions ofequations (7) and (8) in terms of Dhm /Dt with C m provided
the transfer flux Gw , which was compared to the solution ofequation (5).
[9] To evaluate the second-order term, different textures,
initial conditions, matrix sizes, and boundary conditions
were assumed, initially limiting the analysis to rectangularslab-type geometry of matrix blocks. Geometry was represented by a factor b in equations (7) and (8). Values for b inequation (7) could be derived, for example, by comparing
the solution of the first-order transfer term, assuming acertain geometry and a step-function boundary condition,with a Laplace transform of the linearized horizontal waterflow equation (e.g., b equals 3 for slab-type matrix blocks)[Gerke and van Genuchten, 1993a]. By analogy, the sameb values were used in equation (8).
[10] For the standard scenario, rectangular slab-type(b = 3) matrix blocks with a half width a of 5 cm wereconsidered. Two values of hi, equal to 100 and 1000 cm,
were selected for initial conditions. To represent the matrix
soil hydraulic properties, van Genuchten [1980] parameters
representing textures ‘‘silty clay,’’ ‘‘silt,’’ and ‘‘loamy sand’’
were chosen from a database of Carsel and Parrish [1988](Table 1). The following expressions for calculating theeffective hydraulic conductivity, K m , in equation (8) wereused:Matrix

K m ¼ Km ðhm Þ

ð9aÞFracture 
K m ¼ Km hf

ð9bÞ

ð6aÞ
Arithmetic 
K m ¼ Km ðhm Þ þ Km hf =2

ð9cÞ

Geometric 0:5
K m ¼ Km ðhm ÞKm hf

ð9dÞ

ð6bÞ

0

where qm,r and qm,s are the local residual and saturated watercontents of the matrix, respectively, and a and n areempirical constants for the matrix (here without index m forbetter legibility of equation (6a)).

Integral3 of 141Km ¼
hf  hm

Zhf
Km ðhÞdh
hm

ð9eÞ

KÖHNE ET AL.: EVALUATION OF A SECOND-ORDER WATER TRANSFER TERM

W07409

Table 1. Van Genuchten [1980] Parameters Used for
Horizontal, Two-Dimensional, and Vertical Dual-Permeability
Flow Simulationsa

Silty clay
Silt
Sandy loam
Sand

qr

qs

a, 1/cmn

Ks, cm/d

0.070
0.340
0.065
0.045

0.36
0.46
0.41
0.36

0.005
0.016
0.075
0.145

1.09
1.37
1.89
2.68

0.48
6.0
106.1
712.8

a

Parameters from Carsel and Parrish [1988].

When the van Genuchten [1980] model is used to calculateKm, the integral approach (9e) requires numerical integration. Hence, for each time step within the numerical solutionof equation (8), equation (9e) was itself numerically solved
utilizing an integral, adaptive Simpson quadrature as
implemented in the Quad routine of Matlab (Mathworks,
Inc.). Note that for consistency, C m was always evaluatedwith the same respective scheme as K m . Evaluation of C m
affects the prediction of hm but not of transfer flux, Gw .
[11] Statistical analysis was performed for results ofdifferent approaches subject to the step-increase boundary
condition in order to select the most appropriate approachfor further analyses using transient boundary condition.
The root-mean-square error (RMSE), quantifying deviations between observations Oi (horizontal flow equation
solution) and estimates Ei (mass transfer term solution),
was calculated and normalized with the arithmetic meanof observations, O, to yield a normalized RMSE orcoefficient of variation (CV) [Green and Stephenson,
1986]:100
CV ½% ¼
OPni¼1 ðOi  Ei Þ

2

!1=2

n

¼

100
RMSE
O

ð10Þ

2.2. Results for the Fracture-to-Matrix Transfer Usinga Step-Increase Boundary Condition[12] The following initial condition and step-increaseboundary conditions were assumed for the horizontal flow
equation (5):hm ðx; t ¼ 0Þ ¼ hi

ð11aÞhm ðx ¼ 0; t > 0Þ ¼ h0

ð11bÞqðx ¼ a; tð11cÞ0Þ ¼ 0;conceptual boundary condition. Practically, equation (12b)means that hf in equations (7) and (8) was kept constantwhen solving the initial value problem.
[13] Simulated water transfer at the initial pressure headof 1000 cm is shown in Figure 2. Note the decrease intimescale and increase in total water transfer from silty clay,
silt, to sandy loam. Using the matrix approach (8a), thesecond-order model significantly underestimated cumulative transfer for all textures (Figure 2). The estimation
somewhat improved with finer texture, characterized bysmaller hydraulic conductivity. This may explain why
Zimmerman et al. [1996] found an excellent agreement
between the second-order model and the numerical reference simulation of unsaturated horizontal imbibition into arock slab. Zimmerman et al. [1996] considered the saturated
hydraulic conductivity of rock matrix 4 orders of magnitude
smaller than Km of our silty clay.
[14] The second-order term with the fracture approach (9b)overestimated the water transfer (Figure 2) with CV values
smaller than those for the matrix approach (Table 2). Only
minor improvements were gained with the arithmetic approach (9c) that provided predictions relatively close to thefracture approach (Figure 2). The geometric approach (9d)revealed the strongest sensitivity to a change in texture(Figure 2). The best overall approximation, relatively
independent of texture, was achieved with the integral
approach (9e) (Figure 2, Table 2). The scaled first-ordertransfer model (2) always overestimated water transfer(Figure 2). Only for sandy loam (Figure 2c) was theagreement relatively good with a comparably low CV value(Table 2). No approach satisfactorily described the reference
solution for all three textures.
[15] Figure 3 shows results for the initial pressure head of
100 cm. The reference curve was somewhat better matchedby all approaches and particularly by the first-order term (2).
The CV values of most approaches were smaller for the
initial pressure head of 100 cm than for 1000 cm, whichshows that the variably saturated transfer term predictions
improve for smaller pressure head gradients between thedomains. Still better results were found for the initialpressure head of 10 cm (not shown).
[16] Using scaling factors as described by Gerke and van
Genuchten [1993a] to match the dimensionless cumulativetransfer at 0.5 with the second-order term yielded results
only slightly better than the first-order term (2) for all threesoils and both initial conditions of 100 and 1000 cm (notshown). Therefore the scaling approach was not consideredfurther. Instead, in order to better match both early- and latetime water transfers, we tested the following weighted
arithmetic averaging scheme for evaluating K m :Weighted arithmetic

where equation (11c) implies zero flow across the midplaneof rectangular slabs. For water transfer terms (7) and (8), thecorresponding initial and boundary conditions werehm ðt ¼ 0Þ ¼ hf ðt ¼ 0Þ ¼ hi

ð12aÞhf ðt > 0Þ ¼ h0 :

ð12bÞSince the water transfer terms are spatially nondimensional,
ordinary differential equations, equation (12b) is only a

W07409 
K m ¼ pKm ðhm Þ þ Km hf =ð p þ 1Þ; ð13Þ

where p is a weighting factor calibrated to achieve the best
fit to the reference model (5). Excellent agreement between
the results obtained with the second-order term using theweighted arithmetic scheme (13) and the cumulativereference water transfer (Table 2) could be achieved during
the entire time range for both initial conditions of 1000
and 100 cm and all textures (Figures 4 and 5). Note thatno considerable improvement over the scaling approach of
Gerke and van Genuchten [1993a] could be obtained by

4 of 14

KÖHNE ET AL.: EVALUATION OF A SECOND-ORDER WATER TRANSFER TERM

W07409W07409Figure 2. Cumulative water transfer into silty clay, silt, and sandy loam matrix slabs of 5-cm half width
at the initial pressure head of 1000 cm for a step increase of the boundary pressure head to 0 cm.
Simulation results are for the horizontal flow equation (5) (denoted ‘‘reference’’ in the legend), the firstorder-term (7) (denoted ‘‘first-order’’), and the second-order term (8) with different evaluation schemes
for K m (matrix, fracture, arithmetic, geometric, and integral, corresponding to equations (9a), (9b), (9c),(9d), and (9e), respectively).
implementing the weighted arithmetic scheme (13) into thefirst-order transfer term (7) (results not shown).
[17] For the second-order term, the weighting factor pvaried between 8 (sandy loam) and 59 (silty clay) (Table 2).
This variation of p with texture can be explained to somedegree by comparing the progression of the horizontal distributions of pressure heads and hydraulic conductivities
with time. For the initial pressure head of 1000 cm,
Figure 6 shows horizontal profiles of pressure heads andrelative hydraulic conductivity calculated with the reference
approach (5) for different fractions of the maximum time tmax,
defined as the time when simulated pressure heads exceeded0.005 cm in the entire profile. Pressure head increase
became more abrupt in the textural order of silty clay, silt,
and sandy loam, leading to sharp hydraulic gradients in the
sandy loam. For similar positions of the pressure front, thehydraulic conductivity dropped faster from its maximum
value at the (fracture) boundary for silty clay than for silt
and sandy loam. Hence, for horizontal infiltration into a finetextured porous medium assuming a step-increase boundary condition, ‘‘effective’’ hydraulic conductivity can be

expected to depend more on pressure heads inside the domain
than at the (constant head) fracture boundary. This phenomenon is particularly significant for fine textured soils and lessso for coarse textured soils, which is reflected by decreasing pvalues in the order silty clay, silt, and sandy loam (Table 2).
2.3. Results for the Fracture-to-Matrix Transfer Usinga Variable-Head Boundary Condition[18] The step-type boundary conditions (11b) and (12b)represented only an approximation of the real process, sincepressure heads in capillary fractures do not change in steps.
In our subsequent simulations, the imposed pressure head at
the boundary, h0, in the boundary conditions (11b) and(12b) was assumed to be a time-dependent function, h0(t).
The values for h0(t) were generated using the HYDRUS-1D
simulations of infiltration at the zero pressure head upperboundary condition into a vertical, 80-cm-deep highly
permeable porous medium, conceptually representing a‘‘fracture system,’’ being initially at either 1000 or
100 cm pressure head. A 30-min time series of thesimulated pressure heads at 40-cm depth was taken as

Table 2. Coefficients of Variation Between Numerical Solutions of Horizontal Flow Equation (5) and the Second-Order Term (7) With
Different Estimation Methods for K m as Given by Equations (9a) – (9e) and (12), and Between Equation (5) and the First-Order Term (6)Coefficient of Variation, %MatrixFracture

Arithmetic

Silty clay
Silt
Sandy loam

1899
5238
533300

84.4
77.0
71.5

73.4
71.2
57.8

Silty clay
Silt
Sandy loam

20
85
8460

46.9
63.0
64.7

44.5
49.5
81.7

Integral

Weighted
Arithmeticp[]FirstOrderhi = 1000 cm
42.3
19.3
2948.0

16.9
22.0
24.9

3.1
4.1
3.7

59
29
8

45.5
44.4
19.8

hi = 100 cm
30.5
35.1
58.5

15.9
27.0
26.5

1.0
2.4
3.6

45
19
8

11.5
6.5
6.4

Geometric

5 of 14

W07409

KÖHNE ET AL.: EVALUATION OF A SECOND-ORDER WATER TRANSFER TERM

W07409Figure 3. Same as Figure 2, but for the initial pressure head of 100 cm.
h0(t) in equations (11b) and (12b) for water transfer into thematrix blocks. This procedure gives values for h0(t) similar
to those that would be obtained for the fracture region of adual-permeability model with zero or negligible transfer
into the matrix. In this way, the h0(t) values represent themost abrupt, ‘‘worst-case’’ increase of fracture pressure
heads to be expected under natural conditions. Figure 7
shows the h0(t) boundary condition values rising from1000 to 0 cm in less than 10 min and the matrix pressure
head response calculated with the reference model (5)(reference-h). The figure also shows the relative watertransfer estimated with equations (5) (reference-q), (7),and (8) in terms of the cumulative transfer normalized withthe total transfer necessary to obtain full saturation.
[19] The reference matrix pressure heads rose faster in silt
than in silty clay (Figure 7) due to faster water transfer into
the more conductive silt. The first-order approach (7) did
not match well with the relative water transfer for either
texture (Figure 7). Similarly, the second-order term with thematrix, fracture, geometric, and arithmetic approaches failed
to describe the reference curves, while the integral approachgave a good approximation for silty clay but not for silt (notshown). On the average, the second-order term with theweighted arithmetic scheme and p = 17 most accurately
described water transfer into silty clay and silt (Figure 7).
For the initial pressure head of 100 cm, best agreement
was also found for p = 17 (Figure 8). Note that the h0(t)values increase rather suddenly close to saturation between
3 and 0 cm, since no matrix transfer was considered in thederivation of h0(t). This sudden increase causes deviationsof the first-order term results from the reference curve forsilty clay but not so much for silt (Figure 8) since between
3 and 0 cm, hydraulic conductivity increases tenfold forsilty clay but only twofold for silt. Only minor deviations of
the second-order term results were observed (Figure 8).
[20] The value of the empirical weighting factor p most
likely cannot be derived theoretically. However, for different textures and initial conditions, values of p around 17
were obtained for step-increase boundary condition in thereference model, provided that the second-order term wassolved with a variable-head boundary condition employing
the simulated pressure head values in the reference model

Figure 4. Same as Figure 2, but only simulation results for the horizontal flow equation (5)(‘‘reference’’), and the second-order term using K m -evaluation scheme (13) (‘‘weighted arithmetic’’).
6 of 14

W07409

KÖHNE ET AL.: EVALUATION OF A SECOND-ORDER WATER TRANSFER TERM

W07409Figure 5. Same as Figure 4, but for the initial pressure head of 100 cm.
domain at the first numerical node just inside of the(fracture) boundary (results not shown). Hence a continuous
instead of sudden pressure head increase at the boundary
strongly reduces the dependence of p on texture and initialcondition.
2.4. Results for Matrix-to-Fracture Transfer[21] To date, studies on soil water transfer have focused
only on fracture-to-matrix transfer. However, during drain-

age and evaporation periods, matrix-to-fracture transfer may
be more important. We analyzed water transfer, or horizontal drainage, from an initially saturated matrix block subject
to a variable-head boundary condition at the fracture face.
On the basis of the fracture continuum concept [e.g.,
Dykhuizen, 1990; Zimmerman et al., 1993], our approachassumes continuous water retention and hydraulic conductivity functions in a fracture continuum sustaining negative
pressures during matrix-to-fracture transfer. Such a fracture

Figure 6. Simulation results for the horizontal flow equation (5) of pressure heads and relative
hydraulic conductivities for cumulative water transfer into silty clay, silt, and sandy loam matrix slabs of
5-cm half width at the initial pressure head of 1000 cm for a step increase of the boundary pressure headto 0 cm. Graphs represent results at 1/16, 1/8, 1/4, 1/2, and 1 of the respective simulation time tmax.
7 of 14

W07409

KÖHNE ET AL.: EVALUATION OF A SECOND-ORDER WATER TRANSFER TERM

W07409Figure 7. Cumulative water transfer into silty clay and silt matrix slabs of 5-cm half width with the initialpressure head of 1000 cm and increasing boundary pressure heads: the boundary pressure head (BC)and the average matrix pressure head calculated using the horizontal flow equation (5) (‘‘reference-h’’);and water transfer calculated using the horizontal flow equation (5) (‘‘reference-q’’), the firstorder-term (7) (‘‘first-order’’), and the second-order term (8) using the weighted arithmetic scheme (13).
continuum may physically represent a fracture network or asoil fracture pore system (or macropore) partially refilledwith loose soil materials. We once again assumed rectangular slab geometry and 5-cm half width. The h0(t) values
for boundary conditions (11b) and (12b) were generated bysimulating 10 days of drainage at a constant evaporation
rate of 0.3 cm/d from a vertical, 80-cm-deep, initiallysaturated fracture system using HYDRUS-1D. Simulatedpressure heads at 2-cm depth below the evaporation boundary were taken as h0(t) in equations (11b) and (12b). Again,h0(t) values represent the fastest likely reaction of fracture
pressure heads under conditions of drainage and evaporation. Figure 9 shows h0(t) values applied as a boundary
condition for water transfer out of initially saturated matrix
slabs, the average matrix pressure head response calculatedwith the reference model (5) using equation (6), and thecumulative matrix-to-fracture transfer as calculated withequations (5), (7), and (8) and normalized with the reference
solution at 10 days. Average matrix pressure heads wererelatively close to the fracture pressure heads, especially for

Figure 8. Same as Figure 7, but with the initial pressure head of 100 cm.
8 of 14

W07409

KÖHNE ET AL.: EVALUATION OF A SECOND-ORDER WATER TRANSFER TERM

W07409Figure 9. Matrix-to-fracture water transfer for initially saturated silty clay and silt matrix slabs of 5-cm
half width: the boundary pressure head (BC) and the average matrix pressure heads calculated using thehorizontal flow equation (5) (‘‘reference-h’’); and matrix-to-fracture water transfer calculated using thehorizontal flow equation (5) (‘‘reference-q’’), the second-order term (8) using K m -evaluation schemes (9a)(‘‘second-order matrix’’) and (9b) (‘‘second-order fracture’’), and the first-order term (1) (‘‘first-order (1)’’)and (7) using K m -evaluation schemes (9a) (‘‘first-order matrix’’) and (9b) (‘‘first-order fracture’’).
silt (Figure 9). During matrix-to-fracture transfer, there
was only mild hydraulic nonequilibrium as compared withfracture-to-matrix transfer, where differences between
matrix and fracture pressure heads were much larger (seeFigure 8). For both first- and second-order terms, only
results for matrix (9a) and fracture (9b) based Km evaluation
schemes are shown, as all other approaches yielded results
in between these two. The second-order term with thematrix and fracture schemes, as well as the original firstorder term (1) without scaling coefficient, gave excellent
approximations of the reference solution of equation (5) for
both silty clay and silt (Figure 9). The first-order term (7)with gw equal to 0.4 [Gerke and van Genuchten, 1993a]underestimated matrix-to-fracture transfer (Figure 9). Thismeans that the scaling factor gw only improves the firstorder term description of fracture-to-matrix transfer, while
it is not applicable to matrix-to-fracture transfer. The direction-dependent use of gw is similar to a pseudohysteresis in
the first-order model.
2.5. Results for Reversible Water Transfer Between the
Fracture and the Matrix Domain[22] First- and second-order transfer terms were evaluatedfor ‘‘field-like’’ conditions with changing water transfer(matrix-to-fracture or fracture-to-matrix) directions.
The boundary pressure heads were again calculated using
HYDRUS-1D. The scenario (Figure 10) included 30-min
water transfer into the silty clay matrix slab initially at100 cm, followed by matrix-to-fracture transfer for
14400 min (10 days), and later water transfer back into
the matrix for the last 60 min. The prediction using the firstorder term underestimated, while the second-order termwith the weighted arithmetic scheme ( p = 17) matched
the first fracture-to-matrix transfer (Figure 10a). Althoughboundary condition or fracture pressure heads decreased

after 30 min, they remained larger than the averagematrix heads for some time, and water transfer back into
the fracture was delayed (Figure 10a). For the second
fracture-to-matrix water transfer event that started after
10 days at the pressure head of 230 cm (Figure 10b),
the second-order term (8) with the weighted arithmetic
scheme (13) matched the reference solution poorly because
the initial pressure head hi was fixed at the original value of
100 cm. When hi was adjusted to the new initial value(230 cm), the result was in excellent agreement with thereference solution (Figure 10b). Hence the assumption
underlying the second-order term of pressure head equilibrium between matrix and fracture continuum at the time oftransfer reversal, while simplified with respect to the actual
matrix pressure head profile (not shown), evidently allows aclose approximation of the real process.
2.6. Effect of Matrix Block Size on
Water Transfer Results[23] We analyzed the sensitivity of the second-ordertransfer term with the weighted arithmetic scheme to thesize of rectangular slabs. Figure 11 shows results for watertransfer into silty clay matrix blocks with half widths of 1, 5,
and 10 cm. In general, the second-order transfer term
approach was more accurate than the first-order transfer
term approach. Results for the second-order term with theweighted arithmetic scheme (13) came close to the reference
solution for the smaller half widths, whereas agreement
deteriorated for large matrix blocks of 20-cm diameter(Figure 11). The poor agreement was caused mainly by
the boundary pressure heads increasing steeply from 3 cm
to zero. Changing the weighting factor p from 17 to 59 at
this intermediate step of pressure increase as discussed in
the previous section (the value of 59 being calibrated forwater transfer subject to a step-increase boundary condition)

9 of 14

W07409

KÖHNE ET AL.: EVALUATION OF A SECOND-ORDER WATER TRANSFER TERM

W07409improved weighted arithmetic simulation (Figure 11) but is
not a feasible modeling scheme. However, for most fieldsoils, the half width of aggregates will be smaller than
10 cm, justifying the use of a constant weighting factor in
the weighted arithmetic scheme.

3. Evaluation of the Second-Order Transfer Termin a Dual-Permeability Model
3.1. Procedure[24] The second-order transfer term with the weighted
arithmetic approach (13) with p = 17 was evaluated withregard to its effect on dual-permeability model simulations ofpreferential flow. The dual-permeability model of Gerke andvan Genuchten [1993b] as implemented into HYDRUS-1D[Simunek et al., 2003] uses two mixed-type Richards’
equations to describe water flow in the fractures (14a) andthe matrix (14b):


@qf
@hf
@
Gw
¼
þ Kf 
Kf
@t
@z
wf
@z

ð14aÞ



@qm
@
@hm
Gw
Km
:
¼
þ Km þ@z@t@z
1  wf

ð14bÞ

A positive mass transfer term Gw denotes the fracture-tomatrix transfer, while a negative Gw signifies the matrix-tofracture transfer. The system of two dual-permeability flow
equations for matrix and fracture was solved sequentially,first for the matrix and then for the fracture pore system,
using the method of finite differences. The volume fractionof the fracture pore system, wf, was calculated depending on
geometry as follows:
wf ¼

wf ¼

Figure 10. Cumulative water transfer for silty clay matrix
slabs of 5-cm half width for a scenario with variableboundary pressure heads: the boundary pressure head (BC)and the average matrix pressure head calculated using thehorizontal flow equation (5) (‘‘reference-h’’); and watertransfer calculated using the horizontal flow equation (5)(‘‘reference-q’’), the first-order term (7) (‘‘first-order’’), andthe second-order term (8) using evaluation scheme (13) with
p = 17 (‘‘weighted arithmetic’’ and ‘‘weighted arithmetic, hi
variable’’). In the hi variable scheme, the initial head inequation (8) was reset at 14,436 min to the simulatedpressure head at that time. (a) Entire experiment duration.
(b) Results for time interval from 14,420 to 14,500 min.

b
aþb

b2
ða þ bÞ2

parallel plates;

ð15aÞ

hollow cylindrical macropore;

ð15bÞ

where a [L] is the half width of the matrix and b [L] is thehalf width of the fracture or the radius of a cylindricalmacropore. For a dual-permeability medium consisting of
cylindrical macropores surrounded by soil matrix, thefollowing expression was adopted to calculate B [Gerkeand van Genuchten, 1996] in transfer terms (2) and (4):b¼

1
½0:19 lnð16V0 Þ2;

1 < V0 < 100

V0 ¼ ða þ bÞ=b;

ð16aÞ

ð16bÞ

where z0 [ ] is the outer-to-inner-radius ratio of the hollow
cylinder. The transfer terms (2) and (4) were implemented
into the dual-permeability model. As a reference model, weapplied the two-dimensional numerical model HYDRUS-2D[Simunek et al., 1999] to a transport domain with twovertical layers representing the matrix and fracture regions

10 of 14

W07409

KÖHNE ET AL.: EVALUATION OF A SECOND-ORDER WATER TRANSFER TERM

W07409Figure 11. Cumulative water transfer into silty clay matrix slabs with half width a = 1 cm, a = 5 cm,
and a = 10 cm for a scenario of variable fracture pressure heads (BC) calculated using the horizontal flow
equation (5) (‘‘reference’’), the first-order term (7) (‘‘first-order’’), and the second-order term (8) using
evaluation scheme (13) with p = 17 (‘‘weighted arithmetic’’) and with p changing from 17 to 59 at 10 min(‘‘weighted arithmetic, p(17, 59)’’).

(Figure 12). HYDRUS-2D solves the mixed-type Richards
equation in two dimensions using the Galerkin finite element
method. We modified HYDRUS-2D to calculate watertransfer across the interfaces between two vertical layers
using a procedure similar to that for fluxes through theboundary nodes with prescribed Dirichlet boundary conditions [Simunek et al., 1999]. Nodal fluxes are calculated
from a finite element matrix equation that is obtained bydiscretization of the Richards equation and assembled for all
elements on one side of the interface. The finite element
matrix equation for internal fluxes is calculated at the last
iteration after convergence for given time level had beenachieved. Fluxes calculated directly from the Richards
equation are much more accurate than fluxes calculatedusing Darcy’s law and have the same accuracy as the overallsolution.
[25] Five scenarios as defined in Table 3 were selectedassuming an 80-cm-deep vertically structured soil profile.
The fracture system was represented by ‘‘sand’’ with van
Genuchten parameters given in Table 1 and was assigned ahalf width b of 0.25 cm. A seepage face lower boundary
condition was assumed.
3.2. Results
3.2.1. Infiltration Into a Parallel-Slab-Type
Dual-Permeability Medium With the Silty Clay Matrix(Scenario 1)[26] Figure 13a shows the depth profiles of water transfer
Gw at 15 and 60 min during infiltration as simulated using
the dual-permeability model (DPM) with first- and secondorder mass transfer term and HYDRUS-2D assuming adual-permeability medium consisting of a silty clay matrix
with slab geometry adjacent to a fracture (see Figures 1
and 12). The Gw profiles at both 15 and 60 min were similar
for HYDRUS-2D and the second-order DPM, whereas thefirst-order DPM failed to describe high peak values of
Gw near the infiltration front (Figure 13a).
3.2.2. Effect of the Silt Matrix (Scenario 2)[27] The assumption that the matrix consisted of silt
instead of silty clay (in scenario 1) resulted into a significantly slower vertical movement of the moisture front, since
the larger hydraulic conductivity of silt caused larger watertransfer rates into the matrix (Figure 13b). However, watertransfer in both two-dimensional and dual-permeability
model simulations was reduced to zero about 20 cm above
the infiltration front (Figure 13b). In this region, complete
equilibration between pressure heads in fracture and matrix
was obtained due to large transfer rates. Both first- andsecond-order DPMs gave satisfactory approximations of the

Figure 12. Finite element grid used in the twodimensional reference simulations.

11 of 14

KÖHNE ET AL.: EVALUATION OF A SECOND-ORDER WATER TRANSFER TERM

W07409W07409

Table 3. Scenarios Used in HYDRUS-2D and the Dual-Permeability Model for Evaluating First- and Second-Order Water Transfer
TermsaScenario

Flow
Process

Geometry

Matrix
Texturebhi,c cm

a, cm

b

wf

Upper Boundary
Condition

1
2
3
4
5

infiltration
infiltration
infiltration
infiltration
drainagerectangular
rectangularhollow cylindricalhollow cylindrical
rectangularsilty clay
silt
silty clay
silty clay
silty clay

100
100
100
100
0

5
5
5
1
5

3
3
0.82
1.44
3

0.04760.04760.00230.04000.0476

h=0
h=0
h=0
h=0
q=0

a
In all scenarios, profile depth = 80 cm, fracture pore system = sand (see Table 1 for van Genuchten parameters), b = 0.25 cm, seepage face lower
boundary condition.
b
See Table 1 for van Genuchten parameters.
cInitial pressure head.

reference transfer rate, with the second-order DPM matching closer to the transfer peaks. The good approximation of
the first-order term is due to the relatively fast transfer
reaching quickly its ‘‘late-time’’ phase for which the firstorder term is accurate [Zimmerman et al., 1993].
3.2.3. Effect of a Cylindrical Macropore With theLarge Matrix Mantle (Scenario 3)[28] The cylindrical macropore occupied only 0.23% of
the dual-permeability medium (wf = 0.0023, Table 3) in this
scenario as compared with the first scenario where for thesame values of a (5 cm) and b (0.25 cm), the fracture
between plate-like matrix blocks occupied 4.76% of thedomain (wf = 0.476, Table 3). Figure 13c shows depth
profiles of simulated water transfer, Gw . In scenario 3 thefirst-order DPM gave a slightly closer representation of thereference transfer profiles than the second-order DPM. Notethat the first-order term is specifically suited for hollow
cylindrical systems. For instance, it described solute diffusion in a hollow cylindrical system better as compared withother geometries [van Genuchten and Dalton, 1986].
3.2.4. Effect of a Cylindrical Macropore With the
Small Matrix Mantle (Scenario 4)[29] A matrix mantle radius a of only 1 cm resulted in amacropore volume fraction wf of 0.04, thus being 17 times
larger than in scenario 3. Figure 13d shows depth profiles ofwater transfer Gw at 60 and 180 min as simulated using
HYDRUS-2D, and the first- and second-order DPMs.
Reference transfer profiles revealed faster infiltration inspite of larger transfer than was found for scenario 3 witha equal to 5 cm (compare Figures 13c and 13d). This canbe explained by the larger macroporosity wf for a equal to
1 cm. The second-order DPM approximated the sharply
peaked reference transfer profiles considerably better than
the first-order DPM.
3.2.5. Effect of Flow Direction: Drainage (Scenario 5)[30] Finally, a draining scenario is analyzed assuming asoil system similar to scenario 1. The initially saturated
profile was subject to a seepage face lower boundary
condition and a flux upper boundary condition representing
evaporation from the matrix at 0.3 cm/d. Figure 13e showsthat the fracture-to-matrix (positive) transfer took place after
60 min in the upper part of the soil profile where evaporation resulted in more negative pressure heads in the matrix
than in the fracture. Simultaneously, the matrix-to-fracturetransfer occurred in the lower part of the soil profile where
fast draining of the macropore caused pressure heads in themacropore to drop below the corresponding values in
the adjacent matrix. The overall transfer rates are anorder-of-magnitude lower than in the first scenario (compare Figures 13e and 13a). After 600 min, pressure heads
in two regions equilibrated near the soil surface, and thus
water transfer between the domains approached zero
(Figure 13f ). Both the first- and second-order DPMs
approximated the reference transfer profiles satisfactorily.
[31] Scenarios 1 and 2 demonstrated the significantly
higher accuracy of the second-order model for infiltration
into a dual-permeability medium of rectangular matrix slabs
separated by parallel fractures. As mentioned above, thesecond-order term was derived for spherical blocks and has
been demonstrated to be also applicable to relatively similargeometries as cubes [Zimmerman et al., 1993, 1996]. For
infiltration into a cylindrical system with large outer-toinner-radius ratio, z0 (scenario 3), the second-order term did
not give satisfactory results. It was documented in theliterature that the exact analytical solution for spherical
diffusion fails to describe diffusion in a cylindrical macropore system, even when the geometry is embedded infactor b [van Genuchten and Dalton, 1986]. The geometry
factor b accounts for the smaller surface-to-volume ratio of
the cylindrical macropore system as compared with spheres[Gerke and van Genuchten, 1996], but it cannot always
compensate for the fundamentally different flow field in
the cylindrical macropore system. At early times, thenormalized water uptake is proportional to the surface-tovolume ratio. The hollow cylinder has a small surface-tovolume ratio (reflected by a small value of b) particularly
for a large outer-to-inner-radius ratio, z0 [Gerke and van
Genuchten, 1996], and so it will initially have a muchsmaller water uptake than a sphere (cube, solid cylinder).
Only at late times will the flux be nearly the same forthe hollow cylinder and the sphere. However, the fourth
scenario indicated that the second-order term derived forspherical geometry can be applied to a cylindrical macropore system with small z0, i.e., increasing values forsurface-to-volume ratio and b.
[32] For soil systems with small macroporosity, a different
transfer equation may be needed. Unfortunately, analyticalsolutions for diffusion in hollow cylindrical systems [e.g.,van Genuchten et al., 1984] are mathematically rather
complicated, such that a procedure as that presented by
Vermeulen [1953] to derive a second-order expression for
diffusion into spheres may not be easily applicable for thehollow cylindrical system.
[33] The primary significance of the water transfer
calculation might be its effect on advective solute transfer
which in dual-permeability models is calculated as

12 of 14

W07409

KÖHNE ET AL.: EVALUATION OF A SECOND-ORDER WATER TRANSFER TERM

W07409Figure 13. Water transfer rate versus depth as calculated with the two-dimensional (HYDRUS-2D)
model, and the dual-permeability model with the first- (DPM first-order) and second-order (DPM secondorder) water transfer terms using evaluation scheme (13) with p = 17 (weighted arithmetic). Infiltration
into dual-permeability medium with (a) silty clay and (b) silt matrix slabs of a half width of 5 cm;
infiltration into cylindrical dual-permeability medium with silty clay matrix mantle of (c) 5 cm and(d) 1 cm; drainage and evaporation out of dual-permeability medium with silty clay matrix slabs of 5-cm
half width at (e) t = 60 min and (f ) t = 600 min.

13 of 14

W07409

KÖHNE ET AL.: EVALUATION OF A SECOND-ORDER WATER TRANSFER TERM

water transfer times concentration [e.g., Gerke and van
Genuchten, 1993b]. A detailed analysis of the impact of
advective solute transfer on the simulation of preferential
solute transport, however, is beyond the scope of thisstudy and will be addressed in a separate paper.

4. Concluding Remarks[34] A second-order term originally derived for watertransfer between fractures and matrix in fractured rock
was adapted for variably saturated structured soils. The
hydraulic conductivity in the second-order term was evaluated as a weighted average function of pressure heads inthe matrix and fracture. This modification significantly
improved agreement between water transfer calculated withthe second-order term and a numerical solution of a reference equation for horizontal flow. A single weighting factorof 17 proved successful to describe water transfer for arange of initial and boundary conditions, as well as fordifferent soil matrix block sizes and textures. Furthermore,
the initial pressure head in the second-order term wastreated as a dynamic parameter to be reset upon change offlow direction between matrix and fracture domains. With
this adjustment the second-order term does not show
pseudohysteresis for varying directions of water transfer,
as opposed to the second-order term with constant initialpressure head and the first-order term with scaling factor.
Implemented in a dual-permeability model, the modified
second-order transfer term improves the accuracy ofpreferential soil water flow simulations as compared witha first-order term, particularly for geometries other thanhollow cylinders. The modified second-order term is asversatile as first-order terms with regard to its application indual-permeability models for studying variably saturated
flow in field soils.
[35] Acknowledgment. This work was supported by the National
Science Foundation (grant EAR 0296158 and STC Program, EAR
9876800).

W07409Gerke, H. H., and M. T. van Genuchten (1996), Macroscopic representationof structural geometry for simulating water and solute mass transfer indual-porosity media, Adv. Water Resour., 19, 343 – 357.
Green, I. R. A., and D. Stephenson (1986), Criteria for comparison of single
event models, Hydrol. Sci., 31(3), 395 – 411.
Gwo, J. P., P. M. Jardine, G. V. Wilson, and G. T. Yeh (1995), A multiplepore-region concept to modeling mass transfer in subsurface media,
J. Hydrol., 164, 217 – 237.
Ray, C., T. R. Ellsworth, A. J. Valocchi, and C. W. Boast (1997), An
improved dual porosity model for chemical transport in macroporous
soils, J. Hydrol., 193, 270 – 292.
Simunek, J., M. Sejna, and M. T. van Genuchten (1998), The HYDRUS-1D
software package for simulating the one-dimensional movement of water,
heat, and multiple solutes in variably-saturated media, version 2.0, Rep.
IGWMC-TPS-70, 202 pp., Int. Ground Water Modell. Cent., Colo.
School of Mines, Golden, Colo.
Simunek, J., M. Sejna, and M. T. van Genuchten (1999), The HYDRUS-2D
software package for simulating two-dimensional movement of water,
heat, and multiple solutes in variably-saturated media, version 2.0,
Rep. IGWMC-TPS-53, 251 pp., Int. Ground Water Modell. Cent., Colo.
School of Mines, Golden, Colo.
Simunek, J., N. J. Jarvis, M. T. van Genuchten, and A. Gärdenäs (2003),
Review and comparison of models for describing nonequilibrium andpreferential flow and transport in the vadose zone, J. Hydrol., 272,
14 – 35.
van Genuchten, M. T. (1980), A closed-form equation for predicting thehydraulic conductivity of unsaturated soils, Soil Sci. Soc. Am. J., 44,892 – 898.
van Genuchten, M. T., and F. N. Dalton (1986), Models for simulating salt
movement in aggregated field soils, Geoderma, 38, 165 – 183.
van Genuchten, M. T., D. H. Tang, and R. Guennelon (1984), Some exact
solutions for solute transport through soils containing large cylindricalmacropores, Water Resour. Res., 20, 335 – 346.
Vermeulen, T. (1953), Theory of irreversible and constant-pattern solid
diffusion, Ind. Eng. Chem., 45, 1664 – 1670.
Warren, J. E., and P. J. Root (1963), The behavior of naturally fractured
reservoirs, Soc. Pet. Eng. J., 3, 245 – 255.
Young, D. F., and W. P. Ball (1997), Effects of column conditions on thefirst-order rate modeling of nonequilibrium solute breakthrough: Cylindrical macropores versus spherical media, Water Resour. Res., 33, 1149 –1156.
Zimmerman, R. W., G. Chen, T. Hadgu, and G. S. Bodvarsson (1993), A
numerical dual-porosity model with semianalytical treatment of fracture/matrix flow, Water Resour. Res., 29, 2127 – 2137.
Zimmerman, R. W., T. Hadgu, and G. S. Bodvarsson (1996), A new
lumped-parameter model for flow in unsaturated dual-porosity media,
Adv. Water Resour., 19(5), 317 – 327.

References
Carsel, R. F., and R. S. Parrish (1988), Developing joint probability
distributions of soil water retention characteristics, Water Resour. Res.,24, 755 – 769.
Dykhuizen, R. C. (1990), A new coupling term for dual-porosity models,
Water Resour. Res., 26, 351 – 356.
Gerke, H. H., and M. T. van Genuchten (1993a), Evaluation of a first-orderwater transfer term for variably saturated dual-porosity flow models,
Water Resour. Res., 29, 1225 – 1238.
Gerke, H. H., and M. T. van Genuchten (1993b), A dual-porosity model for
simulating the preferential movement of water and solutes in structured
porous media, Water Resour. Res., 29, 305 – 319.




H. H. Gerke, Center for Agricultural Landscape and Land Use Research(ZALF), Institute of Soil Landscape Research, Eberswalder Strasse 84,
D-15374 Müncheberg, Germany. (hgerke@zalf.de)J. M. Köhne and B. P. Mohanty, Department of Biological andAgricultural Engineering, Texas A&M University, 140 Scoates Hall,
College Station, TX 77843-2117, USA. (mkoehne@cora.tamu.edu;
bmohanty@tamu.edu)
J. Simunek, Department of Environmental Sciences, University ofCalifornia, Riverside, 900 University Avenue, A135 Bourns Hall, Riverside, CA 92521, USA. ( jiri.simunek@ucr.edu)

14 of 14