Geoderma 281 (2016) 39–48

Contents lists available at ScienceDirect

Geodermajournal homepage: www.elsevier.com/locate/geoderma

The effects of rock fragment shapes and positions on modeled hydraulicconductivities of stony soils
Hana Hlaváčiková a,⁎, Viliam Novák a, Jirka Šimůnek babInstitute of Hydrology, Slovak Academy of Sciences, Dúbravská cesta 9, 84104 Bratislava, Slovakia
Dept. of Environmental Sciences, Univ. of California Riverside, Riverside, CA 92521, USAa r t i c l e

i n f oArticle history:
Received 21 December 2015
Received in revised form 6 June 2016
Accepted 27 June 2016
Available online xxxx
Keywords:
Soil matrix
Rock fragments
Effective saturated hydraulic conductivity
HYDRUS-2D modelinga b s t r a c t
Mountainous soils usually contain a large number of rock fragments, particles with a diameter larger than 2 mm,
which can inﬂuence soil hydraulic properties that are required to quantitatively describe soil water movement instony soils. The objective of this study was to numerically estimate both the saturated hydraulic conductivity of astony soil and its dependence on a relative content of rock fragments (stoniness), and the shape, position and distribution of rock fragments in a soil matrix. The assessment method was based on a numerical version of Darcy's
classic experiment that involved steady-state ﬂow through a porous material under a unit hydraulic gradient. Our
experiments, involving hypothetical stony soils in this particular case, were simulated using mainly the two-dimensional (2D) numerical model, HYDRUS-2D. A limited number of simulations were carried out using a threedimensional HYDRUS model. Three different shapes of hypothetical rock fragments were used in the study: asphere, an ellipsoid with two different positions, and a pyramid, all represented by their 2D cross-sections (i.e.,a circle, an ellipse, and a triangle, respectively). The mean relative effective saturated hydraulic conductivity(Krs) for the same stoniness was almost the same for all simulated scenarios and ﬁne soil textures. A stoninessbetween 0.07 and 0.5 cm3 cm−3 can cause a decrease of Krs in the range of 0.17–0.70. Numerical experiments
were divided into 3 scenarios. The largest and the smallest values of Krs were different for different shapes of
RFs (scenario A), different orientations of the slab-sided elliptical RFs (scenario B), and regular or irregular distributions of spherical RFs (scenario C). The largest difference between Krs values (0.26) was found in scenario B
when the slab-sided elliptical RFs were oriented either horizontally or vertically for stoniness of 0.24 or 0.31
cm3cm−3. Simulated Krs values were underestimated in all scenarios as compared to the Ravina and Magier(1984) function. The smallest differences (−1.1%–2.5%) between numerically simulated and calculated (theCorring and Churchill (1961) method for a cylindrical shape of RFs) Krs values were found for scenario A withits 2D representation of spherical rock fragments. Calculated (the Corring and Churchill (1961) method for aspherical shape of RFs) Krs values corresponded well with those simulated using a 3D representation of sphericalrock fragments. Numerical models provide a unique opportunity to evaluate the effects of different factors on thesaturated hydraulic conductivity of stony soils that may be nearly impossible to measure in practice.
© 2016 Elsevier B.V. All rights reserved.

1. Introduction
Soils containing a signiﬁcant fraction of rock fragments (RF) are generally denoted as stony soils and can be present in many forested andmountainous areas. The shape, size, degree of weathering, and geological origin of rock fragments can all strongly inﬂuence the soil's
Abbreviations: RF, rock fragments; REV, representative elementary volume; RM, theRavina and Magier (1984) function; CCS, the Corring and Churchill (1961) function for aspherical shape of rock fragments; CCC, the Corring and Churchill (1961) function for acylindrical shape of rock fragments.
⁎ Corresponding author.
E-mail address: hlavacikova@uh.savba.sk (H. Hlaváčiková).

http://dx.doi.org/10.1016/j.geoderma.2016.06.034
0016-7061/© 2016 Elsevier B.V. All rights reserved.

hydrophysical properties, especially the soil's water retention and hydraulic conductivity (Brouwer and Anderson, 2000; Cousin et al.,2003). According to Poesen and Lavee (1994), about 30% of soils in
Western Europe and about 60% of soils in the Mediterranean region
are stony soils. According to Šály (1978), up to 80% of Slovak forest
soils contain stones and their stone content generally increases withdepth. Furthermore, about 47% of Slovak agricultural soils are referred
to as stony soils (Hraško and Bedrna, 1988). The spatial distribution ofrock fragment in hillslopes is mostly controlled by slope gradient andtopographic position (Chen et al., 2012). It is expected that rock fragments, and their size, shape, position, and spatial distribution in thesoil, can strongly inﬂuence the stony soils' properties and can affect40

H. Hlaváčiková et al. / Geoderma 281 (2016) 39–48

soil water movement, inﬁltration (Al-Qinna et al., 2014; Chen et al.,2012), and the occurrence of runoff (Hlaváčiková et al., 2015).
Stony soils are composed of a soil matrix, small particles with a diameter of less than 2 mm, and larger rock fragments, e.g., gravel, cobbles, stones, and boulders. The most important characteristics of astony soil are its stone content (stoniness), water retention curves,and the hydraulic conductivity functions of both the soil matrix androck fragments, as well as bulk characteristics representing the stonysoil as a whole. Stoniness (Rv) is the ratio of the volume of rock fragments to the total volume of the soil.
The presence of stones can affect the hydraulic conductivity of a soil
in several ways. On the one hand, stones reduce the effective cross-sectional area through which water ﬂows, and combined with the fact thatan increase in stoniness results in greater curvatures of ﬂow paths, thiscan result in lower hydraulic conductivities (Bouwer and Rice, 1984;
Childs and Flint, 1990; Ma et al., 2010; Novák et al., 2011; Ravina andMagier, 1984). In contrast, shrink-swell phenomena may create temporal lacunar pores (i.e., voids along the soil/stone interface) that cancause preferential ﬂow and thus an increase in the saturated hydraulicconductivity (Sauer and Logsdon, 2002; Shi et al., 2008; Verbist et al.,2009; Zhou et al., 2009).
Rock fragments are relatively large compared to ﬁne soil particles ofthe soil matrix. Therefore, it is necessary to evaluate the bulk soil characteristics for a “representative elementary volume” (REV) (Bear,
1972), the size of which depends mostly on the size and spatial distribution of the rock fragments. The larger the rock fragments are in the stonysoil, the larger the REV needed. Buchter et al. (1994) recommended thatthe dry mass of a stony soil sample should be at least 100 times the massof the largest particle. However, there is no rule for how large the REV ofa stony soil should be for measuring its hydraulic characteristics. The dimension of the REV can vary from decimeters to meters, with its volume
extending to about 1 m3 when the rock fragments have a diameter of
10 cm or larger. The presence of rock fragments further presents problems for measuring the bulk soil hydraulic properties, water contents,
water potential, or ﬂow regime in general, due to such practical issues
as the difﬁculty of inserting probes in stony soils (e.g., TDR probes andtensiometers) or installing lysimeters (Cousin et al., 2003; Ma et al.,2010).
Since it is technically difﬁcult to perform hydraulic conductivity
measurements on large samples with different stoniness, Novák et al.
(2011) proposed the use of numerical models to simulate the classic
Darcian ﬂow experiment and to calculate the corresponding saturated
hydraulic conductivity. This was done by embedding spherically shaped
stones of different sizes (5, 10, and 20 cm in diameter) into a soil matrixof known hydraulic conductivity and then calculating the effective saturated hydraulic conductivity of the bulk sample with stones (Kbs ).
Novák et al. (2011) showed that the effective hydraulic conductivityof a soil with a given stoniness is smaller when it contains a single“large” stone than when it contains multiple smaller stones. However,they only considered circular stones in the stoniness range of 0.07–
0.31 cm3 cm−3. From their study, it is not clear if the ﬁne soil texture
markedly affects simulated saturated hydraulic conductivities. However, their study indicated that one can similarly assess other factors thatcan potentially inﬂuence the hydraulic resistance to water ﬂow suchas different shapes, orientations, positions, and spatial distributions ofrock fragments.
Therefore, the objective of this study is to describe and quantify theinﬂuence of stoniness and different shapes and positions of rock fragments on the bulk (effective) saturated hydraulic conductivity ofstony soils using numerical modeling. This goal is achieved by answering the following questions:
1. How will A) the shape of rock fragments, B) the orientation of rockfragments, and C) the regular and irregular distributions of rock fragments affect the effective saturated hydraulic conductivity of stonysoils?

2. How much will the effective saturated hydraulic conductivities ofstony soils simulated using the numerical model differ from those
calculated using existing empirical equations?
3. Is it possible to assess the effects of soil texture of the soil matrix on
the effective saturated hydraulic conductivities of stony soils using
a numerical model?
4. From simulated results, is it possible to propose a relationship between stoniness and the effective saturated hydraulic conductivity
for different shapes, orientations, or distributions of rock fragments
in a stony soil?
2. Theory
Only a few empirical equations for evaluating an effective saturated
hydraulic conductivity of stony soils exist (e.g., Brakensiek et al., 1986;
Ma et al., 2010; Ravina and Magier, 1984), and most are derived fromlaboratory experiments. The equation of Ravina and Magier (1984) isoften used either in the absolute form to get the effective saturated hydraulic conductivity:
K s b ¼ ð1−Rv ÞK s f

ð1Þ

or in the relative form to get the relative effective saturated hydraulicconductivity, Krs:
K rs ¼

Ksb
Ks f¼ 1−Rv

ð2Þ

where Kbs and Kfs are the effective saturated hydraulic conductivities ofthe bulk soil (denoted by index “b”) and soil matrix (denoted by index“f”), respectively, and Rv is the relative volume of rock fragments in astony soil (stoniness). According to Eq. (2), Krs linearly decreases withan increasing relative volume of rock fragments. The use of this approach is recommended for sandy soils because their deformation due
to water ﬂow is not signiﬁcant. The preliminary assumption is thatproperties of the soil matrix such as porosity are invariant with respect
to the content of rock fragments.
The effective saturated hydraulic conductivity of a stony soil can also
be calculated using the Corring and Churchill (1961) equations. Theirequations were obtained from analytical solutions of heat transfer equations in dispersive environments. The equation for a stony soil withspherical rock fragments with a negligible or zero retention capacitythat are dispersed in a soil matrix can be expressed as follows:
Ksb ¼ Ks f
2ð1−Rv Þ
2 þ Rv

ð3Þ

For cylindrical shaped rock fragments, it can be written as:Ksb ¼ Ks f

1−Rv
1 þ Rv


ð4Þ

Peck and Watson (1979) also derived an equation for spherical rockfragments identical to Eq. (3).
3. Methods
3.1. Numerical Darcy experiments
In this study, numerical Darcy experiments were performed in asimilar fashion to Novák et al. (2011). A virtual stony soil was createdby distributing rock fragments of different shapes and positions in thesoil matrix, the hydraulic properties (i.e., the saturated hydraulic conductivity and the retention curve) of which were assumed to be
known. Using HYDRUS (2D/3D) (Šimůnek et al., 2008), vertical waterﬂow was then simulated in a vertical cross–section (1 × 1 m2) with

H. Hlaváčiková et al. / Geoderma 281 (2016) 39–48

41

Fig. 1. Isometric view (a) and views in the x-(or y-) (b) and z-direction (c) of the three-dimensional transport domain with 403 regularly distributed spheres representing rock fragments(Rv of 0.211).

impermeable side walls. The ﬁnal steady-state conditions were considered to represent steady-state water ﬂow through a stony soil. Since thefull saturation of the soil proﬁle was considered as the initial condition,
the steady-state conditions were reached very quickly; simulations had
to be run for only about 90–150 min. To establish a unit pressure headgradient, the pressure head of 1 cm was applied at both the surface
and bottom boundaries of a 1-m long soil proﬁle. Consequently, the average water ﬂux across the two boundaries was equal to the effective
saturated hydraulic conductivity of a virtual stony soil.
Water movement in the virtual stony soil was simulated using the
HYDRUS model, which numerically solves the Richards equation describing water ﬂow in unsaturated and/or saturated porous media.
This type of model is needed because both saturated and unsaturated
regions are present in the ﬂow domain, with unsaturated regions locallydeveloping below rock fragments. The saturated soil water content,
which corresponds to a zero pressure head, was used as the initialcondition.
Several simplifying assumptions were made. It was assumed thatthe retention capacity and permeability of rock fragments were negligible, and that there was a close contact between rock fragments and thesoil matrix [i.e., no lacunar pores (Fiès et al., 2002; Verbist et al., 2009)were present]. It was further assumed that water ﬂow occurs only inthe soil matrix (i.e., no ﬂow through rock fragments) and thatmacropores and processes of soliﬂuction and sealing can be neglected.
Finally, it was assumed that ﬂow through a fully three-dimensional system can be approximated using two-dimensional (2D) simulations. A
limited number of fully three-dimensional simulations was carried out
to evaluate the correspondence between results of two- and three-dimensional simulations.
Numerical experiments were divided into 3 scenarios. Scenario A
represents regularly distributed RFs of different shapes (in cross-section
circles, triangles, and ellipses), scenario B represents different orientations (horizontal or vertical) of slab-sided ellipses, and scenario C represents regularly and irregularly distributed spherical rock fragments.
Dimensions of all rock fragment shapes are discussed below.

and 8 cm) were positioned both horizontally and vertically. The baseand height of triangles were 13 and 12.07 cm long, respectively. All circles, representing stones, were the same size, with a diameter of 10 cm.
The cross-sectional area of 78.5 cm2 was the same for all rock fragment
shapes (i.e., circles, ellipses, and triangles). In three-dimensional simulations with regularly distributed spherical rock fragments, all spheres
had a diameter of 10 cm and a volume of 523.6 cm3.
In our numerical two-dimensional simulations, stoniness (Rv) is deﬁned as the ratio of the cross-sectional area of rock fragments and thetotal cross-sectional area. It is assumed to be the same as the relative
volume of rock fragments in the soil volume and was considered withvalues of 0.07, 0.16, 0.24, and 0.31 cm3 cm−3. It was not technically possible to arrange the regular distribution of rock fragments for stoninesshigher than 0.31 cm3 cm−3, especially for triangular and ellipsoidal RFs.
In the three-dimensional simulations, stoniness of 0.157 and
0.211 cm3 cm−3 (corresponding with 300 and 403 rock fragments, respectively) was considered (Fig. 1). Simulations with zero stoniness
represented soil samples without rock fragments. Hydraulic propertiesof the soil matrix, assumed to have the texture of either sandy loam,
loam, or clay loam (Table 1), were taken from the soil catalog implemented in the HYDRUS software. These three textures were selected
to represent a range of textures from coarse to ﬁne.
3.3. Regular and irregular rock fragment distributions (scenario C)
Additional numerical experiments were performed to evaluate theinﬂuence of regular and irregular distributions of circular rock fragments in a stony soil on the effective saturated hydraulic conductivity
Kbs (and Krs). Six different values of stoniness Rv (0.07, 0.16, 0.24, 0.31,
0.42, and 0.5 cm3 cm−3) and eleven different spatial distributions ofrock fragments were considered. While the ﬁrst simulation assumed aregular distribution of rock fragments, the next ten simulations foreach amount of stoniness considered irregular distributions of rock fragments. Irregular distributions of rock fragment positions were obtainedusing a random numbers pair generator. Generated coordinates of centers of circular rock fragments were adjusted so that rock fragments

3.2. Shapes and orientations of rock fragments (scenarios A and B)
To assess the inﬂuence of rock fragment shapes on the saturated hydraulic conductivity of a virtual stony soil, rock fragments of different
shapes, where 2D cross-sections are circles, ellipses, and triangles,
were regularly distributed in a two-dimensional cross-section of thesoil matrix. Two types of elliptical shapes were considered using the ratios of major and minor axes of 1.56 and 2.56. While the former ellipticalshapes (ellipse No.1; the lengths of main axes were 4 and 6.25 cm) were
only positioned horizontally (i.e., the major axis was horizontal), the latter elliptical shapes (ellipse No.2; the lengths of main axes were 3.125

Table 1
Soil hydraulic parameters for the analytical model of van Genuchten (1980) for three textural classes of the USDA soil textural triangle according to Carsel and Parrish (1988) (taken from the HYDRUS soil catalog).
θs

θr
3

Textural class(cm cmSandy Loam
Loam
Clay Loam

0.065
0.078
0.095

−3

)

α
3

(cm cm
0.41
0.43
0.41

−3

)

(cm

−1

0.075
0.036
0.019

)

n

Kfs

(−)

(cm h−1)

1.89
1.56
1.31

4.42
1.04
0.26

42

H. Hlaváčiková et al. / Geoderma 281 (2016) 39–48

Fig. 2. Examples of regular and irregular distributions of rock fragments in a stony soil with a stoniness Rv of 0.31 cm3cm−3.

were not overlapping or in contact with each other. Examples of regularand irregular distributions of rock fragments are shown in Fig. 2. Sinceno signiﬁcant differences were observed between results of simulations
evaluating the effect of the spatial distribution of rock fragments whendifferent ﬁne soil textures were considered to represent the soil matrix,
only results for sandy loam are reported below.
Results of numerical experiments with respect to effective hydraulicconductivities for corresponding stoniness were evaluated in relativevalues. Mean values (for all scenarios, with different shapes, orientations, and distributions of RFs for a given stoniness) and their standard
deviations (for scenario C), as well as the range of values (min andmax) for scenarios A and B, together with 95% conﬁdence intervals for
scenario C, were estimated. Conﬁdence intervals (±2 standard deviations) were calculated from eleven discrete values for every considered
stoniness, assuming Gaussian distribution of simulated Krs values. Furthermore, the relative effective saturated hydraulic conductivities

simulated for all scenarios were compared with values calculated according to both Ravina and Magier (1984) (Eq. (2); RM), and Corringand Churchill (1961) for spherical (Eq. (3); CCS) and cylindrical (Eq.
(4); CCC) shapes of RFs, divided by a corresponding saturated hydraulicconductivity of the soil matrix. Decreases or increases of simulated
values of relative effective saturated hydraulic conductivities wereexpressed in percentages of calculated values. The effects of shape, orientation, or distribution of rock fragments and stoniness on the relative
saturated hydraulic conductivity were assessed using a linear or nonlinear regression.

4. Results
The results presented in the following sections refer to two-dimensional simulations, unless speciﬁed otherwise.

Fig. 3. Spatial distributions of pressure heads (cm) in soils with rock fragments regularly distributed in the sandy loam matrix and a stoniness Rv of 0.24 cm3cm−3. Water ﬂow is from top tobottom.

H. Hlaváčiková et al. / Geoderma 281 (2016) 39–48

Fig. 4. The spatial distribution of pressure heads (cm) in a soil with rock fragments
irregularly distributed in the sandy loam matrix and a stoniness Rv of 0.24 cm3cm−3.
Water ﬂow is from top to bottom.

4.1. Spatial distributions of pressure heads, water contents, and water ﬂuxes
Spatial distributions of pressure heads, volumetric soil water contents, and water ﬂuxes in soils with rock fragments either regularly or
irregularly distributed in the soil matrix are presented in Figs. 3–9. Spatial distributions of pressure heads (Figs. 3 and 4) show that there areindeed regions immediately above rock fragments with positive pressure heads greater than the boundary pressure head of 1 cm. Immediately below rock fragments, there are regions with negative pressure
heads, indicating unsaturated conditions. This phenomenon was observed in all simulations with RFs. The minimum and maximum differences in pressure heads for regularly distributed rock fragments were43

observed for elliptic rock fragments oriented vertically and horizontally,
respectively. The maximum differences in pressure heads were observed for conditions with an irregular distribution of rock fragments.
The maximum and minimum pressure heads of 14.6 cm and −4.9 cm,
respectively, were observed in the soil with stoniness of
0.24 cm3 cm−3 and the sandy loam soil matrix (Fig. 4).
Almost entire soil systems, with the exception of small regions immediately below the rock fragments, were fully saturated. The extent of theunsaturated regions depended on the shape of rock fragments. While theunsaturated regions were very small, or sometimes even did not develop,
for spherical and elliptical (oriented vertically) rock fragments, larger unsaturated regions developed below elliptical, horizontally oriented rockfragments or triangular rock fragments. Similarly as for pressure heads,
the greatest differences in soil water contents were observed for conditions with an irregular distribution of rock fragments (Fig. 6).
Non–uniform water ﬂuxes (Figs. 7–8) were the consequence of nonuniform pressure head distributions (Figs. 3–4). Water ﬂuxes were
smallest (approaching zero) immediately above and below the rockfragments. They were signiﬁcantly larger (by up to 50%) than the saturated hydraulic conductivity between the rock fragments. The largest
water ﬂuxes were observed for conditions with an irregular distributionof rock fragments (Fig. 8).
Fig. 9 compares the spatial distribution of water ﬂuxes in a soil withregularly distributed rock fragments in a two- and three-dimensional
simulations. A cross-section through the middle of the domain isshown for a three-dimensional simulation. While a two-dimensional
domain contains 8 rows of rock fragments, a cross-section through athree-dimensional domains contains only 4 rows, since the other 4
rows were moved in the y-direction. A similar patterns of water ﬂuxes
is apparent around rock fragments present in the cross-section. Notethat the largest ﬂuxes in a three-dimensional simulation were calculated for regions with ‘missing’ (moved) rock fragments, where ﬂow was
accelerated at the surface of these fragments.

Fig. 5. Spatial distributions of volumetric soil water contents in soils with rock fragments regularly distributed in the sandy loam matrix and a stoniness Rv of 0.24 cm3cm−3. Water ﬂow isfrom top to bottom.

44

H. Hlaváčiková et al. / Geoderma 281 (2016) 39–48

Fig. 6. The spatial distribution of volumetric soil water contents in a soil with rockfragments irregularly distributed in the sandy loam matrix and a stoniness Rv of
0.24 cm3cm−3. Water ﬂow is from top to bottom.

All simulated scenarios indicate that the effective or relative saturated hydraulic conductivities always decrease with increasing stoniness.
The decrease in hydraulic conductivities depends not only on stoniness,but also on the shape, orientation, and distribution of rock fragments inthe soil, as well as on the texture of the soil matrix, as discussed in moredetail below.
4.2. Relative effective saturated hydraulic conductivities
4.2.1. The effect of shapes and orientations of rock fragments (scenarios A
and B)The largest Krs values for all simulated stoniness were obtained forslab-sided ellipses (No. 2) oriented vertically, followed by circles, triangles, and No.1 ellipses. The smallest Krs values were obtained for slab-

Fig. 8. The spatial distribution of water ﬂuxes (cm/min) in a soil with rock fragments
irregularly distributed in the sandy loam matrix and a stoniness Rv of 0.24 cm3cm−3.
Water ﬂow is from top to bottom.

sided ellipses (No. 2) horizontally oriented. These trends in Krs are
displayed in Fig. 10a, together with linear regression equations expressing the dependence of Krs on stoniness Rv.
While mean Krs values for scenarios A and B were almost identical,
their smallest and largest values were different (Table 2). The effect ofsimilarly shaped rock fragments in scenario A on the smallest and largest Krs values was smaller than the effect of the orientation of slab-sided
ellipses in scenario B. The difference between the largest and smallest
Krs values due to different orientations of the same elliptical rock fragments is quite large. For example, for a stoniness of 0.24 cm3 cm−3,
the difference between Krs of a system with elliptical rock fragments either horizontally or vertically oriented was 0.26 (the range of values for
scenario B) (Table 2, Fig. 10a, ellipses denoted as “2-H″ or “2-V″). In

Fig. 7. Spatial distributions of water ﬂuxes (cm/min) in soils with rock fragments regularly distributed in the sandy loam matrix and a stoniness Rv of 0.24 cm3cm−3. Water ﬂow is from top
to bottom.

H. Hlaváčiková et al. / Geoderma 281 (2016) 39–48

45

Fig. 9. Comparison of the spatial distributions of water ﬂuxes (cm/min) in a soil with rock fragments regularly distributed in the sandy loam matrix in two-dimensional (a stoniness Rv of
0.31 cm3cm−3) (left) and three-dimensional (a stoniness Rv of 0.211 cm3cm−3) (right) systems. Water ﬂow is from top to bottom.

contrast, the difference between the smallest and largest Krs values forsoils with different rock fragment shapes in scenario A (circles, triangles,
and ellipses with small axes ratio, denoted in Fig. 10a as “1-H″) for astoniness of 0.24 cm3 cm−3 was only 0.11.
Numerical experiments indicated that the reduction in the relative
effective saturated hydraulic conductivity Krs as a function of stoniness
Rv does not depend on hydraulic properties of the soil matrix. Although
the results displayed in Fig. 10a are only for a sandy loam matrix, almost
identical results were obtained for the other two soil matrix textures.
Numerical experiments further indicated that the relative effective saturated hydraulic conductivity Krs decreases much faster as a function ofstoniness Rv than predicted by the equation of Ravina and Magier(1984) and that this decrease depends on the shape and orientation ofrock fragments. Therefore, we generalized the equation of Ravina andMagier (1984) for different shapes and orientations of rock fragments
as follows:
K rs ¼Ksb
Ks f¼ b−a:Rv

ð5Þ

and the following parameters a and b (Fig. 10a) were ﬁtted using linear
regression for:1. spherical rock fragments with a diameter of 10 cm: a = 1.45, b =0.97 with r2 = 0.98,
2. ellipsoidal rock fragments with a long axis of 4 and 6.25 cm and triangular rock fragments with a height of 12.07 cm and a base of
13 cm: a = 1.75, b = 0.96 with r2 = 0.98,
3. ellipsoidal rock fragments with a long axis of 3.125 and 8 cm vertically oriented: a = 1.23, b = 0.99 with r2 = 0.99, and4. ellipsoidal rock fragments with a long axis of 3.125 and 8 cm horizontally oriented: a = 1.98, b = 0.92 with r2 = 0.93.
Simulated relative effective saturated hydraulic conductivities were
either underestimated or overestimated in comparison with calculated
values. This under- or over-estimation depends on a) the calculation
method, b) the modeled scenario, and c) stoniness. The relative deviations in percentages between simulated and calculated values are presented in Table 3.
Simulated values of Krs were underestimated in all scenarios compared with the Ravina and Magier (1984) function (RM), which also
predicts a smaller decrease of Krs as a function of stoniness than both
equations of Corring and Churchill (1961). With respect of the effectof the shapes of RFs (scenario A), the smallest deviations between
modeled and calculated values were obtained for the Corring and

Fig. 10. The relative saturated hydraulic conductivity (Krs = Kbs /Kfs) as a function of stoniness Rv and the shape of rock fragments (Fig. 10a, symbols V and H denote vertical and horizontal
positioning of ellipses, respectively) and the spatial distribution of circular rock fragments (Fig. 10b). All results were obtained using the two-dimensional simulations, except for “3D
sphere”, which were obtained using the three-dimensional simulations. Relationships denoted as RM, CCS, and CCC were calculated according to Ravina and Magier (1984); Corringand Churchill (1961) for spheres, and Corring and Churchill (1961) for cylinders, respectively. Results are displayed for a sandy loam matrix.

46

Table 2
Mean, smallest and largest Krs values, as well as the range of relative effective saturated hydraulic conductivities of stony soils for different modeling scenarios.
Mean valuesa
Stoniness, Rv

0

0.07

0.16

0.24

Smallest and largest valuesb, /range/
0.31

0.42

0.5

0.07

0.16

0.24

Krs, relative values

0.31

0.42

0.5

Krs, relative values

Regularly dist. RF. Scenario A

1

0.83

0.68

0.56

0.48

–

–

0.81–0.86/0.05/

0.65–0.72/0.07/

0.51–0.62/0.11/

0.45–0.54/0.09/––

Regularly dist. RFs Scenario B

1

0.82

0.67

0.56

0.49

–

–

0.74–0.9/0.16/

0.54–0.8/0.26/

0.43–0.69/0.26/

0.37–0.61/0.24/––

Reg. & irr. dist. RFs Scenario C

1
(0)

0.832
(0.023)0.698(0.030)0.570(0.036)0.460(0.038)0.365(0.019)0.297(0.020)0.785–0.878/0.093/0.638–0.758/0.12/

95% Conﬁdence interval
0.497–0.643
0.384–0.537/0.146//0.153/0.327–0.402/0.075/0.257–0.337/0.08/

Table 3Differences (%) between simulated relative effective saturated hydraulic conductivities of stony soils and their corresponding values calculated from an empirical equation for stoniness from 0.07 to 0.5 cm cm−3. Positive values indicate that simulated
values were larger than calculated values. RM denotes to Ravina and Magier (1984) function, CCS denotes to Corring and Churchill (1961) function (for spheres), CCC denotes to Corring and Churchill (1961) function (for cylinders).
Mean values
Stoniness, Rv

0

0.07

0.16

0.24

0.31

Smallest and largest values
0.42

0.5

0.07

0.16

0.24

0.31

0.42

0.5

Relative differences (in %) between simulated and calculated Krs values according to particular equations (RM, CCS or CCC)

Regularly dist. RFs Scenario A

RM
CCS
CCC

0.0
0.0
0.0

−10.8
−7.6
−4.5

−19.0
−12.6
−6.1

−26.3
−17.5
−8.6

−30.4
−19.7
−8.9

–
–
–

–
–
–

−12.9 (−7.5)−9.9 (−4.3)−6.8 (−1.1)−22.6 (−14.3)−16.4 (−7.4)−10.2 (−0.6)

Min (max)−32.9 (−18.4)−34.8 (−21.7)−24.8 (−8.6)−24.7 (−9.6)−16.8 (1.2)−14.6 (2.5)

–
–
–

–
–
–

Regularly dist. RFs Scenario B

RM
CCS
CCC

0.0
0.0
0.0
0.0

−11.8
−8.7
−5.7

−20.2
−13.9
−7.5

−26.3
−17.5
−8.6

−29.0
−18.0
−7.0

–
–
–

–
–
–

−20.4 (−3.2)−17.6 (0.2)−14.9 (3.5)−35.7 (−4.8)−30.6 (2.9)−25.4 (10.5)

Min (max)a−43.4 (−9.2)−46.4 (−11.6)−36.6 (1.7)−38.1 (2.1)−29.8 (12.6)−29.8 (15.8)

–
–
–

–
–
–

RM
CCS
CCC

0.0
0.0
0.0
0.0

−10.5
−7.4
−4.3

−16.9
−10.3
−3.6

−25.0
−16.0
−7.0

−33.3
−23.0
−12.7

−37.1
−23.9
−10.6

−40.6
−25.8
−10.9

−15.6 (−5.6)−12.6 (−2.3)−9.7 (1.0)−24.0 (−9.8)−18.0 (−2.5)−11.9 (4.7)

95% Conﬁdence interval, min (max)−34.6 (−15.4)−44.3 (−22.2)−26.8 (−5.2)−35.7 (−10.1)−18.9 (4.9)−27.1 (2.0)−43.6 (−30.7)−31.8 (−16.1)−19.9 (−1.6)−48.6 (−32.6)−35.8 (−15.8)−22.9 (1.1)

Rreg. & irr. dist. RFs Scenario Ca

The smallest values correspond to RFs oriented horizontally, the largest values correspond to RFs oriented vertically.

H. Hlaváčiková et al. / Geoderma 281 (2016) 39–48

Scenario A represents regularly distributed RFs of different shapes: circles (diameter 10 cm), triangles (base 13 cm, height 12.07 cm), ellipses (main axes 4 cm and 6.25 cm).
Scenario B represents different orientations of the ellipses (with main axes 3.125 cm and 8 cm); the minimum values are for RFs oriented horizontally, the maximum values are for RFs oriented vertically.
Scenario C represents regularly and irregularly distributed spherical RF (diameter 10 cm).
a
In parentheses are standard deviations.
b
Smallest and largest values in scenarios A and B, in scenario C they correspond to 95% conﬁdence intervals.

H. Hlaváčiková et al. / Geoderma 281 (2016) 39–48

Churchill (1961) equation for cylindrical RFs (CCC). Compared with the
RM method, the Krs values for soils with spherical rock fragments were
underestimated by 7–21% in the stoniness range of 0.07–
0.31 cm3 cm−3. Whereas the Krs values were only underestimated by
4–9% compared to the Corring and Churchill (1961) equation for spherical RFs (CCS). Compared to the Corring and Churchill (1961) equation
for cylindrical RFs (CCC), the Krs values were either underestimated by
1.1–0.6% in the range of stoniness of 0.07–0.16 cm3 cm− 3, oroverestimated by 1.2–2% in the range of stoniness of 0.24–
0.31 cm3 cm−3.
The good agreement between Krs values predicted by the Corringand Churchill (1961) equation for cylindrical RFs (CCC) and the two-dimensional simulations with spherical rock fragments is not surprising.
While the goal of these simulations was to evaluate the effects of spherical rock fragments, the two-dimensional simulations can only consider
a circular projections of these fragments and neglect the effects of thethird direction. These calculations thus indeed simulate the effects of cylindrical, rather than spherical, rock fragments. On the other hand, theresults of the three-dimensional simulations, which could fully accountfor spherical shapes of the rock fragments, corresponded closely withvalues predicted by the Corring and Churchill (1961) equation for
spherical RFs (CCS) (Fig. 10a). The good agreement between the resultsof both two- and three-dimensional simulations with cylindrical andspherical rock fragments, respectively, with the Corring and Churchill(1961) equations conﬁrms the robustness of these equations, as wellas the correctness of our numerical simulations.
The Krs values of stony soils with slab-sided elliptical RFs were
underestimated for the entire range of stoniness of 0.07–
0.31 cm3 cm−3 compared with all calculation methods when RFs werehorizontally oriented. These values represented the minimum Krs values
for scenario B. When RFs were vertically oriented, the Krs values
corresponded to the maximum Krs values for scenario B. In this case,
the Krs values were only underestimated when compared to the RM
equation (by 3–11%), and overestimated when compared to the CCS
or CCC methods (by 0.2–2.1% and 3.5–15.8%, respectively).
4.2.2. The effect of regular and irregular distributions of rock fragments(scenario C)
The largest Krs values for scenarios with irregularly distributed RFs
were for all amounts of stoniness close to the Krs values when rock fragments were regularly distributed. The mean Krs values for scenario C in
the range of stoniness of 0.07–0.31 cm3 cm−3 were almost the same as
the mean values for scenarios A and B. The conﬁdence intervals show an
overlap between the minimum Krs values for a particular stoniness andthe maximum Krs values for the next higher stoniness.
The mean and largest and smallest Krs values for scenario C wereﬁtted using non-linear regression polynomial functions (Fig. 10b). The
largest Krs values were similar to those for scenario A for regularly distributed spherical shaped RFs (Table 2). The difference between thelargest and smallest Krs values due to regular and irregular distributionsof RFs varied from 0.075 to 0.153 in the range of stoniness of 0.07–
0.5 cm3 cm−3, with the largest difference identiﬁed at a stoniness of
0.31 cm3 cm−3. Table 2 summarizes the mean Krs values, their standard
deviations, and 95% conﬁdence intervals (the mean value ±2 times thestandard deviation) for different classes of stoniness.
The mean and smallest Krs values were always underestimated inthis scenario compared to all calculation methods (Table 3). The largest
Krs values were underestimated compared to the RM and CCS methods.
Compared to the Corring and Churchill (1961) equation for cylindricalRFs (CCC), the largest Krs values were almost always overestimated(by 1–4.9%) in the range of stoniness of 0.07–0.5 cm3 cm−3, except fora stoniness of 0.42 cm3 cm−3 when a small underestimation (by 1.6%)
was found.
Results in Figs. 10a and b are shown only for stoniness smaller than
0.31 or 0.5 cm3 cm−3, respectively, since it was difﬁcult to regularly or
irregularly distribute rock fragments of all shapes for a larger stoniness.

47

The presence of rock fragments decreased the cross-sectional area of thesoil through which water ﬂows and correspondingly, the relative saturated hydraulic conductivity, as expressed by Eq. (2) of Ravina andMagier (1984). The shape, orientation and spatial distribution of rockfragments provided an additional inﬂuence on the curvatures of waterﬂow paths and consequently, the resistance to water ﬂow through thesoil.
5. Discussion
Due to many factors such as the required size of the REV for stonysoils, as discussed by Novák et al. (2011), it is difﬁcult to compare values
calculated using numerical experiments with measured values. Empirical evidence for the effect of the size of rock fragments on stony soil saturated hydraulic conductivity was presented by Novák et al. (2011). For
example, Novák et al. (2011) showed that large rock fragments with adiameter of 20 cm decreased the relative effective saturated hydraulicconductivity more than smaller rock fragments with diameters of 10
or 5 cm with the same total volume. Similarly, Ma and Shao (2008),
for RFs with diameters in the range of 5–40 mm, showed a more significant decrease in the inﬁltration rate due to larger rock fragments thandue to a larger number of smaller stones of the same volume. They
also pointed out that rock fragments of a spherical shape decrease theinﬁltration rate less than rock fragments of cylindrical or prismatic
shapes. Ma and Shao (2008) explained this phenomenon by different
resistances to water ﬂow for differently shaped rock fragments.
Using ﬁeld measurements, Verbist et al. (2009) observed both a decrease and an increase in Kbs with stoniness. Stoniness was expressed intheir study using different measures such as a relative volume, a relative
mass, or a mean weighted diameter of rock fragments. Increasing or decreasing trends, with moderate (correlation coefﬁcients equal to 0.51–
0.63) and weak (correlation coefﬁcients equal to 0.14–0.19) correlations, respectively, were the result of two opposing phenomena: 1)
rock fragments causing greater curvature of ﬂow paths and a decreasein Kbs , and 2) the presence of interconnected lacunar pores or othermacropores causing an increase in Kbs . Verbist et al. (2013) further
showed that the saturated hydraulic conductivity of a stony soil determined in the ﬁeld can be highly variable, not only because of soil heterogeneity, but also due to the measurement method, the number of
replications, and the calculation method. Furthermore, it is not clear
how much different measurement results are affected by the stony
character of the soil or different measurement methods.
Since the possible presence of lacunar pores (or macropores) was
not accounted for in our numerical experiments, calculated Kbs values
are decreasing with increasing stoniness. Calculated values of Kbs are always lower than those calculated according to Ravina and Magier(1984). One reason for this may be that Eq. (1) or (2) of Ravina andMagier (1984) do not consider the inﬂuence of additional resistances
caused by different shapes, orientations, and positions of rockfragments.
Our numerical simulations that evaluated irregular distributions ofrock fragments did not consider situations when rock fragments werein contact with each other. In such a situation, not only the inﬂuenceof the position of rock fragments on Kbs would be evaluated, but also
their size effect, since such rock fragments would behave as larger
stones, which had already been evaluated in an earlier study (Novák
et al., 2011). In addition, situations when all rock fragments are concentrated in one part of the ﬂow domain were not considered.
Shapes and positions of rock fragments can have a signiﬁcant inﬂuence on the hydrophysical properties of stony soils. This can be especially manifested for saturated conditions. The distribution and orientationof rock fragments in mountainous areas are strongly affected by the process of pedogenesis, slope gradient, and topographic position (Chen etal., 2011, 2012). Rock fragments can be oriented parallel with theslope. Consequently, such soils can have different hydraulic conductivities in horizontal and vertical directions, producing macroscopic

48

H. Hlaváčiková et al. / Geoderma 281 (2016) 39–48

anisotropy. This anisotropic behavior can, for example, accelerate hypodermic ﬂow at the hillslope scale and speed up the formation of subsurface ﬂow (Capuliak et al., 2010).
Water ﬂow in stony soils is a complex, three-dimensional phenomenon. Therefore, since majority of numerical experiments, the results ofwhich were discussed above, have been carried out using a two-dimensional model, the results are only approximate. Nevertheless, we anticipate, as was partially conﬁrmed by a limited number of threedimensional simulations reported above, that differences between ourestimated hydraulic conductivities and those that would be estimatedusing a fully three-dimensional model would be small, and that the general relationships and various trends discussed above would bepreserved.
6. Conclusions
Numerical experiments showed that the effective hydraulic conductivity of saturated stony soils Kbs decreased as stoniness increased. However, the Kbs estimates depended not only on stoniness (Rv) as previously
shown by others such as Ravina and Magier (1984), but they were also
sensitive to other factors such as the shape, distribution, and orientationof the rock fragments in stony soils. For example, the largest (0.9, 0.8,
0.69, 0.61) and smallest (0.74, 0.54, 0.43, 0.37) values of the relative saturated hydraulic conductivities Krs for stoniness of 0.07, 0.16, 0.24, and0.31 cm3 cm−3 were estimated for rock fragments of elliptical shapethat were either vertically or horizontally oriented, respectively. The
Ravina and Magier (1984) function (RM) unpredicted a decrease of Krs
as a function of stoniness compared to the results of numerical simulations for all shapes, orientations, and distributions of rock fragments, aswell as compared to both equations of Corring and Churchill (1961).
Therefore, we have generalized the equation of Ravina and Magier(1984) to account for the shape and orientation of rock fragments andidentiﬁed the values of additionally required parameters a and b.
The results of the two-dimensional simulations with circular inclusions agreed well with the Krs values predicted by the Corring and
Churchill (1961) equation for cylindrical RFs (CCC). Similarly, the resultsof the three-dimensional simulations with spherical shapes of the rockfragments corresponded closely with values predicted by the Corringand Churchill (1961) equation for spherical RFs (CCS). The good agreement between the results of both two- and three-dimensional simulations with cylindrical and spherical rock fragments with the Corringand Churchill (1961) equations conﬁrms the robustness of these equations in describing the effects of rock fragments of speciﬁc shapes on thehydraulic conductivity of stony soils, as well as the correctness of ournumerical simulations.
The relative effective saturated hydraulic conductivities Krs also
depended on whether the rock fragments were regularly or irregularly
distributed in a stony soil. The largest difference of 0.153 in Krs estimates
due to irregular distribution of rock fragments in the stony soil was obtained for a stoniness of 0.31 cm3 cm−3. Differences in Krs estimates due
to irregular distributions of rock fragments decreased for both large(0.5 cm3 cm−3) and small (0.07 cm3 cm−3) stoniness.
Finally, results of steady-state water ﬂow simulations indicated thatthe relative saturated hydraulic conductivities Krs do not depend on thehydraulic properties of the soil matrix of stony soils for a particularshape and orientation of rock fragments.
Acknowledgements
We are grateful to two anonymous reviewers for their careful reviews of our manuscript, all their editorial suggestions and requests

for clariﬁcations, and for their guidance in helping us to signiﬁcantly improve our manuscript.
This contribution was partially supported by a grant from the Slovak
Academy of Sciences VEGA (project No. 2/0055/15, No. 2/0013/15), andis the result of the project implementation ITMS 26240120004 from theCentre of Excellence for Integrated Flood Protection of Land, as supported by the Research & Development Operational Program funded by theEuropean Regional Development Fund.
References
Al-Qinna, M., Scott, H.D., Brye, K.R., Brahana, J.V., Sauer, T.J., Sharpley, A., 2014. Coarse fragments affect soil properties in a mantled-karst landscape of the Ozark Highlands. Soil
Sci. 179, 42–50.
Bear, J., 1972. Dynamics of Fluids in Porous Media. American Elsevier, New York.
Bouwer, H., Rice, R.C., 1984. Hydraulic properties of stony vadose zones. Ground Water
22, 696–705.
Brakensiek, D.L., Rawls, W.J., Stephenson, G.R., 1986. Determining the saturated hydraulicconductivity of a soil containing rock fragments. Soil Sci. Soc. Am. J. 50, 834–835.
Brouwer, J., Anderson, H., 2000. Water holding capacity of ironstone gravel in typic
Phlintoxeralf in Southeast Australia. Soil Sci. Soc. Am. J. 64, 1603–1608.
Buchter, B., Hinz, C., Flühler, H., 1994. Sample size for determination of coarse fragment
content in a stony soil. Geoderma 63, 265–275.
Capuliak, J., Pichler, V., Flühler, H., Pichlerová, M., Homolák, M., 2010. Beech forest density
control on the dominant water ﬂow types in andic soils. Vadose Zone J. 9, 755–769.
Carsel, R.F., Parrish, R.S., 1988. Developing joint probability distributions of soil water retention characteristics. Water Resour. Res. 24, 755–769.
Chen, H., Liu, J., Wang, K., Zhang, W., 2011. Spatial distribution of rock fragments on steep
hillslopes in karst region of northwest Guangxi, China. Catena 84, 21–28.
Chen, H., Liu, J., Zhang, W., Wang, K., 2012. Soil hydraulic properties on the steep karst
hillslopes in northwest Guangxi, China. Environ. Earth Sci. 66, 371–379.
Childs, S.W., Flint, A.L., 1990. Physical properties of forest soils containing rock fragments.
In: Gessel, S.P., Lacate, D.S., Weetman, G.F., Powers, R.F. (Eds.), Sustained Productivityof Forest Soils, Proceedings of the 7th North American Forest Soils Conference. University of British Columbia, Faculty of Forestry Publication, Vancouver, pp. 95–121.
Corring, R.L., Churchill, S.W., 1961. Chem. Eng. Prog. 57, 53–59.
Cousin, I., Nicollaud, B., Coutadeur, C., 2003. Inﬂuence of rock fragments on the water retention sand and water percolation in a calcareous soil. Catena 53, 97–114.
Fiès, J.C., De Louvigny, N., Chanzy, A., 2002. The role of stones in soil water retention. Eur.
J. Soil Sci. 53, 95–404.
Hlaváčiková, H., Novák, V., Holko, L., 2015. On the role of rock fragments and initial soil
water content in the potential subsurface runoff formation. J. Hydrol. Hydromech.
63, 71–91.
Hraško, J., Bedrna, Z., 1988. Applied Soil Science. Priroda Publ. House, Bratislava (inSlovak).
Ma, D.H., Shao, M., 2008. Simulating inﬁltration into stony soils with a dual-porosity
model. Eur. J. Soil Sci. 59, 950–959.
Ma, D.H., Shao, M.A., Zhang, J.B., Wang, Q.J., 2010. Validation of an analytical method for
determining soil hydraulic properties of stony soils using experimental data.
Geoderma 159, 262–269.
Novák, V., Kňava, K., Šimůnek, J., 2011. Determining the inﬂuence of stones on hydraulicconductivity of saturated soils using numerical method. Geoderma 161, 177–181.
Peck, A.J., Watson, J.D., 1979. Hydraulic Conductivity and Flow in Non-Uniform Soil.
Worshop on Soil Physics and Soil Heterogeneity. CSIRO Division of Environmental
Mechanics, Canberra, Australia.
Poesen, J., Lavee, H., 1994. Rock fragments in top soils: signiﬁcance and processes. Catena
23, 1–28.
Ravina, I., Magier, J., 1984. Hydraulic conductivity and water retention of clay soils containing coarse fragments. Soil Sci. Soc. Am. J. 48, 736–740.
Šály, R., 1978. Soil - Basic Component of Forest Production. Priroda Publ. House, Bratislava(in Slovak).
Sauer, T.J., Logsdon, S.D., 2002. Hydraulic and physical properties of stony soils in a small
watershed. Soil Sci. Soc. Am. J. 66, 1947–1956.
Shi, Z., Wang, Y., Yu, P., Xu, L., Xiong, W., Guo, H., 2008. Effect of rock fragments on the
percolation and evaporation of forest soil in Liupan Mountains, China. Acta Ecol.
Sin. 28, 6090–6098.
Šimůnek, J., van Genuchten, M.T., Šejna, M., 2008. Development and applications of the
HYDRUS and STANMOD software packages and related codes. Vadose Zone J. 7,
587–600. http://dx.doi.org/10.2136/vzj2007.0077.
van Genuchten, M.T., 1980. A closed – form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. A. J. 44, 987–996.
Verbist, K., Baetens, J., Cornelis, W.M., Gabriels, D., Torres, C., Soto, G., 2009. Hydraulic conductivity as inﬂuenced by stoniness in degraded drylands of Chile. Soil Sci. Soc. Am. J.
73, 471–484.
Verbist, K.M.J., Cornelis, W.M., Torfs, S., Gabriels, D., 2013. Comparing methods to determine hydraulic conductivities on stony soils. Soil Sci. Soc. Am. J. 77, 25–42.
Zhou, B., Shao, M., Shao, H., 2009. Effects of rock fragments on water movement and solute transport in a loess plateau soil. C.R. Geoscience 341, 462–472.