Table of Contents Author Guidelines Submit a Manuscript
Advances in Materials Science and Engineering
Volume 2016, Article ID 2684297, 13 pages
Research Article

Unsaturated Seepage Analysis of Cracked Soil including Development Process of Cracks

1School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiaotong University, Shanghai 200240, China
2College of Civil Engineering & Architecture, China Three Gorges University, Yichang 443002, China

Received 5 November 2015; Revised 17 February 2016; Accepted 21 February 2016

Academic Editor: Giorgio Pia

Copyright © 2016 Ling Cao et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


Cracks in soil provide preferential pathways for water flow and their morphological parameters significantly affect the hydraulic conductivity of the soil. To study the hydraulic properties of cracks, the dynamic development of cracks in the expansive soil during drying and wetting has been measured in the laboratory. The test results enable the development of the relationships between the cracks morphological parameters and the water content. In this study, the fractal model has been used to predict the soil-water characteristic curve (SWCC) of the cracked soil, including the developmental process of the cracks. The cracked expansive soil has been considered as a crack-pore medium. A dual media flow model has been developed to simulate the seepage characteristics of the cracked expansive soil. The variations in pore water pressure at different part of the model are quite different due to the impact of the cracks. This study proves that seepage characteristics can be better predicted if the impact of cracks is taken into account.

1. Introduction

Cracked soils are common in nature, particularly with expansive soils. The cracks break the integrality of the soil, thereby reducing its strength. The pathways provided by cracks significantly increase the hydraulic conductivity of the soil. As such, it is easier for cracked soil slopes to slip or collapse during rainy seasons [1, 2]. Therefore, in geotechnical engineering, much attention has been given to cracked soil seepage analysis [3, 4]. Most slopes are not affected by the infiltration of rainwater if there are no cracks. However, for slopes with cracks, it is necessary to determine the stability by including the effects of cracks when carrying out saturated-unsaturated seepage analysis. This has been done in previous studies [36].

Over the last few decades, many researchers studied the SWCC and used it to predict the soil hydraulic conductivity. However, for cracked soil, it is a crack-pore dual medium. Zhang and Fredlund [5] showed that the SWCCs for unsaturated cracked rocks depend on the pore distribution. They assumed that both the crack networks system and the pore network system follow lognormal distribution. The SWCC for the cracked rock can then be obtained by combining the estimated SWCCs of the two systems. The double-peak specificity of this SWCC is similar to certain mathematical models developed by other researchers [6]. However, this theory is based on rock mechanics, which means that they assumed that the crack width is invariable. This is a complete contrast to the characteristic of soil cracks. The morphological character of cracks changes with the water content during the repetitive rainfall-evaporation process. Therefore, the hydraulic conductivity of this deformable medium is not a constant. The conductivity changes with the morphological character of the cracks and the suction of the soil. Hence, to determine the hydraulic conductivity of cracked soil, it is necessary to consider the change in crack volume in the analysis of unsaturated seepage. The main objectives of this research are to determine the mechanisms of crack initiation and development under natural conditions, to develop a method to determine the SWCCs of an unsaturated cracked soil, and to use the SWCCs to analyze the unsaturated seepage based on a proper flow model.

The cracked soil can be represented as an overlapping continuum of pore media and crack media. As such, in this study, the pore-crack dual media model has been used to simulate the unsaturated seepage in the cracked soil. The fractal model has been used to describe the SWCC of the dual media, and normal SWCC has been used for the pore media. However, for the crack media, its morphological parameters are not constants. To solve this problem, laboratory experiment has been designed to simulate the dynamic changes of the soil cracks during the rainfall-evaporation cyclic processes. The morphological characteristics of the cracks during the drying-wetting processes can then be obtained from the experiment. Thereafter, the relationships between the morphological parameters, such as fractal dimension, surface crack ratio, and the water content, can be determined. With fractal theory, the SWCC for the crack media has been developed on the basis of the experiment. Finally, the pore-crack dual media seepage model has been used to simulate the seepage process due to ponding in the cracked soil.

2. Dynamic Development of Cracks in Drying and Wetting Cycles

Many laboratory and field experiments have been conducted to study the development of desiccation cracks [79]. The experimental results showed that the cracks development process can be divided into three stages. In Stage I, the water content decreases slowly with only a few cracks developing. In Stage II, the development of the cracks accelerates, while the water content continues to decrease. In Stage III, the development of the cracks tends towards equilibrium. On the other hand, with the increase in water content, the crack width decreases during the wetting process. In highly plastic soil, the cracks may even close up. In lowly plastic soil, the cracks may remain open.

To study the mechanisms of crack initiation and development during repetitive drying-wetting cycles, a laboratory experiment has been conducted. This experiment focuses on the development of surface cracks. The relationships between the crack morphological parameters and the water content during the drying-wetting cycle have been developed.

2.1. Experimental Material

The soil sample has been taken from a highway construction site in Nanjing. The basic physical properties of the soil are summarized in Table 1. The soil has weak swell ability.

Table 1: Basic physical properties of soil sample.
2.2. Experimental Equipment

A set of experimental equipment has been designed to simulate the cracks evolution process of expansive soil under the iterative rainfall-evaporation cycles, as shown in Figure 1. In this experiment, two identical samples with the same initial water content, dry density, and dimensions (29.5 cm × 39.5 cm × 3 cm) have been prepared, one for capturing images of the cracks and one for the measurement of water content.

Figure 1: Sketch of the experiment equipment for cracking test.

Using a fixed camera, photographs were taken at the appointed time so as to record the evolution progress of the soil cracks. During this photographic session, all the curtains in the laboratory were drawn and only the indoor light was on. In this way, the brightness of the light was the same throughout the session.

2.3. Experimental Method

The dry soil was first blended with water with initial water content of 32%, which is slightly greater than the plastic limit. The mixture was then sealed in a plastic bag and put inside a moisturizing cylinder for 24 hours so as to ensure the moisture was evenly distributed.

Based on the designed dimensions of the soil sample and the dry density ( g/cm3 which is close to the natural dry density), the weight of the soil mixture was calculated. Half of the mixture was then poured into a sample mould (29.5 cm × 39.5 cm × 10 cm). Using a jack, the soil mixture was compressed into the designed dry density. The soil sample surface was then roughened with a brush and the second layer was compacted using the same method.

In order to simulate the evaporation process and to accelerate the development of cracks, electric fan was used to generate a steady wind. The drying process was ceased when the weight of the soil sample was constant for three hours. This is the end of the initial drying process.

Then, the electric fan was turned off. To simulate the rainfall process, water was sprayed onto the dry and cracked soil samples. The wetting process was ceased when the soil surface was covered with water and the cracks were stable. This is the end of the first wetting-drying cycle, and the cycle was repeated with the procedures as described in this and the previous paragraphs.

The wetting-drying cycle process has been repeated six times, so as to simulate the continuous rainfall-evaporation process under the natural condition.

2.4. Experimental Result

Figure 2 shows the images of the cracks at the end of each drying process. After the first drying process, only a few major cracks appeared. With subsequent cycles, more cracks appeared but the crack width became smaller. The number of cracks did not change substantially after five drying-wetting cycles. Thus, it may be inferred that the effect of drying-wetting cycles on soil cracking is limited during the first few cycles. Hence, in this study, only the results from the fifth cycle are considered to represent the long-term effect of evaporation-rainfall cycle on the soil cracks.

Figure 2: Cracks images at each end of drying process.

Figures 3 and 4 show the cracks development during the sixth drying-wetting cycle. The crack porosity changes with the water content obviously.

Figure 3: Crack development during a drying process.
Figure 4: Crack development during a wetting process.

To quantitatively analyze the morphological parameters of the cracks at different stages, a MATLAB-based photoprocessing program was used. After cutting, binarization, and bridging, the photographs were converted into binarization image of the same size as the original cracks, as shown in Figure 5. Thereafter, the program quantified the parameters by counting the number of black pixels (crack area) and white pixels (soil area).

Figure 5: Binarization image of cracks.

Using the crack porosity, , as an indicator of the total distribution characteristics, the equation is where is the area of the th crack in the soil sample surface, which can be obtained by the “bwarea” function of MATLAB; is the number of cracks; is the total area of the soil sample surface; is the number of black pixels, which is the cracking part of the sample; is the number of white pixels, which is the noncracking part; and is the total number of pixels.

The width indicates the opening degree of the cracks. It is not uniform along the length of every crack. So, in this study, is the average width, which can be determined from the area and the length of the crack.

Box-counting dimension method [10] is the most popular method for estimating the crack fractal dimension . The procedure is as follows: using a grid consisting of square elements with sides of dimension to cover the binarization image. The number of elements that have black pixel (crack part), , changes with . The relationship between and iswhere is a constant. Hence, values can be derived from the slop of versus plots.

Figure 6 shows the changes in crack porosity, crack fractal dimension, and the maximum crack width with the water content. The functional equations between the morphological parameters and the water content can be derived.

Figure 6: Changes of crack morphological parameters with water content. (a) Crack porosity . (b) Crack fractal dimension . (c) The maximum crack width .

3. The SWCCs and Permeability Function for Cracked Soil

3.1. The SWCCs and Permeability Function for Soil

As the water flow in cracked soil is dependent on the pore-crack structure, the properties of the pore and the crack structure systems can be quantitatively described by the fractal model. Studies [1113] have shown that the spatial distribution of the pore and the crack structure is self-similar at certain scales. Hence, based on the fractal model, the SWCCs and permeability functions for the pore medium are derived from their respective volume-size distribution, as follows:where is the air-entry value; is the effective degree of saturation; is the volumetric water content of unsaturated soil at suction ; is the volumetric water content of saturated soil; is the soil fractal dimension; is the relative coefficient of permeability of unsaturated soil; is the saturated coefficient of permeability; is the coefficient of permeability at suction . For a two-dimensional space, change to in (3) and to in (4).

The drying water retention curve of a soil can be obtained through the SWCC test using a pressure plate extractor. To get a SWCC in a larger range of pressure, the experimental studies of the drying water retention curve have been carried out with the pressure plate extractor (5 Bar1D-SDSWCC) and the temper pressure membrane apparatus (15 Bar LAB523), as shown in Figure 7. The first stage experiment between 0 and 1 Bar was carried out with the pressure plate extractor. The soil sample was then moved to the temper pressure membrane apparatus after the suction became stable. Finally, the experiment was completed between 1 and 8 Bar.

Figure 7: Apparatus for the measurement of SWCC.

By fitting (3) to the measured drying SWCC gives the drying SWCC over the entire suction range. Figure 7 shows the fitted drying curve. Pham et al. (2005) [14] showed that the distance between the wetting and the drying curves is 20% log-cycle, and the slope of the wetting water retention curve is the same as that of the drying curve. Based on these findings, the wetting water retention curve for the crack network has also been predicted, as shown in Figure 8.

Figure 8: SWCC for soil.
3.2. The SWCCs and Permeability Function for Crack

For cracked soil, the hydraulic properties are closely related to the pore-size distribution. A cracked soil can be considered as a dual medium of porosity consisting of crack sand soil pores. For the cracks part, the pore-size distribution refers to the distribution of crack volume with respect to the crack aperture. In this study, the cracks are assumed not to vary with depth. So, the pore-size distribution of the cracks can be obtained by the surface crack morphological parameters, such as crack porosity , crack fractal dimension , and the crack width .

Although (3) and (4) were developed for soil, they should also be valid for any porous medium that behaves in accordance with the capillary law. Reitsma and Kueper (1994) [15] showed that the width distribution in a rough-walled crack is homologous to the distribution of pore throat radii in a porous medium. Thus, in this study, (3) and (4) are considered applicable to the water retention characteristics of cracks in soils.

According to the Laplace Equation, the relationship between the suction and the width of the cracks can be written as where is the surface tension, is the contact angle, is the suction, and is the crack width. Substituting (5) into (3), the equation becomeswhere is the maximum crack width, is the volumetric water content of crack at suction , is the volumetric water content of saturated soil, which equals the crack porosity numerically, and is the crack fractal dimension. For a two-dimensional space, (6) becomesGiven that  N/m at 20°C and °, (7) can be rewritten aswhere , , and are the function equations between the morphological parameters , , and and the water content, respectively. They can be derived on the basis of the curves shown as in Figure 6.

The crack SWCCs at various water contents can be obtained based on relationships between the water content and the crack morphological parameters such as , , and , as shown in Figure 9.

Figure 9: The SWCC for the crack network at various water contents.

4. The Crack-Pore Dual Media Model

According to the media identity, the crack seepage models fall into three main groups: the continuum models, the discrete medium models, and the dual media model. Among these three groups, the dual media model is the most suitable for studying the changes in a seepage field with time [16]. In this model, cracks may be represented by 2D elements in 3D problem or 1D elements in 2D problem. These two meshes are embedded and the overlapped elements share the same nodes with the same hydraulic head. When soil properties, initial condition, and boundary condition for matrix and cracks are given, the flow can be simulated. Many researchers have investigated the flow in an unsaturated cracked porous medium or cracked rock mass [1722].

4.1. The Control Equation

The saturated-unsaturated seepage in a porous medium is governed by Darcy’s Law. The saturated flow through a fracture is described by Hagen-Poiseuille equation [23, 24]. Under the laminar flow conditions, the water flux within a small gap, , between two parallel plates iswhere is the hydraulic head, is the distance, is the hydraulic conductivity of flow through a fracture, and it is defined aswhere is the fluid density, is the gravitational acceleration, and is the dynamic viscosity of the fluid.

Equation (9) shows that the saturated flow in crack is also in the form as Darcy’s Law. Hence, it can be solved by Richards’ equation. If the compression of water phase can be neglected, the hydraulic head becomes the main variable and the 2D-control equation can be expressed aswhere , and it is the slope of SWCC; is the pore water pressure; is the unit weight of water. In the saturated zone, and are the respective saturated hydraulic conductivity in - and -directions, is the saturated volumetric water content, and equals the coefficient of compressibility of soil for one-dimensional consolidation.

4.2. Finite Element Modeling

In the dual media model, the interpolating functions for the 1D and 2D elements in terms of local coordinates arewhere the local coordinates and at node are −1 and 1, respectively. If node is located on the negative side, the coordinate or is −1; otherwise or is 1.

In the dual medium model, the 1D cracks are embedded into the 2D soil mass. So, the water flow velocity vector in the 1D crack is in the same direction as the crack. However, the velocity vector in the global 2D coordinate system still has two components (i.e., - and -directions). In a crack element, the hydraulic head distribution is defined by the function ofwhere and are the hydraulic heads of the nodes at the ends of a linear crack of length ; and are the nodal locations. In the local coordinate system, and when and , respectively.

According to Darcy’s Law, the derivative of (13) is

Then the gradient in an element can be written as

The flow velocity in crack is

Considering the velocity vector resolution, when the angle between the crack and horizontal direction is , the velocity components at and directions are

Although the overlapped 1D and 2D elements share the same node number with the same hydraulic head, the material parameters are different. The 1D and 2D element stiff matrixes are assembled into a single total stiff matrix. Then, the model is solved by a code written in FORTRAN, which is updated from the unsaturated subsurface flow model by Tung et al. [25].

5. Use of Crack-Pore Dual Media Model to Simulate Seepage

Using a flowerpot of soil with water ponding on the soil surface as an example, the numerical simulation of a crack-pore dual media model is explained in this section.

5.1. Computational Model

The computational model consists of a 0.3 m × 0.3 m area, as shown in Figure 7. Within the area, it contains 400 2D soil elements of 0.015 m × 0.015 m, and 92 1D crack elements of 0.005 m width and 0.015 m length. The width of the crack element changes with water content as shown in Figure 6(c). Two sections (AA′ and BB′) have been selected to show the distribution of pore water pressure as shown in Figure 10(c). Section AA′ is 0.015 m from the left boundary, and Section BB′ is at the centre of the model.

Figure 10: Grid map of computational model of soil containing cracks.

The fractal dimension for the soil has been measured in the laboratory using mercury injection technique, which is 2.78 for 3D calculations and 1.78 for 2D calculations. As for the crack, Figure 6(b) shows the fractal dimension.

5.2. Initial and Boundary Conditions for Model Calculations

Both the left and the right boundaries have been modeled as impervious boundaries. The top boundary has been modeled in two stages. For the first stage, the water head has been set at 0.1 m for the first five minutes. For the second stage, the top water head has been set to zero for the next 180 minutes.

At the bottom of the model, there is a 0.06 m diameter hole in the centre. All the other boundaries have been modeled as impervious except for the hole. Similar to the boundary at the top, the boundary conditions for the hole have also been modeled in two stages. For the first stage, the flux boundary and the water head are the -coordinate. For the second stage,  m.

At the start of the calculations, the water level is at the bottom of the pot, and the maximum negative pore pressure is 8 m. Figures 11(a) and 11(b) show the pore water pressure contours at the end of the first stage without and with cracks in the soil.

Figure 11: Pore water pressure contours after 5 minutes without or with cracks.

A comparison of these two figures shows that the cracks have significant impact on the infiltration through the soil. In Figure 11(a), the water moves very slowly and stays at the top of the pot. In Figure 11(b), the water flows out of the hole at the bottom of the pot. Further, the shape of the pore water pressure contours is similar to that of the cracks. In the cracked soil, the water quickly flows to the bottom because there is not enough time for the soil around the cracks to absorb the water. Hence, the pore water pressure of the soil remains negative. As time goes on, the water is absorbed slowly and then rises due to the capillary forces. The pore water pressure then rises towards equilibrium. Figure 12 shows the pore water pressure contours in the cracked soil at the end of the second stage.

Figure 12: Pore water pressure contours after 185 minutes (with cracks).

Figure 13 shows the variation of the pore water pressure with time at the nodes and at various depths along Section AA′. The initial pore water pressure at the top is −80 kPa, and it quickly increases to 0 kPa as the ponded water starts to infiltrate. After two minutes, the pore water pressure of the soil above 0.018 m from the bottom increases with height and changes from 0 kPa to −78.5 kPa. However, for the soil below 0.018 m, there is little change in pore water pressure because the wetting front has not reached there. After five minutes, the changes in pore water pressure are mainly within the soil above 0.0165 m, as the wetting front continues to move downwards. In the second stage, since there is no hydraulic head at the top, the water moves deep into the soil through the cracks. For this part of the soil, the pore water pressure decreases to negative. At the same time, the soil near the bottom absorbs water under the capillary forces, which results in a balance in pressure. Hence, at about 65 min, the pore water pressure distribution curve is practically a straight line. With further infiltration, the pore water pressure begins to increase. After about 185 min, this pressure stabilizes at around 0 kPa.

Figure 13: Variation of pore water pressure with time along Section AA′.

Figure 14 shows the variation of the pore water pressure with time at the nodes and at various depths along Section BB′. Similar to Section AA′, the pore water pressure at the top also quickly increases to 0 kPa as the ponded water starts to infiltrate. However, the variation of the wetting front in this section is not as regular as that in Section AA′ because Section BB′ is located at the centre of the model. It is therefore affected by the cracks on both sides of the section. Anyway, the general trend is downwards. After the first 5 minutes, the pore water pressure reverses from negative to positive. The pore water pressure then gradually decreases and becomes negative again at about 25 min due to subsequent drainage. As the soil absorbs water under capillary forces, the absolute pressure decreases (e.g., at 105 min) to zero at about 185 min.

Figure 14: Variation of pore water pressure with time along Section BB′.

Figure 15 shows the pore water pressure at Points I, II, III, and IV. These points are used to examine the distribution of the pore water pressure. Point I is at the lower part of the soil and is relatively far away from the cracks. Thus, the pore water pressure at Point I is not seriously affected by the quick infiltration of the water in the cracks during the ponding infiltration process. The pressure at Point I keeps a negative value around −80 kPa. The next 180 mins is drainage and capillary process, so capillary dominates the pore water pressure changes of Point I and the pressure rises with time until the value is close to 0 kPa.

Figure 15: Variations of pore water pressure with time at key points.

The other three points are relatively near the cracks, so the high permeability of crack has a greater influence on the pore water pressure of these points. Point II is in the upper part of the soil, along the same section as Point I. Hence, it responds quicker to the change in pore water pressure as compared to Point I. However, Point II is only affected by the crack on the right side, because the boundary on its left side is impervious. After the first five minutes, the pore water pressure at Point II increases from −78.5 kPa to −3.43 kPa. Point III is at the same depth as Point II. Since Point III lies between two cracks on the left part of the model, after the first five minutes, the pore water pressure increases from −78.5 kPa to 0.441 kPa. Point IV is at the same depth as Point I with cracks on both the left and the right, and the cracks are symmetrical. Hence, its response to the change in pore water pressure is faster than that at Point I. After the first five minutes, the pore water pressure at Point IV increases from −78.5 kPa to −2.21 kPa. However, the response at Point IV is slower than that at Point III because Point IV is at the lower part of the soil.

6. Conclusion

A laboratory experimental study has been conducted to gain insight into the mechanisms of crack development in soil under the dry-wet cycle. The soil specimens were weighed and photographed at a regular interval throughout the experiment. A MATLAB-based photoprocessing program has been used to convert the photographs of the cracks in soil into binarization images. From these images, the parameters, such as crack porosity and crack width, have been calculated.

The relationships between the crack geometric parameters and the water content are presented in this paper. Moreover, using the fractal model, the SWCC for the crack development process during the drying-wetting cycles has been developed.

The cracked soil has been separated into two systems: a crack network system and a soil matrix system. Using the crack-pore dual media model, infiltration of ponded water into the cracked soil has been simulated. The simulation results show that cracks have a large impact on the seepage field. The pathways provided by cracks significantly increase the hydraulic conductivity of the expansive soil. The closer the distance between the soil element and the cracks, the stronger the impact of cracks on the pore water pressure. The water then flows within the crack and reaches the deep part of the cracked soil quickly. As such, the bottom part of the cracks quickly becomes a saturated zone. As a consequence, the pressure head rises and the suction decreases in this zone. This is the cause of most expansive slope failures in the vicinity of cracks.

Competing Interests

The authors declare that there are no competing interests regarding the publication of this paper.


This project has been supported by the National Natural Science Foundation of China (no. 51409149), the Science Research Key Project of Hubei Provincial Department of Education (no. D20151201), and the Science Research Foundation of Yichang Science and Technology Bureau (no. A14-302-a10).


  1. Y. Xu and L. M. Zhang, “Breaching parameters for earth and rockfill dams,” Journal of Geotechnical and Geoenvironmental Engineering, vol. 135, no. 12, pp. 1957–1970, 2009. View at Publisher · View at Google Scholar · View at Scopus
  2. J. R. Talbot and C. E. Deal, “Rehabilitation of cracked embankment dams,” Geotechnical Special Publication, no. 35, pp. 267–283, 1993. View at Google Scholar
  3. H. L. Yao, S. H. Zheng, and S. Y. Chen, “Analysis on the slope stability of expansive soils considering cracks and infiltration of rain,” Chinese Journal of Geotechnical Engineering, vol. 23, no. 5, pp. 606–609, 2001 (Chinese). View at Google Scholar · View at Scopus
  4. J. H. Li, Field experimental study and numerical simulation of seepage in saturated/unsaturated cracked soil [Ph.D. thesis], The Hong Kong University of Science and Technology, Hong Kong, 2009.
  5. L. M. Zhang and D. G. Fredlund, “Characteristics of water characteristic curves for unsaturated fractured rocks,” in Proceedings of the Second Asian Conference on Unsaturated Soils (UNSAT-ASIA '4), pp. 425–428, Osaka, Japan, 2004.
  6. C. A. Burger and C. D. Shackelford, “Soil-water characteristic curves and dual porosity of sand-diatomaceous earth mixtures,” Journal of Geotechnical and Geoenvironmental Engineering, vol. 127, no. 9, pp. 790–800, 2001. View at Publisher · View at Google Scholar · View at Scopus
  7. C. J. Miller, H. Mi, and N. Yesiller, “Experimental analysis of desiccation crack propagation in clay liners,” Journal of the American Water Resources Association, vol. 34, no. 3, pp. 677–686, 1998. View at Publisher · View at Google Scholar · View at Scopus
  8. H.-J. Vogel, H. Hoffmann, and K. Roth, “Studies of crack dynamics in clay soil: I. Experimental methods, results, and morphological quantification,” Geoderma, vol. 125, no. 3-4, pp. 203–211, 2005. View at Publisher · View at Google Scholar · View at Scopus
  9. C. Tang, B. Shi, C. Liu, L. Zhao, and B. Wang, “Influencing factors of geometrical structure of surface shrinkage cracks in clayey soils,” Engineering Geology, vol. 101, no. 3-4, pp. 204–217, 2008. View at Publisher · View at Google Scholar · View at Scopus
  10. H. O. Peitgen, H. Jurgens, and D. Saupe, Chaos and Fractals: New Frontiers of Science, Springer, Berlin, Germany, 1st edition, 1992.
  11. Y. F. Xu and P. Dong, “Fractal approach to hydraulic properties in porous media,” Chaos, Solitons & Fractals, vol. 19, no. 2, pp. 327–337, 2004. View at Google Scholar
  12. M. Rieu and G. Sposito, “Fractal fragmentation, soil porosity, and soil water properties: I. Theory,” Soil Science Society of America Journal, vol. 55, no. 5, pp. 1231–1238, 1991. View at Publisher · View at Google Scholar · View at Scopus
  13. L. Guarracino, “A fractal constitutive model for unsaturated flow in fractured hard rocks,” Journal of Hydrology, vol. 324, no. 1–4, pp. 154–162, 2006. View at Publisher · View at Google Scholar · View at Scopus
  14. H. Q. Pham, D. G. Fredlund, and S. L. Barbour, “A study of hysteresis models for soil-water characteristic curves,” Canadian Geotechnical Journal, vol. 42, no. 6, pp. 1548–1568, 2005. View at Publisher · View at Google Scholar · View at Scopus
  15. S. Reitsma and B. H. Kueper, “Laboratory measurement of capillary pressure-saturation relationships in a rock fracture,” Water Resources Research, vol. 30, no. 4, pp. 865–878, 1994. View at Publisher · View at Google Scholar · View at Scopus
  16. Y.-T. Zhang, Rock Hydraulics and Engineering, China Water Power Press, Beijing, China, 2005 (Chinese).
  17. E. M. Kwicklis and R. W. Healy, “Numerical investigation of steady liquid water flow in a variably saturated fracture network,” Water Resources Research, vol. 29, no. 12, pp. 4091–4102, 1993. View at Publisher · View at Google Scholar · View at Scopus
  18. J. J. Nitao and T. A. Buscheck, “Infiltration of a liquid front in an unsaturated, fractured porous medium,” Water Resources Research, vol. 27, no. 8, pp. 2099–2112, 1991. View at Publisher · View at Google Scholar · View at Scopus
  19. K. Pruess, “On water seepage and fast preferential flow in heterogeneous, unsaturated rock fractures,” Journal of Contaminant Hydrology, vol. 30, no. 3-4, pp. 333–362, 1998. View at Publisher · View at Google Scholar · View at Scopus
  20. K. Pruess, “A mechanistic model for water seepage through thick unsaturated zones in fractured rocks of low matrix permeability,” Water Resources Research, vol. 35, no. 4, pp. 1039–1051, 1999. View at Publisher · View at Google Scholar · View at Scopus
  21. R. Therrien and E. A. Sudicky, “Three-dimensional analysis of variably-saturated flow and solute transport in discretely-fractured porous media,” Journal of Contaminant Hydrology, vol. 23, pp. 1–44, 1996. View at Google Scholar
  22. J. S. Y. Wang and T. N. Narasimhan, “Hydrologic mechanismsgoverning fluid flow in a partially saturated, fractured,porous medium,” Water Resources Research, vol. 21, no. 12, pp. 1861–1874, 1985. View at Publisher · View at Google Scholar · View at Scopus
  23. I. G. Currie, Fundamental Mechanics of Fluids, Marcel-Dekker, New York, NY, USA, 1993.
  24. National Research Council, Conceptual Model of Flow and Transport in the Fractured Vadose Zone, National Academy Press, Washington, DC, USA, 2001.
  25. Y. K. Tung, H. Zhang, C. W. W. Ng et al., “Transient seepage analysis of rainfall infiltration using a new conjunctive surface-subsurface flow model,” in Proceedings of the 57th Canadian Geotechnical Conference and the 5th Joint CGS-IAH Conference, vol. 7C, pp. 17–22, Quebec City, Canada, 2004.