Soil & Water Res., 1, 2006 (1): 16–22

Numerical Study of Macropore Impact on Ponded
Infiltration in Clay Soils
Radka Kodešová1, Josef Kozák1 and Jiří Šimůnek2
1

Department of Soil Science and Geology, Faculty of Agrobiology, Food and Natural
Resources, Czech University of Agriculture in Prague, Prague, Czech Republic;2
Department of Environmental Sciences, University of California Riverside, CA, USA
Abstract: The single-porosity and dual-permeability models in HYDRUS-1D (Šimůnek et al. 1998, 2003) wereused to simulate variably-saturated water movement in clay soils with and without macropores. Numericalsimulations of water flow for several scenarios of probable macropore compositions show a considerable impact
of preferential flow on water infiltration in such soils. Preferential flow must be considered to predict water
recharge in clay soils.
Keywords: ponded infiltration; macropores and fractures; preferential flow; single-porosity system; dual-permeability system; HYDRUS-1D

Estimation of water infiltration in heavy soils
and their subsequent water saturation is very important for assessment of slope stability of newly
designed or existing reclaimed landfills. Usuallyapplied methods and numerical models often
significantly underestimate the total amount of
infiltrated water since they assume equilibrium
water flow, and use the Richards equation andsoil hydraulic properties determined for the entire
flow domain. It has been indicated in many studies(for instance in Valla et al. 2000) that in reality
the infiltration rate in such soils is much higher
due to the presence of macropores and fractures.
Models that allow consideration of non-equilibrium
water flow must be used to obtain more precise
estimates of amounts of infiltrated water and water
saturation of the soil profile.
Mathematical description of the equilibrium
water flow and solute transport assumes a single
continuum approach. To describe the non-equi-

librium water flow and solute transport in soils,
numerical models that assume two continuum
approach have been developed. Soil porous system is in such models divided into two domains.
Each domain is characterized with its own setof transport properties and equations describing transport processes. Dual-porosity approachdefining a water and solute transport in systems
consisting of domains of mobile and immobile
water was presented by Phillip (1968) and Van
Genuchten and Wierenga (1976). The dualporosity formulation for water flow and solutetransport is based on a set of partial differential equations describing water flow and solutetransport in the mobile domain and mass balance equations describing moisture dynamics andsolute content in the immobile domain. On theother side, dual-permeability approach assumes
water flow and solute transport in both domains.
The dual-permeability formulation for water flow

Supported by the Ministry of Education, Youth and Sports of the Czech Republic, Project No. MSM 6046070901,
by the Ministry of Agriculture of the Czech Republic, Project No. QF3250, and by the Czech Science Foundation,
Project No. GACR 103/05/2143.

16

Soil & Water Res., 1, 2006 (1): 16–22
and solute transport is based on a set of equations describing water flow and solute transport
separately in each domain (matrix and macropore
domains). Different equations may be used tosimulate water flow in the mobile and macropore
domains. A kinematic wave approach was used byGermann (1985), Germann and Beven (1985)
and Jarvis (1994) to describe flow in macropores.
The Richards equation was used by Gerke and
Van Genuchten (1993, 1996) to describe flow

in both matrix and macropore domains. Other
approaches can be based on Poiseuille equation(Ahuja & Hebson 1992), and Green-Ampt orPhilip infiltration equations (Ahuja & Hebson
1992; Chen & Wagenet 1992).
Single-porosity and dual-permeability models in
HYDRUS-1D (Šimůnek et al. 1998, 2003) basedon the numerical solution of the Richards equation in all domains, as was proposed by Gerke
and Van Genuchten (1993, 1996), were applied

Figure 1. Examples of macropores and fractures in clay soils

17Soil & Water Res., 1, 2006 (1): 16–22
in this study. The same approach was also used
for instance by Vogel et al. (2000). Different scenarios were designed to simulate one-hour ponded infiltration into heavy soils with and without
macropores. The aim of this study was to show
the effects of different volume, configuration andshape of macropores on the infiltration rate andsaturation of the soil profile.

Macropore origin is usually associated withtransformation of pores between and inside aggregates. While neglecting dynamic changes in thepore structure and assuming capillary behavior inboth matrix and macroporous domains, the dualpermeability model may describe the impact ofmacropores on water flow in the soil profile.
Concepts of single-porosityand dual-permeability models

Material and methods
Impact of macropores and fractures
Water flow in soils is traditionally described using
the Richards equation that solves the equilibrium
water flow in variably saturated rigid soil porous
media. Soil hydraulic properties are defined by soil
water content retention and unsaturated hydraulicconductivity curves for the whole flow domain.
Water regime in the soil profile is described byvalues of soil water contents and pressure heads,
that are defined as average values of water contents
or pressure heads in different pore domains (matrix
soil porous system and macroporous system). The
impact of macropores and fractures, which exist
due to soil particle aggregation, volume changes,
influence of living organisms, non-consolidation of soils and so on, on flow and transport isneglected. The impact of macropores is greaterin clay soils, where the hydraulic conductivity is
expected to be very low and where macropores
may transfer a bulk of water and solute, than inloams and sands. Water in clay soils penetrates
deep into the soil profile mostly due to preferential
flow in macropores and fractures. Shapes and sizesof aggregates, and consequently of macroporesand fractures, are different in different soils. They
depend on the mineralogical composition, stageof disintegration, soil water content, way of soil
deposition, consolidation and so on. Soil structure
in clay soils at four locations in Tušimice and Bílina
is documented on Figure 1:(a) angular blocky structure in consolidated moist
soil,(b) angular blocky structure with appearance of
cubic-shape segregates and fractures in dry
soil,(c) clay disintegration under extreme conditions
into the plate-like structure and following
reclamation,(d) fractures in non-consolidated soil or in soil
disturbed by plant roots.
18

Water flow in the soil profile is in this study
simulated using the single-porosity and dualpermeability models (Gerke & Van Genuchten
1993, 1996) implemented in HYDRUS-1D. While
the single-porosity model describes a uniform flow
in soil porous media, the dual-permeability model
can represent a non-equilibrium/preferential flow
in a complex system of soil aggregates (matrix) andintermediate spaces (fractures or macropores). The
Richards equation, describing the one-dimensional
isothermal Darcian flow in a variably saturated rigid
porous medium is used in both models.
The single Richards equation is used for the single-porosity system:
wΌ
w §
wh
·
		
 K h ¸¸  S¨¨ K h
wt

wz ©

wz

(1)

¹

where:
θ – volumetric soil water content (l3/l3)h – pressure head (L)K 	– soil water content (L/T)
S – sink term (T-1)
t – time (T)z 	 – vertical axes (L)Equation (1) is solved for the entire flow domainusing one set of soil water retention and hydraulicconductivity curves.
In the case of the dual-permeability model, theRichards equation is applied separately to each of
the two pore regions – macropore (fracture) andmatrix domains:
wΌ f
wt
wΌ m
wt

w §¨
Kf
wz ¨©

hf

w §
¨ Km hm
wz ¨©

wh f
wz
wh m
wz

hf

·
̆
¸ S  wf¸
fw¹

 Km hm̆w·¸¸  S m 
1 f w
¹

 Kf

(2)

where:
θf , θm – volumetric soil water contents in the macropore
		 and matrix pore domains (l3/l3), respectively

Soil & Water Res., 1, 2006 (1): 16–22
hf,, hm – pressure heads in both domains (L)Kf, , Km 	– hydraulic conductivities (L/T)
Sf , Sm – sink terms (T–1)
Γw– mass exchange between the matrix and macro		 pore regions (T–1)
fw– ratio between volume of the macropore domain
		 and the total flow domain (–)

Equations (2) are solved using two sets of soil
water retention and hydraulic conductivity curves
that are defined for each domain.
Analytical expressions proposed by Van Genuchten (1980) for the soil water content retention curve, θ(h), and the hydraulic conductivity
function, k(θ), are:
Ό ( h )  Όr
		
=
Όe =
Όs  Όrθe = 11
1 + |΅ h| n

m

h<0

(3)

h≥0

K(θ) = Ks θeλ [1 – (1 – θe1/m)m]2
K(θ) = Ks

h≥0

h<0

(4)

where:
θe– effective soil water content (–)
Ks– saturated hydraulic conductivity (L/T)θr , θs – residual and saturated soil water contents(L3/L3)
λ– pore-connectivity parameter (–) (λ = 0.5 in this		 case)α (L–1), n, m (= 1 – 1/n) – empirical parameters (–)

Parameters describing aggregate shapes are theshape factor β (–) (= 15 for spherical aggregates,= 3 for cubic aggregates), the characteristic lengthof aggregate a (L) (sphere radius or half size of the
cube edge), and the dimensionless scaling factor
γ w (–) (= 0.4). The mass exchange between thematrix and macropore regions, Гw, is calculatedusing following equation:
E
J h f  hm
*w
		
(5)
2 Ka wawhere:
Ka – the effective saturated hydraulic conductivity of the		 interface between the two pore domains (L/T)

Scenarios of one-hour ponded infiltration
in soils without and with macropores
The impact of macropores on water flow in thesoil profile is demonstrated on the simulation of

the one-hour ponded infiltration with the ponded
depth of 1 cm on the top of the soil profile. The
initial conditions were given by the pressure headof –200 cm in all domains. The bottom boundary
condition was defined as free drainage. Sink terms
representing the root water uptake were assumed
to be equal to zero. Parameters characterizing
hydraulic properties were set such to study the
impact of macropore from the following pointsof view:(a) shapes and sizes of aggregates,(b) quality of macropore (described by variable soil
water retention curves and saturated hydraulicconductivities), and(c) the ratio of the macropore domain to the wholepore domain.
The basic assumption was that the soil water
retention and hydraulic conductivity curves obtained for the mono-modal soil porous system also
define hydraulic properties of the matrix domainof the dual-permeability system. Van Genuchten
parameters were as follows: θ r = 0.09 cm 3 /cm 3 ,θ s = 0.408 cm 3/cm 3, α = 0.012 cm –1, n = 1.28 andK s = 0.0043 cm/min.
Soil hydraulic parameters of the macropore
domain were estimated assuming aggregate development and hydro-physical properties. Soil
porosities at observed locations varied between
40% and 55% and maximum capillary adsorptions
varied between 34% and 44%. Experimentally
determined values of the saturated hydraulicconductivity in clay soils may vary by two or moreorders of magnitude (Kočárek et al. 2005). Twolimiting cases of aggregate shapes were studied.
Shapes of aggregates were defined as either spherical or cubic. Spherical aggregates may correspondwith the angular blocky structure and cubic aggregates may represent the cubic segregates. Flowin fractures may be approximated by using thegreater characteristic length of cubic aggregates.
Ratios of the macropore domain to the total soilvolume were assumed to be 0.05, 0.1 and 0.15.
Parameters of the soil water retention curve were
set to ensure that the total saturated soil watercontent of the dual-permeability system washigher than for the single-porosity system andthat macropores were filled with water only forpressure heads close to zero (step-shaped soil
water retention curves). The saturated hydraulicconductivities of the macropore domain were
set to be 100 and 500 times greater than for thematrix domain.
19

Soil & Water Res., 1, 2006 (1): 16–22
Results and Discussion

Cumulative infiltration (cm)Cumulative infiltrations and infiltration rates
versus time for all scenarios (except for B3 andB4) are shown in Figures 3 and 4. Shapes of thewetting front in the soil profile 20, 40 and 60 min
after the beginning of the ponded infiltration for
scenarios A1, B1 and B2 are shown in Figure 5.
While in the soil profile without macropores theTime (min)Figure 3. Cumulative infiltration at the top of the soil
profileInfiltration rate (cm/min)

The following nine scenarios were studied. Scenario A1 represents a soil without macropores.
Scenarios B1, B2, B3 and B4 were designed to asses
the impact of the size and shape of aggregates.
The Van Genuchten parameters in all macropore
domains were defined as follows: θ r = 0 cm 3/cm 3,
θ s = 0.7 cm 3 /cm 3 , α = 0.07 cm –1 , n = 3 and K s =0.43 cm/min. The ratio of the macropore domain,f w , was 0.05, the shape factor, β, was 3 (B1 and B2)
and 15 (B3 and B4) and the characteristic length ofan aggregate, a, was 1 (B1 and B3) and 2.5 cm (B2and B4). Scenarios C1, C2 and C3 were designed
to assess macropore quality. The Van Genuchten
parameters in the macropore domain were define
as follows: θ r = 0 cm 3/cm 3, θ s = 0.7 cm 3/cm 3, α =0.07 cm –1, n = 3 and K s = 0.43 cm/min (C1), θ r = 0
cm3/cm3, θs = 0.9 cm3/cm3, α = 0.15 cm–1, n = 3 andK s = 0.43 cm/min (C2), and θr = 0, θs = 0.7, α = 0.07
cm –1, n = 3 and K s = 2.15 cm/min (C3). The ratioof the macropore domain, f w , was 0.1, the shapefactor, β, was 3 and the characteristic length of anaggregate, a, was 1 cm. Scenario D1 together withscenarios B1 and C1 were designed to assess theimpact of the ratio of the macropore domain. All
parameters in the D1 scenario were the same as in
scenarios B1 and C1, with the exception the ratioof the macropore domain, f w , that was equal to0.15. The effective saturated hydraulic conductivity, K a, was same for all scenarios and equaled to
0.0001 cm/min. The total soil hydraulic propertiesof the dual-permeability system are given by thesum of the properties for each domain multipliedby its fraction in the system. Soil water retention
curves, θ(h), for the single-porosity system (A1)and the dual-permeability systems (C1 or C3) are
shown in Figure 2. Hydraulic conductivity curves,
K(h), may be plotted similarly.

Time (min)Figure 4. Infiltration rate at the top of the soil profile

Single-porosity––

Dual-permeability
bothFigure 2. Soil water retention curves for single-porosity
model (scenario A1) and dual-permeability models (scenarios C1 and C3)

Depth (cm)Pressure head (cm)- - - Dual-permeability
macropores

Soil water content (cm3/cm3)

20

Soil water content (cm3/cm3)

Dual-permeability
matrix pores

A1 – 20 min
A1 – 40 min
A1 – 60 min
B1 – 20 min
B1 – 40 min
B1 – 60 min
B2 – 20 min
B2 – 40 min
B2 – 60 min

Figure 5. Soil water content distributions in the soil profile– wetting fronts 20, 40 and 60 minutes after the beginningof the ponded infiltration for scenarios A1, B1 and B2

Soil & Water Res., 1, 2006 (1): 16–22
Soil water content (cm3/cm3)

C1 – 20 min
C1 – 40 min
C1 – 60 min
C2 – 20 min
C2 – 40 min
C2 – 60 min
C3 – 20 min
C3 – 40 min
C3 – 60 min

B1 – 20 min
B1 – 40 min
B1 – 60 min
C1 – 20 min
C1 – 40 min
C1 – 60 min
D1 – 20 min
D1 – 40 min
D1 – 60 min

Depth (cm)Depth (cm)Soil water content (cm3/cm3)

Figure 6. Soil water content distributions in the soil profile– wetting fronts 20, 40 and 60 minutes after the beginningof the ponded infiltration for scenarios C1, C2 and C3

Figure 7. Soil water content distributions in the soil profile– wetting fronts 20, 40 and 60 minutes after the beginningof the ponded infiltration for scenarios B1, C1 and D1

wetting front after 60 min was only in the depthof 12 cm, in the case with macropores the wetting
front was in the depth of 30 and 43 cm for scenarios
B1 and B2, respectively. The total cumulative infiltration was 0.70 cm in the soil profile without
macropores, and 2.72 (B1) and 2.61 (B2) cm in thesoil profiles with macropores. The wetting front
in the scenario with bigger aggregates moved to
the greater depth than in the scenario with smaller
aggregates. However, the total cumulative infiltration was slightly lower due to the irregularity of
the wetting front. Simulated results correspondwith the probable type of water flow in soils withthe cubic segregates and fracture system. Shapesof the wetting front for scenarios B3 and B4 (notshown) exhibited the same trends as those for B1
and B2 scenarios. The wetting front moved only
to the depth of 20 and 32 cm for scenarios B3
and B4, respectively, due to round shape of soil
aggregates that represent the soil porous systemwith angular blocky structure. The total cumulative infiltration was similar to those in scenarios
B1 and B2, i.e., 2.71 and 2.71 cm for scenarios B3
and B4, respectively.
Shapes of the wetting front in the soil profile 20,
40 and 60 min after the beginning of the ponded
infiltration for scenarios C1, C2 and C3 are shownin Figure 6. A comparison of scenarios C1 and C2
shows that, in spite of the higher total saturated
water content in scenario C2, the wetting front
moved to smaller depths than in scenario C1 dueto the less gradual shape of the soil water retention
curve. The total cumulative infiltration (4.54 and4.12 cm for scenarios C1 and C2, respectively) was
also smaller. A comparison of scenarios C1 and C3

shows that the saturated hydraulic conductivity
has the greatest impact on water flow in the soil
profile. The soil profile for scenario C3 is at theend of the simulated infiltration almost completely
saturated. The total cumulative infiltration was
16.01 cm.
Shapes of the wetting front in the soil profile
20, 40 and 60 min after the beginning of the ponded infiltration for scenarios B1, C1 and D1 areshown in Figure 7. Higher ratios of the macropore
domain caused progress of the wetting front into
greater depths and increased total cumulative
infiltrations as follows: 2.72 (B1), 4.54 (C1) and6.31 (D1) cm.
The presented results show that water flow inclay soils may be predominantly affected by characteristics of macropores and fractures. Water istransported downwards mainly through macropores and fractures and then redistributed into
the matrix pores. As a result, at the end of theinfiltration event, the saturation of a soil withmacropores is higher than would be expected ina soil without macropores and fractures.
Conclusions
Macropores and fractures play an importantrole in water transport. The impact of macropores is greater in clay soils where the hydraulicconductivity is expected to be very low. Resultsof numerical simulations showed that the cumulative infiltration into the soil with and without
macropores may differ by two or more orders ofmagnitude. The appropriate models must thenbe used to describe non-equilibrium flow in such
21

Soil & Water Res., 1, 2006 (1): 16–22
soils. Since the water in dual-permeability systems is predominantly influenced by propertiesof macropores and fractures, suitable methods
must be used to determine characteristics of themacropore domain.
Acknowledgement. Thanks are due to L. Borůvka,K. Němeček, M. Rohošková and V. Penízek for their
help with the field and laboratory work and help withphotographic documentation.

References
Ahuja L.R., Hebson C. (1992): Root Zone Water Quality
Model. GPSR Technical Report No. 2, USDA, ARS,
Fort Collins, CO.
Chen C., Wagenet R.J. (1992): Simulation of water andchemicals in macropore soils. Part 1. Representation of
the equivalent macropore influence and its effect on soil
water flow. Journal of Hydrology, 130: 105–126.
Gerke H.H., Van Genuchten M.Th. (1993): A dual-porosity model for simulating the preferential movementof water and solutes in structured porous media. Water
Resources Research, 29: 305–319.
Gerke H.H., Van Genuchten M.Th. (1996): Macroscopic
representation of structural geometry for simulatingwater and solute movement in dual-porosity media.
Advances in Water Resources, 19: 343–357.
Germann P.F. (1985): Kinematic wave approach to infiltration and drainage into and from soil macropores.
Transactions of the ASAE, 28: 745–749.
Germann P.F., Beven K. (1985): Kinematic wave approximation to infiltration into soils with sorbing macropores.
Water Resources Research, 21 (7): 990–996.
Jarvis J. (1994): The MACRO model. Technical description and sample simulation. Reports and dissertations
19. Department of Soil Science, Swedish University ofAgricultural Science, Uppsala, Sweden, 51.

Kočárek M., Kodešová R., Kozák J., Drábek O., Vacek
O. (2005): Chlorotoluron behaviour in five varying soil
types. Plant, Soil and Environment, 51: 304–309.
Phillip J.R. (1968): The theory of adsorption in aggregatedmedia. Australian Journal of Soil Research, 6: 1–19.
Šimůnek J., Šejna M., Van Genuchten M.Th. (1998):
The HYDRUS-1D Software Package for Simulating theOne-dimensional Movement of Water, Heat and Multiple Solutes in Variably-saturated Media. Version 2.0.
IGWMC-TPS-53. International Ground Water ModelingCenter, Colorado School of Mines, Golden, CO.
Šimůnek J., Jarvis N.J., Van Genuchten M.Th., Gärdenäs A. (2003): Review and comparison of models
for describing non-equilibrium and preferential flow
and transport in the vadose zone. Journal of Hydrology,
272: 14–35.
Valla M., Kozák J., Ondráček V. (2000): Vulnerability of aggregates separated from selected anthrosols
developed on reclaimed dumpsites. Rostlinná Výroba,46 (12): 563–568.
Van Genuchten M.Th. (1980): A closed-form equation
for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal, 44:
892–898.
Van Genuchten M.Th., Wierenga P.T. (1976): Masstransfer studies in sorbing porous media: I Analyticalsolutions. Soil Science Society of America Journal,
40: 473–480.
Vogel T., Gerke H.H., Zhang R., Van Genuchten
M.Th. (2000): Modeling flow and transport in a twodimensional dual-permeability system with spatial
variable hydraulic properties. Journal of Hydrology,
238: 78–89.
Received for publication November 14, 2005
Accepted after corrections February 15, 2006

Corresponding author:
Doc. Ing. Radka Kodešová, CSc., Česká zemědělská univerzita v Praze, Fakulta agrobiologie, potravinových
a přírodních zdrojů, katedra pedologie a geologie, Kamýcká 129, 165 21 Praha 6-Suchdol, Česká republika
tel.: + 420 224 382 592, fax: + 420 220 921 643, e-mail: kodesova@af.czu.cz

22