INFILTRATION OF WATER INTO SOIL WITH CRACKS
By V. Novák,1 J. Šimůnek,2 and M. Th. van Genuchten3
ABSTRACT: This paper presents the physical basis of the FRACTURE submodel for simulating infiltration of
precipitation/irrigation water into relatively dry, cracked, fine-textured soils. The FRACTURE submodel forms
part of the HYDRUS-ET variably saturated flow/transport model. Infiltration into the soil matrix is formallydivided into two components: (1) Vertical infiltration through the soil surface; and (2) lateral infiltration via soil
cracks. The first component is described and solved using the 1D Richards’ equation. Excess water that does
not infiltrate through the soil surface is either considered to be runoff, if no soil cracks are present, or routed
into soil cracks from where it may laterally infiltrate into the soil matrix. Horizontal infiltration from soil cracks
into the soil matrix is calculated using the Green-Ampt approach and incorporated as a positive source/sink term
S f in the Richards’ equation describing flow in the matrix. In addition to the hydraulic properties of the soil
matrix, the FRACTURE submodel requires parameters characterizing the soil cracks, notably the specific crack
length per surface area lc and the relationship between crack porosity Pc and the gravimetric soil water content
w. An example problem shows that infiltration from soil cracks can be an important process affecting the soil
water regime of cracked soils. A comparison with the more traditional approach, involving surface infiltration
only, indicates important differences in the soil water content distribution during a rainfall/irrigation event. Thisextension of the classical approach to include crack infiltration significantly improves the identification andprediction of the soil water regime.

INTRODUCTION
Soils often exhibit a variety of heterogeneities such as fractures, fissures, cracks, macropores of biotic origin, and interaggregate pores (Gerke and van Genuchten 1993a; Jarvis
1998). These heterogeneities can significantly affect water andsolute movement in soils by creating nonuniform velocity
fields with spatially variable flows. The resulting nonuniform
flow process is often referred to as preferential flow. Preferential flow may be particularly significant in fine-textured,shrinking/swelling soils containing cracks. Infiltration of precipitation and/or irrigation water into such low-permeable soils
is normally very slow and frequently accompanied with surface runoff. The presence of drying cracks often decreases theamount of surface runoff by increasing the total infiltration
rates and, concomitantly, the soil water content in deeper partsof the soil profile. Hence, cracks greatly affect the irrigation
efficiency, as well as minimize the amount of irrigation water
lost to surface runoff. At the same time crack infiltration is
also associated with accelerated solute transport, with surfaceapplied solutes generally penetrating much deeper into the soil
profile, thus posing risks for soil and ground-water pollution,
including nutrients moving quickly below the root zone of
agricultural crops.
Adequate descriptions of nonuniform water infiltration into
initially dry, cracked, fine-textured soils are lacking in mostsoil-water-plant-atmosphere models. Ignoring the infiltrationof water via soil cracks into the soil matrix usually leads toseverely underestimated infiltration rates, too high predictionsof water accumulating at or near the soil surface, overestimation of surface runoff, and, consequently, unrealistic descriptions of the soil water regime.
Soils generally contain two types of pores (Doležal andKutı́lek 1972): (1) Micropores that form a relatively homo1
Sr. Sci., Inst. of Hydro., Slovak Acad. of Sci., P.O. Box 94, Račianská
75, 830 08 Bratislava, Slovak Republic.
2
Asst. Res., U.S. Salinity Lab., US Dept. of Agr., ARS, 450 West Big
Springs Rd., Riverside, CA, 92507-4617.
3
Res. Leader, U.S. Salinity Lab., US Dept. of Agr., ARS, 450 West
Big Springs Rd., Riverside, CA.
Note. Discussion open until July 1, 2000. To extend the closing date
one month, a written request must be filed with the ASCE Manager ofJournals. The manuscript for this paper was submitted for review andpossible publication on March 22, 1999. This paper is part of the Journalof Irrigation and Drainage Engineering, Vol. 126, No. 1, January/February, 2000. 䉷ASCE, ISSN 0733-9437/00/0001-0041–0047/$8.00 ⫹$.50 per page. Paper No. 20554.

geneous system, depending mostly upon soil texture; and (2)macropores that can be divided into major groups. One groupconsists of relatively stable macropores created by plant roots
and soil fauna. Such macropores are more or less independentof the soil water content. A second group consists of soil
cracks whose dimensions are dependent upon the soil watercontent and on the soil mechanical properties. If the soil matrix
reflects a bimodal porosity, micropores can be further subdivided into interpedal and intrapedal micropores (Kutı́lek andNielsen 1994).
When modeling a soil containing relatively stable macropores, its saturated hydraulic conductivity Ks can be characterized by a lumped parameter to be estimated at the scale ofa representative elementary volume of soil (e.g., using the
Guelph permeameter) (Reynolds 1993). Soils containing drying cracks cannot be characterized in such a way as the cracks
are relatively unstable. Their geometry changes with timedepending upon the soil water content, and their hydraulicconductivity is extremely high as compared with the soil
matrix Ks.
The importance of soil cracks, particularly for infiltration
into fine-textured (clay) soils during rainfall or irrigation is
widely recognized [e.g., Kutı́lek and Novák (1976) and Mitchell and van Genuchten (1993)]. Fig. 1 shows a typical crack
net developed under a spring barley canopy on June 6, 1995,
at the Trnava experimental site in southern Slovakia. The vol-FIG. 1. System of Soil Cracks in 1-m Square Area of Soil Surface at Trnava Experimental Site in Southern Slovakia, June 6,
1995. A Real Crack Porosity at Soil Surface Pc(0) = 0.046 m2/m2
and the Specific Length of the Cracks lc = 11.1 m/m2

JOURNAL OF IRRIGATION AND DRAINAGE ENGINEERING / JANUARY/FEBRUARY 2000 / 41

ume of cracks under a 1-m2 area of soil surface at the time of
the measurement was about 0.005 m3 (i.e., a water layer of
5 mm was needed to fully fill the crack system). The surfaceof the crack walls, again under a 1-m2 area of the soil surface,
was 12.1 m2, which implies that the maximum potential infiltration surface is 12.1 times larger than the soil topographic
surface. These numbers illustrate the potential importance of
cracks in affecting the soil water regime of fine-textured,swelling soils.
A variety of models describing the movement of water anddissolved solutes in soils containing cracks have been developed, most of them based on schematic representations of thecrack network (van Genuchten 1991; Ma and Selim 1998; van
Genuchten and Sudicky 1999). Results obtained with thesemodels are still relatively poor, in part because of conceptual
weaknesses in the models, and in part because of difficulties
in accurately estimating crack input parameters. Flow in structured porous media is frequently also described using dualpermeability models (Pruess and Wang 1987; Gerke and van
Genuchten 1993a; Jarvis 1994, 1998). Approaches of this typeassume that the soil consists of two regions, one associatedwith macropores (the crack network), and the other with theless-permeable matrix region. The difficulty in applying thisapproach to cracked soils is that flow in both regions is described using the Richards’ equation, an assumption that is
likely not valid for flow in large drying cracks. For example,
the MACRO dual-permeability model of Jarvis (1994) uses
the Richards’ equation and the convection-dispersion equation
to model soil water flow and solute transport in the matrix(soil micropores), whereas flow and transport in the macropores is described using a simplified capacitance-type approach. The exchange between matrix and macropore regions
is simulated with an approximate quasi-empirical first-order
rate expression (Jarvis 1994).
Another class of models uses the statistical properties ofcracks (e.g., depth and width). This approach does not requiredetailed knowledge of the crack system and/or the spatial distribution of macropores (Slawinski et al. 1996; van Dam et al.
1997). For example, the SWIMv2.1 model (Verburg et al.
1996) assumes that water and solutes from runoff, after the
infiltration capacity of the soil is exceeded and a critical depthof water is formed at the soil surface, can move directly to
specified bypass depths. Downward bypass will not occur if
the matrix absorbs all water through the soil surface. The latest
version of the SWATRE model, SWAP (van Dam et al. 1997),
also adds infiltrated water from cracks as a source term toRichards’ equation. This model, however, assumes lateral infiltration from the cracks into the soil matrix to be constant in
time, as long as water is present in the cracks, with the infiltration rate changing only as a function of the active crack
area. The infiltration rate itself is calculated using Darcy’s law.
In this study the FRACTURE submodel that quantitatively
describes the infiltration of rain or irrigation water into an
initially relatively dry, cracked, fine-textured soil is presented.
Variably saturated flow in the soil matrix is described using
the Richards’ equation. Excess water at the soil surface thatcannot infiltrate because of a low infiltration capacity is either
removed by surface runoff, or allowed to fill the soil cracks
from where it infiltrates laterally into the soil matrix. Infiltration from the soil cracks into the soil matrix is locally assumed
to be horizontal and described using the Green-Ampt approach.
CONCEPTUAL MODEL
A schematic of the conceptual model forming the basis of
the FRACTURE submodel is shown in Fig. 2. Precipitation
or irrigation q0(t) falls on the soil surface and, as long as thesoil surface is unsaturated, is equal to the actual infiltration

FIG. 2. Schematic of FRACTURE Submodel. Applied Surface
Flux q0(t ) Is Divided between Soil Surface Infiltration Rate q (t )and Flow into Cracks qf (t ). Sf (z, t ) Represents Horizontal Infiltration Rate from Cracks into Soil Matrix

rate q(t). The model assumes that all precipitation falls directlyon the soil surface. After the soil surface becomes saturated(at ponding time tp), excess water is first used to form a surface
layer, the maximum thickness of which is a function of thesurface roughness. The soil surface is hence assumed to havea certain roughness that causes the detention of some waterand prevents it from running off or flowing into cracks. Although the surface roughness can change with time due toimpact of the kinetic energy of rain (Verburg et al. 1996), theroughness is considered to be constant with time. It is further
assumed that surface runoff or flow into cracks can start only
after the surface water layer reaches a critical thickness hs. At
that time water starts flowing into the cracks as long as theprecipitation or irrigation rate is higher than the actual surface
infiltration flux. When the applied surface flux becomes
smaller than the actual infiltration rate, flow into cracks stops
and water in the soil surface layer infiltrates directly into thesoil until it is completely used.
The infiltration process can thus be divided into several
stages:1. Unsaturated infiltration when the soil surface is still unsaturated: q(t) = q0(t), tp < t
2. The formation (or disappearance) of a surface layer ofwater after the soil surface becomes saturated: q(t) < q0(t)[or q(t) > q0(t)]3. Flow into cracks when the soil surface is saturated andthe surface water layer has reached a certain criticalheight hs: q(t) < q0(t)4. Runoff when cracks are either full or not considered: q(t)< q0(t)5. Horizontal infiltration from cracks into the soil matrix
It is further assumed that soil cracks during a certain precipitation event do not change their dimensions. The physicalprocesses of crack formation and swelling, and particularly its
quantification into mathematical models, still needs further investigation.
Horizontal infiltration of water into the soil matrix is considered only through crack surfaces that are in direct contact
with water standing in the cracks. Hence, the model does notconsider infiltration of film water flowing along crack walls.
It is further assumed that excess surface water moves directly
to the prevailing water level in the cracks.
MATHEMATICAL MODEL
The FRACTURE submodel described here is part of the
HYDRUS-ET code (Šimůnek et al. 1997). The 1D Richards’

42 / JOURNAL OF IRRIGATION AND DRAINAGE ENGINEERING / JANUARY/FEBRUARY 2000

equation is assumed to describe water flow in the soil matrix.
Matrix and preferential flow are mutually linked using an extension of the Richards’ equation as follows (Feddes et al.
1988):
⭸␪ ⭸
=
⭸t ⭸z

冋

冉

⭸h
⫹1
⭸z

K(h, z)冊册

⫺ Sr (z) ⫹ Sf (z)(1)

where h = pressure head; ␪ = volumetric water content; t =time; z = vertical coordinate (positive upward); K = unsaturated hydraulic conductivity; Sr (z) = sink term quantifying thevolume of water extracted from soil by roots; and Sf (z) = horizontal infiltration rate of water from the water-filled part ofsoil cracks into the soil matrix. This term is calculated using
the Green-Ampt approach (Green and Ampt 1911)
Sf =

冉

Kh(z)

冊

h0 ⫺ hf
lf

Ac

(2)

where Kh = hydraulic conductivity of the crack-matrix interface; h0 = positive pressure head at the point of infiltration; hf= pressure head at the wetting front at horizontal distance lf
away from the crack surface; and Ac = specific surface of thecracks. The hydrostatic pressure is used as the pressure headh0 at nodes where lateral fracture flow occurs. The wetting
front distance lf can be calculated according to the Green-Ampt
approach as follows:
lf =

冑

2Kh

h0 ⫺ h f
tf
␪s ⫺ ␪ i(3)

where ␪i and ␪s = initial and saturated volumetric water contents, respectively; and tf = time interval since the start ofinfiltration.
The Richards’ equation for flow in the matrix is solved numerically subject to a set of initial and boundary conditions
as summarized below.

⫺Ks

冉

冊

⭸h
⫹1
⭸z

h(z) = hi (z)fort = t0(4)

where hi = initial pressure head; and t0 = initial time.
Upper Boundary Conditions
The upper boundary condition during an infiltration event
depends on the actual status of the soil surface. Three different
boundary conditions at the soil surface are used to characterizedifferent stages of the infiltration process.
Stage 1. When the soil surface is unsaturated a flux boundary condition is used
⫺K

冉

⭸h
⫹1
⭸z

冊

= q0(t)

for

h<0

(5)

where q0 = applied surface flux (i.e., the difference between
evaporation and precipitation/irrigation). Because of the orientation of the z-coordinate, q0 is positive for evaporation andnegative for infiltration. Although HYDRUS-ET considers
evaporation from the soil surface (not across crack walls),
evaporation is negligibly small for a timescale used in theexample given below.
Stage 2. When the surface is saturated and the surface water layer is either being formed or depleted, a ‘‘surface reservoir’’ boundary condition is used (Mls 1982):

dh
dt

for

0 < h < hs(6)where Ks = saturated hydraulic conductivity of the soil matrix;
and hs the critical thickness of the water layer on the soil
surface when surface runoff is initiated, or when water starts
flowing into cracks. The surface reservoir boundary condition[(6)] permits a water layer to either build up at the soil surface,or be used up for infiltration. The left-hand side represents theactual infiltration rate into the soil profile through the soil surface. The first term on the right side, q0, represents the applied
surface flux, and the second term is the change in the thicknessof the water layer at the soil surface.
Boundary condition (6) applies when the soil surface is already saturated, but the value of the critical head hs is not yet
reached. During this time the surface water layer is being
formed, depending upon the actual infiltration rate and theapplied surface flux. Boundary condition (6) is also used after
the precipitation/irrigation event has ended (or when the precipitation/irrigation rate is or has become lower than the sumof the actual infiltration rate and evaporation), and surface water is being used for infiltration and evaporation.
Stage 3. When the surface layer reaches the critical thickness hs, a boundary condition is used that either specifies theflow rate into the soil cracks, or the rate of runoff. In that caseh = hsfor兩q0(t)兩 > 兩q(t)兩

(7)

where q = actual infiltration rate. Flow into the cracks, or
surface runoff qf , is calculated as follows:
qf (t) = q0(t) ⫹ Ks

冉

⭸h
⫹1
⭸z

冊

(8)

Boundary condition (7) is used when the surface layer
reaches the critical thickness hs; the applied surface flux q0 is
still larger that the actual infiltration rate q; and all excesswater either flows into cracks or is removed by surface runoff.
The amount of water in the cracks Vf is calculated as

Initial Condition
The initial condition is the result of hydrologic processes
prior to the infiltration event and given by the soil water pressure head profile= q0(t) ⫺

Vf =

冉冕

qf (t) dt ⫺

t

冕
t

冊

qfs(t) dt(9)where qfs = water flux rate from the cracks into the soil matrix.
Lower Boundary Conditions
A variety of boundary conditions, the same as those usedin the HYDRUS-ET model (Šimůnek et al. 1997), can be prescribed at the lower boundary. The conditions include constantor variable pressure heads or fluxes, seepage faces, and deepor free drainage conditions.
Input Data
In addition to parameters needed by HYDRUS-ET, the coupled HYDRUS-ET FRACTURE model needs parameterscharacterizing flow into and from the soil cracks. First, themodel requires the crack porosity Pc as a function of the soil
gravimetric water content w, Pc(w) (i.e., a soil shrinkage characteristic curve) (Mitchell 1992). This relationship is a soil
characteristic and can be estimated in the laboratory on undisturbed soil samples. As the soil sample slowly dries, changes
in its height and diameter are measured simultaneously withchanging soil water content w. The shrinkage curve Pc(w) canbe easily estimated assuming that the soil sample deformation
is equal to the soil crack porosity Pc (Novák, unpublished paper, 1999). The crack porosity profile Pc(z) can then be estimated using the measured or calculated soil water content profile w(z) and the relationship Pc(w). Second, the specific lengthof cracks lc per unit soil surface area is needed. This length

JOURNAL OF IRRIGATION AND DRAINAGE ENGINEERING / JANUARY/FEBRUARY 2000 / 43

is assumed to be constant over the depth of crack formation(i.e, constant between the soil surface and the bottom part ofthe cracks). The length can be estimated in the field by direct
measurement or by image analysis of the site under consideration.

Horizontal infiltration from the cracks into the soil matrix
is simulated using the hydraulic conductivity of the crack-matrix interface Kh [see (2) and (3)], which is related to the matrix
saturated hydraulic conductivity Ks as follows:

MODEL APPLICATION

where rk = reduction factor accounting for hydraulic resistances across the crack-soil matrix interface. Thomas et al.
(1992) conducted experiments showing that due to coatings
the hydraulic conductivity of the fracture-matrix interface of
fracture rock may be much smaller than the matrix conductivity. Gerke and van Genuchten (1993b) and van Genuchten and
Sudicky (1999) reviewed experimental evidence indicatingthat the interface hydraulic conductivity can be smaller by several orders of magnitude than the conductivity of the matrix
interior. They reasoned that this is probably due to the changed
physical and chemical properties of the interface as a result of
repeated opening and closing of fractures and of coating of
the interface by relicts of roots and other organic matter. Techniques for directly measuring the hydraulic conductivity offracture-matrix interfaces (or of soil macropore walls) are only
recently being perfected (Gerke and Köhne 1999). Field observations have shown that soil cracks form along the sameinternal soil surfaces that represent areas of minimum mechanical strength.
The following crack characteristics were used in our numerical experimental: zc = 40 cm, af = 0.178, Pc 0 = 0.429, andrk = 0.1, where zc is the depth of soil cracks. Simulations werecarried out assuming an irrigation rate q0 of 25 cm/day for aperiod of t = 2 h. Surface runoff, or flow into cracks, wasassumed to start when the thickness of the soil surface waterlayer reached a critical value hs of 0.1 cm.

An example is now presented to illustrate differences between the traditional approach of calculating infiltration andour method for simulating infiltration into a soil containing
drying cracks. The model was applied to a fine-textured (clay)
soil from the Trnava area of southern Slovakia. The hydraulicparameters of the soil are as follows: residual soil water content ␪r = 0.03, saturated soil water content ␪s = 0.407, parameter ␣ = 0.774 cm⫺1, parameter n = 1.235, and saturated hydraulic conductivity of the soil matrix Ks = 5 cm/day (vanGenuchten 1980). The soil cracks were estimated in the field
to have a specific length of the soil cracks lc = 0.046 m/m2.
The shrinkage curve Pc(w) (i.e., the relationship between crack
porosity and the soil water content) was measured in the laboratory. The Pc(w) data, and its linear approximation for theupper soil layer at the Trnava experimental site, are shown in
Fig. 3. The linear approximation is given by (Novák, unpublished paper, 1999)
Pc = ⫺af w ⫹ Pc 0

(10)

where af = slope; and Pc 0 = maximum crack porosity corresponding to a zero soil water content w.

Kh = rk Ks(11)

RESULTS AND DISCUSSION

FIG. 3. Crack Porosity Pc versus Gravimetric Soil Water Content w for Trnava Experimental Site in Southern Slovakia

Fig. 4 shows calculated soil water content profiles for different times during infiltration into the soil for three scenarios:(1) A soil without crack with water accumulating at the soil
surface [Fig. 4(a)]; (2) a soil without cracks but with surface
runoff [Fig. 4(b)]; and (3) a soil with cracks [Fig. 4(c)]. The
figure shows a distinct area above the maximum crack depth(40 cm) through which water and solutes infiltrate into thesoil profile when crack infiltration is being considered (Scenario 3).
The infiltration rates corresponding to the different scenarios
described above are shown in Fig. 5. The vertical infiltration
rate with surface runoff present (Scenario 2) is the same as
the infiltration rate through the soil surface of the cracked soil

FIG. 4. Soil Water Content Profiles during Infiltration into Soil: (a) without Cracks with Water Accumulating on Soil Surface; (b) without Cracks with Surface Runoff; (c) with Cracks
44 / JOURNAL OF IRRIGATION AND DRAINAGE ENGINEERING / JANUARY/FEBRUARY 2000

(Scenario 3). This rate (Curve 3) is slightly lower than theinfiltration rate into the soil when water is allowed to accumulate on top of the soil surface (Scenario 1; no surface runoff, and no cracks) (Curve 2, Fig. 5). The instantaneous infiltration rate, as well as the cumulative infiltration rate, into thesoil is higher when ponding at the soil surface is allowed
above the critical surface layer (Scenario 1) as compared to
Scenarios 2 and 3 where the surface water layer cannot exceed
the critical thickness hs. This is caused by water accumulating
on top of the soil surface that increases the pressure head and,
consequently, the pressure gradient across the soil surface.
Ponded infiltration will continue after the end of the precipitation event until all water accumulated on the soil surface has
infiltrated (Curve 2, Fig. 5).
The important differences between the various infiltration

events involving different boundary conditions (i.e., with andwithout a surface layer, runoff, or cracks) are further illustrated
in Fig. 6, which shows plots of the cumulative infiltration rateof water for all simulated cases. Depending upon the scenarios
involved, applied (rain or irrigation) water exceeding the infiltration capacity of the soil is either temporarily stored in asurface layer or in the soil cracks (Fig. 7), or is lost due tosurface runoff.
The ponded water layer on the soil surface for the first simulation (Scenario 1) starts forming at 0.01 day, reaches a maximum of about 1.3 cm at the end of the irrigation event (0.1day), and is fully depleted at about 0.305 day. Fig. 7 presents
the amount of water stored in the cracks for the third simulation (Scenario 3). Although water flowing into the cracks(Curve 5, Fig. 5) can significantly increase the total actual

FIG. 5. Infiltration Rates versus Time for Three Scenarios Considered in This Study: Curve 1, Irrigation Rate q0 (t ); Curve 2, Infiltration q(t ) into Soil without Cracks with Water Accumulating on Soil Surface; Curve 3, Infiltration q(t ) through Soil Surface into Soil with
Surface Runoff (with or without Cracks Present); Curve 4, Infiltration qf (t ) from Cracks into Soil Matrix; Curve 5, Flow into Cracks[q0 (t ) ⴚ q(t )]FIG. 6. Cumulative Infiltration versus Time for Three Scenarios Considered in This Study: Curve 1, Cumulative Irrigation Rate; Curve
2, Cumulative Infiltration into Soil without Cracks with Water Accumulating on Soil Surface; Curve 3, Cumulative Infiltration throughSoil Surface into Soil with Surface Runoff (with or without Cracks Present); Curve 4, Cumulative Infiltration from Cracks into Soil Matrix; Curve 5, Cumulative Flow into Cracks; Curve 6, Total Cumulative Infiltration into Soil Matrix for Soil with Cracks
JOURNAL OF IRRIGATION AND DRAINAGE ENGINEERING / JANUARY/FEBRUARY 2000 / 45

infiltration rate into the soil matrix, applied water accumulated
on top of the soil surface prolongs the infiltration process until
after the end of precipitation/irrigation (Curve 2, Fig. 5). Notice that the total cumulative infiltration rate into the cracked
soil [i.e., infiltration through both the soil surface and the crack
walls (Curve 6 in Fig. 6)] most closely resembles the applied
surface flux (Curve 1). Thus, the cracks can significantly increase the total infiltration capacity of the soil in comparison
with other simulated scenarios when the applied irrigation rate
is higher than the saturated hydraulic conductivity.
The unique contribution of cracks to the infiltration process
is clearly visible in Figs. 5 and 6. Infiltration via the cracks(Curve 4) is initially zero and remains so until after 0.01 day
when ponding develops and water flows into the cracks. Thisdelay includes the time interval during which the thin waterlayer developed on the soil surface. Only when the water layer

FIG. 7.

reached the critical thickness hs did water start to flow into thecracks and from there infiltrated laterally into the soil matrix.
The infiltration rate from the cracks subsequently increased
quickly as more water flowed into the cracks and a larger
surface area of the cracks became involved in the infiltration
process. After reaching a maximum, the infiltration rate fromthe cracks started to decrease before the end of the precipitation/irrigation event, mostly because of the decreasing rate ofhorizontal infiltration (in response to decreasing lateral pressure gradients). However, the water level as well as the water
volume in the crack kept increasing up to the end of precipitation (Fig. 7). The highest infiltration rates from the cracks
into the soil occurred across crack walls that just became active for infiltration, due to locally high gradients in the soil
water pressure head associated with infiltration into relatively
dry soil. As shown in Fig. 8, infiltration rates from the crack

Water Level (Depth below Soil Surface) and Amount of Water in Cracks during Infiltration

FIG. 8.

Local Infiltration Rates through Crack Walls

46 / JOURNAL OF IRRIGATION AND DRAINAGE ENGINEERING / JANUARY/FEBRUARY 2000

to the soil slowly decrease as the infiltration process proceeds.
Notice also that, as expected, the infiltration maximum slowly
shifts toward the soil surface where relatively dry soil is continually being encountered as the water level rises.
CONCLUSIONSThe example presented in the paper demonstrates the importance of soil cracks in determining rates of infiltration intoa soil during irrigation. In our example the infiltration capacityof the soil without cracks is less than half (34%) of the infiltration capacity of the soil with cracks. This corresponds very
closely with values measured by Mitchell and van Genuchten(1993) for two fallow irrigations (36.1 and 35.4%) of acracked soil in a large weighing lysimeter. A comparison withthe more traditional approach involving only surface infiltration indicates important differences in the soil water content
distribution during an irrigation event. An extension of theclassical Richard’s equation approach to include crack infiltration can lead to significantly improved predictions of the water
regime of swelling clay soils.
APPENDIX I.

REFERENCES

Doležal, F., and Kutı́lek, M. (1972). ‘‘Flow of water in swelling soils.’’
Proc., 2nd Symp. Fundamentals of Transport Phenomena in PorousMedia, IAHR-ISSS, Vol. 1, 292–305.
Feddes, R. A., Kabát, P., van Bakel, J., Bronswijk, J. J. B., and Halbertsma, J. (1988). ‘‘Modeling soil water dynamics in the unsaturated
zone—state of art.’’ J. Hydro., Amsterdam, 100, 69–111.
Gerke, H. H., and Köhne, J. M. (1999). ‘‘Impact of clayey coatings on
hydraulic properties of soil aggregates.’’ Geophys. Res. Abstracts, 1(2),324.
Gerke, H. H., and van Genuchten, M. Th. (1993a). ‘‘A dual-porosity
model for simulating the preferential movement of water and solutes
in structured porous media.’’ Water Resour. Res., 29(2), 305–319.
Gerke, H. H., and van Genuchten, M. Th. (1993b). ‘‘Evaluation of a firstorder water transfer term for variably saturated dual-porosity flow models.’’ Water Resour. Res., 29(4), 1225–1238.
Green, W. H., and Ampt, G. A. (1911). ‘‘Studies on soil physics: 1. Flowof air and water through soils.’’ J. Agr. Sci., 4, 1–24.
Jarvis, N. J. (1994). ‘‘The MACRO model (version 3.1)—Technical description and sample simulations.’’ Rep. and Dissertations 19, Dept.
of Soil Sci., Swedish University of Agricultural Sciences, Uppsala,Sweden.
Jarvis, N. J. (1998). ‘‘Modeling the impact of preferential flow on nonpoint source pollution.’’ Physical nonequilibrium in soils: modelingand application, J. M. Selim and L. Ma, eds., Ann Arbor Press, Chelsea, Mich., 195–221.
Kutı́lek, M., and Nielsen, D. R. (1994). Soil hydrology. Catena Verlag,
Cremlingen, Germany.
Kutı́lek, M., and Novák, V. (1976). ‘‘The influence of soil cracks upon
infiltration and ponding time.’’ Proc., Symp. Water and Heavy Soils,
M. Kutı́lek and J. Šutor, eds., Vol. 1, Institute of Hydrology, Slovak
Academy of Sciences, Bratislava, Slovakia, 26–134.
Ma, L., and Selim, H. M. (1998). ‘‘Coupling of retention approaches tophysical nonequilibrium models.’’ Physical nonequilibrium in soils:
Modeling and application, J. M. Selim and L. Ma, eds., Ann Arbor
Press, Chelsea, Mich., 83–115.
Mitchell, A. R. (1992). ‘‘Shrinkage terminology: Escape from ‘normalcy’.’’ Soil Sci. Soc. Am. J., 56, 993–994.
Mitchell, A. R., and van Genuchten, M. Th. (1993). ‘‘Flood irrigation ofa cracked soil.’’ Soil Sci. Soc. Am. J., 57, 490–497.
Mls, J. (1982). ‘‘Formation and solution of fundamental problems of vertical infiltration.’’ Vodohosp. Čas, 30, 304–313 (in Czech).
Pruess, K., and Wang, J. S. Y. (1986). ‘‘Numerical modeling of isothermal
and nonisothermal flow in unsaturated-saturated fractured rock—a review.’’ Flow and transport through unsaturated fracture rock, D. D.
Evans and T. J. Nicholson, eds., American Geophysics Union, Washington, D.C., 11–22.
Reynolds, W. D. (1993). ‘‘Unsaturated hydraulic conductivity: Field measurement.’’ Soil sampling and methods of analysis, M. R. Carter, ed.,Lewis, Boca Raton, Fla., 633–644.
Šimůnek, J., Huang, K., Šejna, M., van Genuchten, M. Th., Majerčák, J.,
Novák, V., and Šutor, J. (1997). The HYDRUS-ET software package
for simulating the one-dimensional movement of water, heat and mul-

tiple solutes in variably—Saturated media. Version 1.1. Institute ofHydrology, Slovak Academy of Sciences, Bratislava, Slovakia.
Slawinski, C., Sobczuk, H. A., and Walczak, R. T. (1996). ‘‘Submodel of
bypass flow in cracking soil. Part 1. Theory.’’ Int. Agrophysics, Lublin,
Poland, 10, 189–196.
Thoma, S. G., Gallegos, D. P., and Smith, D. M. (1992). ‘‘Impact offracture coating on fracture/matrix flow interaction in unsaturated porous media.’’ Water Resour. Res., 28, 1357–1367.
van Dam, J. C., et al. (1997). ‘‘Theory of SWAP, version 2.0. Simulationof water flow, solute transport and plant growth in the soil-water-atmosphere-plant environment.’’ Rep. 71, Tech. Document 45, Dept. Water Resources, Wageningen Agricultural University, The Netherlands,
Dienst Landbouwkundig Onderzoek (Agricultural Research Department) Winand Staring Centre, Wageningen.
van Genuchten, M. Th. (1980). ‘‘A closed-form equation for predicting
the hydraulic conductivity of unsaturated soils.’’ Soil Sci. Soc. Am.
Proc., 44, 892–898.
van Genuchten, M. Th. (1991). ‘‘Progress and opportunities in hydrologic
research.’’ Rev. of Geophys., Supplement, 189–192.
van Genuchten, M. Th., and Sudicky, E. A. (1999). ‘‘Recent advances invadose zone flow and transport modeling.’’ Vadose zone hydrology:
cutting across disciplines, M. B. Parlange and J. W. Hopmans, eds.,
Oxford University Press, New York, 155–193.
Verburg, K., Ross, P. J., and Bristow, K. L. (1996). ‘‘SWIMv2.1 user
manual.’’ Divisional Rep. No. 130, Commonwealth Scientific and Industrial Research Organization, Australia.

APPENDIX II.

NOTATION

The following symbols are used in this paper:
Ac = specific surface of cracks (L2 L⫺3);
af = slope of Pc(w) function (nondimensional);
h = pressure head (L);hf = pressure head (negative) at wetting front at distance lf
away from crack surface (L);hi = initial pressure head (L):hs = critical head (L) (i.e., thickness of water layer at soil surface when surface runoff is initiated) or when water starts
flowing into soil cracks (L);h0 = positive pressure head at point of infiltration (L);K = unsaturated hydraulic conductivity (L T⫺1);
Kh = hydraulic conductivity of crack-matrix interface (L T⫺1);
Ks = saturated hydraulic conductivity of soil matrix (L T⫺1);
lc = specific length of cracks per unit soil surface area (L L⫺2);
lf = distance of wetting front from infiltration surface (L);n = parameter of van Genuchten’s equation (1980) (nondimensional);
Pc = crack porosity (L3 L⫺3);
Pc 0 = maximum crack porosity corresponding to zero soil watercontent w (L3 L⫺3);
q = actual infiltration rate (L3 L⫺2 T⫺1);
qf = flow into soil cracks or surface runoff (L3 L⫺2 T⫺1);
qf s = water flux rate from cracks into soil matrix (L3 L⫺2 T⫺1);
q0 = potential infiltration rate (L3 L⫺2 T⫺1);
rk = reduction factor for saturated hydraulic conductivity ofcrack-soil matrix interface (nondimensional);
Sf = horizontal infiltration rate of water from cracks into soil
matrix (L3 L⫺3 T⫺1);
Sr = sink term (L3 L⫺3 T⫺1) quantifying volume of water extracted from soil by roots (root extraction term);t = time (T);tf = time interval since the start of infiltration (T);tp = ponding time (T);t0 = initial time (T);Vf = volume of water in cracks (L3 L⫺2) under unit square areaof soil surface;w = soil gravimetric water content (M3 M⫺3);z = vertical coordinate (L) (positive upward);
␣ = parameter of van Genuchten’s equation (1980) (L⫺1);
␪ = volumetric water content (L3 L⫺3);
␪i = initial volumetric water contents (L3 L⫺3);
␪r = residual soil water content (L3 L⫺3); and␪s = saturated soil water content (L3 L⫺3).

JOURNAL OF IRRIGATION AND DRAINAGE ENGINEERING / JANUARY/FEBRUARY 2000 / 47