Parameters obtained from the fractal model of PSD

## 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 *B* of *cover* of *B* is a family *U* of sets

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, *U*_{i}, where i ranges over the natural numbers. The diameter of a subset of *δ*-cover of *U*_{i})≤*δ*, where *δ* is a positive real number and where diam(.) is the diameter function on subsets of *δ*-cover are indicated by *U*_{i}^{δ}. The *s*-dimensional Hausdorff measure *H*_{s}(*B*) of *B* is defined as:

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 *B* is defined as:

The Hausdorff dimension of a subset of *n*. If dim_{H} *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 *S*: *λ*>0 such that, for all *x*, *y*

where |.| denotes the Euclidean distance in *s*, under a similarity transformation *S* with scale factor *λ* the *s*-dimensional Hausdorff measure transforms according to the formula:

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:

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

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:

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 *E*_{0} be the interval [0,1], *E*_{1} is the set obtained by deleting the middle third of *E*_{0} so *E*_{1} consists of the two intervals *E*_{1} gives *E*_{2}, thus *E*_{2} comprises the four intervals *E*_{i} is obtained by deleting the middle third of each interval in *E*_{i-1} and *E*i consists of 2^{i} intervals with length of 3^{-i}. Using Eq. (8), the fractal dimension of the middle-third Cantor set is

## 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/*a*^{3} 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/*a*^{3}-*N* subcubes are then each divided into 1/*a*^{3} 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 *a*^{3} needed to cover the pores equal to or larger than *a*^{3} 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

When we covered a soil pore surface using balls with the same radius, the soil pore surface is covered by *N*_{1} balls with radius *r*_{1}, the same surface is covered by *N*_{i} balls with radius *r*_{i}, and it is covered by *N*_{j} balls with radius *r*_{j}, and usually *r*_{1}>*r*_{i}>*r*_{j}, it is found that the surface area of the soil pore surface contents has the following relationship:

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:

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

where in 11 a) *A*_{p}(≤*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 *V*_{p}(≤*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:

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.,

where *A*_{p}(*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 *A*_{p}(*p*), the surface fractal dimension of soil pores is given by

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 log*r* vs. log*V*_{p}. If the slope of the fitting line in the log*r*-log*V*_{p} plane is *κ*, the fractal dimension of the pore surface is written as

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 *V*_{v} 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 *V*_{p}(≤*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 log*r* vs. log(*V*_{p}(≤*r*)/*V*_{v}) 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 *V*_{p}(≤*r*)/*V*_{v}=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-02 | 0.34 | 2.66 | -0.34 | 0.0125 | 12 | Watabe et al., 2000 |

S-03 | 0.37 | 2.63 | -0.37 | 0.06 | 2.5 | |

S-04 | 0.49 | 2.51 | -0.49 | 0.006 | 25 | |

A horizon | 0.22 | 2.78 | -0.22 | 0.2 | 0.75 | Vogel and Roth, 1998 |

B horizon | 0.13 | 2.87 | -0.13 | 0.25 | 0.6 | |

Toyoura sand | 1.33 | 1.67 | -1.33 | 0.09 | 1.67 | Uno et al., 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 2*r* were removed. The cumulative volume *V*_{p}(>*r*)/*V*_{T} was obtained using the application of successively larger structuring elements by Vogel and Roth (1998). The following relationship is obtained between *V*_{p}(>*r*)/*V*_{T} and *V*_{p}(≤*r*)/*V*_{T}.

where *V*_{T} is the total volume of soil, *V*_{v} 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*)/*V*_{T} can be obtained from the cumulative volume *V*(>*r*)/*V*_{T}. The cumulative volume of *V*(≤*r*)/*V*_{T} can be translated into *V*(≤*r*)/*V*_{v} according to Eq.(16). The cumulative volume of *V*(≤*r*)/*V*_{v} 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*)/*V*_{v}=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., *V*_{p}/*V*_{v}=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.

### 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

where *V*_{w} is the water volume absorbed by clay, *V*_{m} 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 *P*_{0} 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%.

where *p*_{s} is swelling pressure, the maximum axial pressure which is needed to maintain the original sample height, *T* is the Kelvin temperature, *M*_{w} 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

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, *C*_{b} 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 (*e*_{c}) and vertical pressure is shown in Fig. 8. Relationship between the clay void ratio (*e*_{c}) and normalized water volume is given by

where *C*_{m} is the montmorillonite content, and *C*_{m}=75%. It is seen that the slope of log*e*_{c}-log*p* equals to that of log(*V*_{w}/*V*_{m})-log*p* from Eq. (20) for *C*_{m} being constant. Thus, the surface fractal dimension of Wyoming bentonite can be obtained from the linear relationship of log*e*_{c}-log*p*. The surface fractal dimensions of Wyoming bentonite are 2.64, which obtained from the swelling deformation tests conducted by different authors.

### 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 *N*_{g}, the pore radius *r* is related to the radius *r*_{g} of the smallest grain as follows (Watabe et al., 2000)

where *f*_{p} is a constant less than 1.0. If the particles are randomly arranged, the probability of having a grain of radius less than *r*_{g} is d*P*_{g}, the probability of having the *N*_{g}-1 grains with a radius greater than *r*_{g} is (1-*P*_{g})^{Ng-1}. The probability *P* of having a pore with a radius less than *r* is given by (Watabe et al., 2000)

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

where *D*_{g} is the fractal dimension of the grain-size distribution (GSD), *R*_{g} is the maximum radius of soil grains. The probability *P* is given by (Watabe et al., 2000)

If *r*/(*f*_{g}*R*_{g})<<1, Eq.(24) can be written as

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*)/*V*_{v} and *P*_{g}=*M*(≤*r*)/*M*_{T}, *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.

## 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*→d*r* to the water content is given by

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

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

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

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

where *θ*, *θ*_{r} and *θ*_{s} are the volumetric water content, residual volumetric water content and saturated volumetric water content, *S* and *S*_{r} 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 *S*_{r}=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.

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.

Sample | Method | Bulk density (g/cm3) | θs | θr | D | R (mm) | ψe (kPa) |

FH31 | Pressure plate | 1.58 | 0.32 | 0.02 | 0.85 | 0.065 | 2.31 |

W3 | MSO | 1.42 | 0.33 | 0.058 | 1.51 | 0.0175 | 8.6 |

Mix8 | Rosetta | 1.45 | 0.41 | 0.061 | 1.52 | 0.00053 | 290 |

MZ | Pressure plate | 1.60 | 0.44 | 0 | 2.29 | 0.00035 | 429 |

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.

## 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):

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.

*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.

*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):

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=*(*L*_{e}/*L*)^{2}, in which *L*_{e} 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):

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

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

*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

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):

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)

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 d*q* flowing through the connections between pores of radius *r*_{1}→*r*_{1}+d*r*_{1} on one side of the slab (at *x*) and pores of radius *r*_{2}→*r*_{2}+d*r*_{2} on the other side of the slab (at *x*+d*x*) can be expressed by (Mualem, 1976)

where *β* is a constant which incorporates the fluid with the matrix properties; *r*_{e}(*r*_{1}, *r*_{2}, *q*) is the effective radius; *A*_{e}(*r*_{1}, *r*_{2}, *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 *k*_{w} is given by

The relatively hydraulic conductivity (*k*_{r}) is the ratio of the hydraulic conductivity at any volumetric water content (*k*_{w}) to the hydraulic conductivity at saturation (*k*_{s}). The relative hydraulic conductivity (RHC) is related to the radius of the soil pores, and is written as (Mualem, 1976)

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 *r*_{1}→*r*_{1}+d*r*_{1} at *x* to encounter pores of radius *r*_{2}→*r*_{2}+d*r*_{2} at *x*+d*x* is given by

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)

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

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

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:

where *d* is the soil-water diffusivity (SWD), *k*_{w} is the saturated hydraulic conductivity. Substituting Eqs. (44) and (45) into Eq. (46), the soil-water diffusivity (SWD) is written as

where *ς*=2*D*-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 32^{3}=32768 nodes. The hydraulic conductivity, *k*_{w}, is determined at each step of desaturation by imposing a pressure gradient ∆*p* across the ends of the network. Water flow *q*_{ij} 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*/(*A*∆*p*), 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.

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.

**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 *k*_{s} 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 type | D | δ | ξ | ξ/δ | ς | ψe (kPa) |

Haploxeroll loam | 2.63 | -0.37 | -3.11 | 8.41 | -1.74 | 0.14 |

Toyoura sand (wetting) | 1.60 | -1.40 | -6.2 | 4.43 | -3.8 | / |

Toyoura sand (drying) | 1.26 | -1.74 | -7.22 | 4.15 | -4.48 | / |

Loamy sand | 2.68 | -0.32 | -2.96 | 9.25 | -1.64 | 2.0 |

McGee Ranch soil | 2.51(GSD) | -0.49 | -3.47 | 7.08 | -1.98 | 4.6 |

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).

Soil | θs | Sr (%) | D | ψe (kPa) | Reference |

Lakeland | 0.375 | 30-100 | 1.75 | 2.5 | Elzeftawy and Cartwright, 1981 |

Superstition | 0.5 | 30-100 | 2.7 | 1.2 | Richards, 1952 |

Columbia sandy loam | 0.458 | 50-100 | 2.29 | 4.5 | Brooks and Corey, 1964 |

Touchet silt loam | 0.43 | 20-100 | 2.09 | 7 | Brooks and Corey, 1964 |

Silt loam | 0.396 | 50-100 | 2.59 | 15 | Reisenauer, 1963 |

Guelph loam | 0.52 | 45-100 | 2.75 | 3 | Elrick and Bowmann, 1964 |

Yolo light clay | 0.375 | 45-100 | 2.87 | 1.5 | Moore, 1939 |

Speswhite kaolin | 0.56 | 55-100 | 2.85 | 15 | Peroni et al., 2003 |

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).

## 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):

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)

where *k*_{w} is the unsaturated hydraulic conductivity at any water content, *k*_{s} 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,

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

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

where *m*=1-2/*n*.

Note that for both Eqs. (52) and (54), because *m* is less than 1, the derivative d*k*/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

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

where *a*_{1}=*m*+(2+*c*)/*n*, *a*_{2}=2*m*+(2+*c*)/*n*, *b*=1-(2+*c*)/*n*, *I*_{x}(*a*, *b*) and *B*(*a*, *b*) are the Incomplete Beta and Beta functions of positive arguments *a* and *b*, respectively, and given by:

The results presented by Touma (2009) were obtained with *c*=0.5 and *a*_{1}=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 *k*_{r} 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

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

when BC is applied with the Burdine condition, and

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.

Soil type | Fractal model | G-M model | ||||

δ | D | ξ | ξ/δ | α (kPa-1) | n | |

Beit Netofa clay | -0.11 | 2.89 | -2.33 | 21.2 | 0.004 | 1.17 |

Hanford glacial soil 0-099 | -0.62 | 2.38 | -3.86 | 6.23 | 0.015 | 1.71 |

Hanford glacial soil 2-1637 | -0.82 | 2.18 | -4.46 | 5.44 | 0.076 | 1.89 |

Grenoble sand | -0.64 | 2.36 | -3.92 | 6.13 | 0.4 | 2.17 |

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.

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.

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.

## 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:

where *H*_{w} is the total water energy potential comprised of both pressure and elevation potentials; *P*_{a} is the pore air pressure; *k*_{w} is hydraulic conductivity; *y* is y coordinate; *m*_{w} is the slope of SWCC; *t* is time; and *Q*_{w} 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:

where *k*_{a} is air permeability; *θ*_{a} is volumetric air content; *ρ*_{a} is air density; *γ*_{oa} is initial air unit weight; *ρ*_{oa} is initial air density; *T* is temperature; *H*_{w} is total head and *P*_{a} 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 *k*_{s}=1×10^{-5}m/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).

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 *k*_{s} 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 *k*_{s}, 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):

where *τ* is the shear strength; *c*' is the effective cohesion; (*σ*-*u*_{a}) 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):

where *d*_{cr} is the critical depth; *β* is slope angle; *h*_{p} is the pore passive pressure and *h*_{c} 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:

The critical depth for infinite slope failure is now a function of the suction, *ψ*, pressure head *h*_{p}, 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/m^{3}, *γ*_{w}=10kN/m^{3}, 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.

### 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 (*k*_{s}=1×10^{-5}m/s, 1×10^{-6}m/s, 1×10^{-7}m/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.

### 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 *k*_{s} equals to 1×10^{-5}m/s, and for the situation that *k*_{s} equals to 1×10^{-6}m/s, the time can range from 10 to 100 hours. While much more time is needed for the lower saturated hydraulic conductivity soil type.

### 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*=*I*_{r}(rainfall intensity); ah, de, fg=*q*=0m^{3}/s (i.e. no flow boundary); and ef, gh=*h*_{t}(total head at the side). Four slope heights, *H*_{s},(3m, 4m, 5m and 6m), three slope angles, *α*,(26.6, 33.7 and 45.0), eight initial depths of groundwater table (GWT), *H*_{w},(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, *k*_{s} (10^{-4}m/s, 10^{-5}m/s, 10^{-6}m/s, 10^{-7}m/s and 10^{-8}m/s) and five rainfall intensities *I*_{r}(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/m^{3}.

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:

where *H* is the total head; *k*_{x} are *k*_{y} ars hydraulic conductivity in x, y direction respectively；*Q* is the applied boundary flux; *m*_{w} 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, *F*_{s}.

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, *u*_{w}, 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 *H*_{s}=6m, initial groundwater table depth *H*_{w}*=*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.

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 *F*_{s} and *k*_{s}, shown in Fig. 46. The saturated hydraulic conductivity, *k*_{s}, was varied with five different values of 10^{-8}m/s, 10^{-7}m/s, 10^{-6}m/s, 10^{-5}m/s and 10^{-4}m/s for a homogeneous soil slope of constant soil type S10, 2.35 (*ψ*_{e}=10kPa, *D*=2.35), *H*_{w}=2m, *H*_{s}=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 *k*_{s} values of 10^{-8}m/s and 10^{-7}m/s, respectively, are less affected by rainfall. Contrarily, soil with *k*_{s} values of 10^{-5}m/s and 10^{-4}m/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 *F*_{s}(min)-*k*_{s} critical curve was presented in the broken line in Fig. 47 as follows:

where *F*_{s}(min) is minimum factor of safety; *k*_{s} 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×10^{5}, and the corresponding coefficient of correlation, r^{2} is 0.9998.

## 8.2. Effect of rainfall intensity

The relationship between minimum factor of safety, *F*_{s}(min), versus logarithmic of rainfall intensity, *I*_{r} are plotted in Fig. 47. The semi log plot shows that generally the *F*_{s}(min) and *I*_{r} relationships follow a sigmoid shape in Fig. 47. There are two inflect points for the sigmoid shape line. The *F*_{s}(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 *F*_{s}(min) versus *I*_{r} relationship observed in Fig. 47 can be described by a sigmoid equation (MMF Line) as the form of

where *F*_{s}(min) is minimum factor of safety; *I*_{r} is rainfall intensity; and *a, b, c, d* are fitting parameters. The values for the fitting parameters and the corresponding coefficient of correlation, r^{2}, for the sigmoid line in Fig.47.

## 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 *F*_{s}(ini) and *F*_{s}(min) with *H*_{w}, shown in Fig. 48. The initial depth of water table, *H*_{w}, 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, *k*_{s}=10^{-5}m/s), *H*_{s}=6m, *α*=33.7°, subjected to rainfall for 24h with five rainfall intensities of 3, 9, 15, 36, and 100mm/h. The relationship between *F*_{s}(ini), and *H*_{w} shown in Fig. 48 appears to be linear up to a depth of 7.5m beyond which *F*_{s}(ini) remains constant. This is because the initial pore-water pressure profiles generated for the slope at *H*_{w}=7.5m as same as those generated when *H*_{w}=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, *F*s. And the initial depth of water table, *H*_{w}, mainly determines the value of initial factor of safety, *F*_{s}(ini). The *F*_{s}(ini) is smaller for slopes with a shallower *H*_{w} which means that the soil slope stability is much lower. Therefore, slopes with a shallow *H*_{w} are more likely to fail due to a rainfall compared with slopes which has a deep *H*_{w}.

## 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 (*H*_{s}). The effects of slope angel on the stability of a homogenous soil slope are reflected through the relationship between *F*_{s}(ini) and *F*_{s}(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, *k*_{s}=10^{-5}m/s), *H*_{s}=6m, *H*_{w}=2m, subjected to rainfall for 24h with five rainfall intensities of 3, 9, 15 and 36mm/h. The relationship between *F*_{s}(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 *F*_{s}(ini) and *F*_{s}(min) with *H*_{s}, shown in Fig. 50. The slope height, *H*_{s}, 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, *k*_{s}=10^{-5}m/s), *H*_{w}=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.

## 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.