Next Article in Journal
Multiple-Frequency Force Estimation of Controlled Vibrating Systems with Generalized Nonlinear Stiffness
Next Article in Special Issue
An Analytical Framework for Assessing the Unsaturated Bearing Capacity of Strip Footings under Transient Infiltration
Previous Article in Journal
Fuzzy Domination Graphs in Decision Support Tasks
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Global Sensitivity Analysis Method for Embankment Dam Slope Stability Considering Seepage–Stress Coupling under Changing Reservoir Water Levels

1
College of Water Conservancy and Hydropower Engineering, Hohai University, Nanjing 210098, China
2
State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering, Hohai University, Nanjing 210098, China
3
College of Ocean Science and Engineering, Shanghai Maritime University, Shanghai 201306, China
*
Authors to whom correspondence should be addressed.
Mathematics 2023, 11(13), 2836; https://doi.org/10.3390/math11132836
Submission received: 12 May 2023 / Revised: 21 June 2023 / Accepted: 22 June 2023 / Published: 24 June 2023

Abstract

:
Ensuring the long-term, efficient, and safe operation of reservoir dams relies on the slope stability of embankment dams. Periodic fluctuations of the reservoir water level due to reservoir scheduling operations make the slope of the reservoir bank vulnerable to instability. To investigate the influence of various factors and their interactions with embankment dam slope stability under changing reservoir water levels, a global sensitivity analysis method is proposed that accounts for seepage–stress coupling. An embankment dam in Shaanxi Province, China, is studied as an example, with COMSOL Multiphysics software simulating the seepage and slope stability of the dam under fluctuating reservoir water level conditions and seepage–stress coupling. The global sensitivity analysis of factors affecting dam slope stability is accomplished by combining Plackett–Burman and Box–Behnken experimental designs, with ANOVA determining the sensitivity of each factor and interaction term. The results demonstrate that during the impoundment period of the reservoir, the saturation line is concave, and the overall stability safety of the dam slope increases first and then tends to be stable, according to the coefficient. The internal friction angle φ, cohesion c, and soil density ρs represent the three most sensitive factors affecting the stability and safety of the dam slope, while c × ρs is a second-order interaction term with significant sensitivity to the stability and safety coefficient of the dam slope. The reservoir drainage period infiltration line is convex, and dam slope stability first reduced and then increased. The magnitude of water level change H, internal friction angle φ, cohesion c, and soil density ρs are the four most sensitive factors for the coefficient of safety of dam slope stability, while c × ρs, H × ρs, and φ × ρs are the second-order interaction terms with significant sensitivity to the coefficient of safety of dam slope stability. These research findings and methods can offer valuable technical support and reference for the investigation and evaluation of the stability of embankment dam slopes.

1. Introduction

Water conservancy engineering represents a crucial means of improving water resource regulation capabilities and promoting the development of green energy. Dams constitute the central pillar of such engineering, playing an integral role in power generation, flood control, and water supply as water resources continue to be developed. China boasts a vast array of reservoirs, with initial statistics pointing to the construction of over 98,000 such facilities, of which embankment dams comprise over 90% of the total number due to the advantages of locally sourced materials and ease of construction [1]. Sliding instability is one of the most common failure modes among embankment dams constructed [2]. Therefore, the slope stability of embankment dams in service has been a major focus of scholarly attention and research [3,4,5,6,7,8,9]. However, the stability of a dam slope is determined by numerous factors, in addition to internal factors such as the physical parameters of the soil body. During normal operation of embankment dams, the reservoir level is often in a state of flux due to water storage or discharge, and the mutual replenishment of water inside and outside the waterfront of the dam slope causes the internal seepage field of the dam to constantly change which, in turn, affects the overall or local stability of the dam [10,11,12]. Hence, slope stability is often determined by a combination of internal and external factors. However, the degree of influence of these factors on the slope is complex, ambiguous, and variable [13,14,15,16]. Sensitivity analysis of slope stability factors involves quantitatively analyzing the correlation between the factors that affect slope stability and slope safety coefficient, studying the influence of changes in each factor on slope safety coefficient, identifying the dominant factor of slope instability, and providing targeted guidance for the prevention and control of slope instability disasters. Therefore, studying the slope stability sensitivity of embankment dams under conditions of changing reservoir water levels is of great significance in ensuring the long-term, efficient, and safe operation of reservoir dams.
In order to accurately determine the sensitivity of each factor to embankment dam slope stability under varying reservoir water levels, it is essential to establish a reasonable intrinsic structure model. To simplify the calculation, seepage and stability problems in slope stability analysis are often considered separately for many embankment dams [17,18,19,20]. That is, seepage calculation is conducted first, and dam stress and slope stability analysis are performed based on the results of the seepage calculation. However, engineering practice has shown that seepage and stress fields interact with each other [21,22,23], and changes in the seepage field due to reservoir water level fluctuations affect the distribution of pore water pressure acting on the soil; in addition, the effective stress transferred between the soil skeleton changes, which directly impacts the stress and deformation of the dam. Consequently, this leads to changes in the soil’s pore and permeability properties which, in turn, affect the seepage field inside the dam [24]. Therefore, the problem of dam slope stability caused by the fluctuation of reservoir water level is a complex nonlinear problem. If the seepage field and stress field of the dam are considered separately, the simulation results are likely to be inaccurate. On the other hand, studying the mutual coupling of the seepage and stress fields can accurately describe the dynamic changes of the internal pore pressure field of the dam body during the water level fluctuation and its influence on the stability of the dam slope. This approach can reveal the mechanism of embankment dam landslide disasters and provide calculation results that are more consistent with actual situations.
In recent years, sensitivity analysis has become increasingly popular in slope stability studies. Hu [25] assessed the sensitivity of cohesion, internal friction angle, and gravity to the stability safety factor using partial Spearman rank correlation coefficients. Ahmad [26] used a tree-augmented naive Bayes (TAN) model to study the sensitivity of cohesion, internal friction angle, pore pressure ratio, slope angle, heaviness, and slope angle on slope stability. Cross [27] investigated the sensitivity of physical and mechanical parameters of soil to the coefficient of safety of shallow plane sliding of residual soil based on the ISSA model. Ghadrdan [28] analyzed the sensitivity of material properties, including internal friction angle, cohesion, permeability, and creep, to slope stability in a deep foundation pit using both the limit equilibrium method (LEM) and the finite element method (FEM). Karthik [29] explored the sensitivity of soil strength parameters, soil heaviness, Young’s modulus, Poisson’s ratio, slope angle, slope height, and depth of subgrade below the slope to the stability factor of safety of homogeneous soil slopes by utilizing the finite element method. However, most of the previous studies have only conducted local sensitivity analyses, focusing on the effects of individual soil property parameters on the analyzed indices, while ignoring the interactions between these factors. Although this approach can save researchers a lot of time and resources, it has limitations in terms of linearity, normality, and local variability and cannot fully capture the nonlinear mapping and interaction relationships between sensitive factors and response indicators. However, in the actual operation of dams, uncertainty problems arise due to complex changes in water level environments and the synergistic effects of multiple factors, which cannot be overlooked. Therefore, global sensitivity analysis that considers the interaction effects among factors is crucial in slope stability analysis.
In summary, the current literature primarily focuses on local sensitivity analysis of the impact of soil parameters on slope stability, while global sensitivity analysis of dam slope stability considering the coupling effect of reservoir water level changes and seepage–stress has been rarely reported. To investigate the influence of various factors and their interactions on dam slope stability under the conditions of reservoir water level fluctuation in a more comprehensive manner, this paper proposes a global sensitivity analysis method that considers the coupling of seepage and flow stress. The method combines Plackett–Burman and Box–Behnken experimental designs to conduct global sensitivity analysis of factors affecting dam slope stability, and the sensitivity of each factor and its interaction term is determined by analysis of variance (ANOVA). The proposed method is expected to provide significant technical support and theoretical reference for studying and evaluating the slope stability of embankment dams under variable water level conditions.

2. Methodology and Theory

2.1. Theory

2.1.1. Coupled Model of Seepage and Stress

Water flow within an embankment dam presents a complex saturated–unsaturated seepage challenge. To tackle this intricate issue, the saturated–unsaturated Richards equation is employed in this study to precisely describe the seepage field of the dam [30]:
Q m ρ ( C m ρ g + S e S ) p t ρ ( k s k r ( θ ) μ ( p + ρ g z ) ) = 0
where ρ is the density of water, kg/m3; Cm is the specific moisture absorption, m−1; g is the acceleration of gravity, m/s2; Se is relative saturation; S is the elastic water storage rate, Pa−1; p is pressure, Pa;   is the Laplace operator; θ is water content, m3/m3; ks is the permeability coefficient of saturated soil, m/s; kr(θ) is the relative hydraulic conductivity of unsaturated soil, which is a function of water content θ, m/s; μ is the dynamic viscosity of water, Pa s; z is the base elevation, m; and Qm is the source and sink term of the seepage field, kg/(m3s).
The key to addressing unsaturated seepage field challenges lies in predicting the hydraulic conductivity of unsaturated soils using mathematical models [31,32,33,34]. The van Genuchten model [35,36,37] is frequently utilized for studying unsaturated infiltration in different porous media due to its strong agreement between the fitting of the hydraulic characteristic curve and measured data [38]. In this study, the hydraulic conductivity coefficients are derived by adopting the van Genuchten model.
θ = { θ r + S e ( θ s θ r ) , p < 0 θ s , p 0 }
S e = { 1 [ 1 + | α H p | n ] m , p < 0 1 , p 0 }
C m = { α m 1 m ( θ s θ r ) S e 1 m ( 1 S e 1 m ) m , p < 0 0 , p 0 }
K r = { S e l [ 1 ( 1 S e 1 m ) m ] 2 , p < 0 1 , p 0 }
Here, θ, θs, and θr are the water content, saturated water content, and residual water content, m3/m3, respectively; α is the reciprocal of the water characteristic curve taking water, 1/m; Hp is the pressure head; and m and n are indicators of the slope of the water characteristic curve, obtained by fitting SWCC, m = 1 − 1/n.
This study establishes the intrinsic and equilibrium equations for coupled seepage–stress analysis based on the effective stress principle. Drawing on the Bishop unsaturated effective stress principle and assuming that pore gas pressure equals atmospheric pressure, the Bishop factor represents soil saturation, enabling the expression of saturated–unsaturated effective stress:
σ = σ + p w S e
where σ is the total stress, Pa; σ′ is the effective stress, Pa; pw is the pore water pressure, Pa; and Se is the soil saturation (for saturated soil Se = 1, for dry soil Se = 0; usually, 0 < Se < 1).
An elastoplastic intrinsic model is used to describe the stress-deformation behavior of the rock mass, and the Mohr–Coulomb criterion is used for the elastic–plastic analysis of the dam [39]. The Mohr–Coulomb yield function and the corresponding plastic potential are
f = ( cos Θ 1 3 sin φ sin Θ ) J 2 C cos φ 1 3 I 1 sin φ
Here,
Θ = 1 3 sin 1 ( 3 3 2 J 3 ( J 2 ) 3 )
where f is the yield function, φ is the internal friction angle,  Θ [ π 6 , π 6 ] , and J1, J2, and J3 represent the first invariant of effective stress tensor and the second and third invariants of effective stress partial tensor, respectively.
The coefficient of safety for dam slope stability is calculated using the finite element strength reduction method [40], which divides the shear strength of the geotechnical material by the reduction factor for finite element elastoplastic calculations and uses the reduction factor when the structure reaches damage as the safety factor [41]. The expressions for the reduced strength parameters c and φ of the geotechnical body are calculated as:
c = C F O S , φ = tan 1 ( tan φ F O S )
where FOS is the factor of safety, c is the cohesion, and φ is the angle of internal friction of the soil.

2.1.2. Plackett–Burman Experimental Design

The Plackett–Burman design [42] is a highly effective screening method for evaluating variables and identifying important factors by filtering out the factors that have a significant impact on the studied subject from a large number of factors through a small number of experiments. In Plackett–Burman testing, the results being examined are referred to as indicators; the parameters that may influence the test indicators are referred to as factors; and the specific test conditions to be compared in the test for each factor are referred to as levels. The Plackett–Burman method designs the test using two levels of high (+1) and low (−1) for each factor, and the significance of the factors is determined by comparing the difference between the two levels of each factor with the overall difference. The Plackett–Burman experimental design is generated according to a set of rules, which can result in nonunique arrangements. In actual studies, it is necessary to retain more than three dummy variables. Specifically, if M represents the number of trials, the number of factors studied should be less than or equal to M-4. Moreover, M should be a multiple of four, excluding the power of two values. Commonly used values for M include 12, 20, 24, 28, and 36.

2.1.3. Box–Behnken Experimental Design

The Box–Behnken experimental design [43] is a commonly used statistical approach in response surface methodology that can solve multivariate problems with three to seven factors through a three-level experimental design. With this method, the functional relationship between factors and response values can be fitted using a multivariate quadratic variance analysis with fewer trials. The Box–Behnken design requires at least three factors, and Figure 1 provides a visualization of the three-factor design. This spherical design conforms to rotatability or almost rotatability, meaning that any point within the test region is equidistant from the design center point. Furthermore, all test points are located on equidistant endpoints, and the tests generated at the vertices of the cube by the high and low levels of each variable are excluded.

2.1.4. Analysis of Variance

To distinguish the impact of changes in factor levels from test errors and assess the sensitivity of factors, this study utilized analysis of variance (ANOVA) [44] to analyze the test results obtained from the two experimental design methods. The ANOVA technique decomposes the total sum of squared deviations in test data into two components: the sum of squared deviations caused by each factor and the sum of squared deviations caused by the error. It then computes the F-statistic and performs an F-test to determine the degree of influence of each factor on the test index. The ANOVA method follows the principles described below:
S T = S T i = j = 1 r g = 1 k ( x j g i ) 2 ( j = 1 r g = 1 k x j g i ) 2 r g , ( i = 1 , 2 , , m )
S i = j = 1 r ( g = 1 k x j g i ) 2 g ( j = 1 r g = 1 k x j g i ) 2 r g , ( i = 1 , 2 , , m )
S e = S e i = j = 1 r g = 1 k ( x j g i ) 2 j = 1 r ( g = 1 k x j g i ) 2 g , ( i = 1 , 2 , , m )
where ST is the total sum of squared deviations; Si is the sum of squared deviations of the factors; Se is the sum of squared deviations from the error;  S T = i = 1 m S i + S e ; k is the total number of trials of factor i at level j x j g i  is the test result corresponding to the gth trial of factor i at level j; and  { x j g i | i ( 1 , m ) , j ( 1 , r ) , g ( 1 , k ) } . df is the degree of freedom of each factor, where dfT is the degree of freedom of ST d f T = r g 1 ; dfi is the degree of freedom of Si d f i = r 1 ; and dfe is the degree of freedom of Se d f e = r ( g 1 ) .
The mean square of each factor is  M S i = S i / d f i . The mean square of the error is  M S e = S e / d f e . The value of the F-test sensitivity statistic Fi is  F i = M S i / M S e . The corresponding F0.01 and F0.05 values were verified in the F distribution table and compared with the calculated F values to determine the significance level of the factors. If  F i > F 0.01 , it indicates that the factor has a highly significant effect on the test results, which is noted as “**”. If  F 0.01 > F i > F 0.05 , it indicates that the factor has a significant effect on the test results, which is noted as “*”. If  F i < F 0.05 , it indicates that the effect of the factor on the test results is not significant, and “*” is not used.

2.2. Global Sensitivity Analysis Method Considering Seepage–Stress Coupling

This paper proposes a global sensitivity analysis method that considers seepage–stress coupling to investigate the effects of various factors and their interactions on the stability of dam slopes under variable reservoir water levels. To illustrate the approach, an embankment dam project in Shaanxi Province, China, is employed as an example, and the method is depicted in Figure 2. It is important to note that due to the large number of factors being examined in the global sensitivity analysis, it is not feasible to consider the interaction of all factors simultaneously. Therefore, the proposed method utilizes the Plackett–Burman test design to screen the significance of all factors, followed by the Box–Behnken test design in the response surface method to conduct a global sensitivity analysis on the screened significant factors. This approach aims to obtain the sensitivity of each significant factor and the interaction term between them on the coefficient of safety of dam slope stability. The specific operational steps are outlined below:
1. Geotechnical partitioning and physical parameters of the dam and the dam site were determined based on engineering design data and hydrogeological conditions. A typical design section of the dam was obtained, and the actual operating conditions of the project’s water level were determined.
2. A three-dimensional finite element model of the dam and the dam site was constructed. Material parameters and boundary conditions of the model were determined, and COMSOL Multiphysics software was utilized to analyze the seepage behavior and slope stability of the embankment dam under the coupling effect of reservoir water level fluctuation and seepage–stress.
3. The Plackett–Burman experimental design was utilized to preliminarily screen the eight factors affecting dam slope stability. Combining ANOVA, factors with insignificant sensitivity to the safety coefficient of dam slope stability were filtered out, and significant factors were retained.
4. The Box–Behnken test design in the response surface method was utilized to conduct main effect analysis and interaction effect analysis on the significant factors screened by the Plackett–Burman test design. Contour plots were drawn, and factors with significant sensitivity to the safety coefficient of dam slope stability and their interaction terms were found by combining the ANOVA method.

3. Study Area

3.1. Project Overview

The embankment dam project examined in this study is situated in Shaanxi Province, China, with a medium scale and engineering grade III. The reservoir’s typical operating level is 1082.00 m, with a corresponding storage capacity of 6796.5 million m3, a dead water level of 1066.05 m, and a storage capacity of 996.5 million m3. The design flood level and check flood level are 1082.11 m and 1082.94 m, respectively. The total storage capacity of the reservoir is 72.81 million m3, and the regulation capacity is 58 million m3. This reservoir spans across a 5 km2 area and serves as a multiyear regulating reservoir, fulfilling several purposes including water supply, downstream irrigation, and water for production, everyday life, and ecological conservation. The embankment dam in question is a homogeneous structure with a crest length of 668 m, a maximum height of 33.8 m, a crest elevation of 1085.80 m, and a width of 8 m at its crest. At its base, the maximum width of the dam is 205 m. The upper and lower slopes of the dam are protected by concrete and turf, respectively. Both the upstream and downstream slopes of the dam are compound sections. The upstream slope is defined by an elevation of 1070.8 m, with a slope ratio of 1:3.0 above and 1:3.5 below, while the downstream slope is also defined by an elevation of 1070.8 m, with a slope ratio of 1:2.5 above and 1:2.75 below. At both the upstream and downstream change slopes, there is a 2 m wide berm, and prism drainage is installed at the foot of the dam. The dam foundation is designed to prevent seepage by using a 0.6 m thick concrete cutoff wall, with the lower part of the wall penetrating 5 m into the relatively impermeable layer. A typical cross-section of the dam is illustrated in Figure 3.

3.2. Three-Dimensional Model

Based on the engineering geological conditions and design data of the reservoir dam, a 3D model of the dam was constructed using COMSOL Multiphysics (Figure 4). To enhance the accuracy of the model calculation, the entire model was meshed and refined, resulting in a total of 192,806 domain cells (Figure 5).

3.3. Boundary Conditions and Material Parameters

3.3.1. Boundary Conditions

The boundary conditions of the coupled model were determined based on the actual conditions of the project [45]. For the seepage field calculation, the constant head boundary was set to the water level below the upstream and downstream areas:
p = ρ g ( H 0 z )
where H0 is the given head, m, and z is the elevation, m.
In the model, the seepage-free surface is incorporated using the permeable layer boundary provided in COMSOL Multiphysics [46]. This boundary is essentially a hybrid boundary and is utilized when conducting unsaturated seepage calculations. It can be expressed as
n ρ u = ρ R b ( H b H )
where Hb is the external pressure head, m; H is the internal boundary head, m; n is the normal vector outside the boundary; u is the flow velocity, m/s; Rb is the external resistance, 1/s, Rb = ks/l; ks is the saturation hydraulic conductivity, m/s; and l is the coupling length scale.
As the seepage overflow point is dependent on the upstream and downstream head boundaries as well as the permeability of the dam material, it remains unknown. This study aims to identify the seepage overflow point on the downstream slope of the dam using the calculated pore pressure p. When p ≥ 0, the corresponding boundary is transformed into the position head boundary H = z (where z represents the position head), and for p < 0, it is transformed into the flow boundary q = 0.
The initial condition for the seepage field is set to zero. As for the remaining boundaries, they are set as impermeable boundaries.
n ρ u = 0
When calculating the stress field, the boundary below the upstream water level is defined as the pressure boundary. The bottom boundary of the model is considered a fixed constraint, while the surrounding truncated boundary is set as a roller support, and all other boundaries are free. The initial displacement field and velocity field are both set to zero.

3.3.2. Material Parameters

The calculation parameters of the coupled model are determined by the engineering design data, hydrogeological data, and field tests. The details are shown in Table 1.

4. Case Study Results

To perform transient seepage analysis, it is necessary to specify the initial conditions (initial pore pressure field), and steady-state seepage analysis results are often used as the initial conditions for transient seepage analysis. In this study, the actual situation of reservoir storage was considered. During the impoundment phase, the reservoir level rises, and the steady seepage field under the dead water level of 1066.05 m was used as the initial condition. The storage rate of 1 m/d rises to the restricted storage level of 1074.05 m and stays at this level for 8 days. During the drainage phase, the reservoir level decreases, and the stable seepage field under the restricted storage level of 1074.05 m was used as the initial condition. The reservoir level then decreases to the dead water level of 1066.05 m with a storage rate of 1 m/d and stays at this level for 8 days. The total time required for storage and discharge is 16 days. Rainfall infiltration and other dynamic water conditions during this period were not considered.

4.1. Seepage Characterization

For ease of analyzing the seepage characteristics of the reservoir dam, the elevation range of the water level slope between 1066.05 m and 1074.05 m was horizontally shifted 0.1 m downstream to serve as the monitoring cross-section, as illustrated in Figure 6.
Figure 7 illustrates the changes in saturation observed in the monitoring section as the water level gradually rises to a constant level during the impoundment period. As shown, the saturation value of the monitoring section increases over time, with the magnitude and response time of the increase being directly proportional to the elevation. Specifically, higher elevations of the monitoring section exhibit a longer response time and a greater increase in saturation. This phenomenon can be attributed to the fact that points located at higher elevations above the initial reservoir water level experience a longer response time during the impoundment period. Moreover, the magnitude of saturation rise is found to be inversely proportional to the initial saturation, with smaller initial saturation values resulting in a greater rise in saturation.
Figure 8 shows the saturation changes of the monitoring section during the drainage period as the water level drops to a constant level. The saturation value of the monitoring section decreases over time during the drainage period. The elevation of the monitoring section is inversely proportional to the response time and directly proportional to the decrease in saturation. This means that the higher the elevation of the monitoring section, the shorter the response time, and the greater the decrease in saturation. This phenomenon can be explained by the fact that during the drainage period, the initial reservoir level is lower for points located at higher elevations, leading to an earlier response and shorter response time. Moreover, a higher initial saturation leads to a greater magnitude of the saturation drop.
Due to the dynamic distribution of head contours and streamlines within the dam during fluctuations in water levels, only four significant moments (0 d, 4 d, 8 d, and 16 d) have been chosen to present the dynamic changes in the seepage field of a representative section of the dam. These changes are illustrated in Figure 9 and Figure 10, with the limited space available necessitating the presentation of only these four characteristic moments.
Figure 9 presents the seepage field of a representative section of the dam at various significant moments during the impoundment period. As the water level rises, the groundwater level near the waterfront rim of the dam slope increases due to reservoir water recharge, but the rate of rise decreases with elevation. This phenomenon can be attributed to the relationship between the unsaturated soil’s permeability coefficient and saturation. Within the dam, the higher the elevation, the smaller the saturation and permeability coefficient, resulting in less significant groundwater level uplift. Due to the soil material and the permeability of the concrete cutoff wall, the groundwater level in the aquifer within the dam did not rapidly respond to the rising water level in front of the slope, and the groundwater level lifting effect was weaker as it approached downstream, lagging behind. Generally, the infiltration line displays an upward concave trend. As the reservoir level continues to rise, the hysteresis effect becomes more pronounced, leading to a more significant upward concave trend. When the reservoir reaches a constant water level retention period, the water level in front of the slope stabilizes, and the influence of gravity causes the groundwater level of the aquifer in the dam to decline. As a result, the infiltration line gradually returns to a flat state.
Figure 10 illustrates the seepage field of a typical dam section at various significant moments during the drainage period. As the water level drops during the drainage period, the groundwater level in the dam decreases, leading to a reduction in the pressure head. As a result, the saturated zone contracts while the unsaturated zone gradually expands, causing the infiltration line to change from a gentle curve to an upwardly convex curve. However, due to soil material and concrete cutoff wall permeability constraints, the groundwater level response of the aquifer inside the dam produces an apparent hysteresis effect. The downstream side exhibits a more significant hysteresis effect than the upstream side. Therefore, when the reservoir water level falls to the corresponding elevation, the groundwater on the upstream side discharges slowly, resulting in a higher peak head than that of the two sides in the upstream dam aquifer. During the constant water level retention period, the hydraulic gradient from the peak head high-point to both sides of the flow leads to the dam aquifer recharge. In the process, the peak head gradually dissipates, and the pore water is discharged, resulting in the gradual stabilization of the saturated and unsaturated areas. The infiltration line is restored to a smooth curve, and the pore water seepage from within the slope is redirected from outside the slope to inside the slope.
Based on the aforementioned analysis, it is evident that the infiltration line at the riverside closely tracks the fluctuations of the reservoir water level whether it is increasing or decreasing. The changes in the internal seepage field, or infiltration line, are mainly concentrated in the upstream section of the dam body, while the middle section of the dam body exhibits a noticeable delay in the changes of its infiltration line. Furthermore, the downstream part of the concrete cutoff wall shows a more significant delay in its infiltration line changes.

4.2. Dam Slope Stability Analysis

In order to study the influence of changes in the seepage field of the dam on the stability of the dam slope during the rise and fall of the reservoir water level, this section analyzes the stability of the dam slope at four characteristic moments under both reservoir storage and discharge conditions by means of a multifield coupled model with an embedded finite element strength reduction method. The equivalent plastic strain is typically used as the criterion for slope instability in finite element strength reduction analysis. If the equivalent plastic strain pervades from the slope’s foot to the top, the slope is deemed destabilized. Additionally, since geotechnical body deformation is more pronounced in areas with significant equivalent plastic strain, this criterion can effectively identify the potential sliding surface of the geotechnical body. Figure 11 and Figure 12 illustrate the equivalent plastic strain and the overall stability safety factor of the dam slope for each characteristic moment of the embankment dams during the impoundment periods, respectively, while Figure 13 and Figure 14 illustrate the equivalent plastic strain and the overall stability factor of safety of the dam slope for each characteristic moment of the embankment dams during the drainage periods.
As depicted in Figure 11 and Figure 12, prior to the impoundment, the water level in front of the slope is at the reservoir dead level, and the initial seepage field is the stable seepage field. With the reduction of the dam material strength, a plastic failure surface emerges in the upstream slope, and the computed overall stability safety factor of the slope is 1.434, indicating that this embankment dam is more susceptible to unstable failure than the downstream slope when operating under the dead water level. After 4 days of water impoundment, the water level in front of the slope rises to 1070.05 m, and the overall stability factor of safety of the slope is 1.548. The sliding surface of the dam slope instability appears downstream. The water level in front of the slope further increases to 1074.05 m after 8 days of water impoundment, and the overall stability factor of safety of the slope is 1.55. The sliding surface of the dam slope instability also occurs downstream. Following 8 days of constant water level, the overall stability factor of safety of the slope is 1.543, and the sliding surface of the dam slope instability emerges downstream.
According to the principles of hydrostatics and coupling theory, the stability of a dam slope under impoundment conditions depends primarily on two factors. Firstly, the increase in groundwater level leads to an increase in pore water pressure and a decrease in effective stress of the soil body in the upstream dam, which adversely affects the stability of the upstream slope. Secondly, the elevation of the water level in front of the slope increases the infiltration pressure of reservoir water on the front edge of the slope body, thereby enhancing the anti-slip force on the landslide body, which has a positive impact on maintaining the stability of the dam slope. In this study, the dam fill material has a high compaction and a low permeability coefficient. As the reservoir level rises, the groundwater level within the dam body rises more slowly than the reservoir water outside the slope. Therefore, the infiltration pressure of the reservoir water towards the slope is high, while the reduction in effective stress of the soil body is minimal. At this time, the overall stability safety factor of the dam slope increases, making the upstream slope less prone to sliding than the downstream slope of the dam. With the continuous reduction of the strength parameter of the dam material, a plastic failure surface appears in the downstream slope. The seepage analysis shows that due to the limited permeability of the soil material and the concrete cutoff wall, the reservoir impoundment has a minimal impact on the seepage field of the downstream dam, and the overall stability factor of safety of the dam slope remains essentially unchanged.
Based on Figure 13 and Figure 14, it is evident that the water level in front of the slope is restricted before water release, and the initial seepage field is stable. As the strength parameters of the dam material continue to decrease, a plastic failure surface appears in the downstream slope. At this point, the overall stability safety factor of the dam slope is calculated to be 1.528, indicating that the embankment dam is more prone to unstable failure than the upstream slope when operating under restricted water level conditions. After 4 days of water discharge, the water level in front of the slope rose to 1070.05 m, and the sliding surface of the dam slope instability appeared upstream, resulting in a coefficient of safety of 1.508 for the dam slope stability. After 8 days of water discharge, the water level in front of the slope further decreased to 1066.05 m, causing a significant reduction in stability. The overall coefficient of safety of the dam slope stability was 1.36, and the sliding surface of the dam slope instability appeared upstream. Following 8 days of constant water level, the water level in front of the slope remained at 1066.05 m, and the sliding surface of the dam slope instability appeared upstream. After an additional constant water level period of 8 days, the stability rebounded, resulting in an overall stability safety factor of 1.389 for the slope, and the sliding surface when the slope was destabilized appeared upstream.
During the initial stage of water release (0–4 days), the high reservoir water level leads to greater water thrust in front of the slope, which enhances the stability safety of the upstream dam slope. As the dam parameters are continuously reduced, a plastic failure surface appears in the downstream slope of the dam. However, the drop in water level does not significantly affect the downstream seepage field, and the overall stability safety factor of the slope remains steady. In the later stage of water release (4–8 days), the water pressure generated by the reservoir water on the upstream slope decreases significantly. The poor permeability of the soil itself hinders timely water discharge, resulting in seepage from the inside of the slope to the outside of the slope under the hydraulic gradient. This creates a seepage pressure that points from the inside of the slope to the outside, adversely affecting the stability of the upstream dam slope. As a result, the overall stability safety factor of the dam slope decreases significantly. During the water level retention period (8–16 days), the water in the aquifer of the upstream dam slope is gradually discharged. This causes the unsaturated zone to expand and the pore water pressure to decrease, increasing the effective stress of the soil and the stability of the upstream dam slope. Consequently, the overall stability safety factor of the dam slope shows a small rebound and eventually stabilizes.

4.3. Sensitivity Analysis

4.3.1. Factor Significance Screening

Sensitivity analysis methods include local sensitivity analysis and global sensitivity analysis. The global sensitivity analysis method used in this paper can better reflect the interaction between factors. However, due to the large number of factors, it is not possible to take into account the interaction of all factors, so it is necessary to first screen out a few parameters that are more sensitive to the test index and then conduct a global sensitivity study on the screened parameters. The coupled model has a total of 13 influencing factors: magnitude of water level change (H), rate of water level change (v), infiltration coefficient (ks), cohesion (c), angle of internal friction (φ), modulus of elasticity (E), Poisson’s ratio (ν), soil density (ρs), saturated water content (θs), residual water content (θr), and van Genuchten model parameters (α, m, and n). Among them, the five parameters θs, θr, α, m, and n have minimal influence on the coefficient of safety of dam slope stability and can be neglected [47], so they are not involved in the discussion in this study. The eight factors H, v, ks, c, φ, E, ν, and ρs were selected to fully reflect the difference in sensitivity of each factor to the model output results. Each factor was examined at two levels, with the high level (+1) increasing by 10% according to the basic working condition and the low level (−1) decreasing by 10% according to the basic working condition, resulting in a 20% difference between high and low levels, as shown in Table 2. The test was arranged according to the Plackett–Burman experimental design, and the factors with significant sensitivity to the coefficient of safety of dam slope stability were screened out using ANOVA in order to further investigate their main effects as well as interaction effects; the specific test protocol and test results are shown in Table 3.
As can be seen from Table 4, when the reservoir is stored with water, the sensitivity of each influencing factor is ranked as follows in descending order: φ > ρs > c > E > v > ks > H > ν. Notably, φ, c, and ρs have a highly significant influence on the overall stability and safety coefficient of the dam slope, while E, v, ks, H, and ν have an insignificant effect. Conversely, during the reservoir discharge, the sensitivity of each influencing factor is ranked as follows: φ > c > H > ρs > ks > v > ν > E. Among these, φ, c, H, and ρs have a highly significant impact on the overall stability and safety coefficient of the dam slope. However, the effects of water levels ks, v, ν, and E on the overall stability safety coefficient of the dam slope were not significant. It can be seen that the sensitivity of the soil mechanical parameters (φ, c, and ρs) to the overall stability safety factor of the dam slope is significant for both water storage and water release. The sensitivity of the soil saturation permeability coefficient ks to the overall stability safety factor of the dam slope is not significant, which may be due to the high degree of compaction of the earth and rock dam fill in this study and the small permeability coefficient of the soil. Furthermore, the sensitivity of the seepage field influences (H, v, and ks) on the overall stability factor of the dam slope is less pronounced during the reservoir storage phase in comparison to the discharge period. This disparity in sensitivity can be attributed to the varying effects of reservoir storage and release on the overall stability of the dam slope. During the storage period, it is the downstream slope of the dam that is more likely to experience landslides. Through seepage and stability analysis, it becomes evident that reservoir water level fluctuations primarily affect the stability of the dam slope by modifying the waterfront and the internal seepage field of the upstream slope. Conversely, these fluctuations have minimal impact on the seepage field of the downstream slope. Therefore, the influence of seepage field factors (H, v, and ks) on the stability of the downstream slope should be minor and not significantly sensitive.

4.3.2. Analysis of Main Effects and Interaction Effects

A Box–Behnken test design was carried out to explore the primary effects of each significant factor and the interaction effects between them for the factors identified by the Plackett–Burman test to have significant effects on the overall stability and safety coefficient of the dam slope. For each factor, three levels were chosen, with the high level (+1) increasing by 10% from the base operating conditions, the low level (−1) decreasing by 10% from the base operating conditions, and the base level (0) being maintained according to the base operating conditions. The test factors and their corresponding levels for the Box–Behnken design are presented in Table 5. The specific test procedures and results can be found in Table 6 and Table 7.
The effects of influencing factors on the overall stability and safety coefficients of dam slopes during the impoundment and drainage period are presented in Figure 15. The main effect of each factor on the overall stability safety coefficient of the dam slope during impoundment is ranked as φ > ρs > c, while during drainage, the ranking is φ > c > H > ρs, consistent with the Plackett–Burman test results According to Figure 15a, as values of φ and c increase, the overall stability safety factor of the dam slope during the reservoir period gradually increases. In contrast, as the value of ρs increases, the overall stability safety factor of the dam slope during the reservoir period gradually decreases. Additionally, based on Figure 15b, the factors φ and c cause significant differences in the overall stability safety factor of the dam slope during the drainage period. As the values of φ and c increase, the overall stability safety factor of the dam slope during drainage significantly increases, while the factors of H and ρs leads to a gradual decrease.
Based on the Box–Behnken experiment results, the sensitivity of significant factors and their second-order interactions on the overall stability safety factor of the dam slope were analyzed using variance analysis. Table 8 shows the results of the variance analysis. Third-order and fourth-order interactions were not analyzed in this study as they were relatively small. From Table 8, it can be observed that for the overall stability safety factor of the dam slope during impoundment, the sensitivity ranking of the three factors φ, c, and ρs and their interaction effects were φ > ρs > c > c × ρs > φ × ρs > c × φ. Among these, the effects of φ, c, and ρs were extremely significant, while the effect of c × ρs was significant. For the overall stability safety factor of the dam slope during drainage, the sensitivity ranking of the four factors H, c, φ, and ρs and their interaction effects were φ > c > H > ρs > c × ρs > H × ρs > φ × ρs > H × φ > H × c > c × φ. Among them, the effects of φ, c, H, and ρs were extremely significant, while the effects of c × ρs, H × ρs and φ × ρs were significant. The sensitivity ranking of the factors obtained through the Box–Behnken experiment is consistent with the Plackett–Burman experiment discussed earlier. Furthermore, factors c, φ, and ρs and their second-order interaction, c × ρs, have significant effects on the overall stability of the dam slope during both water storage and drainage periods and should receive priority attention in dam slope stability calculations.
To further investigate the influence of interaction terms on the stability safety factor of the dam slope, noninteraction factors were fixed at their optimal levels, and four sets of contour maps depicting the overall stability safety factor of the dam slope under the four interaction effects were generated (Figure 16). As shown in Figure 16, when φ is held constant at 15°, the overall stability safety factor of the dam slope during impoundment gradually declines with decreasing c and increasing ρs. When H is kept at 8 m and φ is maintained at 15°, the overall stability safety factor of the dam slope during drainage gradually declines with decreasing c and increasing ρs. Additionally, when c is maintained at 33 kPa and φ is kept at 15°, the overall stability safety factor of the dam slope during drainage gradually declines with increasing H and ρs. Lastly, when H is fixed at 8 m and c is held constant at 33 kPa, the overall stability safety factor of the dam slope during drainage gradually decreases with decreasing φ and increasing ρs.

5. Conclusions and Future Work

In order to address the need for sensitivity analysis on the stability of dam slopes considering seepage–stress coupling during reservoir water level fluctuations, a novel global sensitivity analysis method is proposed. To demonstrate its applicability, a case study was conducted on an earth and rock dam project located in Shaanxi Province, China. A three-dimensional finite element model was developed, and the seepage and slope stability of the dam under the influence of seepage–stress coupling and reservoir water level fluctuations were analyzed using COMSOL Multiphysics software. Furthermore, the global sensitivity of the earth and rock dam stability under reservoir water level fluctuations were explored through a combination of Plackett–Burman and Box–Behnken experimental designs. The main findings and conclusions drawn from this analysis are summarized as follows:
1. During the rise and fall of the reservoir level, the infiltration line on the waterfront of the dam slope is synchronized with the change of the reservoir level and the change of the internal seepage field of the dam body is mainly concentrated upstream, clearly lagging behind the slope surface, while there is no obvious change of the seepage field downstream of the dam body. When the reservoir is filled with water, the infiltration line of the dam body is depressed; the groundwater level rises; the saturation of the soil body rises; and the rise rate and response time are proportional to the elevation. The overall stability safety factor of the dam slope rises first and then tends to stabilize. The dam infiltration line is convex when the reservoir is discharged. The groundwater level decreases; the soil saturation decreases; the decline is proportional to the elevation; and the response time is inversely proportional to the elevation. The overall stability factor of safety of the dam slope decreases and then increases, with the minimum value occurring at the moment of the end of the release, when the upstream slope is more prone to sliding compared to the downstream slope of the dam. Therefore, monitoring of upstream dam slope landslides should be strengthened during reservoir releases, and safety warnings should be conducted when necessary.
2. Based on the Plackett–Burman test, the factors that have significant sensitivity to the overall stability coefficient of safety of the dam slope during the storage and release periods were determined separately. During reservoir storage, φ, c, and ρs have significant effects on the overall slope stability factor, while H, v, ks, ν, and E have insignificant effects on the overall slope stability factor. In the case of reservoir discharge, the effects of H, φ, c, and ρs on the overall stability coefficient of safety of the dam slope are very significant, while the effects of v, ks, ν, and E on the overall stability coefficient of safety of the dam slope are not significant. Therefore, in practical application, the observation accuracy of H, φ, c, and ρs should be improved and determined by test and inversion, if necessary, to ensure the accuracy of the simulation results, while factors with insignificant sensitivity can be determined by referring to the relevant literature or by engineering analogy.
3. Based on the Box–Behnken test, the main effects and interaction effects of the selected significant factors were analyzed, and it was concluded that the main effect of φ was the largest for the overall stability and safety factor of the dam slope during the storage period, and the interaction effect of c × ρs was significant. The main effects of φ and c are larger for the overall stability factor of the dam slope during the discharge period; the interaction effects of c × ρs, H × ρs, and φ × ρs are significant. The second-order interaction term c × ρs has a significant effect on the overall stability of the dam slope during both the storage and discharge periods, and it cannot be ignored in the calculation of dam slope stability.
In this paper, a great amount of work has been carried out on the sensitivity analysis of seepage and slope stability of earth and rock dams under fluctuating reservoir water level conditions, but there are still research limitations, and many aspects need further study. For example, the combined effect of rainfall and reservoir level rise and fall is not considered. Under some conditions, especially in rainy seasons, reservoir scheduling often faces extreme rainfall conditions. On the one hand, these rains will raise the reservoir level and affect the change of seepage field at the waterfront of the dam slope. On the other hand, rainfall infiltration will lead to a change of seepage field at the surface of the dam body, and the combined effect of these two aspects will further change the stress conditions inside the dam body and increase the risk of dam slope instability. The risk of dam slope instability increases. In addition, the effects of different combinations of reservoir level rise and fall rate, reservoir level rise and fall magnitude, rainfall intensity, rainfall duration, and rainfall type on seepage and slope stability of earth and rock dams need to be further explored.

Author Contributions

C.Z.: Conceptualization, Investigation, Formal analysis, Visualization, Data curation, Writing—original draft, and Writing—review and editing. Z.S.: Supervision, Formal analysis, and Writing—review and editing. L.X.: Formal analysis, Methodology, Funding acquisition, Supervision, and Writing—review and editing. Y.S.: Validation, Supervision, Funding acquisition, and Writing—review and editing. W.Z.: Data curation, Supervision, and Writing—review and editing. H.Z.: Funding acquisition, Supervision, and Writing—review and editing. J.P.: Visualization, Investigation, and Writing—review and editing. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation of China (grant No. 52179130), the Fundamental Research Funds for the Center Universities (grant No. 2019B70814), and the Project of Natural Science Foundation of Jiangsu Province (grant No. BK20201312).

Data Availability Statement

Not applicable.

Acknowledgments

We sincerely acknowledge all reviewers. Their helpful comments have undoubtedly contributed to improving the article.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Fu, X.; Gu, C.S.; Su, H.Z.; Qin, X.N. Risk analysis of earth-rock dam failures based on fuzzy event tree method. Int. J. Environ. Res. Public Health 2018, 15, 886. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  2. Jiang, X.; Zhu, Z.; Chen, H.; Deng, M.; Niu, Z.; Deng, H.; Zou, Z. Natural dam failure in slope failure mode triggered by seepage. Geomat. Nat. Hazards Risk 2020, 11, 698–723. [Google Scholar] [CrossRef] [Green Version]
  3. Jadid, R.; Montoya, B.M.; Bennett, V.; Gabr, M.A. Effect of repeated rise and fall of water level on seepage-induced deformation and related stability analysis of Princeville levee. Eng. Geol. 2020, 266, 105458. [Google Scholar] [CrossRef]
  4. Wang, L.; Wu, C.; Tang, L.; Zhang, W.; Lacasse, S.; Liu, H.; Gao, L. Efficient reliability analysis of earth dam slope stability using extreme gradient boosting method. Acta Geotech. 2020, 15, 3135–3150. [Google Scholar] [CrossRef]
  5. Cen, W.J.; Li, D.J.; Wang, H. Impact of transient seepage on slope stability of earth-rock dams with geomembrane barrier defects. Environ. Geotech. 2020, 7, 581–590. [Google Scholar] [CrossRef]
  6. Guo, X.; Dias, D.; Carvajal, C.; Peyras, L.; Breul, P. Three-dimensional probabilistic stability analysis of an earth dam using an active learning metamodeling approach. Bull. Eng. Geol. Environ. 2022, 81, 40. [Google Scholar] [CrossRef]
  7. Wu, Z.; Chen, C.; Lu, X.; Pei, L.; Zhang, L. Discussion on the allowable safety factor of slope stability for high rockfill dams in China. Eng. Geol. 2020, 272, 105666. [Google Scholar] [CrossRef]
  8. Guo, X.; Dias, D.; Carvajal, C.; Peyras, L. Reliability analysis of embankment dam sliding stability using the sparse polynomial chaos expansion. Eng. Struct. 2018, 174, 295–307. [Google Scholar] [CrossRef]
  9. Liang, R.Y.; Nusier, B.O.; Malkawi, A.H. A reliability based approach for evaluating the slope stability of embankment dams. Eng. Geol. 1999, 54, 271–285. [Google Scholar] [CrossRef]
  10. Elia, G.; Falcone, G.; Cotecchia, F.; Rouainia, M. Analysis of the Effects of Seasonal Pore Pressure Variations on the Slope Stability Through Advanced Numerical Modelling. In Geotechnical Research for Land Protection and Development: Proceedings of CNRIG 2019 7; Springer International Publishing: Berlin/Heidelberg, Germany, 2020; pp. 184–194. [Google Scholar] [CrossRef]
  11. Zhang, W.; Shi, D.; Shen, Z.; Wang, X.; Gan, L.; Shao, W.; Tang, P.; Zhang, H.; Yu, S. Effect of Calcium Leaching on the Fracture Properties of Concrete. Constr. Build. Mater. 2023, 365, 130018. [Google Scholar] [CrossRef]
  12. Zhang, W.; Shi, D.; Shen, Z.; Shao, W.; Gan, L.; Yuan, Y.; Tang, P.; Zhao, S.; Chen, Y. Reduction of the Calcium Leaching Effect on the Physical and Mechanical Properties of Concrete by Adding Chopped Basalt Fibers. Constr. Build. Mater. 2023, 365, 130080. [Google Scholar] [CrossRef]
  13. Kala, Z. New importance measures based on failure probability in global sensitivity analysis of reliability. Mathematics 2021, 9, 2425. [Google Scholar] [CrossRef]
  14. Bai, Z.; Wei, H.; Xiao, Y.; Song, S.; Kucherenko, S. A Vine Copula-Based Global Sensitivity Analysis Method for Structures with Multidimensional Dependent Variables. Mathematics 2021, 9, 2489. [Google Scholar] [CrossRef]
  15. Peng, X.; Xu, X.; Li, J.; Jiang, S. A Sampling-based sensitivity analysis method considering the uncertainties of input variables and their distribution parameters. Mathematics 2021, 9, 1095. [Google Scholar] [CrossRef]
  16. Song, K.; Chang, I.; Pham, H. A testing coverage model based on NHPP software reliability considering the software operating environment and the sensitivity analysis. Mathematics 2019, 7, 450. [Google Scholar] [CrossRef] [Green Version]
  17. Liu, D.; Lin, T.; Gao, J.; Xue, B.; Yang, J.; Chen, C.; Zhang, W.; Sun, W. Study on the Influence of Water Level on Earth Dam Reinforced by Cut-Off Wall: A Case Study in Wujing Reservoir. Water 2023, 15, 140. [Google Scholar] [CrossRef]
  18. Wang, L.; Wu, C.; Gu, X.; Liu, H.; Zhang, W. Probabilistic stability analysis of earth dam slope under transient seepage using multivariate adaptive regression splines. Bull. Eng. Geol. Environ. 2020, 79, 2763–2775. [Google Scholar] [CrossRef]
  19. Xie, Q.; Liu, J.; Han, B.; Li, H.; Li, X. Experimental and numerical investigation of bottom outlet leakage in earth-fill dams. Perform. Constr. Facil. 2019, 33, 04019037. [Google Scholar] [CrossRef]
  20. Wang, R.; Wan, J.; Cheng, R.; Wang, Y.; Wang, Z. Physical and Numerical Simulation of the Mechanism Underpinning Accumulation Layer Deformation, Instability, and Movement Caused by Changing Reservoir Water Levels. Water 2023, 15, 1289. [Google Scholar] [CrossRef]
  21. Jin, J.; Song, C.; Chen, Y. Investigation of a fluid–solid coupling model for a tailings dam with infiltration of suspended particles. Environ. Earth Sci. 2017, 76, 758. [Google Scholar] [CrossRef]
  22. Su, Z.; Chen, G.; Meng, Y. Study on seepage characteristics and stability of core dam under the combined action of the variation of reservoir water level and rainfall. Geotech. Geol. Eng. 2021, 39, 193–211. [Google Scholar] [CrossRef]
  23. Meng, Q.; Qian, K.; Zhong, L.; Gu, J.; Li, Y.; Fan, K.; Yan, L. Numerical Analysis of Slope Stability under Reservoir Water Level Fluctuations Using a FEM-LEM-Combined Method. Geofluids 2020, 2020, 1–12. [Google Scholar] [CrossRef]
  24. Jie, Y.X.; Wen, Y.F.; Deng, G.; Chen, R.; Xu, Z.P. Impact of soil deformation on phreatic line in earth-fill dams. Comput. Geosci. 2012, 46, 44–50. [Google Scholar] [CrossRef]
  25. Hu, L.; Takahashi, A.; Kasama, K. Effect of spatial variability on stability and failure mechanisms of 3D slope using random limit equilibrium method. Soils Found. 2022, 62, 101225. [Google Scholar] [CrossRef]
  26. Ahmad, F.; Tang, X.W.; Qiu, J.N.; Wróblewski, P.; Ahmad, M.; Jamil, I. Prediction of slope stability using Tree Augmented Naive-Bayes classifier: Modeling and performance evaluation. Math. Biosci. Eng. 2022, 19, 4526–4546. [Google Scholar] [CrossRef]
  27. Cross, M. Sensitivity analysis of shallow planar landslides in residual soils on south Pennine hillslopes, Derbyshire, UK. Bull. Eng. Geol. Environ. 2019, 78, 1855–1872. [Google Scholar] [CrossRef]
  28. Ghadrdan, M.; Shaghaghi, T.; Tolooiyan, A. Sensitivity of the stability assessment of a deep excavation to the material characterisations and analysis methods. Geomech. Geophys. Geo-Energy Geo-Resour. 2020, 6, 59. [Google Scholar] [CrossRef]
  29. Karthik, A.V.R.; Manideep, R.; Chavda, J.T. Sensitivity analysis of slope stability using finite element method. Innov. Infrastruct. Solut. 2022, 7, 184. [Google Scholar] [CrossRef]
  30. Richards, L.A. Capillary conduction of liquids through porous mediums. Physics 1931, 1, 318–333. [Google Scholar] [CrossRef]
  31. Kirkland, M.R.; Hills, R.G.; Wierenga, P.J. Algorithms for solving Richards’ equation for variably saturated soils. Water Resour. Res. 1992, 28, 2049–2058. [Google Scholar] [CrossRef]
  32. Barry, D.A.; Parlange, J.Y.; Sander, G.C.; Sivaplan, M. A class of exact solutions for Richards’ equation. J. Hydrol. 1993, 142, 29–46. [Google Scholar] [CrossRef]
  33. Tocci, M.D.; Kelley, C.T.; Miller, C.T. Accurate and economical solution of the pressure-head form of Richards’ equation by the method of lines. Adv. Water Res. 1997, 20, 1–14. [Google Scholar] [CrossRef]
  34. Paniconi, C.; Troch, P.A.; vanLoon, E.E.; Hilberts, A.G.J. Hillslope-storage Boussinesq model for subsurface flow and variable source areas along complex hillslopes: 2. Intercomparison with a three-dimensional Richards equation model. Water Resour. Res. 2003, 39, 1317. [Google Scholar] [CrossRef]
  35. 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. [Google Scholar] [CrossRef] [Green Version]
  36. Gardner, W.R.; Hillel, D.; Benyamini, Y. Post-irrigation move-ment of soil water: 1. Redistribution. Water Resour. Res. 1970, 6, 851–861. [Google Scholar] [CrossRef]
  37. Brooks, R.H.; Corey, A.T. Hydraulic properties of porous media and their relation to drainage design. Trans. ASAE 1964, 7, 0026–0028. [Google Scholar] [CrossRef]
  38. Zhang, W.B.; Shen, Z.Z.; Ren, J.; Gan, L.; Wang, F.; Yu, B.; Li, C. Modeling and comparative analysis of a flow and heat cou-pling model of the riparian zone based on thermal conductivity empirical models. J. Hydrol. 2020, 582, 124539. [Google Scholar] [CrossRef]
  39. Thomas, H.R.; He, Y. Modelling the behaviour of unsaturated soil using an elastoplastic constitutive model. Géotechnique 1998, 48, 589–603. [Google Scholar] [CrossRef]
  40. Kumar, J.; Samui, P. Stability determination for layered soil slopes using the upper bound limit analysis. Geotech. Geol. Eng. 2006, 24, 1803–1819. [Google Scholar] [CrossRef]
  41. Chen, D.H.; Du, C.B. Application of strength reduction FEM to dynamic anti-sliding stability analysis of high gravity dam with complex dam foundation. Water Sci. Eng. 2011, 4, 212–224. [Google Scholar] [CrossRef]
  42. Plackett, R.L.; Burman, J.P. The Design of Optimum Multifactorial Experiments. Biometrika 1946, 33, 305–325. [Google Scholar] [CrossRef]
  43. Box, G.E.P.; Behnken, D.W. Some new three level designs for the study of quantitative variables. Technometrics 1960, 2, 455–475. [Google Scholar] [CrossRef]
  44. Zhang, W.; Shen, Z.; Ren, J.; Gan, L.; Xu, L.; Sun, Y. Phase-field simulation of crack propagation in quasi-brittle materials: COMSOL implementation and parameter sensitivity analysis. Model. Simul. Mater. Sci. Eng. 2021, 29, 055020. [Google Scholar] [CrossRef]
  45. di Lernia, A.; Cotecchia, F.; Elia, G.; Tagarelli, V.; Santaloia, F.; Palladino, G. Assessing the influence of the hydraulic boundary conditions on clay slope stability: The Fontana Monte case study. Eng. Geol. 2022, 297, 106509. [Google Scholar] [CrossRef]
  46. Chui, T.F.M.; Freyberg, D.L. The use of COMSOL for integrated hydrological modeling. In Proceedings of the COMSOL Conference, Boston, MA, USA, 2–4 October 2007; pp. 217–223. [Google Scholar]
  47. Zhang, W.; Shen, Z.; Ren, J.; Bian, J.; Xu, L.; Chen, G. Multifield Coupling Numerical Simulation of the Seepage and Stability of Embankment Dams on Deep Overburden Layers. Arab. J. Sci. Eng. 2022, 47, 7293–7308. [Google Scholar] [CrossRef]
Figure 1. Schematic layout of the 3-factor Box–Behnken design.
Figure 1. Schematic layout of the 3-factor Box–Behnken design.
Mathematics 11 02836 g001
Figure 2. Flow chart of global sensitivity analysis method considering seepage–stress coupling.
Figure 2. Flow chart of global sensitivity analysis method considering seepage–stress coupling.
Mathematics 11 02836 g002
Figure 3. Typical cross-section of the reservoir dam.
Figure 3. Typical cross-section of the reservoir dam.
Mathematics 11 02836 g003
Figure 4. 3D model of the reservoir dam.
Figure 4. 3D model of the reservoir dam.
Mathematics 11 02836 g004
Figure 5. 3D finite element mesh of the reservoir dam.
Figure 5. 3D finite element mesh of the reservoir dam.
Mathematics 11 02836 g005
Figure 6. Schematic diagram of the location of the monitoring section.
Figure 6. Schematic diagram of the location of the monitoring section.
Mathematics 11 02836 g006
Figure 7. Surface of saturation variation in the monitoring section during the impoundment period.
Figure 7. Surface of saturation variation in the monitoring section during the impoundment period.
Mathematics 11 02836 g007
Figure 8. Surface of saturation change in the monitoring section during the water drainage period.
Figure 8. Surface of saturation change in the monitoring section during the water drainage period.
Mathematics 11 02836 g008
Figure 9. Seepage field of a typical section of the dam at each characteristic moment during the impoundment period. (a) 0 d, (b) 4 d, (c) 8 d, (d) 16 d.
Figure 9. Seepage field of a typical section of the dam at each characteristic moment during the impoundment period. (a) 0 d, (b) 4 d, (c) 8 d, (d) 16 d.
Mathematics 11 02836 g009
Figure 10. Seepage field of a typical section of the dam at each characteristic moment during the drainage period. (a) 0 d, (b) 4 d, (c) 8 d, (d) 16 d.
Figure 10. Seepage field of a typical section of the dam at each characteristic moment during the drainage period. (a) 0 d, (b) 4 d, (c) 8 d, (d) 16 d.
Mathematics 11 02836 g010
Figure 11. Equivalent plastic strain before dam destabilization at each characteristic moment during the impoundment period. (a) 0 d, (b) 4 d, (c) 8 d, (d) 16 d. The arrows represent the trend of soil sliding.
Figure 11. Equivalent plastic strain before dam destabilization at each characteristic moment during the impoundment period. (a) 0 d, (b) 4 d, (c) 8 d, (d) 16 d. The arrows represent the trend of soil sliding.
Mathematics 11 02836 g011
Figure 12. Variation curve of the overall stability safety factor of the dam slope during the impoundment period.
Figure 12. Variation curve of the overall stability safety factor of the dam slope during the impoundment period.
Mathematics 11 02836 g012
Figure 13. Equivalent plastic strain before dam destabilization at each characteristic moment during the drainage period. (a) 0 d, (b) 4 d, (c) 8 d, (d) 16 d. The arrows represent the trend of soil sliding.
Figure 13. Equivalent plastic strain before dam destabilization at each characteristic moment during the drainage period. (a) 0 d, (b) 4 d, (c) 8 d, (d) 16 d. The arrows represent the trend of soil sliding.
Mathematics 11 02836 g013
Figure 14. Variation curve of the overall stability safety coefficient of the dam slope during the drainage period.
Figure 14. Variation curve of the overall stability safety coefficient of the dam slope during the drainage period.
Mathematics 11 02836 g014
Figure 15. Main effect of each significant factor on the overall stability and safety coefficient of the dam slope. (a) Coefficient of safety for the overall stability of the dam slope during the impoundment. (b) Coefficient of safety for the overall stability of the dam slope during the drainage period.
Figure 15. Main effect of each significant factor on the overall stability and safety coefficient of the dam slope. (a) Coefficient of safety for the overall stability of the dam slope during the impoundment. (b) Coefficient of safety for the overall stability of the dam slope during the drainage period.
Mathematics 11 02836 g015
Figure 16. Contour map of the overall stability safety factor of the dam slope. (a) Coefficient of safety for the overall stability of the dam slope during the impoundment. (b) Coefficient of safety for the overall stability of the dam slope during the drainage period.
Figure 16. Contour map of the overall stability safety factor of the dam slope. (a) Coefficient of safety for the overall stability of the dam slope during the impoundment. (b) Coefficient of safety for the overall stability of the dam slope during the drainage period.
Mathematics 11 02836 g016
Table 1. Calculated parameters for each part of the dam.
Table 1. Calculated parameters for each part of the dam.
ks
(m/s)
c
(kPa)
φ
(°)
νE
(MPa)
ρ
(kg/m3)
θs
(m3/m3)
θr
(m3/m3)
α
(1/m)
mn
Dam body4.7 × 10−633.015.00.35802006.120.500.051.120.542.17
Layer14.4 × 10−59.028.00.334002115.19/// /
Layer22.0 × 10−616.831.00.328002214.07/// /
Concrete cutoff wall1.0 × 10−9//0.326002344.55/// /
Drainage prism2.5 × 10−3027.30.341202064.220.480.021.080.512.03
Table 2. Plackett–Burman experiment levels of each factor.
Table 2. Plackett–Burman experiment levels of each factor.
Level of
the Factor
H (m)v (m/s)ks (m/s)c (kPa)φ (°)νE (MPa)ρs (kg/m3)
−17.20.94.23−629.713.50.315721805.508
+18.81.15.17−636.316.50.385882206.732
Table 3. Plackett–Burman experimental design scenarios and calculation results of each scenario during impoundment and drainage.
Table 3. Plackett–Burman experimental design scenarios and calculation results of each scenario during impoundment and drainage.
SchemeH
(m)
v
(m/s)
ks
(m/s)
c
(kPa)
φ
(°)
νE
(MPa)
ρs
(kg/m3)
FOS
(Impoundment)
FOS
(Drainage)
1+1−1−1−1+1+1+1−11.6651.360
2−1−1+1+1+1+1+1+11.6401.519
3−1+1−1−1−1+1+1+11.3351.234
4+1−1+1+1−1+1−1−11.5751.354
5+1+1−1+1+1+1+1−11.7891.494
6+1+1+1−1+1+1−1+11.5401.334
7−1+1+1+1−1+1+1−11.5791.424
8−1+1+1−1+1+1−1−11.6651.439
9+1−1+1−1−1+1+1+11.3351.184
10−1−1−1−1−1+1−1−11.4551.289
11+1+1−1+1−1+1−1+11.4301.290
12−1−1−1+1+1+1−1+11.6401.519
Table 4. Results of sensitivity analysis of Plackett–Burman experiment.
Table 4. Results of sensitivity analysis of Plackett–Burman experiment.
IndicatorsStatistical
Quantities
HvkscφνEρs
FOS
(impoundment)
Fi0.561.180.56608.52 **2119.90 **0.562.13917.07 **
FOS
(drainage)
Fi198.62 **0.145.60687.46 **944.96 **0.130.1193.76 **
If marked with “**”, it indicates that the factor has a highly significant effect on the test results, and if marked with “*”, it indicates that the factor has a significant effect on the test results. If the factor does not have a significant effect on the test result, the “*” mark is not used.
Table 5. Box–Behnken test levels of each factor.
Table 5. Box–Behnken test levels of each factor.
Level of the FactorH
(m)
v
(m/s)
ks
(m/s)
c
(kPa)
φ
(°)
νE
(MPa)
ρs
(kg/m3)
−17.20.94.23−629.713.50.315721805.508
08.01.04.70−633.015.00.350802006.120
+18.81.15.17−636.316.50.385882206.732
Table 6. Box–Behnken experimental design scheme and calculation results of each scheme (impoundment).
Table 6. Box–Behnken experimental design scheme and calculation results of each scheme (impoundment).
Schemec (kPa)φ (°)ρs (kg/m3)FOS (Impoundment)
1−1+101.600
2+10+11.535
3+1+101.709
4+10−11.685
50001.549
60001.549
70−1−11.515
8−10−11.560
90−1+11.385
10−10+11.439
110+1+11.590
120+1−11.729
13−1−101.390
140001.549
15+1−101.499
Table 7. Box–Behnken experimental design scheme and calculation results of each scheme (drainage).
Table 7. Box–Behnken experimental design scheme and calculation results of each scheme (drainage).
SchemeH (m)c (kPa)φ (°)ρs (kg/m3)FOS (Drainage)
10−10−11.329
200−1−11.320
30−1+101.384
400−1+11.264
500001.369
60−1−101.229
7−10−101.324
800001.369
9+1+1001.398
100+1−101.354
1100+1−11.470
120+1+101.509
13+100−11.359
14−1−1001.340
15−1+1001.464
16−100−11.434
170−10+11.289
180+10+11.404
19+10+101.410
200+10−11.464
21+1−1001.274
2200+1+11.428
23−100+11.374
24−10+101.484
25+10−101.260
26+100+11.314
2700001.369
Table 8. Results of sensitivity analysis of Box–Behnken test.
Table 8. Results of sensitivity analysis of Box–Behnken test.
Factors and Their Second-Order Interaction TermsFi
FOS (Impoundment)FOS (Drainage)
H/5841.74 **
c1941.36 **19,929.14 **
φ7122.04 **31,009.64 **
ρs2948.98 **3254.17 **
H × c/0.04
H × φ/9.23
H × ρs/23.37 *
c × φ00.02
c × ρs17.35 *42.67 *
φ × ρs1.8220.02 *
If marked with “**”, it indicates that the factor has a highly significant effect on the test results, and if marked with “*”, it indicates that the factor has a significant effect on the test results. If the factor does not have a significant effect on the test result, the “*” mark is not used.
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Zhou, C.; Shen, Z.; Xu, L.; Sun, Y.; Zhang, W.; Zhang, H.; Peng, J. Global Sensitivity Analysis Method for Embankment Dam Slope Stability Considering Seepage–Stress Coupling under Changing Reservoir Water Levels. Mathematics 2023, 11, 2836. https://doi.org/10.3390/math11132836

AMA Style

Zhou C, Shen Z, Xu L, Sun Y, Zhang W, Zhang H, Peng J. Global Sensitivity Analysis Method for Embankment Dam Slope Stability Considering Seepage–Stress Coupling under Changing Reservoir Water Levels. Mathematics. 2023; 11(13):2836. https://doi.org/10.3390/math11132836

Chicago/Turabian Style

Zhou, Congcong, Zhenzhong Shen, Liqun Xu, Yiqing Sun, Wenbing Zhang, Hongwei Zhang, and Jiayi Peng. 2023. "Global Sensitivity Analysis Method for Embankment Dam Slope Stability Considering Seepage–Stress Coupling under Changing Reservoir Water Levels" Mathematics 11, no. 13: 2836. https://doi.org/10.3390/math11132836

APA Style

Zhou, C., Shen, Z., Xu, L., Sun, Y., Zhang, W., Zhang, H., & Peng, J. (2023). Global Sensitivity Analysis Method for Embankment Dam Slope Stability Considering Seepage–Stress Coupling under Changing Reservoir Water Levels. Mathematics, 11(13), 2836. https://doi.org/10.3390/math11132836

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop