Experimental fitting models of the SWCC.
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
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) .
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 .
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 [3–5]; 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 , 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 , 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).
The SWCC represents a fundamental relationship that can be used to describe the behaviour of unsaturated soils . 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 .
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 . 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 [15–17].
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:
Saturated capillary zone. The soil zone that is maintained in a saturated state, whose defining limit coincides with the air-entry value , 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 . 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:
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 , 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 .
|Van Genuchten (1980)|
|Van Genuchten and Mualem (1976)|
|Van Genuchten (1980)|
|ws||saturated gravimetric water content|
|wr||residual gravimetric water content|
|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|
|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|
|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|
|mf||adjustment parameter related to the residual water content of soil|
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:
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  and Van Genuchten and Mualem .
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 .
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 .
A summary of the main estimation methods expressed as a function of hydraulic conductivity is presented in Table 2.
|Van Genuchten and Burdine (1953)|
|Brooks and Corey (1964)|
|Van Genuchten and Mualem (1976)|
|Fredlund and Xing (1994)|
|Leong and Rahardjo (1997)|
|hydraulic conductivity function|
|ks||saturated hydraulic conductivity coefficient|
|ψaev||soil suction at the air-entry value|
|pore size distribution index|
|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|
|upper integration limit (1,000,000 kPa)|
|fictitious variable that describes soil suction (logarithm scale)|
|variable derived from the function of soil storage capacity (characteristic curve)|
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
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 , the one-dimensional water flow through an unsaturated soil can be expressed as
is the hydraulic gradient
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 , 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:
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).
is the suction head
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 . 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:
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 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 . 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:
Eq. (13) is known as the Richards  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 :
is the shear force between slices
is the normal force
percentage (in decimal form) of the expressed function
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
However, the factor of safety with respect to the balance of forces is determined by
is the effective cohesion
is the effective friction angle
is the pore pressure
is the normal force at base of slice
is the weight of slice
is the load line
, , , , , are the geometric parameters
is the slope.
One important factor that is relevant for the previous equations is the normal force , which is a term defined as
5. Practical application: coupled water flow-slope stability analyses of a tailings dam
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.
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 [37–40]. 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.
|Fine material||Silt with sand||ML||1.582 × 10−8||0||25|
|Coarse material||Silty sand||SM||1.000 × 10−6||0||35|
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).
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  SWCC equation was applied.
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.
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 . 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.
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).
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  and from considering the results of the steady-state groundwater flow analyses.
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.
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.
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 . 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.
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 .
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.
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.
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 , and in Table 4, a comparison is made of the number of elements required for each model.
|3D realistic topography||2889||1580|
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.
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.
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.
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.
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.