Open access peer-reviewed chapter

Study of Unsaturated Soils by Coupled Numerical Analyses of Water Flow-Slope Stability

Written By

Norma Patricia López-Acosta and José Alfredo Mendoza-Promotor

Submitted: 21 January 2016 Reviewed: 22 April 2016 Published: 27 July 2016

DOI: 10.5772/63903

From the Edited Volume

Groundwater - Contaminant and Resource Management

Edited by Muhammad Salik Javaid

Chapter metrics overview

2,398 Chapter Downloads

View Full Metrics


The geotechnical engineering, among the problems related to water flow, is specifically interested in soil and water that it contains, and also on the movement of water through their pores, in addition to the laws governing this phenomenon. A very important subject is to quantify the retention and filtration of water within the soil structure; however, the emphasis should be not only on how much water flows through the soil but also on the state of pore water pressures because this pressure, either positive or negative, has a direct influence on the stress state and changes in volume of soil. Several publications address the issue of water flow in saturated state; however, only some of them consider the flow under unsaturated conditions. In this chapter, the main emphasis is focused on the study of water flow in unsaturated soils.


  • water flow
  • unsaturated soils
  • numerical analysis
  • coupled water flow-slope stability analysis
  • soil-water characteristic curve
  • hydraulic conductivity function

1. Introduction

Unsaturated soils are characterized by negative pore-water pressure. In practical terms, partially saturated and unsaturated are synonymous: both terms indicate a degree of saturation lower than one; however, in specific terms, unsaturated implies the introduction of a third phase (gaseous) to the two-phase system already present in the saturated soils (liquid and solid phases only) [1].

Currently, computer programs are helpful for solving problems of water flow and facilitate both the study of transient-state flow and the characterization of unsaturated soils, which, from an analytical standpoint, are complicated and laborious tasks [2].

To demonstrate the application of the theoretical foundations presented in this chapter, an analysis of a tailings dam is presented as a case study. These structures are generally found in an unsaturated state. Thus, the primary goals of this analysis are to describe the applied methodology and to establish criteria and recommendations that can be followed to solve related problem types.


2. General theoretical foundations of unsaturated soils

2.1. Soil suction

Several definitions exist that define soil suction and its significance [35]; however, for practical engineering applications, proposed definitions are either unsuitable or complex because they are largely based on thermodynamic concepts. In simpler terms, suction can be defined as a state of negative water pressure in the soil, which is influenced by several factors [6], including temperature, gravity, capillary effects, salt content, and electrical charge (van der Waals forces), among others.

Total suction, or simply suction, is composed of two variables: (a) matric suction and (b) osmotic suction. Matric suction describes the difference between air pressure and water pressure [6], whereas osmotic suction is defined as the negative pressure resulting from the effect of the dissolved salts in the water of the soil matrix. Osmotic suction is commonly disregarded due to its lesser influence on total suction in addition to the difficulty of separating these two variables.

2.2. Water storage function or SWCC

The storage function, which is also defined by the soil-water characteristic curve (SWCC), describes the relationship between the water content of the soil (degree of saturation, gravimetric, or volumetric water content) and the soil suction. The nature of the SWCC is directly associated with grain-size distribution and soil structure. Therefore, the water content-suction ratio varies as a function of the material type (Figure 1).

Figure 1.

Soil-water characteristic curve (SWCC) for different soil types [7].

The SWCC represents a fundamental relationship that can be used to describe the behaviour of unsaturated soils [8]. Similarly, the SWCC can be used to determine other properties, such as the permeability and shear strength of soil and the volumetric changes of soil [9].

Several methods and empirical relationships have been proposed to describe the SWCC. Several of these models use relationships that depend on the air entry value, the residual water content, or the slope of the curve [10, 11]. Other procedures are based on the grain-size distribution curve in addition to other soil properties (volume-mass properties) or the use of statistical correlations between soil data and soil water content [9, 12, 13]. The reliability of the models depends on the quality and quantity of data used for the statistical correlations [14]. Other methods consider the pore-size distribution within the soil, which, in certain cases, can be determined from the grain-size distribution curve of the material [1517].

2.2.1. Interpretation of the SWCC

The SWCC can be associated with three zones that describe the desaturation process of a soil (Figure 2). In addition, the SWCC allows for the saturated and residual water content to be determined, as well as their respective suction values, which are the main parameters that represent the SWCC as an input into one of the fitting methods mentioned in this chapter:

Figure 2.

Zones corresponding to the soil-water characteristic curve (SWCC).

  • Saturated capillary zone. The soil zone that is maintained in a saturated state, whose defining limit coincides with the air-entry value [18], which can be described as the value that the matric suction must exceed before air enters into the soil macropores.

  • Desaturation zone. The zone where water is displaced due to the air entry in the pores. The defining limit of this zone is determined by the residual water content or where water in the pores becomes discontinuous and permeability considerably decreases.

  • Residual zone. This section of the curve represents the zone where increases in suction do not produce significant changes in water content. As water is so scarce, it does not flow between pores, and its removal only occurs by evaporation [19]. This region is characterized by extremely high suction values.

  • Air-entry value. Represents the suction value at which air begins to displace water in pores, which begins with those of greatest size.

  • Residual value. Represents the suction value at which the liquid phase of the soil becomes discontinuous and surrounds soil particles as a thin film.

  • Residual water content. The content of water at which high suction values are required to remove any additional water from the soil mass.

  • Saturated water content. The water content in soil in a saturated state.

2.2.2. Estimation methods of SWCC

When laboratory data on the relationship between suction and water content are unavailable or insufficient, estimation models are one method used to determine the SWCC as a function of the index properties of soil (volume-mass and grain-size distribution relationships). Currently, there are several models that allow for the SWCC to be determined, including the following:

  • Fredlund and Wilson method [9]

  • Scheinost, Sinowski, and Auerswald method [20]

  • Aubertin, Mbonimpa, Bussière, and Chapuis method (modified from Kovacs) [21]

  • Zapata’s correlations [22]

  • Arya and Paris method [15]

2.2.3. Fitting methods of SWCC

Fitting methods are based on experimental or empirical equations that aim to define the SWCC according to data obtained from laboratory tests. Fitting methods are used when laboratory data are available, yet given the dispersion of the values, it is necessary to apply models to adjust or define the data trend to generate a representative curve for the material considered in the study.

Several models have been developed to define the SWCC [23], which consider fitting parameters that provide a range of flexibility to represent distinct materials. In Table 1, several of the most common fitting methods for SWCC are listed. These equations consider gravimetric water content as the primary variable; however, the equations can also use volumetric water content or the degree of soil saturation as inputs [24].

Model Equation
Van Genuchten (1980)  ww=wr+(wswr){ 1[ 1+(avbψ)nvb ]mvb }
Gardner (1958)  ww=wr+(wswr){ 1agψng }
Brooks and
Corey (1964) 
ww=wr+(wswr)[ acψ ]nc
Van Genuchten and Mualem (1976)  ww=wr+(wswr){ 1[ 1+(avmψ)nvm ]mvm }
Van Genuchten (1980)  ww=wr+(wswr){ 1[ 1+(avgψ)nvg ]mvg }
Fredlund and
Xing (1994) 
ww=ws[ 1ln(1+ψψr)ln(1+106ψr) ]{ 1[ ln(e+(ψaf)nf) ]mf }
ws saturated gravimetric water content
wr residual gravimetric water content
ψ soil suction
avb adjustment parameter that depends on the air-entry value of the soil
nvb adjustment parameter that depends on the desaturation velocity of the soil once the air-entry value has been exceeded
mvb adjustment parameter related to the residual water content of soil, which is considered to be mvb=1-2/nvb
ag adjustment parameter derived from the air-entry value of the soil
ng adjustment parameter that depends on the desaturation velocity of the soil once the air-entry value has been exceeded
ac air pressure
nc soil pore size index
avm adjustment parameter that depends on the air-entry value of the soil
nvm adjustment parameter that depends on the desaturation velocity of the soil once the air-entry value has been exceeded
mvm adjustment parameter related to the residual water content of soil, which is considered to be mvm=1-1/nvm
avg adjustment parameter that depends on the air-entry value of the soil
nvg adjustment parameter that depends on the desaturation velocity of the soil once the air-entry value has been exceeded
mvb adjustment parameter related to the residual water content of soil, is considered to be mvg=1-1/n o mvm=1-1/2n
af adjustment parameter that depends on the air-entry value of the soil
nf adjustment parameter that depends on the desaturation velocity of the soil once the air-entry value has been exceeded
mf adjustment parameter related to the residual water content of soil
ψr residual suction
e irrational number

Table 1.

Experimental fitting models of the SWCC.

2.3. Hydraulic conductivity function

The hydraulic conductivity function represents the relationship between hydraulic conductivity and soil suction and can be expressed as a function of the degree of saturation or volumetric water content of soil.

Several estimation methods are based on the SWCC and the saturated hydraulic conductivity at distinct suction intervals. These techniques can be classified into the following categories:

  • Empirical or experimental equations. These equations outline the relationship between the SWCC and a hydraulic conductivity function, including the models of Brooks and Corey [25] and Gardner [26].

  • Statistical equations. These equations consist of a physical model that represents the trajectories of water flow through soil pores of different sizes. Examples of these include the models of Van Genuchten and Burdine [27] and Van Genuchten and Mualem [28].

  • Correlation equations. These models are based on a proposed relationship between the SWCC and the hydraulic conductivity function, for example, the model of Leong and Rahardjo [29].

  • Regression equations. These models use values of hydraulic conductivity obtained from laboratory testing or other estimation method, for example, the model of Fredlund and Xing [30].

A summary of the main estimation methods expressed as a function of hydraulic conductivity is presented in Table 2.

Model Equation
Van Genuchten and Burdine (1953)  k(ψ)=ks{ 1(avbψ)nvb2[ 1+(avbψ)nvb ]mvb }2[ 1+(avbψ)nvb ]2mvb
Gardner (1958)  k(ψ)=ks1+ag(ψρwg)ng
Brooks and Corey (1964)  k(ψ)={ ksψψaevks(ψaevψ)2+3λψψaev
Van Genuchten and Mualem (1976)  k(ψ)=ks{ 1(avmψ)nvm1[ 1+(avmψ)nvm ]mvm }2[ 1+(avmψ)nvm ]mvm/2
Fredlund and Xing (1994)  k(ψ)=kslnψbθ(ey)θ(ψ)eyθ'(ey)dyln(ψaev)bθ(ey)θseyθ'(ey)dy
Leong and Rahardjo (1997)  k(ψ)=ks[ Θd(ψ) ]q
k(ψ) hydraulic conductivity function
ks saturated hydraulic conductivity coefficient
ψ soil suction
ψaev soil suction at the air-entry value
λ pore size distribution index
ρw water density
g acceleration of gravity
ag adjustment parameter related to the air-entry value
avb, avm adjustment parameter related to the inverse of the air-entry value
ng soil parameter that depends on the desaturation process once the air-entry value is exceeded
nvb, nvm adjustment parameter of the characteristic curve obtained from Van Genuchten’s model (1980) 
mvb adjustment parameter (1-2/nvb); value oscillates between zero and one
mvm adjustment parameter (1-1/nvm); value oscillates between zero and one
b upper integration limit (1,000,000 kPa) 
y fictitious variable that describes soil suction (logarithm scale) 
θ' variable derived from the function of soil storage capacity (characteristic curve) 

Table 2.

Estimation models for determining hydraulic conductivity function.


3. Unsaturated seepage theory

In geotechnical engineering, water flow through soil can be categorized in several ways: laminar or turbulent flow (determined by the Reynolds number); one-dimensional (1D), two-dimensional (2D), or three-dimensional (3D) flows (depending on the number of planes); steady-state and transient seepage (constant and variable across time, respectively); and confined and unconfined flows (depending on the limits defining them).

3.1. Steady-state seepage

Figure 3.

One- and two-dimensional water flow through an unsaturated element [24].

Considering that Figure 3(a) represents a soil sample subject to an upward flow and taking into consideration the law of continuity, the following can be proposed:


After applying Darcy’s law [31], the one-dimensional water flow through an unsaturated soil can be expressed as



dhw/dy is the hydraulic gradient

dhw is the hydraulic head

For the example of Figure 3(b), which demonstrates a two-directional water flow, the following expression can be deduced from the continuity equation:


After applying Darcy’s law [31], the following equation can be proposed to describe steady-state seepage in two directions for an unsaturated anisotropic soil:


Under isotropic conditions, the previous equation can be expressed as follows:


Figure 4.

Three-dimensional water flow through an unsaturated element [24].

For three-directional water flow, considering the unsaturated soil sample subjected to the water flow conditions indicated in Figure 4, where the hydraulic conductivities vary in all directions, the following expression can be deduced based on flow continuity:


Once again, after applying Darcy’s law, the equation that describes the steady-state seepage in three directions in anisotropic conditions is detailed as follows:


For isotropic soils, the previous equation can be simplified to


3.2. Transient seepage (the Richards equation)

In the transient seepage analysis and in contrast to steady state seepage, a variable hydraulic head exists over time. Variation occurs due to the changes in the boundaries of the system (due to variation in water levels over time).

For practical applications, Darcy’s law [31] can be generalized to unsaturated water flow problems by considering hydraulic conductivity to be a function of the soil suction or suction head [32, 33]:



hm is the suction head

k(hm) is the hydraulic conductivity function.

If the osmotic pressure head is disregarded, then the total head of an unsaturated soil can be expressed as the sum of the matric suction head and the elevation head (h=hm+z). Thus, if this consideration is substituted in the equation of the law of conservation of matter and a constant water density is assumed, then the following expression can be obtained:

x[ kx(hm)hmx ]+y[ ky(hm)hmy ]+z[ kz(hm)(hmz+1) ]=θtE10

where the additional term in the direction of the z-axis is due to the elevation head.

The term on the right side of Eq. (10) can also be expressed as a function of the matric suction head:


where θ/hm is the slope of the relationship between the volumetric water content and the suction head, which can be directly determined from the SWCC. The slope denotes the specific moisture capacity, which is typically denoted as C. As the soil storage function is not linear, it is necessary to describe the specific moisture capacity as a function of the suction or suction head:


If Eqs. (11) and (12) are substituted in Eq. 10, the expression that describes transient seepage in unsaturated soils can be expressed as follows:

x[ kx(hm)hmx ]+y[ ky(hm)hmy ]+z[ kz(hm)(hmz+1) ]=C(hm)θtE13

Eq. (13) is known as the Richards [33] equation, where given the boundary limits and initial conditions specific to a system, the equation provides the values for suction across space and time. It is highlighted that in using this equation, it is necessary to have data on the SWCC and the hydraulic conductivity function that are specific to the material being studied.


4. Slope stability analyses

One of the most common methods for evaluating slope stability is the general limit equilibrium (GLE) method. A series of equations has been proposed by several authors, who agree in dividing the slide zone into slices. The primary differences are related to the equations that these seek to satisfy and to the differential forces influencing each slice, including the existing relationship between shear and normal forces. The foundation of this method is based on two equations that determine the factor of safety and evaluate the relationship between normal and shear forces [34, 35]. Thus, one equation provides the factor of safety with respect to the equilibrium of moments (Fm) and the other with respect to horizontal forces (Ff).

For the GLE, shear forces are determined according to the equation proposed by Morgenstern and Price [36]:



X is the shear force between slices

E is the normal force

λ percentage (in decimal form) of the expressed function

f(x) is the function used for study.

For the GLE, the factor of safety is determined with respect to the method of moments and is determined by

Fm= [ c'βR+(Nuβ)Rtanφ' ] Wx NfDdE15

However, the factor of safety with respect to the balance of forces is determined by

Ff= [ c'βcosα+(Nuβ)tanφ'cosα ] Nsinα DcosωE16


c' is the effective cohesion

φ' is the effective friction angle

u is the pore pressure

N is the normal force at base of slice

W is the weight of slice

D is the load line

β, R, x, f, d, ω are the geometric parameters

α is the slope.

One important factor that is relevant for the previous equations is the normal force N, which is a term defined as


5. Practical application: coupled water flow-slope stability analyses of a tailings dam

5.1. Introduction

Mining waste (tailings) is disposed of in structures called tailings dams. The most important difference between a tailings dam and a typical water storage dam in the conventional water-retaining sense is that tailings dams do not contain an engineered water barrier and they are built and used simultaneously. For this reason, it is common for the originally conceived project to undergo constant modifications. Current technological advances allow for the numerical modelling of steady and transient state flow analyses for these structures, and along with adequate monitoring, their stability can be verified with coupled water flow-slope stability analyses before the occurrence of any potential errors. In the following sections, a cross section of a tailings dam is analysed (Figure 5), and the importance and benefits that result from the numerical model are highlighted; in addition, further recommendations are given for this type of analysis.

Figure 5.

Maximum tailings dam cross section for two-dimensional model.

5.2. Material properties

The properties of tailings dams are variable because they depend largely on the origin of the materials used for their construction. Several authors have performed tests with different materials and have specified typical values for these [3740]. Based on these references, in Table 3, the properties of the materials from the case study considered in the numerical model of this section are provided.

Material USCS classification k (m/s) c' (kPa) φ' (°)
Name Symbol
Fine material Silt with sand ML 1.582 × 10−8 0 25
Coarse material Silty sand SM 1.000 × 10−6 0 35

Table 3.

Tailings properties for the tailings dam considered.

Assuming that the materials of the tailings dam are found in unsaturated state, for this analysis, the SWCC and hydraulic conductivity functions must be determined. If laboratory data are unavailable, it is possible to use estimates with the aid of material index properties, such as grain-size distribution curves, which are considered in the calculations performed here (Figure 6).

Figure 6.

Grain-size distribution curves for tailings considered in the numerical model.

5.3. Unsaturated soil properties functions

For the dam materials analysed in the present case study, the SWCC (that represent the relationship between water content and soil suction) of Figure 7 were determined by the following method:

  • For coarse material, based on data generated from a laboratory experiment with a pressure plate, the fitting Fredlund and Xing [30] SWCC equation was applied.

  • For fine material, the estimation method of Fredlund and Wilson [9] was initially used, and afterwards, the fitting Fredlund and Xing [30] SWCC equation was applied.

Figure 7.

Soil-water characteristic curves of tailings dam materials.

In Figure 8, the hydraulic conductivity functions considered in the present dam case study are presented. A great similarity between the laboratory data and the estimated data can be observed. Thus, the use of estimation methods when laboratory data are unavailable represents a feasible solution and potential advantage; however, such models should be used with caution and rationality.

Figure 8.

Hydraulic conductivity functions of the tailings dam materials.

5.4. Two-dimensional modelling for steady-state seepage conditions

A steady-state seepage analysis was performed with Eq. (5), which can be solved numerically by the finite element method using the Seep/W algorithm [41]. For two-dimensional models, boundary conditions should be adequately defined in addition to discretizing the flow regions the best possible. Greater emphasis must be placed on areas that require greater detail in the expression of results, where numerical difficulties may be presented, or for areas with high contrast in the permeability of materials by orders of magnitude. In Figure 9, the boundary conditions and the discretization of the flow regions are detailed.

Figure 9.

Boundary conditions and discretization model for the tailings dam with Seep/W [41].

Based on the steady-state seepage analysis, it is possible to estimate the position of the phreatic surface of the tailings dam and for the different pond levels (Figure 10).

Figure 10.

Variation in the phreatic surface due to reduction of the freeboard with Seep/W [41].

Based on the previous results, the stability of the structures of the tailings dam can be evaluated by a GLE method, such as the Morgenstern-Price method, to define the slip surface for each level of the pond. In Figure 11, the slip surfaces are indicated, which were obtained from the Slope/W algorithm [42] and from considering the results of the steady-state groundwater flow analyses.

Figure 11.

Slip surfaces due to phreatic surfaces for different pond levels with Slope/W [42].

5.5. Two-dimensional modelling for transient conditions

A transient model (variable over time) was performed to evaluate the influence of rainfall on the tailings dam defined in Figure 5 by the numerical solution of Eq. 13. In addition, the behaviour of the pond levels is studied when the level varies 3.00 m with respect to the crest. In this case, the boundary conditions are modified to be a time-dependent function while the discretization of the model was maintained without modifications with respect to the mesh used for the steady-state analysis.

The aforementioned model allows for the distinct stages of the analysis to be defined. The first set of calculations corresponds to the period of 0–48 h, after which, a storm would cause an immediate increase in water level of the pond (Figure 12). An evaluation of the infiltration due to rainfall during this period is also shown in Figure 13.

Figure 12.

Variation in the water surface within the tailings structure due to rainfall and loss of freeboard for 48-h period with Seep/W [41].

Figure 13.

Behaviour of the pore-water pressure near the ground surface due to rainfall and loss of freeboard.

Short analysis periods for materials of low permeability are often inadequate because they might not be able to clearly determine the influence of water. For example, in the previous analysis during the water filling of pond, the water surfaces variation was unable to be clearly distinguished. Similarly, rainfall was shown to have low significance because it was unable to infiltrate to a considerable enough depth to affect the slope stability.

Therefore, long-term analyses for these types of materials are more representative because they allow remark the water surfaces variation over a longer period of time. Thus, the results of these analyses are shown in Figure 14.

Figure 14.

Water surfaces variation within the tailings structure for a transient analysis of 20 years with Seep/W algorithm [41].

Following the previous procedure and considering the pond levels indicated in Figure 12, the slip surfaces were defined by the Morgenstern-Price method with the Slope/W algorithm [42]. It is worth mentioning that a stability analysis for a period of 48 h did not provide significant results because during this time, the structure was not considerably influenced by the water. However, for long-term conditions, such stability analyses could play an important role in the study.

Figure 15.

Variation in the factor of safety in the tailings dam over a period of 20 years.

In Figure 15, the variation of factor of safety over time is presented, wherein the variation in the position of water surface over time is considered (Figure 14). A decrease in stability over time is also observed. For this type of evaluations, numerical difficulties may be presented that are subsequently reflected in the result, which leads to erroneous behaviour. Therefore, it is recommended to manipulate the variables that may affect the model, such as the calculated time intervals (smaller time increments may lead to a more detailed response), the finite elements mesh, and the convergence parameters. An adequate manipulation of these variables can significantly improve the results of the numerical model [43].

5.6. Three-dimensional numerical modelling

5.6.1. Model geometry

In current numerical modelling, it is common for 3D domains to be extruded, or in other words, two-dimensional sections are assigned a certain thickness. Then, the model is extended along a horizontal axis, resulting in a continuous, three-dimensional model.

Extrusions are recommended depending on the type of problem to be solved. For example, extrusions are suitable for analyses of protection levees or embankments with regular geometries that extend over long distances. In the case of structures with irregular geometries, such as dams, it is more suitable that calculations consider the specific topography of the site for these to be more representative. This latter method requires a greater amount and detail of information for the numerical model to be successfully resolved. However, it is convenient in this case to standardize certain surface areas of the model to avoid an excessive discretization of each region. In Figure 16, the geometries assumed for both cases performed here are shown.

Figure 16.

Geometries: (a) 3D model extrusion and (b) 3D realistic model.

5.6.2. Boundary conditions

In assigning boundary conditions, sufficient regions must be created to allow the site-specific conditions to be adequately represented and well defined, which thus leads to optimum modelling results. In Figure 17, the conditions assigned in both scenarios of the case study analysis are shown.

Figure 17.

Boundary conditions of the 3D models considered in the calculations with SVFlux algorithm [44].

5.6.3. Discretization of the model

The mesh generation for 3D models is perhaps the most complicated part of this process and even more so when considering site-specific topographic conditions. During these processes, the benefits of the extruded 2D models are evident because the generation of the mesh, as well as its distribution, is more regular. Moreover, 3D models of realistic topographies make mesh generation more difficult; in addition, 3D models require a greater number of elements to adapt to the model. In Figure 18, the distribution of the meshes obtained with the SVFlux algorithm is shown [44], and in Table 4, a comparison is made of the number of elements required for each model.

Figure 18.

Finite element mesh used for calculation of the 3D models with SVFlux algorithm [44].

Model Nodes Elements
2D 220 91
3D extruded 1342 758
3D realistic topography 2889 1580

Table 4.

Comparison of the number of finite elements for different analysis conditions.

5.6.4. Water flow analyses

Important differences are found when 3D model extrusions are compared to 3D models that consider site-specific topography. In Figure 19, it can be observed that the distribution of the hydraulic heads tends to vary in both cases. The 3D model extrusion shows a constant dissipation in the hydraulic head; however, this is not very representative of this type of structure. On the other hand, the 3D realistic model that considers site-specific topography shows a more variable behaviour in the distribution of the hydraulic head, which can be considered to be more representative of the analysed case study.

Figure 19.

Distribution of hydraulic heads (m) with SVFlux algorithm [44].

The previous results can be verified by comparing the distributions of the hydraulic heads at the maximum cross section of the structure. Theoretically, the distributions of the 2D model, 3D model extrusion, and 3D realistic model with site-specific topography should be nearly identical or extremely similar considering the similar conditions and inputs for the analysis. In Figure 20, this comparison is shown.

Figure 20.

Variation in the hydraulic heads (m) for different analysis criteria.

5.6.5. Three-dimensional slope stability analysis

Several studies have demonstrated the importance of defining the relationship between the two-dimensional and three-dimensional models. The majority of studies agree in that the values for the factors of safety obtained from the 2D slope stability analyses are conservative, and these values tend to increase upon considering 3D models, which occurs because these calculations consider steady-state water conditions and disregard the filtration forces generated within the structure [45, 46]. This phenomenon is observed in Figure 21, where the results of a slope stability analysis during dry conditions are shown. In this case, the 3D model extrusion and 3D realistic model (considering site-specific topography) did not present significant differences.

On the contrary, in Figure 22, the safety factors obtained for each scenario from a slope stability analysis with the linear phi-b model are shown, which consider water flow in the tailings structure. For these cases, the differences between each of the analyses and their criteria are distinguishable. According to the previous results, it is important to consider the water flow through the specific medium under study because it can significantly affect the stability of the earth structures. In this case, the 2D and 3D model extrusions show a high value for the factor of safety in comparison with the 3D realistic model. These differences can be attributed to the irregular topography, which causes the distribution of the hydraulic head to exhibit non-linear behaviour.

Figure 21.

Factor of safety for distinct slope stability analyses in dry conditions (without water flow) with SVSlope algorithm [47].

Figure 22.

Factor of safety for distinct slope stability analyses considering water flow through the tailings structure with SVSlope algorithm [47].

It is important to highlight that this type of analysis requires more exhaustive evaluations. In this case, only the importance of these calculations is reinforced, and methodological suggestions are proposed. However, additional numerical analyses and their respective validation in the field are imperative to adequately define the behaviour of these types of numerical models.


6. Conclusions

The primary focus of this study was unsaturated soils. The most important specific concepts related to this subject were detailed with the goal of enabling their application with the aid of specialized computer programs based on the finite element method.

Several of the concluding comments and recommendations derived from the analysis of the case study in this chapter are detailed as follows:

  • When specialized laboratory results are unavailable (to describe variations in water content and permeability with respect to suction), the unsaturated soil property functions required for the analysis can be determined by estimation methods.

  • For the correct estimation of unsaturated soil property functions (soil-water characteristic curve and hydraulic conductivity function), a large amount of laboratory data can be required, specifically, results from different index properties, grain-size distribution curve, and soil permeability experiments.

  • To determine the appropriate unsaturated soil property functions (soil-water characteristic curve and hydraulic conductivity function) for the calculations in analyses performed here, the form and tendency of the soil curves were considered. However, it is possible to use statistical methods to define these functions. Certain programs will perform these calculations (such as SoilVision Database).

  • For numerical modelling that considers materials with high contrast in permeability by orders of magnitude, for example, it is common for numerical difficulties to be presented. To resolve this type of problems, highly permeable materials that do not contribute to the dissipation of the hydraulic head can be omitted; however, as a consequence, these will not have an effect in the final resolution of the model. If the material must be considered in calculations, then it can be substituted by gravel to facilitate the convergence of the model. To obtain satisfactory results, an adequate definition of the SWCC is necessary in addition to the hydraulic conductivity functions, which will be necessary to obtain representative solutions. Despite this, numerical complications may be present, which can only be reduced with convergence parameters, or in this case, by step-by-step analysis until an optimal solution is achieved.

  • Depending on the problem to be analysed and its complexity, more detailed discretizations may be required to decrease the numerical error of calculations. This approach should not be confused with the automatic generation of a dense finite element mesh because for many water flow models, the use of such discretizations is not justified and will only affect calculation times without leading to improvements in the solution of the problem.

  • Water infiltration plays an important role in the stability of deposits of mining wastes (tailing dams), and it is appropriate to consider water movement in these calculations to better understand its behaviour. Therefore, the benefits of numerical modelling are numerous because these methods allow for diverse situations, given relatively short periods of analysis, to be considered, thereby informing and enabling appropriate decision-making.

  • Currently, 3D models have wide advantages over 2D models because the latter do not consider all the variables that influence the behaviour of structures. Some of the main features that can be included in three-dimensional models are irregular geometry of the structure under study, topographic configuration of soil, unsaturated flow that represent more realistically the physical conditions of the problem, among others.

  • In 3D analysis, extrusions are recommended depending on the type of problem to be solved. For example, protection levees or embankments with regular geometries extending in great lengths are suitable extrusions. In the case of structures with irregular geometries, such as dams, it agrees that calculations are performed considering the topography of the site, to give greater representation to them.

Finally, for more representative numerical analyses, it is recommended to apply the theory of unsaturated soils in all those situations where the material is in this state. This chapter demonstrated that the considerations for modelling, estimating, and fitting methods of the soil properties functions provide the necessary elements to carry out these analyses in a simple way. Illustrated examples also evaluated the potential of numerical methods to increase the degree of realism in unsaturated soils analysis. Computer programs facilitate the study of transient flow and unsaturated soil conditions, cases that attempt to analytically solve are complicated and laborious. However, it should be recognized that the computer programs now replaces the judgement of an engineer.


  1. 1. Lu N, Likos W. Unsaturated soil mechanics. 1st ed. Hoboken, New Jersey: John Wiley and Sons; 2004.
  2. 2. López-Acosta N P. Numerical modelling of water flow problems. In: Proceedings of the 27th National Meeting on Geotechnical Engineering. SMIG (November 20–21, 2014), Puerto Vallarta, Jalisco, Mexico; 2014. (in Spanish)
  3. 3. Aitchison G D, Richards B G. A broadscale study of moisture conditions in pavements subgrade throughout Australia. In: Moisture Equilibria and Moisture Changes in Soils Beneath Covered Areas, A Symposium in Print; Butterworths, Sydney, Australia. pp. 184–232.
  4. 4. Aslyng H C et al. Soil physics terminology, draft report. Int. Soc. of Soil Sci. 1962; Bull. 23. pp. 7.
  5. 5. Lee H C, Wray W K. Techniques to evaluate soil suction—a vital unsaturated soil water variable. In: Proceedings of the First International Conference on Unsaturated Soils; Paris, France: A.A. Balkema/Presses de l’ Ecole Nationale des Ponts et Chaussés; 1995.
  6. 6. Zepeda-Garrido J A, editor. Unsaturated soil mechanics. 1st ed. México D.F.: Mexican Society for Soil Mechanics, Autonomous University of Queretaro2; 2004. 322 p. (in Spanish)
  7. 7. Pérez-García N, Determination of Soil-Water Characteristic Curve with pressure plate test. Safandila, Querétaro, México: Secretaria de Comunicaciones y Transportes, Mexican Transportation Institute; 2008. 54 p. Publicación Técnica No. 313 (in Spanish)
  8. 8. Fredlund D G, Rahardjo H. Soil mechanics for unsaturated soils. 1st ed. New York: John Wiley and Sons; 1993. 521 p. DOI: 10.1002/9780470172759
  9. 9. Fredlund M D, Fredlund D G, Wilson G W. Prediction of the soil-water characteristic curve from grain-size distribution and volume-mass properties. In: 3rd Brazilian Symposium on Unsaturated Soils; April 22–25; Rio de Janeiro, Brazil; 1997.
  10. 10. Van Genuchten M T. A closed form-equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 1980;44:892–898.
  11. 11. Fredlund D G, Xing A, Fredlund M D, Barbour S L. The relationship of the unsaturated soil shear to the soil-water characteristic curve. Can. Geotech. J. 1996;33(3):440–448. DOI: 10.1139/t96-065
  12. 12. Aubertin M, Ricard J F, Chapuis R P. A predictive model for the water retention curve: application to tailings from hard-rock mines. Can. Geotech. J. 1998;35(1):55–69. DOI: 10.1139/t97-080
  13. 13. Mbagwu J S C, Mbah C N. Estimating water retention and availability in Nigerian soils from their saturation percentage. Comm. Soil Sci. Plant Anal. 1998;29(7–8):913–922.
  14. 14. Rojas E. Towards a Unified Soil Mechanics Theory. The Use of Effective Stresses in Unsaturated Soils. 1st ed. Sharjah, United Arab Emirates: Bentham eBooks; 2013. 190 p.
  15. 15. Arya L M, Paris J F. Physicoempirical model to predict the soil moisture characteristic from particle-size distribution and bulk density data. Soil Sci. Soc. Am. J. 1981;45:1023–1030.
  16. 16. Arya L M, Dierolf T S. Predicting soil moisture characteristics from particle-size distributions: And improved method to calculate pore radii from particle radii. In: van Genuchten M.Th., Leij F.J., Lund L.J., editors. Indirect methods for estimating the hydraulic properties of unsaturated soils; University of California, Riverside, CA. 1992. p. 115–124.
  17. 17. Basile A, D’Urso G. Experimental corrections of simplified methods for predicting water retention curves in clay-loamy soils from particle-size determination. Soil Technol., Elsevier. 1997;10(3):261–272. DOI: 10.1016/S0933-3630(96)00020-7
  18. 18. Fredlund M D. The role of unsaturated soil property functions in the practice of unsaturated soil mechanics [thesis]. Saskatoon, Saskatchewan, Canada: Department of Civil Engineering, University of Saskatchewan; 1999. 292 p.
  19. 19. Hosagasi-Fuselier T. Evaluation of soil water characteristic curves and permeability functions for modelling of seepage in unsaturated soils [thesis]. Somerville/Medford, Massachusetts, USA: Tufts University; 2006.
  20. 20. Scheinost A C, Sinowski W, Auerswald K. Regionalization of soil water retention curves in a highly variable soilscape, I. Developing a new pedotransfer function. Munich, Germay: Geoderma, 78, Elsevier Science B. V. 1997. p 129–143.
  21. 21. Aubertin M, Mbonimpa M, Bussière B, Chapuis R P. A model to predict the water retention curve from basic geotechnical properties. Can. Geotech. J. 2003;40(6):1104–1122. DOI: 10.1139/t03-054
  22. 22. Zapata C E. Uncertainty in soil-water characteristic curve and impacts on unsaturated shear strength predictions [thesis]. Tempe, Arizona: Arizona State University; 1999.
  23. 23. Siller W S. The mathematical representation of the soil-water characteristic curve [thesis]. Sakatoon, Saskatchewan, Canada: University of Saskatchewan; 1997.
  24. 24. Fredlund D G, Rahardjo H, Fredlund M D. Unsaturated soil mechanics in engineering practice. 1st ed. Hoboken, New Jersey: John Wiley and Sons; 2012. 857 p.
  25. 25. Brooks R H, Corey A T. Hydraulic properties of porous media. Fort Collins, Colorado: Colorado State University Hydrology Papers. 1964;(3):27.
  26. 26. Gardner W R. Some steady state solutions of the unsaturated moisture flow equation with application to evaporation from a water-table. Soil Sci. J. 1958;85(4):228–232.
  27. 27. Burdine N T. Relative permeability calculations from pore size distribution data. Trans. Metall. Soc. AIME. 1953;198:71–78.
  28. 28. Mualem Y. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 1976;12:513–522.
  29. 29. Leong E C, Rahardjo H. Permeability functions for unsaturated soils. J Geotech. Geoenviron. Eng. 1997;123(12):1118–1126.
  30. 30. Fredlund D G, Xing A. Equations for the soil-water characteristic curve. Can. Geotech. J. 1994;31(4):521–532. DOI: 10.1139/t94-061
  31. 31. Darcy H. Histoire des Foundataines Publique de Dijon. Paris: Dalmont; 1856. pp. 590–594.
  32. 32. Buckingham, E. Studies of the movement of soil moisture. Washington, DC: U.S. Department of Agriculture, Bureau of Soils; 1907. Bulletin No.: 38.
  33. 33. Richards L A. Capillary conduction of liquids in porous medium. Physics. 1931;1:318–333.
  34. 34. Fredlund D G, Krahn J. Comparison of slope stability methods of analysis. In: 29th Canadian Geotechnical Conference; Vancouver, British Columbia; 1976. pp. 56–74.
  35. 35. Fredlund D G, Krahn J, Pufahl D E. The relationship between limit equilibrium slope stability methods. In: Proceedings of the International Conference on Soil Mechanics and Foundation Engineering; Stockholm, Sweden; 1981. Vol. 40: pp. 409–416.
  36. 36. Morgenstern N R, Price V E. The analysis of the stability of general slip. Geotechnique. 1965;15(1):77–93.
  37. 37. Manual tailings dams and reservoirs (In Spanish). Mexico D.F.: Association of Mining Engineers, Metallurgists and Geologist of Mexico; 1993. 80 p. (in Spanish)
  38. 38. Mining waste storage. 1st ed. México D.F.: Mexican Society for Soil Mechanics; 2001. 137 p. (in Spanish)
  39. 39. Flores Castrellón O. Dynamic properties of tailings [thesis] (In Spanish). México D.F.: National Autonomous University of Mexico, Faculty of Engineering: 1996. (in Spanish)
  40. 40. Flores Castrellón O. Modulus and Poisson’s ratio dynamic measurements obtained on the fringe of the middle third in granular soil specimens [thesis]. México D.F.: National Autonomous University of Mexico, Faculty of Engineering: 2008. 123 p. (in Spanish). Available from:
  41. 41. GeoSlope International. Geostudio suite, Seep/W groundwater seepage analysis. Calgary, Alberta, Canada: 2007. Version: V7.23.
  42. 42. GeoSlope International. Geostudio suite, Slope/W slope stability analyses. Calgary, Alberta, Canada: 2007. Version: V7.23.
  43. 43. López-Acosta N P, Mendoza-Promotor J A. Water flow in partially saturated soils and its application in geotechnical engineering. 1st ed. Mexico City: Engineering Institute of the UNAM; 2016. 125 p. (In Spanish)
  44. 44. SoilVision Systems LTD. SVFlux 1D/2D/3D Saturated/Unsaturated finite element groundwater seepage modelling. Saskatoon, Saskatchewan, Canada 2009. Version: 2.4.26
  45. 45. Escobedo Rodríguez B. Stability analyses of the tailings dam No. 5, “La Negra” mining unit, Querétaro [thesis]. Mexico D.F.: Universidad Nacional Autónoma de México, Facultad de Ingeniería; 2011. 210 p. (in Spanish). Available from:
  46. 46. Zhang L, Fredlund M D, Fredlund D G, Haihua L. Comparison of 2-D and 3-D slope stability analyses for unsaturated soil slopes. Geo Regina; 2014.
  47. 47. SoilVision Systems LTD. SVSlope 2D/3D Limit equilibrium slope stability analysis; 2009. Version: 2.4.26.

Written By

Norma Patricia López-Acosta and José Alfredo Mendoza-Promotor

Submitted: 21 January 2016 Reviewed: 22 April 2016 Published: 27 July 2016