# Unsaturated Seepage Analysis of Cracked Soil including Development Process of Cracks

Authors: Ling Cao; Zhijian Wang; Yong Chen
Journal: Advances in Materials Science and Engineering
 (2016)

---

## Abstract

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.


---

## Body

## 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 [3–6].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 (
1
) to determine the mechanisms of crack initiation and development under natural conditions, (
2
) to develop a method to determine the SWCCs of an unsaturated cracked soil, and (
3
) 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 [7–9]. 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 Table1. The soil has weak swell ability.Table 1
Basic physical properties of soil sample.
Specific gravity
G
s
Liquid limit
w
L (%)
Plastic limit
w
P (%)
Plasticity index
I
p
Free swelling rate
δ
e
f
Shrinkage limit
w
s (%)
2.72
69.65
29.95
28.66
52%
15.3%


### 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 Figure1. 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 (ρ
d
=
1.44 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

Figure2 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.
(a)
After the initial drying process
(b)
After the firstw
-
d cycle
(c)
After the secondw
-
d cycle
(d)
After the thirdw
-
d cycle
(e)
After the fourthw
-
d cycle
(f)
After the fifthw
-
d cycle
(g)
After the sixthw
-
d cycleFigures3 and 4 show the cracks development during the sixth drying-wetting cycle. The crack porosity δ
c changes with the water content w obviously.Figure 3
Crack development during a drying process.
(a)
w
=
31%, δ
c
=
0
(b)
w
=
30.5%, δ
c
=
0.08%
(c)
w
=
30.16%, δ
c
=
0.57%
(d)
w
=
29.58%, δ
c
=
1.35%
(e)
w
=
28.22%, δ
c
=
2.77%
(f)
w
=
6.11%, δ
c
=
5.002%Figure 4
Crack development during a wetting process.
(a)
w
=
6.11%, δ
c
=
5.002%
(b)
w
=
8.2%, δ
c
=
4.47%
(c)
w
=
10%, δ
c
=
4.09%
(d)
w
=
16%, δ
c
=
2.63%
(e)
w
=
20%, δ
c
=
1.56%
(f)
w
=
31%, δ
c
=
0To 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 Figure5. 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,δ
c, as an indicator of the total distribution characteristics, the equation is(1)
δ
c
=
∑
i
=
1
n
A
i
A
=
n
b
n
b
+
n
w
=
n
b
n
, where A
i is the area of the ith crack in the soil sample surface, which can be obtained by the “bwarea” function of MATLAB; n is the number of cracks; A is the total area of the soil sample surface; n
b is the number of black pixels, which is the cracking part of the sample; n
w is the number of white pixels, which is the noncracking part; and n is the total number of pixels.The widthr indicates the opening degree of the cracks. It is not uniform along the length of every crack. So, in this study, r 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 D. 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), N
(
ε
), changes with ε. The relationship between N
(
ε
) and ε is(2)
log
⁡
N
ε
=
A
-
D
log
⁡
ε
,where A is a constant. Hence, D values can be derived from the slop of log
⁡
N
(
ε
) versus log
⁡
ε plots.Figure6 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δ
c. (b) Crack fractal dimension D. (c) The maximum crack width r
m
a
x.
(a)
(b)
(c)


## 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 Table1. The soil has weak swell ability.Table 1
Basic physical properties of soil sample.
Specific gravity
G
s
Liquid limit
w
L (%)
Plastic limit
w
P (%)
Plasticity index
I
p
Free swelling rate
δ
e
f
Shrinkage limit
w
s (%)
2.72
69.65
29.95
28.66
52%
15.3%


## 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 Figure1. 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 (ρ
d
=
1.44 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

Figure2 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.
(a)
After the initial drying process
(b)
After the firstw
-
d cycle
(c)
After the secondw
-
d cycle
(d)
After the thirdw
-
d cycle
(e)
After the fourthw
-
d cycle
(f)
After the fifthw
-
d cycle
(g)
After the sixthw
-
d cycleFigures3 and 4 show the cracks development during the sixth drying-wetting cycle. The crack porosity δ
c changes with the water content w obviously.Figure 3
Crack development during a drying process.
(a)
w
=
31%, δ
c
=
0
(b)
w
=
30.5%, δ
c
=
0.08%
(c)
w
=
30.16%, δ
c
=
0.57%
(d)
w
=
29.58%, δ
c
=
1.35%
(e)
w
=
28.22%, δ
c
=
2.77%
(f)
w
=
6.11%, δ
c
=
5.002%Figure 4
Crack development during a wetting process.
(a)
w
=
6.11%, δ
c
=
5.002%
(b)
w
=
8.2%, δ
c
=
4.47%
(c)
w
=
10%, δ
c
=
4.09%
(d)
w
=
16%, δ
c
=
2.63%
(e)
w
=
20%, δ
c
=
1.56%
(f)
w
=
31%, δ
c
=
0To 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 Figure5. 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,δ
c, as an indicator of the total distribution characteristics, the equation is(1)
δ
c
=
∑
i
=
1
n
A
i
A
=
n
b
n
b
+
n
w
=
n
b
n
, where A
i is the area of the ith crack in the soil sample surface, which can be obtained by the “bwarea” function of MATLAB; n is the number of cracks; A is the total area of the soil sample surface; n
b is the number of black pixels, which is the cracking part of the sample; n
w is the number of white pixels, which is the noncracking part; and n is the total number of pixels.The widthr indicates the opening degree of the cracks. It is not uniform along the length of every crack. So, in this study, r 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 D. 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), N
(
ε
), changes with ε. The relationship between N
(
ε
) and ε is(2)
log
⁡
N
ε
=
A
-
D
log
⁡
ε
,where A is a constant. Hence, D values can be derived from the slop of log
⁡
N
(
ε
) versus log
⁡
ε plots.Figure6 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δ
c. (b) Crack fractal dimension D. (c) The maximum crack width r
m
a
x.
(a)
(b)
(c)


## 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 [11–13] 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:(3)
S
e
=
θ
θ
s
=
u
s
D
-
3
u
s
e
D
-
3
,
(4)
k
r
=
k
w
k
s
=
u
s
3
D
-
11
u
s
e
3
D
-
11
,where u
s
e is the air-entry value; S
e is the effective degree of saturation; θ is the volumetric water content of unsaturated soil at suction u
s; θ
s is the volumetric water content of saturated soil; D is the soil fractal dimension; k
r is the relative coefficient of permeability of unsaturated soil; k
s is the saturated coefficient of permeability; k
w is the coefficient of permeability at suction u
s. For a two-dimensional space, change D
-
3 to D
-
2 in (3) and 3
D
-
11 to 3
D
-
8 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 Figure7. 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.
(a)
The apparatus in the library
(b)
The pressure plate extractor
(c)
The temper pressure membrane apparatusBy 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δ
c, crack fractal dimension D, and the crack width r.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(5)
u
s
=
2
T
cos
⁡
α
r
,where T is the surface tension, α is the contact angle, u
s is the suction, and r is the crack width. Substituting (5) into (3), the equation becomes(6)
θ
=
θ
s
r
3
-
D
r
m
a
x
3
-
D
=
θ
s
2
T
cos
⁡
α
u
s
·
r
m
a
x
3
-
D
,where r
m
a
x is the maximum crack width, θ is the volumetric water content of crack at suction u
s, θ
s is the volumetric water content of saturated soil, which equals the crack porosity δ
c numerically, and D is the crack fractal dimension. For a two-dimensional space, (6) becomes(7)
θ
=
θ
s
r
2
-
D
r
m
a
x
2
-
D
=
θ
s
2
T
cos
⁡
α
u
s
·
r
m
a
x
2
-
D
.Given that T
=
0.0728 N/m at 20°C and α
=
60°, (7) can be rewritten as(8)
θ
=
θ
s
r
2
-
D
r
m
a
x
2
-
D
=
θ
s
0.0728
u
s
·
r
m
a
x
2
-
D
=
f
3
w
0.0728
u
s
·
f
2
w
2
-
f
1
w
,where f
1
(
w
), f
2
(
w
), and f
3
(
w
) are the function equations between the morphological parameters D, r
m
a
x, and δ
c 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 asD, θ
s, and r
m
a
x, as shown in Figure 9.Figure 9
The SWCC for the crack network at various water contents.
(a)
w
=
6.11%
(b)
w
=
7.08%
(c)
w
=
10.18%
(d)
w
=
14.84%
(e)
w
=
18.72%
(f)
w
=
22.89%
(g)
w
=
25.315%
(h)
w
=
28.22%
(i)
w
=
29.58%
(j)
w
=
30.16%


## 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 [11–13] 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:(3)
S
e
=
θ
θ
s
=
u
s
D
-
3
u
s
e
D
-
3
,
(4)
k
r
=
k
w
k
s
=
u
s
3
D
-
11
u
s
e
3
D
-
11
,where u
s
e is the air-entry value; S
e is the effective degree of saturation; θ is the volumetric water content of unsaturated soil at suction u
s; θ
s is the volumetric water content of saturated soil; D is the soil fractal dimension; k
r is the relative coefficient of permeability of unsaturated soil; k
s is the saturated coefficient of permeability; k
w is the coefficient of permeability at suction u
s. For a two-dimensional space, change D
-
3 to D
-
2 in (3) and 3
D
-
11 to 3
D
-
8 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 Figure7. 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.
(a)
The apparatus in the library
(b)
The pressure plate extractor
(c)
The temper pressure membrane apparatusBy 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δ
c, crack fractal dimension D, and the crack width r.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(5)
u
s
=
2
T
cos
⁡
α
r
,where T is the surface tension, α is the contact angle, u
s is the suction, and r is the crack width. Substituting (5) into (3), the equation becomes(6)
θ
=
θ
s
r
3
-
D
r
m
a
x
3
-
D
=
θ
s
2
T
cos
⁡
α
u
s
·
r
m
a
x
3
-
D
,where r
m
a
x is the maximum crack width, θ is the volumetric water content of crack at suction u
s, θ
s is the volumetric water content of saturated soil, which equals the crack porosity δ
c numerically, and D is the crack fractal dimension. For a two-dimensional space, (6) becomes(7)
θ
=
θ
s
r
2
-
D
r
m
a
x
2
-
D
=
θ
s
2
T
cos
⁡
α
u
s
·
r
m
a
x
2
-
D
.Given that T
=
0.0728 N/m at 20°C and α
=
60°, (7) can be rewritten as(8)
θ
=
θ
s
r
2
-
D
r
m
a
x
2
-
D
=
θ
s
0.0728
u
s
·
r
m
a
x
2
-
D
=
f
3
w
0.0728
u
s
·
f
2
w
2
-
f
1
w
,where f
1
(
w
), f
2
(
w
), and f
3
(
w
) are the function equations between the morphological parameters D, r
m
a
x, and δ
c 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 asD, θ
s, and r
m
a
x, as shown in Figure 9.Figure 9
The SWCC for the crack network at various water contents.
(a)
w
=
6.11%
(b)
w
=
7.08%
(c)
w
=
10.18%
(d)
w
=
14.84%
(e)
w
=
18.72%
(f)
w
=
22.89%
(g)
w
=
25.315%
(h)
w
=
28.22%
(i)
w
=
29.58%
(j)
w
=
30.16%


## 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 [17–22].


### 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 q
f
r within a small gap, b, between two parallel plates is(9)
q
f
r
=
-
K
f
r
d
H
d
l
,where H is the hydraulic head, l is the distance, K
f
r is the hydraulic conductivity of flow through a fracture, and it is defined as(10)
K
f
r
=
ρ
g
μ
b
2
12
,where ρ is the fluid density, g 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 as(11)
∂
∂
x
k
x
∂
H
∂
x
+
∂
∂
y
k
y
∂
H
∂
y
+
Q
=
m
w
γ
w
∂
H
∂
t
,where m
w
=
∂
θ
/
∂
u
w, and it is the slope of SWCC; u
w is the pore water pressure; γ
w is the unit weight of water. In the saturated zone, k
x and k
y are the respective saturated hydraulic conductivity in x- and y-directions, θ is the saturated volumetric water content, and m
w 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 are(12)
N
i
=
1
+
ξ
i
ξ
2
,
i
=
1,2
,
N
i
=
1
+
ξ
i
ξ
1
+
η
i
η
4
,
i
=
1,2
,
3,4
,where the local coordinates ξ
i and η
i at node i are −1 and 1, respectively. If node i is located on the negative side, the coordinate ξ
i or η
i is −1; otherwise ξ
i or η
i 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.,x- and y-directions). In a crack element, the hydraulic head distribution is defined by the function of(13)
h
=
N
→
·
H
→
=
N
1
H
1
+
N
2
H
2
,
l
=
N
→
·
l
→
=
N
1
l
1
+
N
2
l
2
,where H
1 and H
2 are the hydraulic heads of the nodes at the ends of a linear crack of length L; l
1 and l
2 are the nodal locations. In the local coordinate system, l
1
=
-
L
/
2 and l
2
=
L
/
2 when ξ
1
=
-
1 and ξ
2
=
1, respectively.According to Darcy’s Law, the derivative of (13) is(14)
∂
N
→
∂
l
=
∂
N
→
∂
ξ
1
∂
l
/
∂
ξ
=
∂
N
→
∂
ξ
1
∂
N
1
/
∂
ξ
l
1
+
∂
N
2
/
∂
ξ
l
2
=
2
L
∂
N
→
∂
ξ
.Then the gradient in an element can be written as(15)
i
=
∂
h
∂
l
=
∂
N
→
∂
l
·
H
→
=
2
L
∂
N
1
∂
ξ
H
1
+
∂
N
2
∂
ξ
H
2
=
H
2
-
H
1
L
.The flow velocityv in crack is(16)
v
=
k
i
=
k
H
2
-
H
1
L
.Considering the velocity vector resolution, when the angle between the crack and horizontal direction isα, the velocity components at x and y directions are(17)
v
x
=
v
cos
⁡
α
;
v
y
=
v
sin
⁡
α
.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].


## 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 q
f
r within a small gap, b, between two parallel plates is(9)
q
f
r
=
-
K
f
r
d
H
d
l
,where H is the hydraulic head, l is the distance, K
f
r is the hydraulic conductivity of flow through a fracture, and it is defined as(10)
K
f
r
=
ρ
g
μ
b
2
12
,where ρ is the fluid density, g 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 as(11)
∂
∂
x
k
x
∂
H
∂
x
+
∂
∂
y
k
y
∂
H
∂
y
+
Q
=
m
w
γ
w
∂
H
∂
t
,where m
w
=
∂
θ
/
∂
u
w, and it is the slope of SWCC; u
w is the pore water pressure; γ
w is the unit weight of water. In the saturated zone, k
x and k
y are the respective saturated hydraulic conductivity in x- and y-directions, θ is the saturated volumetric water content, and m
w 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 are(12)
N
i
=
1
+
ξ
i
ξ
2
,
i
=
1,2
,
N
i
=
1
+
ξ
i
ξ
1
+
η
i
η
4
,
i
=
1,2
,
3,4
,where the local coordinates ξ
i and η
i at node i are −1 and 1, respectively. If node i is located on the negative side, the coordinate ξ
i or η
i is −1; otherwise ξ
i or η
i 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.,x- and y-directions). In a crack element, the hydraulic head distribution is defined by the function of(13)
h
=
N
→
·
H
→
=
N
1
H
1
+
N
2
H
2
,
l
=
N
→
·
l
→
=
N
1
l
1
+
N
2
l
2
,where H
1 and H
2 are the hydraulic heads of the nodes at the ends of a linear crack of length L; l
1 and l
2 are the nodal locations. In the local coordinate system, l
1
=
-
L
/
2 and l
2
=
L
/
2 when ξ
1
=
-
1 and ξ
2
=
1, respectively.According to Darcy’s Law, the derivative of (13) is(14)
∂
N
→
∂
l
=
∂
N
→
∂
ξ
1
∂
l
/
∂
ξ
=
∂
N
→
∂
ξ
1
∂
N
1
/
∂
ξ
l
1
+
∂
N
2
/
∂
ξ
l
2
=
2
L
∂
N
→
∂
ξ
.Then the gradient in an element can be written as(15)
i
=
∂
h
∂
l
=
∂
N
→
∂
l
·
H
→
=
2
L
∂
N
1
∂
ξ
H
1
+
∂
N
2
∂
ξ
H
2
=
H
2
-
H
1
L
.The flow velocityv in crack is(16)
v
=
k
i
=
k
H
2
-
H
1
L
.Considering the velocity vector resolution, when the angle between the crack and horizontal direction isα, the velocity components at x and y directions are(17)
v
x
=
v
cos
⁡
α
;
v
y
=
v
sin
⁡
α
.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 Figure7. 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 dimensionD 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 headh are the y-coordinate. For the second stage, h
=
0 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. Figures11(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)
Without cracks
(b)
With cracksA comparison of these two figures shows that the cracks have significant impact on the infiltration through the soil. In Figure11(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).Figure13 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′.Figure14 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′.Figure15 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.


## 5.1. Computational Model

The computational model consists of a 0.3 m × 0.3 m area, as shown in Figure7. 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 dimensionD 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 headh are the y-coordinate. For the second stage, h
=
0 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. Figures11(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)
Without cracks
(b)
With cracksA comparison of these two figures shows that the cracks have significant impact on the infiltration through the soil. In Figure11(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).Figure13 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′.Figure14 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′.Figure15 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:(
1
) a crack network system and (
2
) 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.



---

*Source: 2684297-2016-03-15.xml*