Abstract
On the basis of the Boussinesq unsteady seepage differential equation, a new simplified formula for the phreatic line of slopes under the condition of decreasing reservoir water level is derived by means of the Laplacian matrix and its inverse transform. In this context, the expression of normal stress on the slip surface under seepage forces is deduced, and a procedure for obtaining the safety factors under hydrodynamic forces is proposed. A case study of the Three Gorges Reservoir is used to analyze the influences of the water level, decreasing velocity and the permeability coefficient on slope stability.
1. Introduction
During the normal operation of a dam, a decrease in water level is the main cause of bank slope failure [1]. There are many potential landslide sites that exist in the region of the Three Gorges Reservoir, and these are generally triggered by changes in the water level. Landslides cause huge losses both in terms of property and life. For example, the landslide that occurred on 13 July 2003 at Qiangjiangping on the Three Gorges Reservoir claimed many lives and was responsible for considerable infrastructure losses. The accident serves as a profound warning [2, 3] and strongly motivates further study of this type of natural disaster. Many studies have confirmed that changing water levels are the main cause of landslide events [2–9].
The Chinese government has invested considerable resources to ensure slope stability under normal operation and to prevent natural disasters such as landslides. During the landslide survey and design review of the Three Gorges Reservoir, some problems were identified and require further clarification. First, there is no scientific basis for the determination of phreatic lines. A decrease in water level is the most unfavorable situation in terms of slope stability and usually leads to landslides. Furthermore, a decrease in water level can cause seepage problems related to the rate of water level decrease and the slope’s permeability coefficient. Therefore, the correct approach should focus on the above-mentioned factors to determine the phreatic line and on the seepage pressure to analyze the slope stability. However, most of the determination of phreatic lines is based on the designer’s experience, which may lead to incorrect stability calculations and safety risk assessment. The second problem is related to the calculation of water forces, which is ambiguous. For example, both the seepage force and the surrounding hydrostatic pressure are considered in calculations, which cause repetition of the water forces in the equations and hence introduce more uncertainty. The third problem is related to the velocity of the decreasing water level, which is again not clearly defined in the slope stability calculation. In the case of power generation or flood control, the Three Gorges Reservoir releases water in a very short span of time. Such action leads to a rapid decrease in the water level and in some cases triggers landslides.
The determination of the phreatic line constitutes a free surface (unconfined) seepage problem in rock and soil mechanics. In such a case, the key calculation is to determine the free surface that delimits the flow boundaries. The free surface can be found using nonlinear numerical techniques, such as the finite difference method with adaptive mesh [10] and the finite element method with adaptive mesh [11] or fixed mesh [12]. Among the proposed methods, the extended pressure (EP) method proposed by Brezis and Kinderlehrer [13] and further simplified by Bardet and Tobita [14] using finite differences is considered the simplest and most efficient for free surface calculation. Through an extension of Darcy’s law, the EP method substantially reduces the variational inequalities, which can then be applied for computation of the whole domain. By applying variational inequalities, Zheng et al. [15] and Chen et al. [16] made significant contributions in solving the slope and dam free surface seepage problems. To improve the accuracy and computational efficiency of the EP method, an iterative error analysis is used with the finite difference equations [17]. Considerable progress has been made in calculating the slope phreatic line, especially through numerical analysis using the finite difference and finite element methods. However, such numerical methods are not commonly used in engineering practice and are usually ignored in soil mechanics, as they are obtained through complicated derivations and are challenging to implement. Therefore, there is a need for a simple and efficient method for practical engineering applications and educational training.
To explain the effects of water load on the slope stability to engineers and technicians, the surrounding water pressure is introduced by the Swedish Slice method calculated from the drift properties. Based on these experiments, it is concluded that the seepage force, water weight on the soil slices, and the surrounding hydrostatic pressure act as balancing forces for each other [18].
In slope stability analysis under a changing water level, the conventional slice method is the limit equilibrium method. The uncertainty is usually applied to the interslice forces after eliminating the slice bottom normal compressive force by the two equilibrium equations of a single slice [19]. Afterwards, the statically indeterminate problem of the slope limit safety factor can be solved using methods such as the Morgenstern-Price method [20], Spencer’s method [21], the Swedish method [22], the Bishop simplified method [19], or the Janbu simplified method [23]. Many slice methods have demonstrated that the interslice force is important in calculating the slope stability in statically indeterminate problems [24]. Such traditional slice methods are categorized as local analysis methods.
Alternatively, other limit equilibrium methods can be categorized as global analysis methods, including the graphic method [25], the variational method [26], and Bell’s global analysis method [24]. In contrast to the other limit equilibrium methods, Bell’s method considers the entire system instead of using a single slide bar as a research object. Therefore, the interslice force is not required in this method, which offers a new way to solve the strict slice method. Theoretically, the distribution of the slip surface normal stress in global analysis methods should be easier to understand than the distribution of interslice force in the Morgenstern-Price method. However, this method has not received sufficient attention, even in the summary provided by Duncan [27]. Until 2002, researchers in many studies [4, 28–30] have used similar methods. For example, in Bell’s derivation process, Zhu et al. [29] and Zhu and Lee [28] used quadratic interpolation to solve the slip surface normal stress distribution and derive one variable cubic equation with an unknown safety factor. Zheng and Tham [30] used the Green formula to convert the relevant domain integral calculus to boundary integral calculus.
Firstly, this paper applies the Boussinesq unsteady seepage basic differential equations and boundary conditions for derivation of the expression for the phreatic line during decreasing water level conditions. The polynomial fitting method is used to simplify the formula because of its lower complexity. For solving the effects of water load on the slope stability, the surrounding water pressure is calculated from the drift property. Finally, the global analysis method is proposed to analyze the slope stability under seepage forces. A typical landslide in the Three Georges Reservoir is used as a case study to analyze the influential factors for slope stability during conditions of decreasing water level.
2. Phreatic Line Calculations
2.1. Fundamental Assumptions
(1)The aquifer is homogeneous and isotropic with infinite lateral extension;(2)The phreatic flow parallel to slope surface is caused by water level fluctuation and the rainfall infiltration causes the phreatic flow to be perpendicular to slope surface;(3)The reservoir water level is decreasing at a constant speed of ;(4)The reservoir bank is considered as a vertical slope. The reservoir bank within the declining amplitude is much smaller than the ground, and in order to simplify this, it is considered as vertical reservoir bank [5].
As shown in Figure 1, based on the Boussinesq equation and following the above-mentioned assumptions, the differential equation of diving unsteady motion can be expressed aswhere and are two directions of local coordinates, as shown in Figure 1; is the aquifer thickness; is time; is the permeability coefficient; and is the specific yield, which is defined as the amount of water released as a result of gravity in a unit volume of saturated rock or soil and can be expressed as the ratio between the water volume released by the gravity and the rock or soil volume in saturated conditions.
Eq. (1) is a second-order nonlinear partial differential equation with no general analytical solution; however, it can be simplified into a linear problem. The simplified method assumes the aquifer thickness to be a constant, and the value is selected as the mean value of the initial and end diving thickness. Therefore, the simplified one-dimensional unsteady seepage flow equation can be rewritten aswhere ; the value of will be discussed in Section 2.3.
2.2. Models Used in This Study
As shown in Figure 1, at the initial time, , and by assuming the change of water level of the position at time , the water level of each point, , can be calculated aswhere is the reservoir water level change at time and position in local coordinates, is the water level slope in local coordinates, and is the reservoir water level in local coordinates. The change of water level at is .
The reservoir water level decreases at a speed . After the lateral seepage, at , ; as , .
From (3), the above semi-infinite aquifer groundwater unsteady flow during the decreasing water level can be summarized by the following mathematical equations:
By applying the Laplace transform to (4), solving the equations, and applying the inverse transform [31, 32], the solution of the differential equation can be mathematically expressed aswhere and is the residual error function.
By assuming (6) and following Figure 1 and (7), the equation can be rewritten as (8):which is the formula of the phreatic line of the slope when the reservoir water level declines at a constant speed . The function can be calculated by (6), and its curve is shown in Figure 2. According to (6), integral calculation is required, which is difficult to perform and not convenient for engineering applications. To obtain a practical expression for engineering applications, a polynomial fitting of (2) is used, as shown in Figure 2, which is simplified to [18]
Afterwards, the simplified formula of the phreatic line under a constant declining speed can be obtained. The results from this formula are consistent with those from the finite element method under the same conditions, which verifies the accuracy of this formula. The expression can be written aswhere , is the permeability coefficient (m/d), is the aquifer thickness (m), is the specific yield, is the time of the water level decrease (d), and is the reservoir water level before it began to decline.
2.3. Determination of the Aquifer Thickness
For solving (1), the value is selected as the mean value of the initial and final diving flow thickness to replace during the linearization, and is referred to as the aquifer thickness of the diving stream. This assumption is reasonable when the drop height is relatively small but causes larger error when the drop height is relatively large. Thus, the following method is used to determine the aquifer thickness. For the irregular impermeable layer, as shown in Figure 3, the average aquifer thickness is not as simple as above, and this often causes a problem in engineering practice. As a result, it is necessary to propose a method that is convenient for use in engineering applications.
For depicting the geological conditions of the Three Gorges Reservoir, the land surface component usually consists of weathered soils and collapsed accumulated bodies overlying rock. Such geological conditions lead to gravelly soil forming in the upper part and a geological rock structure in the lower part. The permeability coefficient of the rock is very small compared to the soil. Thus, the former can be neglected, and the combined surface between the rock and soil can be used as an impermeable layer. The cross-section is shown in Figure 3. The calculation area is restricted by the rock (impermeable layer); however, when the reservoir water level drops, the influential area cannot be used as the horizontal impermeable layer. In this case, the following method is promising for determining the average thickness of the aquifer.
Actual landslides do not have vertical slopes, but rather slope gradients. When the reservoir water level decreases, the overflow points are usually located on the steep slope surface and the flow follows the slope gradient. In contrast, the phreatic line on the slope is relatively flat. In such situations, the initial water level and slope line can replace the phreatic line in the calculation of the average aquifer thickness, which can be calculated by the following equation:where is the area consisting of the slope surface , the initial water level , the bedrock surface , and the cross point between the water level drop and the slope surface ; and is the intersection of the initial water level with the bedrock surface. The geological section is often drawn using CAD software [5, 18].
2.4. Specific Yield
Specific yield is an important hydrogeological parameter that should be obtained through actual measurements. It is defined as the amount of water released as a result of gravity in a unit volume of saturated rock or soil and can be expressed as the ratio between the water volume released by the gravity and the rock or soil volume in saturated conditions.
In rock and soil, only some pores are filled by water and movement is possible through connected pores. Such pores are generally termed as effective pores in hydrogeology. In hydrology or groundwater dynamics, the effective porosity is the ratio between the pore volume filled by water and the water released from gravitational pull and the total pore volume of these pores. The effective porosity is generally controlled by the physical properties of the rock/soil. It is often expressed in terms of percentage.
According to the test data of gravel and clayey soil [33], the empirical formula of specific yield iswhere is the porosity and is the permeability coefficient (cm/s). In the absence of experimental data, (12) can be used to calculate the specific yield.
3. The Normal Stress in the Slip Surface Induced by Hydrodynamic Forces and Surface Water Forces
3.1. Seepage Force Calculation
The drag force on the soil by water under seepage is known as the seepage force and is often described as a dynamic water force in engineering. Seepage force calculation is the key factor for evaluating the slope stability under seepage. Thus, the accuracy of the seepage force calculation strongly influences the results. It appears that engineers are often conceptually confused about the seepage force and the surrounding hydrostatic force calculations, in which the equation of the water forces is repeated, and the influence of the void ratio on seepage force calculations. To clarify these issues, the calculation of the seepage forces can be studied in terms of the water forces using the boundary of each differential slice; this is the most direct way to analyze the seepage forces [18].
The force diagram of a single soil slice is presented in Figure 4. In the diagram, and are the saturated soil gravity components above and below the phreatic line, respectively. The force is the hydrostatic force at the border; is the hydrostatic force at the border; is the hydrostatic force at the border; is the normal force at the base, which is also known as the effective force; is the angle between the bottom surface and the horizontal direction; and is the angle between the phreatic line and the horizontal direction.
To determine the hydrostatic forces , , and on the , , and borders, the flow grid properties based on the nature of drift and flow lines perpendicular to the equipotential lines are used. In Figure 5, and are perpendicular to the phreatic line, is perpendicular to , and is perpendicular to . Thus, the hydraulic head in and the hydraulic head in C can be expressed according to geometric relationships:The resultant force of the hydrostatic force on the and borders can be represented byThe hydrostatic force on the sliding surface is represented byThe forces in the vertical and horizontal directions can be, respectively, expressed asThe water weight in the soil iswhere is the water gravity below the phreatic line of the soil slice.
By assumingwe obtain
All the water loads below the phreatic line in the soil slice are shown in Figure 6. The -direction component of all the water loads isThe -direction component of all the water loads isFrom the geometrical relationships, as shown in Figure 6, the equation can be represented asTherefore, the forces from all water loads become
The geometrical interpretation of the full water immersion saturated area in the soil slice is that the water gravity and hydraulic gradient are equal to the seepage force and hydrodynamic force, respectively. The direction of this force is the same as the fluid flow direction, and the angle between the force and the horizontal direction is .
This calculation demonstrates that the seepage force is the same as the resultant forces of water gravity and surrounding hydrostatic force of the soil below the phreatic line. Therefore, when the seepage force is used to express the safety factor, the natural gravity is used above the phreatic line, whereas the floating gravity, in this case, is used below the phreatic line. Therefore, the calculation diagram shown in Figure 2 can be replaced by that in Figure 6. Additionally, the water force and gravity of the soil slice can be replaced by the seepage force . As a result, this type of problem is straightforward to solve in practice.
3.2. Normal Stress on the Slip Surface under the Water Loads
The forces along the arc length ds in the antislide direction of the differential slice are shown in Figure 7. The parameters and are the interslice force increments in the horizontal and vertical slice, respectively; and are the bar weight and seismic force, respectively; and and are the components in the horizontal and vertical directions on the outer contour line , respectively. The relationship between the normal face force and tangential force iswhere and are the and directional components along the differential arc length in , as shown in Figure 7.
All the forces are projected on the normal direction of the slip surface; thus, the equation becomeswhere is the bar height above the phreatic line, is the average bar weight, and is the effective average weight.
Substituting (24), (26), and (27) into (25), where , we obtainwhere and are the contributions from the interslice forces and external load to the normal stress in the slip surface, respectively, and are generally expressed as a function of The stress is decomposed intowhere is the contribution from the slider body volume force to the normal stress in the slip surface:Here, is the contribution from the water pressure of the outer contour slope slip surface to the normal stress in the slip surface:where is the outer contour slope surface at point .
Usually the interslice forces and at both ends of the sliding surface are zero; however, their derivatives and may not be zero. Additionally, from (28) and (29), the slip surface normal stress at both ends may not be zero either. Therefore, it is not necessary to satisfy the zero condition at the slip surface ends when determining the normal stress at the slip surface. Hence, the distribution of the slip surface normal stress approximation can be determined as shown in (33):where is the correction function for the slip surface normal stress and and are the two undetermined parameters. The reason for introducing two undetermined parameters is that three equilibrium equations are used to solve for three unknowns, in which the safety factor is one of the unknowns.
In this study, is defined as a linear function [30]:where and are the coordinates of two endpoints on the slip surface in this two-dimensional problem.
4. Slope Stability Analysis Method under Water Loads
As shown in Figure 8, the slip body is defined as a planar region Ω consisting of an outer contour line ACB and a potential slip surface ADB. The background materials are generally formed by a variety of geotechnical movements at the landslide sites. The active forces of the landslide include volume forces (e.g., weight and horizontal seismic forces) and the concentrated force or surface force on the outer contour ACB; the constrained forces include the slip surface normal force and the tangential force .
By taking any three points , which are not in the same straight line as the moment center, and with Ω in equilibrium, the sum of the three moment centers is always equal to zero and can thus be expressed aswhere is the resultant moment of all active forces on the landslide body ; ; and and are the vector components on the slip surface at position ():
If not specifically stated, all subscripts in the equations are free subscripts. In the formula, will take the value of 1, 2, and 3 and provide three moment centers (, ). Additionally, to simplify the expression and by assuming the slope is the right-side slope, it is found that an increase in slope height occurs with an increasing coordinate value.
By assuming that the sliding surface meets the Mohr-Coulomb criterion and the landslide body is in the limit equilibrium state, the equation can be written aswhere is the safety factor and and are the effective shear strength parameters.
By substituting (36) into (35), the unknown function integral equations of the sliding surface normal stress can be obtained aswhereBy substituting (33) into (37), (38) can be derived:where is the third-order nonlinear vector function of and and ~ are the six third-order vectors, which can be defined as follows:Equation (38) can be used as a quasi-Newton method to solve [34]. The required Jacobian matrix of in the solution can be defined as :where the three column vectors are expressed as
In this study, the following condition can be used to terminate iteration during numerical simulations and can be expressed aswhere is the difference between the two consecutive iteration steps and is the selected tolerance level.
In the following examples, is used as a damped Newton’s method [34]. The landslide border is dispersed into several small segment grids. In the smaller grid segments, the slip surface normal stress is assumed to be constant.
5. Influence Factors of Slope Stability
For ease of analysis, the drop height of the reservoir water is assumed to be ; (8) can then be rewritten as follows:
These equations indicate that the factors influencing the phreatic line under a decreasing water level include the permeability coefficient , the specific yield , the decreasing speed of the water level , the aqueous layer thickness , and the decreasing height . From Figure 2, is a decreasing function inversely related to . A higher value causes a reduction in the speed of the free surface water flow and vice versa. Therefore, when , and the speed of the free surface water on the slope and the water level become equal. When , and the free surface water on the slope and the water level are unchanged. It is clear from Figure 2 that when , becomes 0.
The higher water level on the slope causes an unfavorable condition for the slope stability. Based on the above analysis, the smaller is, the faster the free surface water flow on the slope is, which means the free surface water decreases with decreasing water level and vice versa. Similarly, the larger values cause unfavorable conditions for the slope stability. The favorable factors for the slope stability are the decrease in , directly related to the decreasing rate of the water level , the permeability coefficient , the specific yield , and the water thickness . Most engineers generally analyze the effects of decreasing rate of water level and permeability coefficient on the slope stability.
The Woshaxi landslide occurred on 13 July 2003 at the right bank of Qinggan River, a tributary of the Yangtze River. The Qinggan River flows from west to east and passes in front of the landslide site. The landslide was located 6 km from the mouth of the river. The distance from the downstream Qinggan River to the Qianjiangping landslide site is approximately 1.5 km (Yin et al., 2012; Wang and Li, 2007; Wand and Yang, 2006; Wen et al., 2008; Wu et al., 2006; Yin and Peng, 2007; Zhang et al., 2004; Jian et al., 2014); [2–4, 6]. The distance to the Three Gorges dam site is 50 km. The geographical position of the landslide is shown in Figure 9.
The geological cross-section is shown in Figure 10, and the geological materials are listed in Table 1. The Three Gorges Reservoir water level is found to be between 145 m and 175 m. This study evaluates the landslide stability associated with the decreasing water level from the height of 175 m. Below the phreatic line, the soil saturation density is used, and for the natural strength of slip, the soil above the phreatic line is used.
Using different decreasing water level speeds, the other related parameters are shown in Table 1. Figure 11 shows the calculation process for the phreatic line under different decreasing water level rates. Figure 11(a) indicates the change in the phreatic line for different water level conditions and a decreasing water level rate () of 0.1 m/d. Figure 11(b) shows the phreatic line for different water levels at a decreasing water level rate of m/d. Figure 11(c) shows the phreatic line for different water level conditions and a decreasing water level rate of m/d. The safety factors under different water levels at different decreasing water level rates are shown in Figure 12. When the decreasing water level rate is relatively small, there is sufficient time for the water in the slope to discharge from the slope, which causes a slow change in the phreatic line, as shown in Figure 11(a). When the decreasing water level rate is relatively high, there is insufficient time for the water on the slope to discharge, causing a sudden change in the phreatic line, as shown in Figure 11(c). The decreasing water level rate has a profound effect on the phreatic saturation lines and therefore on the slope safety factors. As shown in Figure 12, a more obvious change is observed in the safety factors as the water level slightly decreases; however, no significant changes are observed in the safety factors when the water level decreases by a substantial margin. Overall, as the decreasing water level rate increases, it causes a decrease in the safety factor. It is further observed that at a given decreasing water level rate, the safety factor tends to be more stable with an increase in water level. It is found that when the height of the reservoir water level decreases from 175 m to 145 m, a 5% change is observed in the safety factor at decreasing water level rates of m/d and m/d.
(a) = 0.1 m/d
(b) = 0.4 m/d
(c) = 1.2 m/d
In the current study, the phreatic line is calculated by varying the slope permeability coefficients at a constant water level decreasing rate; that is, m/d. The other parameters required for calculating the phreatic line are given in Table 1. Figures 13(a)–13(c) and 14 show the phreatic line and safety factors, respectively, under different water levels when the permeability coefficient is 0.01 m/d, 0.10 m/d, and 1.00 m/d, respectively. It can be observed that at the same water level decreasing rate, the variation in causes differences in phreatic line behavior. With a relatively small , there is a sudden change in the phreatic line, as sufficient time is not available for water to be transported down the slope (Figure 13(a)). However, with increasing , the change in the phreatic line is gradual because more time is available for the water to be transported down the slope (Figures 13(b)-13(c)). The permeability coefficients have a profound effect on the phreatic lines and the slope safety factors. As shown in Figure 14, significant changes in the safety factors occur when the water level decreases slightly with different slope permeability coefficients. In contrast, the safety factor does not change substantially with large changes in the water level. Overall, the slope permeability coefficient is inversely proportional to the safety factor, which indicates that an increasing permeability coefficient can cause a decrease in the safety factor. A larger slope permeability coefficient implies a higher safety factor. It is also evident from this study that a decreasing water level leads to a more stable safety factor, which can be ascribed to the slope permeability coefficient. A larger slope permeability coefficient is associated with a greater probability of penetration with decreasing water level, which is favorable for slope stability. It is clear from the results that when the height of reservoir water level decreases from 175 m to 145 m, the difference in the safety factor is only approximately 6% at m/d and m/d.
(a) = 0.01 m/d
(b) = 0.10 m/d
(c) = 1.00 m/d
6. Summary
(1)The expression of the phreatic line is based on the Boussinesq unsteady seepage basic differential equations and boundary conditions. The polynomial fitting method is used to obtain a simplified formula that is convenient to use in engineering applications.(2)To simplify the explanations and provide greater understanding to engineers, the flow grid properties are used to determine the water forces, which demonstrates that the seepage and resultant force caused by the water weight in the soil and the surrounding hydrostatic forces balance each other.(3)The expression of the slip surface normal stress under the dynamic water forces is derived, and an equation for solving the safety factor under different dynamic water forces is proposed.(4)Combined with a practical project in the Three Gorges Reservoir, the effects of decreasing water level rate and the slope permeability coefficient are analyzed and presented to describe slope stability.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This study is supported by the National Natural Science Foundation of China [Grant no. 11572009] and the CRSRI Open Research Program [Program SN: CKWV2015212/KY] and the Youth Innovation Promotion Association CAS.