Open access peer-reviewed chapter

Unsaturated Hydraulic Conductivity of Fractal-Textured Soils

By Yongfu Xu

Submitted: March 27th 2012Reviewed: June 3rd 2013Published: December 4th 2013

DOI: 10.5772/56716

Downloaded: 1494

1. Introduction

The increasing concern with groundwater pollution and contamination of soils has stimulated the development of numerous mathematical models of pollutant transport in soils. The most important approaches to model transient water and solute transport in the vadose zone are based on the Richards equation. To solve this equation, the knowledge of the soil hydraulic properties, namely, the soil-water characteristic curve (SWCC) and the unsaturated hydraulic conductivity is required. The laboratory measurements show that the value of unsaturated hydraulic conductivity varies considerably from soil to soil with different water content (Khaleel and Relyea, 1995). Indeed, it is found that the unsaturated hydraulic conductivity decreases by one to three orders of magnitude across a small pressure head range even near saturation (0~10cm pressure head), due to the effects of structural macropores (Jarvis and Messing, 1995). Of all hydraulic properties, the unsaturated hydraulic conductivity is most difficult to measure. Therefore the use of indirect methods has become more and more common to estimate the unsaturated hydraulic conductivity from more easily measured soil properties (van Genuchten et al., 1992).

Modern hydrological models require information on hydraulic conductivity and soil-water retention characteristics. All hydraulic properties, the soil-water characteristics, hydraulic conductivity and soil-water diffusivity (SWD) are closely related to the geometry of a porous media (Brooks and Corey, 1966; Burdine, 1953). Measurements of hydraulic properties are expensive, time-consuming and highly variable (Dirksen, 1991). Prediction of these properties is a viable alternative, especially when the predictive model contains a few parameters sensitive to structural conditions. Porous media (e.g. soils, rocks, etc.) are heterogeneous systems composed of numerous, different and interacting components (van Damme, 1995). The complex nature of these porous media complicates any prediction of their hydraulic properties.

A potentially powerful method results from the pore surface models, in which the soil–water characteristic curve of unsaturated soils is interpreted as statistical measure of its equivalent pore-size distribution (PSD) (van Genuchten, 1980; Corey, 1992). The frequency of different pore radii is related to matric suction by the soil-water characteristic curve (SWCC) and the relation between matric suction and pore radius is described by the Young-Laplace equation. The relative hydraulic conductivity (RHC) of unsaturated soils can be deduced from the soil-water characteristic curve (SWCC) through simplifying assumptions on pore topology and using the Hagen-Poiseuille law as an approximation for water flow. An additional empirical parameter is introduced, which includes all uncertainties. This empirical parameter is often referred to as “tortuosity”, but its physical meaning is unclear (Vogel and Roth, 1998). The choice of the analytical model for the soil-water characteristic curve (SWCC) can significantly affect the predicted function of the unsaturated hydraulic conductivity.

Fractals describe hierarchical systems and are suitable to model the heterogeneous soil structure with tortuous pore space (Rieu and Sposito,1991; Xu and Sun, 2002). Toledo et al. (1990) modeled the soil-water characteristic curve (SWCC) and unsaturated hydraulic conductivity using fractal geometry and thin-film physics. Tyler and Wheatcraft (1990) gave the unsaturated hydraulic conductivity functions based on the fractal model for the soil-water characteristic curve (SWCC) and the relative conductive models developed by Mualem (1976) and Burdine (1953). Crawford (1994) studied the influence of heterogeneity of both the solid matrix and the pore space, as well as the shape of the pore boundary, on the saturated and unsaturated hydraulic conductivity. Fuentes et al. (1996) derived an expression for unsaturated hydraulic conductivity using the fractal dimension obtained from the soil-water characteristic curve (SWCC). Hunt and Gee (2002) applied critical path analysis from percolation theory to calculate the unsaturated hydraulic conductivity of soils with fractal pore space. Xu (2004), Xu and Dong (2004) Xu et al. (2004) derived the unsaturated hydraulic conductivity using the fractal model for the pore surface and gave a simple method to determine the fractal dimension. Models of the unsaturated hydraulic conductivity incorporate fractal dimension characterizing scaling of different properties including parameters representing connectivity and tortuosity. Thus, it is encouraged to derive the functions of the soil–water characteristic curve and unsaturated hydraulic conductivity from the fractal model for the pore surface. In this chapter, the soil-water characteristic curve (SWCC) and relative hydraulic conductivity (RHC) function were derived and expressed by the effective degree of saturation based on the fractal model for the pore surface. The proposed soil–water characteristic curve and unsaturated hydraulic conductivity function were examined in detail by the published experimental data. Comparisons between the prediction using both the fractal model and the van Genuchten-Mualem (G-M) model and the measurement of the relative hydraulic conductivity (RHC) were conducted. Using the fractal model for the soil-water characteristic curve (SWCC) and relatively hydraulic conductivity (RHC), one-dimensional rainfall infiltration and slope stability due to rainfall infiltration were analysized in chapter.

2. Fractals

Some of the most common and useful examples of fractals are: the Koch curve, the Sierpinski gasket and carpet, and the Menger sponge. All of these fractals enjoy the property of self-similarity. Roughly speaking, a subset of nis said to be self-similar if it is a union of a number of smaller similar copies of itself. Some of the essential notions of the theory of fractals as it applies to self-similar sets should be briefly reviewed. A previous acquaintance with at least one of the sets mentioned above is useful but not necessary. Given an arbitrary subset B of n, a cover of B is a family U of sets U0nsuch that:

BnUαE1

If the family U is countable (or finite) the cover is said to be countable (or finite). It is customary in that case to indicate the members of the family with a Latin subscript, namely, Ui, where i ranges over the natural numbers. The diameter of a subset of nis the least upper bound (i.e., the supremum) of the distance between pairs of points in the subset. By convention, the diameter of the empty set is zero. A δ-cover of Bnis a (countable) cover such that, for every i, diam(Ui)≤δ, where δ is a positive real number and where diam(.) is the diameter function on subsets of n. The members of a δ-cover are indicated by Uiδ. The s-dimensional Hausdorff measure Hs(B) of B is defined as:

Hs(B)=lims0infi=1n[diam(Uis]sE2

where the inf extends over all δ-covers. It can be shown that this definition indeed provides an outer measure for all non-negative values of s and for any subset of n. The Hausdorff dimension of B is defined as:

dimHB=inf{s0:Hs(B)=0}E3

The Hausdorff dimension of a subset of ncannot exceed n. If dimH B > 0, then it can be shown that for all values of s strictly smaller (respectively, larger) than the dimension, the s-dimensional Hausdorff measure of B is infinite (respectively, zero). In this sense, the Hausdorff dimension represents a critical value of discontinuity of the s-dimensional Hausdorff measure, thus making Eq. (3) meaningful. Of particular interest is the value of the s-dimensional Hausdorff measure for s equal to the Hausdorff dimension of the set. This is usually called the Hausdorff measure of the set. The Hausdorff measure of a set may turn out to have any value, including zero or infinite.

An important feature of the s-dimensional Hausdorff measure is the scaling property, namely, the way it changes under similarity transformations of n. Recall that a transformation S: nnis called a similarity if there exists a real scale factor λ>0 such that, for all x, y n, the following equation is satisfied:

|S(y)S(x)|=λ|yx|E4

where |.| denotes the Euclidean distance in n. It is not difficult to prove that for all subsets Bnand for all values of s, under a similarity transformation S with scale factor λ the s-dimensional Hausdorff measure transforms according to the formula:

Hs(S(B))=λsHs(B)E5

Another important property of the s-dimensional Hausdorff measure is that it is preserved under Euclidean isometries (translations, rotations, reflections). In other words, congruent sets have the same s-dimensional Hausdorff measure. The behavior of the Hausdorff measure under general affine transformations, on the other hand, cannot be captured under the umbrella of a simple formula.

The properties just described of the s-dimensional Hausdorff measure can be used to obtain a straightforward evaluation of the Hausdorff dimension of self-similar fractals. Indeed, let B be the union of m copies of λ-scaled mutually congruent copies of itself (with λ<1). If these copies are disjoint Borel sets, we have:

Hs(B)=mλsHs(B)E6

by virtue of the scaling property (5). Assume now that the Hausdorff measure of B (namely, the value of Hs(B) for s=dimH(B)) is finite and positive. In that case, it follows Eq. (6) that:

dimHB=logmlogλE7

The fractal dimension for a pore surface can be defined in the following way. Imagine the pores being enclosed by a set of spheres of radius r, the number of spheres N necessary to do this is clearly a function of radius of the spheres. The definition of the fractal dimension D is:

D=limr0[ln(N(A,r))ln(1/r)]E8

Equation (8) is often used to determine the fractal dimension by experiments. A typical fractal set, the middle third Cantor set may be constructed from a unit interval by a sequence of deletion operations (Fig. 1). Let E0 be the interval [0,1], E1 is the set obtained by deleting the middle third of E0 so E1 consists of the two intervals [0,13], [23,1]. Deleting the middle third of E1 gives E2, thus E2 comprises the four intervals [0,19], [29,39], [69,79], [89,99]. Proceeding in this like manner, Ei is obtained by deleting the middle third of each interval in Ei-1 and Ei consists of 2i intervals with length of 3-i. Using Eq. (8), the fractal dimension of the middle-third Cantor set is D=limi{log(2i)log(1/(3i))}=0.63. The Cantor set is often used to model the dust distribution.

Figure 1.

Middle-third Cantor set with fractal dimension D=0.63

3. Fractal model for the pore surface

Many research results show that the soil pore surface is fractal (Avnir and Jaroniec 1989; Xu and Sun 2002). To investigate the impacts of fractal scaling upon hydraulic properties of porous media, a fractal representation of a porous media is developed. A cube of size 1 by 1 is formed by initially subdividing the cube into 1/a3 subcubes, each of size a by a. From the original cube, N subcubes are removed and represent a pore of size a by a. The remaining 1/a3-N subcubes are then each divided into 1/a3 subcubes and N subcubes are removed from each of the original subcubes. Such a recursion algorithm results in a cube everywhere filled with pores of all sizes, with a predominance of small pores. The number of size of a3 needed to cover the pores equal to or larger than a3 is given by N. The Menger sponge (Fig. 2) is often used to model porous materials, such as soils. For the Menger sponge, a=1/3, N=7, and the fractal dimension is given by D=limi{log(20i)log(1/(3i))}=2.73.

Figure 2.

The Menger sponge with fractal dimension D=2.73

When we covered a soil pore surface using balls with the same radius, the soil pore surface is covered by N1 balls with radius r1, the same surface is covered by Ni balls with radius ri, and it is covered by Nj balls with radius rj, and usually r1>ri>rj, it is found that the surface area of the soil pore surface contents has the following relationship: N0r02<Niri2<Njrj2. If the soil pore surface is a self-similar surface, i.e., the pore surface is a fractal surface, the relationship between the number of covered balls and its radius can be written as (Mandelbrot 1982)

N(r)=CrDE9

where C is a constant, r is the size corresponding to the pore radius, D is the fractal dimension of the pore surface. All materials have a surface fractal dimension in the whole range of physically meaningful values, i.e., from 2 to 3. The limit D=2 would correspond to a perfectly regular, smooth (Euclidean) surface, whereas D=3 would correspond to a self-similar surface so intricate that it would be space filling.

The fractal dimension of the pore surface can be determined using mercury intrusion porosimetry and adsorption isotherm. These two methods are introduced in detail as follows.

3.1. Mercury intrusion porosimetry

The d-measure is obtained from Eq. (9), and is expressed as follows:

M(r)=N(r)rd=CrdDE10

For the fractal pore surface, the relationship between the surface area with the radius less than r and the pore radius r can be obtained from Eq. (10). The surface area can written as

Ap(r)=Cr2Da)Vp(r)=Cr3Db)E11

where in 11 a) Ap(≤r) is the surface area of the soil pores with the radius smaller than r. Similarly, the total volume of the pores with the radius less than r is given by 11b), where in Vp(≤r) is the cumulative volume of the soil pores with the radius smaller than r.

Neimark (1992) gave a method to measure the surface fractal dimension of the soil pores using mercury intrusion porosimetry. The soil pore surface can be approximated by the equilibrium interface between mercury and soil particles in the close vicinity of the pore surface. According to the Young–Laplace equation, the mean radius of the equilibrium interface is expressed as follows:

r(p)=2σcosαpE12

where r(p) is the mean radius of the equilibrium interface under pressure p, σ is surface tension. From the thermodynamic viewpoint, the equilibrium interface area can be calculated from the balance between the work of formation of the equilibrium interface and the work of mercury intrusion (Neimark 1992), i.e.,

Ap(p)=0V(p)V(p)TdpE13

where Ap(p) is the equilibrium interface area under pressure p. The surface fractal dimension can be determined theoretically from Eq. (11a). If the slope of the straight line is λ in log r(p) versus log Ap(p), the surface fractal dimension of soil pores is given by

D=2λE14

The pore-size distribution (PSD) is usually obtained using mercury intrusion porosimetry. According to Eq.(11b), the fractal dimension of the pore surface can be obtained from the slope of the regressed linear relation in the plane of logr vs. logVp. If the slope of the fitting line in the logr-logVp plane is κ, the fractal dimension of the pore surface is written as

D=3κE15

The value of fractal dimension spans a large range from 1.0 to 3.0 (Gimenez et al., 1997). Larger fractal dimension is associated with clayey soils (van Damme, 1995).

The accumulative volume of the pores can be measured by many methods. Mercury intrusion porosimetry provides a useful and potentially valuable measurement of the pore volume for the porous medium. The ranges of equivalent pore diameter explored cover almost five orders of magnitude, from several hundred microns down to approximately 16Å (Watabe et al., 2000).

Figure 3 shows the soil pore-size distribution (PSD) obtained from mercury intrusion porosimetry (Watabe et al., 2000). Symbols S-02, S-03 and S-04 are the serial numbers of soil samples, and Vv is the total pore volume in Fig. 3. It is seen from Fig. 3 that the surface of the soil pores can be described by fractal model, and the relationship between Vp(≤r) and r can be expressed by a linear function in log-log plot. The slopes of the regressed linear relation in the plane of logr vs. log(Vp(≤r)/Vv) are 0.34, 0.37 and 0.49, and therefore the fractal dimensions are 2.66, 2.63 and 2.51 for specimens of S-02, S-03 and S-04, respectively. The maximum radius is the radius at which the pore volume reaches the maximum value, and is defined as the intersection between the fitting line and the line of Vp(≤r)/Vv=100%. The maximum radii of soil pores are 0.0125, 0.06 and 0.006 mm for specimens of S-02, S-03 and S-04, respectively. The relationship between the maximum radius R and the air-entry value ψe can be expressed by the Young-Laplace equation, i.e. ψe=2σcosα/R. The parameters α=0 and σ=0.075kPa mm (Watabe et al., 2000), and the air-entry values are 12kPa, 2.5kPa and 25kPa corresponding to the maximum pore radius R of 0.0125, 0.06 and 0.006 mm for specimens of S-02, S-03 and S-04, respectively. The parameters obtained from the fractal model of the pore surface are listed in Table 1. In Table 1, parameters R and κ are obtained from the fractal model of the pore surface in Fig. 3.

Soil typeκDδR (mm)ψe (kPa)Data source
S-020.342.66-0.340.012512Watabe et al., 2000
S-030.372.63-0.370.062.5
S-040.492.51-0.490.00625
A horizon0.222.78-0.220.20.75Vogel and Roth, 1998
B horizon0.132.87-0.130.250.6
Toyoura sand1.331.67-1.330.091.67Uno et al., 1998

Table 1.

Parameters obtained from the fractal model of PSD

Figure 3.

Fractal dimension and the maximum pore radius obtained from PSD for glacial tills (Data from Watabe et al., 2000)

Figure 4.

Fig. 4 Fractal dimension and the maximum radius of the PSD of two silty soils (Data from Vogel and Roth, 1998)

Experimental data of the pore-size distribution (PSD) for agricultural silty soils from two horizons, A and B were also taken from Vogel and Roth (1998). The pore-size distribution (PSD) was determined using the three-dimensional reconstructions of the pore space. The pore space was eroded by a spherical structuring element with a given diameter 2r and subsequently dilated by the same structuring element (Vogel and Roth, 1998). In the resulting image all pores smaller than 2r were removed. The cumulative volume Vp(>r)/VT was obtained using the application of successively larger structuring elements by Vogel and Roth (1998). The following relationship is obtained between Vp(>r)/VT and Vp(≤r)/VT.

V(>r)VT+V(r)VT=VvVT=θsE16

where VT is the total volume of soil, Vv is the total volume of void, θs is the saturated volumetric water content. The saturated volumetric water contents are 0.502 and 0.448 for the soils at A horizon and B horizon, respectively. Hence the cumulative volume V(≤r)/VT can be obtained from the cumulative volume V(>r)/VT. The cumulative volume of V(≤r)/VT can be translated into V(≤r)/Vv according to Eq.(16). The cumulative volume of V(≤r)/Vv is shown in Fig. 4, where the dots denote the results transformed from the experimental data given by Vogel and Roth (1998) and the solid lines denote the results of the fractal model. The relationship between the pore volume and the pore radius satisfies an approximately linear expression in log-log plot. The slopes (κ) of the regression line are 0.24 and 0.13 for the soils at A and B horizons, respectively. According to Eq. (15), the fractal dimensions of the pore surface are 2.76 and 2.87 for A and B horizons, respectively. The maximum radius is defined as the value at which the cumulative volume of pores reaches the maximum, and V(≤r)/Vv=100%. The maximum radii of the soil pores at the A and B horizons are 0.2cm and 0.25cm, respectively. The air-entry values can be calculated using the Young-Laplace equation, and they are 0.75kPa and 0.6kPa for the A and B horizons, respectively. The fractal dimension of the soil pore and the air-entry value of the soils at the A and B horizons are listed in Table 1. It can be seen from Figs. 3-4 that the pore surface can be modeled by fractal theory, and the pore-size distribution (PSD) satisfies with Eq. (11b). The soil structure can be described by the fractal dimension. The larger the fractal dimension, the more tortuous is the pore space. The fractal dimension is varied with the soil structure, and is different for different soils.

Uno et al. (1998) studied the relationship between the pore-size distribution (PSD) and unsaturated hydraulic properties of Toyoura sand. The pore-size distribution (PSD) of Toyoura sand was measured by the air intrusion method, and is shown in Fig. 5. The relationship between the pore volume and pore radius is approximated by a linear function in log-log plot, and satisfies with Eq. (11b). Hence the pore-size distribution (PSD) of Toyoura sand can be described by fractal model, and its fractal dimension is 1.67 obtained from Fig. 5 not increase with radius, i.e., Vp/Vv=1. The maximum pore radius of Toyoura sand is nearly 0.09 mm in Fig. 5. The surface tension (T) of water is 0.075kPa mm. The air-entry value can be calculated from Young-Laplace equation, and is 1.67kPa on the assumption that α=0. The fractal dimension of the soil pore and the air-entry value of Toyoura sand are listed in Table 1.

Figure 5.

The fractal dimension and the maximum radius evaluated from the PSD of Toyoura sand. (Data from Uno et al., 1998.)

Figure 6.

The surface fractal dimension of Tsukinuno bentonite from nitrogen adsorption

3.2. Adsorption isotherm

The fractal dimension of the pore surface serves to characterize the pore surface roughness or irregularities. The magnitude of the fractal dimensions is relevant to many important physicochemical processes namely adsorption, surface diffusion and catalysis. Avnir and Jaroniec (1989) proposed a convenient method to determine the fractal dimension from the nitrogen adsorption. Similarly with the adsorption isotherm equation, the correlation of the water volume absorbed by clay to the vapour pressure is given by

VwVm=k[ln(P0P)]D3E17

where Vw is the water volume absorbed by clay, Vm is the volume of montmorillonite, P is the partial water vapour pressure in equilibrium with clay at some water content w and some temperature T and P0 is the equilibrium water vapour pressure of pure water at temperature T. The surface fractal dimension of Tsukinuno bentonite obtained from nitrogen adsorption is shown in Fig. 6. The surface fractal dimensions of bentonite powder and compacted bentonite are 2.66 and 2.65, respectively.

The adsorption isotherm also allows the calculation of swelling pressure of clay as a function of the water content (Kahr et al., 1990). Satisfactory agreement is found between the calculated and the experimental data for water content bigger than 10%.

ps=R¯TMwνwln(P0P)E18

where ps is swelling pressure, the maximum axial pressure which is needed to maintain the original sample height, R¯is the molar gas constant, T is the Kelvin temperature, Mw is molecular mass of water and νw is partial specific volume of water. Combining Eq. (17) with (18), the relationship between normalized water volume by the clay volume and swelling pressure is written as

VwVmpsD3E19

The surface fractal dimension of clay can be obtained from Eq. (19). The experimental data of the swelling pressure and swelling deformation tests are compiled in Fig.7. It is seen that the relationship between the normalized water content and swelling pressure or vertical overburden pressure can be described by the same linear function which is expressed in Eq. (19). Therefore, Eq. (19) represents a general relationship between the water volume absorbed by clay per unit volume and vertical pressure supplied to specimen in swelling process. The surface fractal dimension of Tsukinuno bentonite is 2.63, which nearly equals to that obtained from Fig. 6. In Fig. 6, Cb is the bentonite content.

Swelling deformation tests of Wyoming bentonite were also conducted by Millins et al. (1996) and Studds et al. (1998). The relationship between the clay void ratio (ec) and vertical pressure is shown in Fig. 8. Relationship between the clay void ratio (ec) and normalized water volume is given by

ec=1CmVwVmE20

where Cm is the montmorillonite content, and Cm=75%. It is seen that the slope of logec-logp equals to that of log(Vw/Vm)-logp from Eq. (20) for Cm being constant. Thus, the surface fractal dimension of Wyoming bentonite can be obtained from the linear relationship of logec-logp. The surface fractal dimensions of Wyoming bentonite are 2.64, which obtained from the swelling deformation tests conducted by different authors.

Figure 7.

Relationship between absobed water volume and swelling pressure or vertical overburden pressure of Tsukinuno bentonite

Figure 8.

Relationship between clay void ratio and vertical overburden pressure of Wyoming bentonite

3.3. Correlation between the fractal dimension of PSD and that of the GSD

The pore-size distribution (PSD) is not often and easily measured experimentally. The grain-size distribution (GSD) is easily measured and the fractal dimension of the pore surface can be conveniently obtained from the grain-size distribution (GSD). The relationship between the pore-size distribution (PSD) and the grain-size distribution (GSD) can be constructed for the soils with a homogeneous fabric. Consider a pore limited by a number of grains Ng, the pore radius r is related to the radius rg of the smallest grain as follows (Watabe et al., 2000)

r=fprgE21

where fp is a constant less than 1.0. If the particles are randomly arranged, the probability of having a grain of radius less than rg is dPg, the probability of having the Ng-1 grains with a radius greater than rg is (1-Pg)Ng-1. The probability P of having a pore with a radius less than r is given by (Watabe et al., 2000)

dP=Ng(1Pg)Ng1dPgE22

It is assumed that the grain-size distribution (GSD) satisfies the fractal model, and is given by

Pg=(rg/Rg)3DgE23

where Dg is the fractal dimension of the grain-size distribution (GSD), Rg is the maximum radius of soil grains. The probability P is given by (Watabe et al., 2000)

P=1[1(rfpRg)3Dg]NgE24

If r/(fgRg)<<1, Eq.(24) can be written as

P=Ng(rfpRg)3DgE25

It is seen that the pore surface has the same fractal dimension as that obtained form the grain-size distribution (GSD) from Eqs. (23) and (24). Comparisons between the pore-size distribution (PSD) and the grain-size distribution (GSD) are shown in Fig. 9, where P=V(≤r)/Vv and Pg=M(≤r)/MT, V and M are the pore volume and grain mass, respectively. It is seen from Fig. 9 that both the pore-size distribution (PSD) and the grain-size distribution (GSD) can be expressed by fractal model, and they have the nearly same fractal dimension.

Figure 9.

Correlation of the pore-size distribution (PSD) to the grain-size distribution (GSD).

Figure 10.

Model test for simulating the water distribution

4. Soil-water Characteristic Curves (SWCC)

The soil-water characteristic curve defines the relationship between water content (gravimetric water content w, volumetric water content θ) or degree of saturation (S) and matric suction (ψ). The soil water characteristic curve (SWCC) has been used extensively in the study of unsaturated soils and relates the water content and soil matric suction. The relationship between the matric suction and the radius of the incurvated surface between the pore-air and pore-water is expressed by the Young-Laplace equation. The hydraulic and mechanical properties of unsaturated soils are correlated to the microstructures of soils through the Young-Laplace equation.

For unsaturated soils, the pore water is usually localized in small micropores, and the micropores are usually fully filled with water. The water content in macropores is very little, sometimes, there is no water in macropores in highly unsaturated soils (Fig. 10). The water distribution (black domain) of model test shown in Fig. 10 proved the above assumptions. In the model test, the soil grains of different sizes were simulated by the aluminum rods with different radius (Matsuoka, 1999). The water mainly distributed in the small pore and no water exists in the large pore for unsaturated soils in the model test.

Three air-water distribution states were identified as saturated funicular state, complete pendular state, and partial pendular state, respectively. All the voids are filled with liquid in the saturated state. Air bubbles are present in the voids in the funicular state, and the liquid phase is continuous. There is no continuity in the water phase in the pendular state. Pore water within unsaturated soils can be divided into three forms (Wheeler and Karube, 1996) (Fig. 11): bulk water within those void spaces that are completely flooded, meniscus water surrounding all inter-particle contact points that are not covered by bulk water and absorbed water (which is tightly bounded to the soil particles and acts as parts of the soil skeleton). The relationship between the soil-water characteristic curve (SWCC) and the pore-water forms is shown in Fig. 11. When the volumetric water content is less than its residual value, the pore-water is tightly absorbed by soil particle and cannot move freely. This absorbed water is seen as a part of soil particles. Thus, the contribution of the filled pores with radius r→dr to the water content is given by

Figure 11.

Relationship between the SWCCs and the pore-water forms of unsaturated soil (Matsuoka, 1999)

Figure 12.

Comparison between the prediction and experiments of the SWCCs (Data from Watabe et al., 2000)

dΛ=N4πr2drVTE26

where Λ is the relative volumetric water content, and Λ=θ-θr, θ and θr are the actual and residual volumetric water content, respectively. The residual volumetric water content is the volumetric water content at which the effectiveness of matric suction to cause further removal of water requires vapour migration. Substituting Eq. (9) into Eq. (26), Λ is given by

Λ=Br3DE27

where B=4πC/[VT(3-D)]. Similarly with Eq. (27), the relatively volumetric water content at saturation is written as follows

Λs=BR3DE28

where R is the maximum radius of soil pores.

The relationship between matric suction ψ and the pore radius is obtained from the Young-Laplace equation. Thus, the soil-water characteristic curve (SWCC) is obtained as follows

Θ=(ψψe)δE29

where δ=D-3, Θ is the normollized volumetric water content, and is written as follows:

Θ=θ-θrθs-θr=SSr100SrE30

where θ, θr and θs are the volumetric water content, residual volumetric water content and saturated volumetric water content, S and Sr are the degree of saturation and residual degree of saturation. The residual volumetric water content is not always made routinely, in which case it has to be estimated by extrapolating available soil-water characteristics data towards lower water content, such as shown in Fig. 11. van Genuchten (1980) defined the residual volumetric water content as the water content for which the gradient dθ/dψ becomes zero at high matric suction. From a practical point of view it seemed sufficient to define θr as the water content at some large matric suction, e.g. at the permanent wilting point ψe=1500kPa.

Equation (29) is the soil-water characteristic curve (SWCC) derived from the fractal model for the pore surface. The normalized volumetric water content is equal to unity at values of matric suction up to the air-entry value and equal to zero after residual saturation. The only two parameters ψe and D, which have obvious physical meaning, are used in Eq. (29) to express the soil-water characteristic curve (SWCC). The fractal dimension of the pore surface and the maximum pore radius can be evaluated from the pore-size distribution (PSD) measured using the mercury intrusion tests. The air-entry value can be calculated from the Young-Laplace equation using the maximum pore radius. Thus, the soil-water characteristic curve (SWCC) can be calculated from Eq. (29) using the fractal dimension and the air-entry value obtained from the mercury intrusion tests.

The soil-water characteristic curves (SWCC) were calculated from Eq. (29) using the fractal dimension of the pore surface and the air-entry value for specimens of S-02, S-03 and S-04, which listed in Table 1. The value of the normalized volumetric water content is calculated from Eq. (30). Here Sr=7.5% (Watabe et al., 2000). Experimental data show in good accord with the calculation of the soil-water characteristic curve (SWCC) in Fig. 12.

The experimental data of the soil-water characteristic curve (SWCC) measured in a multi-step outflow experiment were given by Vogel and Roth (1996) for silty soils. The saturated volumetric water contents are 0.502 and 0.448 for silty soils at the A and B horizons, respectively. The residual volumetric water content is equal 0 for both A and B horizons. Comparisons between the calculations of Eq. (29) and experimental data of the soil-water characteristic curves (SWCC) were shown in Fig. 13. The calculation of the soil-water characteristic curves (SWCCs) for silty soils was obtained from Eq. (29). The fractal dimensions of the soil pore and the air-entry values were listed in Table 1.

Figure 13.

Comparison between prediction and experiments of SWCCs (Data from Vogel and Roth, 1996)

Figure 14.

Comparison of the SWCCs for Toyoura sand (Data from Uno et al., 1998)

Uno et al. (1998) gave the experimental data of the soil water characteristic curves for Toyoura sand. Using the fractal dimension and the air-entry value obtained from the pore-size distribution (PSD), the soil-water characteristics can be simulated from Eq. (29) for Toyoura sand. The parameters obtained from the pore-size distribution (PSD) and used to predict unsaturated hydraulic conductivity are listed in Table 1. θs=0.425, θr=0. Comparison between the prediction of Eq. (29) and experimental data of the soil-water characteristics is shown in Fig. 14 for Toyoura sand.

Stingaciu et al. (2009) evaluated the feasibility of using nuclear magnetic resonance (NMR) relaxometry measurements to characterize pore size distribution and hydraulic properties in four porous samples with different texture and composition. The sandy samples FH31 and W3 were loaded in the penetrometer and packed at the same packing density as for the NMR measurements. The Mix8 and MZ samples were used as solid conglomerates. From the relaxation time distribution functions, the cumulative pore-size distribution functions were calculated with the average surface relaxivity. The normalized cumulative pore-size distributions functions were displayed in Fig. 15. The fractal dimension of pore surface and the maximum pore radius can be obtained from Fig. 15, and were listed in Table 2.

SampleMethodBulk density (g/cm3)θsθrDR (mm)ψe (kPa)
FH31Pressure plate1.580.320.020.850.0652.31
W3MSO1.420.330.0581.510.01758.6
Mix8Rosetta1.450.410.0611.520.00053290
MZPressure plate1.600.4402.290.00035429

Table 2.

Parameters for the prediction of SWCCs

The soil-water characteristic curves (SWCC) were determined using the standard sand bed, pressure cell or multistep outflow method, and were plotted in Fig. 16. SWCC of FH31 and MZ were based on pressure plate measurements. SWCC of W3 was determined using multistep outflow (MSO) and that of Mix8 using ROSETTA software (Schaap et al., 2001). The matric suction as plotted in the abscissa (Fig. 16) was transformed into pore diameter using the Young-Laplace equation. In addition, the ordinate was the normalized volumetric water content. The predictions of SWCC were conducted using Eq. (29) with the parameters listed in Table 2. The parameters were derived from the pore-size distribution (PSD) shown in Fig. 15. Comparisons of the predictions and experimental data were shown in Fig. 16. The predictions of Eq. (29) were in good accord with the experimental data.

In general, the soil-water characteristic curves (SWCC) were usually interpreted as cumulative distribution functions in comparison to pore-size distribution (PSD) functions Using the soil-water characteristic curves (SWCC), pore-size distribution can be extracted for a given porous medium on the basis of an empirical law that related the pore suction to the effective pore radius. The soil-water characteristic curves (SWCC) and the pore-size distributions (PSD) are related by correlations as: (1) the absolute values of the slopes of SWCC and PSD are equal; (2) the air-entry value of SWCC is related to the maximum pore radius of PSD through the Young-Laplace equation.

The following conclusions can be obtained from Figs. 12-16: (1) The surface of the soil pore can be described by fractal model, and the surface fractal dimension of the soil pore and the air-entry value can be evaluated from the pore-size distribution (PSD). The soil-water characteristic curves (SWCCs) can be calculated using fractal model of the soil pore. (2) The fractal dimension and air-entry value obtained from the soil pore-size distribution (PSD) are equivalent to those obtained from the soil-water characteristic curve (SWCC). Hence the fractal dimension and air-entry value can be obtained from the fitting of the soil-water characteristic curve (SWCC) if the pore-size distribution (PSD) were not measured.

Figure 15.

Cumulative pore-size distribution of for the three artificial substrates (FH31, W3, and Mix8) and the natural soil Merzenhausen (MZ), from which the fractal dimension of pore surface and the maximum pore radius were derived (Data from Stingaciu et al. 2009)

Figure 16.

Comparison between the predictions of Eq. (29) and experimental data of SWCCs (Data from Stingaciu et al. 2009)

5. Unsaturated hydraulic conductivity and diffusivity

5.1. Brief review

Brutsaert (2000) made a brief review of the capillary model, which was cited as follows.

It is very difficult to determine unsaturated hydraulic conductivity through experiments. Hydraulic conductivity for unsaturated soils is a function of soil-water content θ. For this reason, over the years many attempts have been made to represent the hydraulic conductivity by simple parametric equations, in terms of other properties of the soil which are easier to determine. One of these was of the following form (Brutsaert, 2000):

kr=ΘωE31

In the past, the power form of Eq. (31) has been derived on the basis of some widely different conceptual models of the pore geometry. There have been some marked differences in the obtained values of ω, depending on the adopted model of the pore geometry.

  1. Uniform pore size models

A simple approach was that porous medium was characterized by some equivalent uniform pore size and the variability of the pore sizes was not considered. This approach invariably resulted in a constant value of ω. The capillary tube model of Averyanov (Polubarinova-Kochina, 1952) led to ω=3.5, whereas the extension of the hydraulic radius model of Kozeny (Kozeny, 1927) to unsaturated soils by Irmay (1954) produced ω=3.

  1. Parallel models

In parallel models, the pore system was assumed to be equivalent with a bundle of uniform capillary tubes of many different sizes. The pore-size distribution was derived from the soil-water characteristics (SWCC), which can be related to the effective pore radius r by the Yong-Laplace equation. The true mean velocity in each pore was described by a Poiseuille-like equation for creeping flow. Purcell (1949) and Gates and Tempelaar-Lietz (1950) first used this approach. Several subsequent applications of the parallel model can be written in the common form (Brutsaert, 1967):

kr=β0Θψ(x)2dxβ001ψ(x)2dxE32

where x is the dummy variable representing Θ. The variable β=β(Θ) is related to the tortuosity, and β0 is its value at satiation, when Θ=1. The tortuosity concept was originally introduced by Carman (1956) as an improvement on the non-uniform hydraulic radius model of Kozeny (1927) and it can be expressed as T=(Le/L)2, in which Le is the actual or microscopic path length of the fluid particles in the pores and L is their apparent or macroscopic path length along the Darcy stream lines. The parameter β was introduced in the parallel model by Wyllie and Spangler (1952), because without it, Eq. (32) yielded values considerably larger than their experimental data. Burdine (1953) proposed on the basis of his experimental data that β/β0=Θ2.

A further development was performed by adopting the following soil-water characteristic curve (SWCC) (Brooks and Corey, 1966):

Θ=1forψψeΘ=(ψψe)bforψψeE33

Brooks and Corey (1966) integrated (32) with Burdine's assumption for β (Burdine, 1953) to obtain Eq. (31), with an exponent

ω=3-2bE34

It should be noted that, without Burdine’s assumption for the tortuosity, the parallel model Eq. (32) applied with Eq. (33) would have yielded

ω=1-2bE35
  1. Series-parallel models

The theoretical construction of this model also started with a bundle of parallel pores, each with a different but uniform size. However, these pores were then cut normally to the direction of flow with two resulting faces, and finally after some random rearrangement of the tubes the faces are joined again. This way account was taken of the random variations of the pore sizes, not only in the plane normal to the direction of flow, but also along the direction of flow. The discharge rate in each single pore, which consisted of two sections in cases, was assumed to be governed by the section with the smaller diameter. The pore-size distribution was derived form of the soil-water characteristics by means of the Yong-Laplace equation. The true velocity in each pore was obtained by means of a Poiseuille-like equation.

This model was originally proposed by Childs and Collis-George (1950) in a finite-difference scheme to calculate hydraulic conductivity from experimental data of the soil-water characteristics (SWCC). It was subsequently reformulated in integral form by Brutsaert (1968), to allow the derivation of more concise analytical expressions for hydraulic conductivity. The integral form can be written as

kw=[(2Tγ)2θ02(1Sr)2G]×[0Θ0xψ(x)2dydx+0Θψ(x)20Θdydx]E36

where G is a geometrical constant. In the case of Poiseuille’s equation, G=8. It was trivial to show by integration by parts that the first double integral on the right was identical to the second. Hydraulic conductivity can be expressed concisely as (Brutsaert, 2000):

kr=0Θ(Θx)ψ(x)2dx01(1x)ψ(x)2dxE37

The capillary model similar with Eq. (37) in form was originally proposed for by Fatt and Dykstra (1951).

5.2. Fractal model for hydraulic conductivity

Soil pore systems are viewed as a collection of interconnected voids. The individual voids are considered to have two equilibrium states: either filled by water or empty. Flow within the soil pores is assumed to be described by (Burdine, 1953; Mualem, 1976)

v=r2gcηdhdxE38

where v is the average velocity within the pore; r is the radius of the pore; g is the gravity acceleration; η is the kinematic viscosity; c is a constant depending on the pore geometry; h is the hydraulic head, x is the position axis. Let us consider a porous slab of thickness ∆x isolated from a homogenous soil column by two parallel cross-sections normal to the x axis. The areal porosity at each face of the slab is the same and it is equal to the volumetric porosity. The major simplifying assumption is that the actual configuration of slab may be replaced by a set of capillary tubes with different radii, parallel to the x axis. The flux dq flowing through the connections between pores of radius r1r1+dr1 on one side of the slab (at x) and pores of radius r2r2+dr2 on the other side of the slab (at x+dx) can be expressed by (Mualem, 1976)

dq=βre2(r1,r2,ρ)Ae(r1,r2,ρ)dhdxdr1dr2E39

where β is a constant which incorporates the fluid with the matrix properties; re(r1, r2, q) is the effective radius; Ae(r1, r2, q) is the effective area; ρ is the maximum radius of the water filled pores at volumetric water content θ. According to the Darcy equation, the hydraulic conductivity kw is given by

kw=qdh/dx=β0ρ0ρre2(r1,r2,ρ)Ae(r1,r2,ρ)dr1dr2E40

The relatively hydraulic conductivity (kr) is the ratio of the hydraulic conductivity at any volumetric water content (kw) to the hydraulic conductivity at saturation (ks). The relative hydraulic conductivity (RHC) is related to the radius of the soil pores, and is written as (Mualem, 1976)

kr=kwks=0r0rre2(r1,r2,r)Ae(r1,r2,r)dr1dr20R0Rre2(r1,r2,r)Ae(r1,r2,r)dr1dr2E41

where r is the maximum radius of the water filled pores at the volumetric water content θ. It is assumed that the pore distribution at the two cross-section is completely random. The probability of pores of radius r1r1+dr1 at x to encounter pores of radius r2r2+dr2 at x+dx is given by

Ae(r1,r2)dr1dr2=f(r1)f(r2)dr1dr2E42

The effective radius is assumed to be equal to the radius r reduced by a factor, which is equal to the ratio between the effective area to flow and the actual pore area. The effective radius is given by (Mualem, 1976)

re=r0r0rAe(r1,r2,r)dr1dr20Rf(r1)dr1=rΛ1/2E43

Substituting Young-Laplace equation and Eqs.(42) and (43) into Eq.(41), the relative hydraulic conductivity (RHC) is obtained as

kr=(ψψe)ξE44

where ξ=3D-11. Substituting Eq. (29) into Eq. (44), the relative hydraulic conductivity (RHC) unction expressed by effective degree of saturation is written as follows

kr=Θξ/δE45

Equations (44) and (45) are the relative hydraulic conductivity (RHC) function derived from the fractal model of the pore surface. The parameters in the proposed hydraulic conductivity function have an obvious physical implication, and can be determined from the pore-size distribution (PSD). Eqs. (44) and (45) offer a powerful, physical-based method to calculate the relative hydraulic conductivity (RHC). Eq. (44) is nearly similar with the Brooks and Corey equation (Brooks and Corey, 1966) in form, and the parameters in Eq. (44) have obvious physical implication. Eqs. (44) and (45) offer a powerful, physical-based method to calculate the relative hydraulic conductivity (RHC).

The expression of the soil-water diffusivity (SWD) can be derived from the unsaturated hydraulic conductivity and soil-water characteristics, and written as follows:

d=kw|dψdθ|E46

where d is the soil-water diffusivity (SWD), kw is the saturated hydraulic conductivity. Substituting Eqs. (44) and (45) into Eq. (46), the soil-water diffusivity (SWD) is written as

d=ksψe(ψψe)ζE47
d=ksψeΘζ/δE48

where ς=2D-7.

As given in Eqs. (29), (44), (45), (47) and (48), methods to evaluate the soil-water characteristics (SWCCs), unsaturated hydraulic conductivity (RHC) and soil-water diffusivity (SWD) are proposed using the fractal model for the pore surface. The fractal dimension and the maximum radius of soil pores can be obtained from the pore-size distribution (PSD). The air-entry value can be calculated from the Young-Laplace equation using the maximum radius of soil pores. The soil-water characteristic curves (SWCC), unsaturated hydraulic conductivity (RHC) and soil-water diffusivity (SWD) can be predicted using the fractal dimension and air-entry value, which can be determined from the pore-size distribution (PSD). According to Eq. (29), the fractal dimension of the pore surface and the air-entry value can also be evaluated from the fitted soil-water characteristic curve (SWCC). Thus, the unsaturated hydraulic conductivity (RHC) and soil-water diffusivity (SWD) can also be estimated from Eqs. (44), (45), (47) and (48) using the fractal dimension and air-entry value evaluated from the soil-water characteristic curve (SWCC).

1. Predictions of SWCC and RHC from PSD

Vogel and Roth (1998) inferred effective hydraulic properties of unsaturated soil from the structure of the pore surface. The water retention characteristic and the hydraulic conductivity were simulated by network models with 323=32768 nodes. The hydraulic conductivity, kw, is determined at each step of desaturation by imposing a pressure gradient ∆p across the ends of the network. Water flow qij in a bond between the nodes i and j through a horizontal surface A is described by Poiseuille’s law. The hydraulic conductivity was calculated as =qL/(Ap), where L is the vertical length of the network. Comparison between the predictions of Eq. (44) and experimental data given by Vogel and Roth (1998) was shown in Fig. 17. The parameters used to do prediction were listed in Table 1, which were obtained from the fractal model for the pore-size distribution (PSD). The predictions were in good accord with the experimental data.

Figure 17.

Comparison between the predictions of Eq. (44) and experimental data of the RHC ( Data from Vogel and Roth, 1998)

Uno et al. (1998) studied the relationship between the pore-size distribution (PSD) and unsaturated hydraulic properties of Toyoura sand. The relationship between the relative hydraulic conductivity (RHC) and the normollized volumetric water content Θ can be calculated from Eq. (45) using the parameters listed in Table 1. The prediction of Eq. (45) satisfactorily agrees with the experimental data of relative hydraulic conductivity (RHC) for Toyoura sand in Fig. 18. The proposed function of relative hydraulic conductivity (RHC) (Eq. (45)) is verified by the experimental data in Fig. 18.

Figure 18.

Comparison between predictions of Eq. (45) and experimental data of the RHC for Toyoura sand. (Data from Uno et al., 1998)

2. Prediction of RHC and SWD from SWCC

The characteristics of the pore-size distribution (PSD) were prerequisite to determine the relative hydraulic conductivity and soil-water diffusivity of unsaturated soils in Eqs. (44), (45), (47) and (48). A close relationship was established between the soil-water characteristic curve (SWCC) and the pore-size distribution by fractal model for the pore surface. The pore-size distribution (PSD) can be estimated from the soil-water characteristic curve because many soils show that the soil-water characteristic curve is equivalent to the pore-size distribution (PSD) function (Watabe et al., 2000). Thus, the fractal dimension of the pore surface and the air-entry value can be estimated from the fitting of the soil-water characteristic curve (SWCC). Thus, unsaturated hydraulic properties, such as relative hydraulic conductivity (RHC) and soil-water diffusivity (SWD) can be predicted using the fractal dimension and air-entry value obtained from the soil-water characteristic curve (SWCC).

Smettem and Kirkby (1990) offered the hydraulic properties of haploxeroll loam. Experimental data of the soil-water characteristic curve (SWCC) were shown in Fig. 19. The fractal dimension of the pore surface can be obtained from the fitted soil-water characteristic curve (SWCC). Through fitting the soil-water characteristic curve (SWCC), we obtained that δ=-0.37 for haploxeroll loam from Fig. 19. Fractal dimension of the pore surface can be calculated from Eq. (29) using the fitting value of δ. The fractal dimension (D) of the pore surface is 2.63 for haploxeroll loam. The air-entry value is 0.14kPa obtained from Fig. 19.

The hydraulic properties of unsaturated soils can be predicted using the fractal dimension of the pore-size distribution (PSD), obtained from fitted soil-water characteristic curve (SWCC). The fundamental parameters for predictions of hydraulic properties were tabulated in Table 3. The comparisons between predicted result and experimental data of the relative hydraulic conductivity (RHC) were shown in Fig. 20. The coefficient of hydraulic conductivity at saturation is 140 mmh-1. Fig. 20 depicted the relationship between the relative hydraulic conductivity (RHC) and water potential for haploxeroll loam. The parameter ξ was -3.11 calculated from Eq. (44) using the fractal dimension D=2.63 for haploxeroll loam. The predictions of the relative hydraulic conductivity (RHC) calculated from the proposed function Eq. (44) are in satisfactory agreement with the experimental data.

The relationship between the soil-water diffusivity and water potential for haploxeroll loam was also shown in Fig. 20. The values of ks and ψe are 140mmh-1 and 0.14kPa, respectively, used in Eq. (47). The parameters ς was -1.74, calculated from Eq. (47) using the fractal dimension D=2.63. The prediction of the soil-water diffusivity was the line passing through the point of (0.14, 0.64) with the slope of -1.74. The prediction of the soil-water diffusivity nearly agreed with the experimental data. The proposed function for the soil-water diffusivity was verified by the comparison between prediction and measured results. The parameters obtained from the soil-water characteristic curve (SWCC) and used to predict the relative hydraulic conductivity (RHC) and soil-water diffusivity (SWD) for haploxeroll loam were listed in Table 3.

Soil typeDδξξ/δςψe (kPa)
Haploxeroll loam2.63-0.37-3.118.41-1.740.14
Toyoura sand (wetting)1.60-1.40-6.24.43-3.8/
Toyoura sand (drying)1.26-1.74-7.224.15-4.48/
Loamy sand2.68-0.32-2.969.25-1.642.0
McGee Ranch soil2.51(GSD)-0.49-3.477.08-1.984.6

Table 3.

Parameters used to predict the unsaturated hydraulic conductivity

Figure 21 shows the experimental data and the fitting of the soil-water characteristic curve (SWCC) for loamy sand. The experimental data are scaled from Simunek et al. (1999). The parameter δ was -0.32, which is obtained from the fitting of the soil-water characteristic curve (SWCC) for loamy sand from Fig. 26. Hence the fractal dimension of the pore surface was 2.68 calculated from Eq. (29) using the fitting value of δ for loamy sand. The air-entry value was given by the intersection between the line expressed by Eq. (29) and the line of Θ=1 in log-log plot. The air-entry value is 2kPa obtained from Fig. 21.

The parameters used to predict the hydraulic properties of loamy sand were listed in Table 3. Substituting D=2.68 and ψe=2kPa in Eq. (44), the predictions of relative hydraulic conductivity (RHC) can be obtained. The comparisons between predictions and experimental data of the relative hydraulic conductivity were shown in Fig. 22. The parameter ξ is -2.96 calculated from Eq. (44) using the fractal dimension D=2.68 for loamy sand. The predictions of the relative hydraulic conductivity were in good accord with the experimental data. The proposed function for the relative hydraulic conductivity was validated by the good agreement between the predictions and measurements.

The parameters of fractal model were first calibrated by matching the measured moisture-suction data. It should be noted that the experimental values of hydraulic conductivity were not matched by adjusting the model parameters. The calibrated SWCC model parameters were, instead, directly used to predict the relative hydraulic conductivity. The soil-water characteristic curves (SWCC) and their calibration of eight soils were shown in Fig. 23. The calibration of SWCC model parameters were listed in Table 4. As seen in these figures, the measured moisture-suction data for these soils were unavailable for the full range (0%-100%) of the degree of saturation. Predicting the suction beyond the available experimental data range was a challenging task since the pattern of variation was unknown (Ravichandran and Krishnapillai, 2011).

Figure 19.

Fractal dimension of the pore surface and air-entry value evaluated from SWCC (Data from Smettem and Kirkby, 1990)

Figure 20.

between predictions and experimental data of RHC and SWD for haploxeroll loam (Smettem and Kirkby, 1990)

Figure 21.

Fractal dimension and the air-entry value obtained from the fitting of SWCC for of loamy sand (Data from Simunek et al., 1999)

Figure 22.

Comparison between measurements of RHC and prediction of Eq. (44) for loamy sand (Data from Simunek et al., 1999)

SoilθsSr (%)Dψe (kPa)Reference
Lakeland0.37530-1001.752.5Elzeftawy and Cartwright, 1981
Superstition0.530-1002.71.2Richards, 1952
Columbia sandy loam0.45850-1002.294.5Brooks and Corey, 1964
Touchet silt loam0.4320-1002.097Brooks and Corey, 1964
Silt loam0.39650-1002.5915Reisenauer, 1963
Guelph loam0.5245-1002.753Elrick and Bowmann, 1964
Yolo light clay0.37545-1002.871.5Moore, 1939
Speswhite kaolin0.5655-1002.8515Peroni et al., 2003

Table 4.

Parameters of the fractal model for eight selected soils

Figure 23.

Fractal dimension and air-entry value obtained from the fitting of SWCC

Figure 24.

Comparison between predictions of Eq. (44) experiments of RHC

The relative hydraulic conductivity of the above mentioned eight soils were predicted using the fractal model with the same fitting parameters listed in Table 4. The fitting parameters were obtained by matching the experimental SWCC. Fig. 24 illustrated the prediction of Eq. (44) and experimental data of relative hydraulic conductivity. It should be noted that the fractal model parameters were not calibrated or adjusted to match the measured hydraulic conductivity values. The predictions of Eq. (44) were comparable for sand and loam, as shown in Fig. 24. The fractal model showed slight deviations from the measured data of clay. In the case of Yolo light clay and Speswhite kaolin, the difference between the predictions and the experimental data increased as the suction increases (Fig. 24). The fractal model showed better predictions at lower suction range. However, the accuracy of fractal model in the higher suction range could not be verified because the experimental results were available only for the lower suction ranges.

It is seen that it is feasible to estimate the relative hydraulic conductivity of unsaturated soils using the fractal dimension of the soil pore surface from Figs. 19-24. A practical method to express the relative hydraulic conductivity was proposed to use the fractal dimension and the air-entry value, which can be obtained from the fitting of the soil-water characteristic curve (SWCC).

3. Prediction of RHC and SWD from GSD

The feasibility to determine the soil hydraulic properties using the fractal dimension of the grain-size distribution (GSD) is examined by the experimental data of McGee Ranch soil. The grain-size distribution (GSD) was measured by Hunt and Gee (2002) and was shown in Fig. 25. The fractal dimension of the grain-size distribution (GSD) is 2.51, obtained from Fig. 25.

Hunt and Gee (2002) gave the values of θs and θr were 0.4 and 0.01, respectively, the air-entry value was 4.6kPa. Using the fractal dimension of the grain-size distribution (GSD), the soil-water characteristic curve (SWCC) and the relative hydraulic conductivity (RHC) can be determined. The parameters of the grain-size distribution (GSD) used to determine the hydraulic properties are listed in Table 3. Comparisons between the predictions of Eq. (29) using the fractal dimension of the grain-size distribution (GSD) and the experimental data of soil-water characteristic curve (SWCC) are shown in Fig. 26. It was seen that a good result is obtained to predict the soil-water characteristic curve (SWCC) using the fractal dimension of the grain-size distribution (GSD).

Hunt and Gee (2002) gave the values of saturated hydraulic conductivity was 0.001cm/s. Comparisons between the predictions of Eqs. (44) and (45) and experimental data were shown in Fig. 27. The prediction of Eqs. (44) and (45) nearly satisfied the experimental data of relative hydraulic conductivity in Fig. 27. In Fig.27, the predictions of Eq. (45) deviated from the experimental data near saturation, and extend to equal to the experimental results at low water content. The soil hydraulic conductivity can be approximately determined using the fractal dimension of the grain-size distribution (GSD).

Figure 25.

Fractal dimension of GSD for McGee Ranch soil (Data from Hunt and Gee, 2002)

Figure 26.

Comparisons between predictions and experiments of SWCC for McGee Ranch soil (Data from Hunt and Gee, 2002)

6. Comparison of the fractal model with the van Genuchten–Mualem (G-M) model

6.1. van Genuchten–Mualem (G-M) model for relative hydraulic conductivity (RHC)

van Genuchten (1980) derived an empirical relationship to describe the soil-water characteristic curve (SWCC):

Θ=1(1+|αψ|n)mE49

where α, n and m are the van Genuchten curve-fitting parameters, and m=1-1/n.

Burdine (1953) and Mualem (1976) presented a similar model to estimate the unsaturated hydraulic conductivity from pore size distributions inferred from soil water retention characteristics. The models share a number of similarities, allowing them to be written in a general form as (Hoffmann-Riem et al., 1999)

kr=kwks=Θl(0ΘψχdΘ01ψχdΘ)γE50

where kw is the unsaturated hydraulic conductivity at any water content, ks is the saturated hydraulic conductivity, Θ is the normalized water content, ψ is matric suction. Generally accepted parameter values for the parameters (l, χ, γ) are (0.5, 1, 2) and (2, 2, 1) for the Mualem (1976) and the Burdine (1953) models, respectively.

Using the Mualem model, van Genuchten (1980) derived a closed-form to determine the relative hydraulic conductivity (RHC) at a degree of saturation,

kr=[1|αψ|n1(1+|αψ|n)m]2(1+|αψ|n)mE51
kr=Θ12[1(1Θ1m)m]2E52

Equations (51) and (52) are the representations of the G-M model, which is the most widely and popularly used to predict the hydraulic conductivity.

The other models for relative hydraulic conductivity (RHC) were also introduced as follows.

6.2. van Genuchten–Burdine (G-B) model for relative hydraulic conductivity (RHC)

Hydraulic conductivity resulting from the Burdine capillary models (Burdine, 1953) was expressed as

kr=Θ2[1(1Θ1m)m]E53

Replacing Eq. (49) in Eq. (53) yielded

kr=[1|αψ|n2(1+|αψ|n)m](1+|αψ|n)2mE54

where m=1-2/n.

Note that for both Eqs. (52) and (54), because m is less than 1, the derivative dk/dθ is infinite for θs, whatever the value of n.

6.3. van Genuchten–Fatt & Dykstra (G-FD) model for relative hydraulic conductivity (RHC)

The capillary model proposed by Fatt and Dykstra (1951) was expressed as

kr=0Θ(Θx)ψ(x)2-cdx01(1x)ψ(x)2cdxE55

Evaluating ψ(x) from Eq. (49), the integrals of Eq. (55) with the condition of m>-(2+c)/2n and n>2+c resulted in:

kr=ΘB(a1,b)IΘ1/m(a1,b)B(a2,b)IΘ1/m(a2,b)B(a1,b)B(a2,b)E56

where a1=m+(2+c)/n, a2=2m+(2+c)/n, b=1-(2+c)/n, Ix(a, b) and B(a, b) are the Incomplete Beta and Beta functions of positive arguments a and b, respectively, and given by:

B(a,b)=01xa1(1x)b1dx,Iy(a,b)=0yxa1(1x)b1dxB(a,b)E57

The results presented by Touma (2009) were obtained with c=0.5 and a1=1.

Even though the conductivity predicted by the capillary model of Fatt and Dykstra (1951) gave the best results compared with both the quasi-analytical solution and observations, the resulting expression was not easy to use because it was necessary to evaluate Incomplete Beta and Beta functions. In order to simplify G-FD model, Θ was expressed by Eq. (49) and kr was according to Eq. (31) with ω=2+2.5/mn. The fitted result in a value of b was close to -mn, and a conductivity curve close to that predicted by Eq. (56).

6.4. Brooks & Corey-Brutsaert model for relative hydraulic conductivity (RHC)

Brutsaert (2000) presented a capillary model expressed as Eq. (37). Combining this capillary model with the soil-water characteristic curve (SWCC) of BC model (Brooks and Corey, 1966), relative hydraulic conductivity was written as Eq. (31) in form, and the parameter ω was given by

ω=2-2bE58

Combining of BC model for the soil-water characteristics with the Mualem and the Burdine capillary models, the relative hydraulic conductivity curves were also expressed as Eq. (31) in form, and the parameter ω was given by

ω=3+2mmE59

when BC is applied with the Burdine condition, and

ω=2.5+2mmE60

when applied with the Mualem condition.

Here we focused on the G-M model, as it was the most widely used in soil science and hydrology.

van Genuchten (1980) found that the predictions of the G–M model were less accurate for Beit Netofa clay. The fractal dimension and the air-entry value evaluated from the fitted soil-water characteristic curve (SWCC) are 2.89 and 60kPa, respectively (Fig. 28) The fitted parameters of the van Genuchten model, α and n are 0.0015 and 1.17, respectively from Fig.28. The parameters for the determination of the relative hydraulic conductivity (RHC) are listed in Table 5.

Figure 27.

Comparisons between predictions and experimental data of RHC for McGee Ranch soil (Data from Hunt and Gee, 2002)

Figure 28.

Fractal dimension and the parameters of G-M model obtained from SWCC for Beit Netofa clay (Data from van Genuchten, 1980)

Soil typeFractal modelG-M model
δDξξ/δα (kPa-1)n
Beit Netofa clay-0.112.89-2.3321.20.0041.17
Hanford glacial soil 0-099-0.622.38-3.866.230.0151.71
Hanford glacial soil 2-1637-0.822.18-4.465.440.0761.89
Grenoble sand-0.642.36-3.926.130.42.17

Table 5.

Parameters obatined from SWCC for prediction of RHC

The saturated and residual water contents of Beit Netofa clay are 0.446 and 0, respectively. Comparisons between the predictions of both the fractal model (Eq. (44)) and the G–M model (Eq. (51)) and the experimental data of the relative hydraulic conductivity (RHC) for Beit Netofa clay are shown in Fig. 29. It is seen that predictions of the fractal model (Eq. (44)) satisfactorily agree with the experimental data, especially at high matric suction, while the predictions of the G–M model (Eq. (51)) are found to deviate from the experimental data. However, the fractal model shows a sharp corner at the air-entry value point, which is disagree with the measured data that usually show smooth transition after the air-entry value point. This default does not reduce the advantages to predict the relative hydraulic conductivity (RHC) using the fractal dimension obtained from the pore-size distribution (PSD), because the hydraulic conductivity near the point of the air-entry value nearly equal to the saturated conductivity.

Figure 29.

Comparison between the predictions of the RHC using fractal model and G-M model (van Genuchten, 1980)

Figure 30.

Fractal dimension and G-M parameters obtained from SWCC for Hanford glacial (Khaleel and Relyea, 1995)

The measured soil–water characteristic curves of Hanford glacial soil by Khaleel and Relyea (1995) were shown in Fig. 30. The saturated water contents (θs) are 0.338 and 0.303 for the samples 0–099 and 2–1637, respectively. The residual water contents (θr) are 0.039 and 0.025 for the samples 0–099 and 2–1637, respectively. From Fig. 30, it is seen that the soil-water characteristic curves (SWCC) of Hanford glacial soil are in good accord with the fractal model. The fractal dimensions (D), according to Eq. (29), are 2.35 and 2.18 obtained from the soil-water characteristic curves for the samples 0–099 and 2–1637 of Hanford glacial soil, respectively.

From the experimental data of the soil-water characteristic curve (SWCC), the parameters for determining relative hydraulic conductivity (RHC) using the fractal dimension are listed in Table 5. The larger the fractal dimension, the larger the complexity of the pore connectivity will be. The larger the fractal dimension, the larger the change in the water content with the same change in matric suction will occur. From the fittings of the soil–water characteristic curve in Fig. 30, the parameters (a and n) in the van Genuchten model (1980) for the soil–water characteristic curve (SWCC) are obtained, and are also listed in Table 5. Comparisons between the predictions of the fractal model (Eq. (45)) and the experimental data of the relative hydraulic conductivity (RHC) for Hanford glacial soil are shown in Fig. 31. The predictions of the fractal model (Eq. (45)) nearly agree with the experimental data of the relative hydraulic conductivity (RHC) for Hanford glacial soil in Fig. 31. The predictions of the G–M model (Eq. (52)) deviate from the experimental data of the relative hydraulic conductivity (RHC) in Fig. 31. It is seen that the predictions of the fractal model are better than those of the G–M model for Hanford glacial soil in Fig. 31.

Figure 31.

Comparison between the predictions using fractal model and G-M model (Data from Khaleel and Relyea, 1995)

Figure 32.

Fitted SWCC curves

The combinations of Eq. (29) with fractal capillary model were tested on Grenoble sand. Data points for the soil-water characteristic curve (SWCC) and hydraulic conductivity were taken from Touma (2009). The infiltration experiment was conducted under a constant head of 2.3cm (Touma & Vauclin, 1986). The solid lines in Figure 37 were the fitted curves for Eq. (29), and the resulting predicted hydraulic conductivity were shown in Fig. 32. Fractal model gave good results for the Grenoble sand. The parameters obtained through the fitting results of soil-water characteristic curves (SWCC) were used to predict hydraulic conductivity curve. Comparison between the fractal model and the G-M model were shown in Figs. 33. It was found that the prediction of hydraulic conductivity curve using Eq. (45) closed to the experiments for Grenoble sand.

Figure 33.

Comparison between fractal model and G-M model

7. Rainfall infiltration of unsaturated soils

Rainfall induced slope failures often occur as relatively shallow failure surfaces orientated parallel to the slope surface and are observed and analyzed by different mechanisms. The effect of seepage on slope stability is introduced in the analyses by calculating the critical depth for an infinite slope with rainfall infiltration, while the effect of both negative and positive pore water pressures on the stabilities of initially unsaturated slopes are explained and coupled with infinite slope analysis and pore-air flow analysis methods in order to present a predictive formulation of slope failures that occurs in rainfall evens and is derived from a fractal model on unsaturated soil. The formulation serves as a baseline analysis method for evaluating potentially unstable slopes.

The pore water pressure pattern that develops in the initially unsaturated soil will occur as a transient process as the infiltration moves downward into the soil profile. Several factors should been taken into account in unsaturated analyses. The shear strength of the soil mass and the development of seepage forces which both depend on the evolution of the pore water pressure profile must been addressed in detail.

Here, an individual soil slice can be treated as a one-dimensional column with vertical infiltration. Suppose that any lateral flow between the adjacent slices will be equal on the up-slope and down-slope boundaries. Then considering vertical infiltration only could meet the flow continuity requirement. These one-dimensional infiltration analyses are addressed by the saturated/unsaturated seepage finite element software SEEP/W (GEO-SLOPE International Ltd.2007) coupled with the air flow software AIR/W (GEO-SLOPE International Ltd.2007).

For a homogenous, isotropic soil, one-dimensional water flow equation (GEO-SLOPE International Ltd.2007) is given by:

mwγwHwt=y(kwHwt)+mwPat+QwE61

where Hw is the total water energy potential comprised of both pressure and elevation potentials; Pa is the pore air pressure; kw is hydraulic conductivity; y is y coordinate; mw is the slope of SWCC; t is time; and Qw has units of length per time.

For the air conversation of mass, we can arrive at the pore air general mass balance equation (GEO-SLOPE International Ltd.2007) as like:

(ρwmw+θaR¯T)Pat=y[ρakaγ0aPay+ρa2kaρ0a]+ρaλwmwHwtE62

where ka is air permeability; θa is volumetric air content; ρa is air density; γoa is initial air unit weight; ρoa is initial air density; R¯is ideal air constant and R¯=287J/(kg K) for dry air; T is temperature; Hw is total head and Pa is pore air pressure. These two governing equations can be used to determine total pressure head profile for the water and the pore air pressure, and then suction is attained.

The relationships relating volumetric water content and hydraulic conductivity to suction must be known to solve the Eq. (61). And similarly to calculate the pore air pressure the pore air permeability function is also wanted. To find the effect of different fractal dimensions and air-entry values of different soil types, series of parameter studies have conducted. Different types of Soil-water characteristic curves and hydraulic conductivity curves were shown in Figs. 34 and 35. In the legends of these pictures, Soil 10, 2.1,-5, for example, means a type of soil with air-entry value, ψe=10kPa, the fractal dimension, D=2.1 and saturated hydraulic conductivity ks=1×10-5m/s. It is assumed that θr=0. The pore air permeability function curve used in these studies was shown in Fig.36, and it remained a constant (GEO-SLOPE International Ltd.2007).

Figure 34.

Soil-water characteristic curves for soils

Figure 35.

Hydraulic conductivity curves f

Figure 36.

Pore air permeability function used in study

Figure 37.

Sketch map of the one-dimensional infiltration model

A 4m-deep soil column was modeled for the one-dimensional analyses, as shown in Fig. 37, and infiltration occurs from the top of the column, while the bottom was set to be a pervious boundary. A linear hydrostatic suction distribution presented for the initial conditions. Analyses for the rainfall infiltration into different types of soils were performed using the unsaturated characteristic curves shown in Figs. 34-36. In applying a top boundary condition to simulate the infiltration of rainfall, it is important to realize that an influx boundary depends on the relationship between the saturated hydraulic conductivity and the rainfall intensity. In the following analyses, if a rainfall intensity greater than or equal to ks with the non-infiltrating rainfall running off the slope, the top boundary condition is set to a total head equal to 4m; while if field results indicate that the rainfall intensity is less than ks, then a flux type boundary condition is more appropriate.

Since initial failures often have small depth-to-length ratios and form failure planes parallel to the slope surface, the use of infinite slope analysis in modeling the infiltration process by vertical one-dimensional analysis makes it justified in describing the physical process of failure initiation. However, the methods used in traditional infinite slope analysis must be modified to take into account the variation of the suction profile that results from the infiltration process. There are two distinct failure mechanisms can be initiated by the infiltration process. The failure takes place due to positive pore pressure and it takes place while the suction still exists.

The shear strength of unsaturated soil can be represented by fractal model, and written as follows (Xu, 2004):

τ=c+(σua)tanϕ+ψe3DψD2tanϕE63

where τ is the shear strength; c' is the effective cohesion; (σ-ua) is the total normal stress; φ' is effective friction angle; ψ is the suction; ψe is the air-entry value. The stability envelope can be represented as (Collins and Znidarcic, 2004):

dcr=c+γwhctanϕ-γwhptanϕγcos2β(tanβtanϕ)E64

where dcr is the critical depth; β is slope angle; hp is the pore passive pressure and hc is the pore negative pressure.

We can get suction with analyses of both pore water pressure and pore air pressure in the soil slope under rainfall infiltration. Combining Eq. (63) with Eq. (64), the new stability envelope can be derived and written as follows:

dcr=c+ψe3DψD2tanϕ-γwhptanϕγcos2β(tanβtanϕ)E65

The critical depth for infinite slope failure is now a function of the suction, ψ, pressure head hp, the given material and slope characteristics, c', φ', γ, ψe, ξ, γw, and β. In this way, the slope stability issues related to the decrease in shear strength from a loss of soil suction and the development of seepage forces from positive pressure head generation can be clearly understood.

Because the infiltration results and the slope stability results were both presented in terms of the pressure head and the suction profile, the two analyses can be coupled to yield a comprehensive method for determining the location and time of failure for a slope if the soil, slope, and rainfall parameters are given. The methodology of the coupled analysis involves plotting a stability envelope as defined by Eq. (65) for specified soil and slope parameters, over a given infiltration profile generated from the particular rainfall and unsaturated characteristic curves. Coulomb failure and the initiation of slope mobilization are defined at the points where the infiltration trace intersects the stability envelope. Intersections of the infiltration trace and the stability envelope indicate points at which the slope is unstable and can be thought of as ‘‘critical depths’’ of failure. It must therefore be assumed that if the slope is initially stable, the point at which the stability envelope inter sects the initial suction distribution line will not contribute any further information about the failure mechanisms from infiltration.

Here, factors controlling the unsaturated soil slope under infiltration are studied and the results which combine infiltration analyses and stability analyses are shown in Figs.43-47. Infiltration analyses were performed by software SEEP/W coupled with AIR/W in order to calculate the suction file and the stability envelope as defined by Eq.(65) for slope parameters, β=40°, γ=19.6kN/m3, γw=10kN/m3, and shear strength parameters c'=5kPa, φ'=15°. Other parameters are studied in a systematic way.

7.1. Effect of fractal dimension

Soil column was saturated due to rainfall infiltration from the initial hydrostatic suction profile. It was shown that infiltration took place in unsaturated soil much faster for the soil with higher fractal dimension. Three fractal dimensions (D=2.1, 2.35 and 2.78) were calculated respectively and the results were shown in Fig.38. It was seen from Fig. 38 that the higher initial hydraulic conductivity existed on the top of soil column which had a higher fractal dimension. On the other hand, the stability envelope for the soil with higher fractal dimension was steeper. There was an obvious inflection point in the stability envelope of the soil with the low fractal dimension. The results also indicated that it was the relative shape of the unsaturated characteristic curves (Fig.34 and Fig. 35) that has a controlling effect on the suction distribution.

7.2. Effect of air-entry value

Three air-entry values (ψe=1kPa, 10kPa and 20kPa) are calculated respectively and the results as shown in Figs.39. Compared with the results of three different air-entry values, both the infiltration trace and the stability envelope have changed with the different air-entry values. When the air-entry value is relatively small (ψe=1kPa), as shown in Fig.39a, passive pore water pressure was generated on the upper soil column at the time t=32h. Due to the low initial permeability on the top of soil column and a state of lower saturation, fewer pores initially filled with water and there are fewer channels available for fluid transport and consequently the flow of water is hindered. With the commencement of infiltration at the top boundary, water is forced into a soil which is not capable of transporting it efficiently, which results in the development of positive pressure heads as the infiltration front progresses downward. While the air-entry value is relatively larger, as shown in Fig.39c, infiltration moves downward much faster and suction has reached to the air-entry value. Passive pore water pressure occurs at the bottom of the soil column which means that there is little unsaturated zone in the middle of soil column during the rainfall infiltration. Besides, the stability envelope is much smooth when the soil has a relatively low air-entry value, and consequently the critical depth is smaller.

Figure 38.

Infiltration results for soils with different fractal dimension and the new stability envelope

Figure 39.

Infiltration results for soils with different air-entry value and the new stability envelope

7.3. Effect of saturated hydraulic conductivity

The dominate factors which control the stability of unsaturated soil slope under rainfall infiltration are the soil saturated hydraulic conductivity and rainfall intensity. While the soil saturated hydraulic conductivity determines the water transportation ability of soil and the infiltration quantity most depends on the rainfall intensity. In this parameter study series, three different saturated hydraulic conductivity values (ks=1×10-5m/s, 1×10-6m/s, 1×10-7m/s) are taken into account and the results are shown in Fig.40.

When the rainfall intensity is much larger than the saturated hydraulic conductivity, soil in the top region of the column has been saturated soon after the infiltration begins and its suction almost reaches to zero at the time t=1.5h. It can be seen that the suction at the base is increasing over time as a positive number, which means that the soil is actually de-saturating at its base while wetting up above. This is caused by the air which becomes trapped when the top has saturated and starts to resist water infiltration. As the situation that the rainfall intensity is much smaller than the saturated hydraulic conductivity as shown in Fig.40a, the suction of the top region decreases with the infiltration developing, however, the top region has remained unsaturated and the infiltration depth is also small due to the small rainfall quantity into the soil.

Figure 40.

Infiltration results for soils due to different rainfall intensities

7.4. Effect of rainfall intensity

The rainfall intensity for a certain saturated hydraulic conductivity has taken into considered to simulate the infiltration process. The range is changed from 0.3mm/h to 100mm/h and the relationship between the time of failure and rainfall intensity for the three soil types of saturated hydraulic conductivity has plotted in Fig.41. These curves decrease as exponent forms similarly while the magnitude of different curves has changed dramatically. The time at which the soil slope failure occurs can be about ten hours for the soil saturated hydraulic conductivity ks equals to 1×10-5m/s, and for the situation that ks equals to 1×10-6m/s, the time can range from 10 to 100 hours. While much more time is needed for the lower saturated hydraulic conductivity soil type.

Figure 41.

The time of failure vs. rainfall intensity

Figure 42.

Infiltration results without air phrase consideration

7.5. Infiltration analyses coupled with the air flow analyses

In the infiltration analyses, air flow movement in the soil pore is also conducted. The effect of the pore air flow movement is especially obvious in the situation when the rainfall intensity is much larger than the soil saturated hydraulic conductivity. Compared with the results of infiltration analyses coupled with air flow analyses, the results without air phrase consideration are shown in Fig.42. In this condition, the negative pore water pressure gradually reaches zero at the top of soil column. While the pore water pressure remains unchanged at the bottom of the column. If we take the pore air flow movement into consideration, the results are different as shown in Fig.40c. Besides, the suction could be calculated more precisely when we attain both the pore water pressure and the pore air pressure by simulating infiltration analyses coupled with air flow analyses.

8. Slope stability analyses due to rainfall infiltration

It was widely recognized that the rainfall infiltration took a great role in causing landslides, while the relative importance of soil properties, rainfall intensity, initial water table depth and slope geometry in inducing instability of a homogenous unsaturated soil slope under different rainfall was investigated through a series of parametric studies (Rahardjo, et al, 2007). Soil properties and rainfall intensity were found to be the primary factors controlling the instability of slopes due to rainfall, while the initial water table depth and slope geometry only played a secondary role.

The factors affecting the stability of a slope were considered to be the soil properties, rainfall intensity, initial depth of the groundwater table, and the slope geometry (i.e., slope angle and slope height). To assess the effects and relative contribution of controlling factors, a series of parametric studies were performed on a typical geometry of a homogeneous soil slope shown in Fig. 43. It had the boundary conditions as follows: ab, bc, cd=q=Ir(rainfall intensity); ah, de, fg=q=0m3/s (i.e. no flow boundary); and ef, gh=ht(total head at the side). Four slope heights, Hs,(3m, 4m, 5m and 6m), three slope angles, α,(26.6, 33.7 and 45.0), eight initial depths of groundwater table (GWT), Hw,(0.5m, 1m, 2m, 3m, 5m, 7.5m, 10m and 15m), three fractal dimensions D (2.1, 2.35 and 2.78),four air-entry values ψe (1kPa, 10kPa, 20kPa and 50kPa), five values of saturated hydraulic conductivity, ks (10-4m/s, 10-5m/s, 10-6m/s, 10-7m/s and 10-8m/s) and five rainfall intensities Ir(3mm/h, 9mm/h, 15mm/h, 36mm/h and 100mm/h each for 24h duration) were used in six cases of parametric studies. The shear strength parameters of the soils used in the parametric study were c’=10kPa, effective angle of internal friction, φ'=26°, and unit weight of soil, γ=20kN/m3.

Figure 43.

Sketch map of the slope in study

The fractal model of the soil-water characteristic curves for unsaturated soil was written as Eq. (29). The function of relative hydraulic conductivity was used as Eqs. (44) and (45). The derived SWCC and RHC for the all six soils were shown in Figs. 34 and 35. In the legends of these pictures, the symbol S in the soil names represent “soil”, the first number means the air-entry value, and the second one means the fractal dimension.

In the seepage analysis, the governing partial differential equation (GEO-SLOPE International Ltd.2007) for a two-dimensional transient water flow used in the finite element seepage model was written as follows:

mwγw(Hy)t=y(kwHwt)+mwPat+QE66

where H is the total head; kx are ky ars hydraulic conductivity in x, y direction respectively;Q is the applied boundary flux; mw is the slope of the SWCC and γw is the unit weight of water. Eq. (66) was solved by SEEP/W software (Geo-Slope 2007). The boundary conditions used in the transient seepage analysis were shown in Fig. 43. And the initial condition for the analyses was a hydrostatic condition with a limiting pore water pressure of -75kPa. Then the SEEP/W software generated the negative pore water pressure as the initial condition. The pore-water pressures obtained from the seepage analysis were then used in the slope stability analyses to calculate the factor of the slope safety, Fs.

The equation for unsaturated shear strength was used as Eq. (63). In the study, we chose c’=10kPa, φ'=26°. The Morgenstern-Price Method was adopted in the slope stability analysis. This method considered both shear and normal inter-slice forces, satisfied both moment and force equilibrium, and allowed for a variety of user-selected inter-slice force function. The slope stability analysis using Morgenstern-Price Method was performed using SLOPE/W software (Geo-Slope 2007). The pore-water pressures, uw, obtained from the transient seepage analyses using SEEP/W were added to SLOPE/W to be incorporated in the slope stability analyses.

The results were presented with attention on the effects of the factors which impacting on the soil slope stability under rainfall for 24h with a combination of various controlling factors.

8.1. Effect of soil properties

The variation in factor of safety with time for a homogeneous soil slope of constant slope height Hs=6m, initial groundwater table depth Hw=2m, subjected to rainfall intensities of 3, 9 and 15mm/h of respective soil for 24h with a combination of various soil types (different fractal dimensions and air-entry values) were analysized as follows.

1. Effect of fractal dimension

In Figs.44, the fractal dimension, D was 2, 2.35 and 2.78, respectively. The air-entry value, ψe remained 10kPa. It is shown that the reduction of slope stability is larger in the condition of larger fractal dimension under the condition of the same rainfall intensities and the saturated hydraulic conductivity. This means that the soil slope is much safe if the fractal dimension is low while the saturated hydraulic conductivity remains constant. The factor of safety decreases in a short time under rainfall infiltration when the fractal dimension is relatively large, as shown in Fig.44c, while the decrease occurs at a longer time for the low fractal dimension soil type, as shown in Fig.44a.

Figure 44.

Effect of fractal dimension on slope stability

If the air-entry value and the saturated hydraulic conductivity remain unchanged, the slope of the hydraulic conductivity function is small when the fractal dimension is large. Then the initial hydraulic conductivity value at the beginning of the rainfall infiltration is much big, which results in the speed of rainfall infiltration is great and the saturated degree of soil inside of slope changes quickly leading unsaturated soil to become saturated. The passive pore water pressure generated by the increase of degree of soil saturation makes a negative impact on the instability of soil slope under rainfall.

2. Effect of air-entry value

The effects of air-entry value on the stability of soil slope due to rainfall infiltration are shown in Figs. 45. The air-entry value, ψe is 1kPa, 10kPa, 20kPa and 50kPa respectively and the fractal dimension, D remains 2.35. The style of the slope safety curve changes a lot with different air-entry values in Figs. 45. In Fig. 45a, the air-entry value, ψe is 1kPa and the decrease of factor of safety occurs during 16~20h after rainfall began. As the air-entry value, ψe is 10kPa, shown in Fig. 45b, the decrease of slope safety factor occurs at 8h after rainfall began. At the beginning of rainfall infiltration, the decrease happens in Figs. 45c-d, in which the air-entry value is 20kPa and 50kPa respectively. This means that a delay time appears which is depended with the air-entry value. And it is obvious that small air-entry value decides a relative long delay time. Besides, the factor of slope safety remains minimum during another 24 hours after rainfall stops in Fig. 45a, which has a low air-entry value (ψe=1kPa). While ψe=50kPa, the recovery of safety factor curve in Fig. 45c occurs and dramatically. This is also due to the style of hydraulic conductivity function, as shown in Fig. 35. The initial hydraulic conductivity value becomes large when the air-entry value increases with fractal dimension, saturated hydraulic conductivity and saturation degree remain the same.

3. Effect of saturated hydraulic conductivity

The effects of saturated hydraulic conductivity on stability of a homogenous soil slope are reflected through the relationship between Fs and ks, shown in Fig. 46. The saturated hydraulic conductivity, ks, was varied with five different values of 10-8m/s, 10-7m/s, 10-6m/s, 10-5m/s and 10-4m/s for a homogeneous soil slope of constant soil type S10, 2.35 (ψe=10kPa, D=2.35), Hw=2m, Hs=6m, α=33.7°, subjected to rainfall for 24h with six rainfall intensities of 3, 9, 15, 36, 100 and 360mm/h. All the plots in Fig. 47 shows that soil with ks values of 10-8m/s and 10-7m/s, respectively, are less affected by rainfall. Contrarily, soil with ks values of 10-5m/s and 10-4m/s are greatly affected by rainfall. It is suggested that soil slopes with a low saturated hydraulic conductivity are relatively safe under short-duration rainfall infiltration and for the soil slopes with a high saturated hydraulic conductivity the stability is affected by the short-duration rainfall greatly. This means that slopes with low saturated hydraulic conductivity need a long-duration rainfall to intrigue the instability.

The Fs(min)-ks critical curve was presented in the broken line in Fig. 47 as follows:

Fs(min)=a1+becksE67

where Fs(min) is minimum factor of safety; ks is saturated hydraulic conductivity; a, b and c are the fitting parameters and e is natural number (i.e., 2.71…). Here, a=1.293, b=-0.492, c=1.2×105, and the corresponding coefficient of correlation, r2 is 0.9998.

Figure 45.

Effect of air-entry value on the slope stability

8.2. Effect of rainfall intensity

The relationship between minimum factor of safety, Fs(min), versus logarithmic of rainfall intensity, Ir are plotted in Fig. 47. The semi log plot shows that generally the Fs(min) and Ir relationships follow a sigmoid shape in Fig. 47. There are two inflect points for the sigmoid shape line. The Fs(min) is almost constant at very low rainfall intensities for all soil types. And it starts to decrease rapidly when the first inflection point is reached. The trend in Fs(min) versus Ir relationship observed in Fig. 47 can be described by a sigmoid equation (MMF Line) as the form of

Fs(min)=ab+cIrdb+IrdE68

where Fs(min) is minimum factor of safety; Ir is rainfall intensity; and a, b, c, d are fitting parameters. The values for the fitting parameters and the corresponding coefficient of correlation, r2, for the sigmoid line in Fig.47.

Figure 46.

Relationship between Ks and minimum factor of safety

Figure 47.

Effect of rainfall intensity on variation of minimum factor of safety

8.3. Effect of initial water table location

The effects of initial groundwater table location on stability of a homogenous soil slope are reflected through the relationship between Fs(ini) and Fs(min) with Hw, shown in Fig. 48. The initial depth of water table, Hw, was varied with five different values of 0.5, 1, 2, 2.5, 5, 7.5, 10 and15m for a homogeneous soil slope of constant soil type S10, 2.35, -5 (ψe=10kPa, D=2.35, ks=10-5m/s), Hs=6m, α=33.7°, subjected to rainfall for 24h with five rainfall intensities of 3, 9, 15, 36, and 100mm/h. The relationship between Fs(ini), and Hw shown in Fig. 48 appears to be linear up to a depth of 7.5m beyond which Fs(ini) remains constant. This is because the initial pore-water pressure profiles generated for the slope at Hw=7.5m as same as those generated when Hw=7.5m. This is due to the limiting pore-water pressure of −75kPa adopted in the analyses. In Fig. 54 it also shows that the rainfall intensity is the dominated factor impacting on the reduction in factor of safety, Fs. And the initial depth of water table, Hw, mainly determines the value of initial factor of safety, Fs(ini). The Fs(ini) is smaller for slopes with a shallower Hw which means that the soil slope stability is much lower. Therefore, slopes with a shallow Hw are more likely to fail due to a rainfall compared with slopes which has a deep Hw.

Figure 48.

Effect of initial groundwater table depth on factor of safety

Figure 49.

Effect of slope angle on factor of safety

8.4. Effect of slope angle

The effect of slope geometry on the stability of a homogenous soil slope is evaluated in terms of slope angle (α) and slope height (Hs). The effects of slope angel on the stability of a homogenous soil slope are reflected through the relationship between Fs(ini) and Fs(min) with α, shown in Fig. 49. The slope angel, α was varied with three different values of 26°, 33.7°, 45°and 63° for a homogeneous soil slope of constant soil type S10, 2.35, -5 (ψe=10kPa, D=2.35, ks=10-5m/s), Hs=6m, Hw=2m, subjected to rainfall for 24h with five rainfall intensities of 3, 9, 15 and 36mm/h. The relationship between Fs(min) and α shown in Fig. 49 appears to be negative linear. In general, the higher the slope angle, the lower the initial factor of safety and the minimum factor of safety. Because a steep slope will yield a lower factor of safety compared with a flat slope.

8.5. Effect of slope height

The effects of slope angel on the stability of a homogenous soil slope are reflected through the relationship between Fs(ini) and Fs(min) with Hs, shown in Fig. 50. The slope height, Hs, was varied with four different values of 3m, 4m, 5m and 6m for a homogeneous soil slope of constant soil type S10, 2.35, -5(ψe=10kPa, D=2.35, ks=10-5m/s), Hw=2m, α=33.7°, subjected to rainfall for 24h with four rainfall intensities of 3mm/h, 9mm/h, 15mm/h and 36mm/h. Fig. 50 shows that initial factor of safety decreases exponentially as the slope height increases. It also suggests that high slopes are generally easier to fail under rainfall due to the low initial factor of safety. The reduction in factor of safety for a high slope is smaller and occurs at a slower rate compared with a low slope.

Figure 50.

Effect of slope height on factor of safety

Acknowledgments

The National Natural Science Foundation of China (Grant No. 41272318) and State Key Laboratory of Ocean Engineering are sincerely acknowledged for their financial support.

© 2013 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Yongfu Xu (December 4th 2013). Unsaturated Hydraulic Conductivity of Fractal-Textured Soils, Hydraulic Conductivity, Vanderlei Rodrigues da Silva, IntechOpen, DOI: 10.5772/56716. Available from:

chapter statistics

1494total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

Calibration of a New Device to Measure Water Content of Rocks

By Maria Clementina Caputo and Rita Masciale

Related Book

First chapter

Biodegradation of Medical Purpose Polymeric Materials and Their Impact on Biocompatibility

By Elisa Tamariz and Ariadna Rios-Ramírez

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us