Assessing the stability of the reservoir rim in moraine deposits for a mega RCC dam

1.   Introduction

Reservoir rim stability is an important but often neglected aspect in the planning and design of dams ^[1]. The Vajont Dam in Italy experienced a well-known example of a mega dam disaster due to landsliding in the reservoir rim. Soon after the dam's construction, as the first filling of the reservoir was completed, a mass of approximately 270 million cubic meters collapsed, generating huge dynamic waves that overtopped the dam and hit the downstream villages ^[2]. The impacts of this event were huge in scale, resulting in around 2, 056 casualties and incalculable economic loss.

Numerical approaches utilizing recent software developments have substantially facilitated the investigation of geohazards related to water interaction in dams and groundwater ^{[3,4,5,6,7,8]}. Despite cataclysmic disasters, considerable research has yet to be done to understand the behavior of reservoir slopes ^[9], especially on the first impounding and alternate wetting and drying cycles. Santosh and Khazanchi (2017) ^[10] performed a reservoir rim stability assessment for the Punatsangchhu Dam in Bhutan. Vulnerable reaches were identified, and a stability assessment was carried out. Subsequently, geoengineering measures for the stabilization of slopes were suggested. MWH Consultants (2013) ^[11] did a stability assessment for the reservoir of the proposed Susitna-Watana Dam in Alaska. Areas of expected landslides were evaluated through the limit equilibrium slope stability analysis. Freeze-thaw and earthquakes are significant factors for the stability and failure of moraine soil slopes in high-altitude, cold, and strong earthquake areas ^[12]. Klimeš et al. (2016) ^[13] studied landslides in moraines as triggers of glacial lake outburst floods, adopting a case study of Palcacocha Lake in Cordillera Blanca, Peru. Ghimire and Gajurel (2020) ^[14] studied moraine dam stability analysis of the Ngozumpa glacier in the Gokyo area in Nepal using the Morgenstern-Price method. Anbalagan and Kumar (2015) ^[15] studied reservoir-induced landslides, adopting a case study of the reservoir rim region of the 260-meter-high Tehri dam in India. Xiaoping and Jingwu (2011) ^[16] performed a stability analysis of bank slope under conditions of reservoir impounding and rapid drawdown using centrifugal model tests, soil laboratory tests, and numerical analysis. Yin et al. (2016) ^[17] published their work on reservoir-induced landslides and risk control at the Three Gorges Dam on the Yangtze River in China.

The most important soil instability is probably the inherent spatial variability of soil properties and its influence on slope safety factors ^[18]. Although the probabilistic stability analysis approach has been widely applied to the safety assessment of geotechnical structures, few studies have been performed to investigate the effects of water level fluctuations on earth dam slope stability, considering uncertainties of soil parameters ^[9]. Researchers performed a probabilistic stability analysis of the Ashigong earth dam under transient seepage, considering the effects of soil parameters uncertainty and water level fluctuation velocity on the earth dam slope failure. Wang et al., in 2020 ^[9], developed an efficient extreme gradient boosting (XGBoost)-based reliability analysis approach for evaluating the earth dam slope failure probability, accounting for the spatial variability of soil properties. Kumar et al. (2023) ^[19] performed a probabilistic slope stability analysis of Mount St. Helens in the USA using Scoops3D and a hybrid intelligence paradigm. In this study, the reliability index and the probability of failure were estimated under seismic and non-seismic situations.

Diamer Basha Dam is a proposed RCC dam to be constructed in the upper reaches of the Indus Valley, which has been glaciated. The dam will be located on the border between the Northern Areas and Khyber Pakhtunkhwa (KP), about 166 km downstream of Gilgit, 40 km downstream of the district headquarters of Chilas, and about 316 km upstream of Tarbela Dam. The location map of the proposed dam is shown in Figure 1. The world's highest RCC gravity dam will create a reservoir about 105 km long, stretching from the dam axis to the Raikot Bridge. The reservoir will be able to store 8.1 million acre-feet (MAF) of water. The Indus River Valley, where the reservoir would be formed, has been glaciated in geological times. It contains thick deposits of glacial origin, mostly in the upper reaches of the reservoir. A number of these deposits stand with high vertical slopes. The thickness of these glacial deposits is quite variable and exceeds hundreds of meters at several locations. These deposits are somewhat stable in their current state. However, there is a serious concern that saturation during reservoir filling, alternate wetting, drawdown during reservoir operation, or a seismic event could cause large mass movement of these slopes into the reservoir. These circumstances could lead to disastrous consequences, including blocking the reservoir, generating dynamic waves that could hit the dam and shoreline facilities, a reduction in the reservoir storage capacity, and a dramatic increase in sediment load.

Figure 1.  Location map of the Diamer Basha Dam (Source: Feasibility Report, Volume I, Diamer Basha Dam Project).

DownLoad: Full-Size Img PowerPoint

Since such stability analysis of the reservoir rim has yet to be performed for a huge national project, this research aimed to identify the potentially unstable deposits around the reservoir rim and perform a stability assessment of these deposits using deterministic, probabilistic, and dynamic analysis approaches. Landslide susceptibility maps were then generated for the reservoir based on the results of the stability analysis and topography of the reservoir area, followed by an assessment of the potential impacts of these landslides.

1.1.   Site description and geology

The general form of the valley in the reservoir area is roughly U-shaped, which is attributed to erosion by major trunk glaciers that have occupied the valley at various times. The valley profile is characterized by steep sides/sub-vertical, deep, and narrow gorges at several locations in the reservoir area. The landform within the valley indicates that it has been eroded and backfilled by repeated flooding and recession events.

The reservoir area is situated in the Kohistan region's northern mountainous area. The Kohistan terrain was made by the tectonic collision of the Eurasian and Indian plates, resulting in the accumulation of a 40-kilometer-thick arrangement of ultramafic, mafic, plutonic, volcanic, sedimentary, and metamorphic rocks ^[20]. The Kohistan terrain is one of the world's best examples of a complete section through island arc crust from the uppermost supra-crustal rocks to the lowermost crust and adjoining mantle. The current literature regards the terrain as a composite that comprises back arc, island arc, intra arc, forearc, and current rift-related intrusions and diapers ^[21]. The Cretaceous-Paleogene Kohistan arc complex is one of the complete sections through a preserved paleo-island arc ^[22]. In addition, several researchers have described in detail the various units of Kohistan ^{[21,23,24,25]}.

The lowermost part of Kohistan comprises several ultramafic-mafic metamorphosed igneous complexes, the largest of which is the Jijal Complex. The Jijal Complex predominantly comprises lower dunites, harzburgites, and pyroxenites ^[25].

Two types of material are found in the reservoir area: quaternary deposits and bedrock. The quaternary deposits occupy about 74% of the total reservoir area. The material ranges from silt to boulder size, and the depth of these deposits changes along the reservoir. The rock outcrops mainly consist of norite and outcrops of amphibolite, hornblendite, gneiss, and diorite. Veins of ultra-basic mafic and ultra-mafic igneous rocks are also present. These outcrops cover about 26% of the reservoir area. Silicic andesitic to rhyolitic igneous rocks are the most prolific in northern Kohistan ^[25].

Over the reservoir length, the river valley is wide and U-shaped, with some narrow sections locally. The valley contains thick deposits of glacial residues through which the river has down-cut its current course. The alluvial slopes are sub-vertical at the river, and the tributary nullahs have deeply eroded the top of the alluvial material. The material ranges from silt to boulder size but is primarily granular.

The Karakoram glaciers, which can reach lengths of up to 60 km, are accompanied by latero-glacial sediment complexes over tens of kilometers. In addition to their large horizontal distribution, they are spread over a considerable vertical range and occur between 2500 m and 5000 m ^[26]. Moraine deposits with steep inclinations are of particular concern because they can potentially slide during reservoir impoundment. Such steeply inclined deposits are present on both sides of the Indus River.

The moraine deposits of Pleistocene to Sub-Recent age are in the reservoir rim area upstream of Chilas. The thickness of these deposits varies from less than a meter to several hundred meters. These extremely heterogeneous deposits can be observed from their exposures along the Karakoram Highway. These moraine deposits contain various materials ranging from silt to poorly graded gravel, large boulders, and cobbles. In the dry state, these slopes could stand in a nearly vertical orientation due to the presence of apparent cohesion due to the silt fraction. The silt loses its apparent cohesion as soon as it is fully saturated with water. Based on their gradation, moraines are further separated into strata as glacial tills, angular boulder gravel material (ABGM), rounded boulder gravel material (RBGM), and glaciofluvial and sand deposits.

2.   Selection of sites

Moraine deposits with a steep inclination and considerable heights are of particular concern because they are anticipated to slide during reservoir impoundment. A three-dimensional computer-generated (3D CG) representation of a terrain's surface can be made through the digital elevation model (DEM). A slope and an elevation map with a cell size of 10 m were generated from the ASTER DEM data in ArcGIS software and are attached as an appendix. A slope map gives a good indication of the slope angles of these deposits, and an elevation map indicates their heights. Using the slope map and elevation map along with Google satellite imagery from Google Earth Pro, fifteen (15) cross-sections from the most vulnerable reaches of the reservoir area were selected for performing the stability analysis, i.e., the sections with the steepest slope angles and maximum elevations. Figure 2 shows the geometries of two typical sections, and Figure 3 shows the location of these selected sections.

Figure 2.  Google Earth images of typical sections (Sections 7 and 10).

DownLoad: Full-Size Img PowerPoint

Figure 3.  Location of the sections chosen for analysis.

DownLoad: Full-Size Img PowerPoint

3.   Laboratory testing

To determine the engineering characteristics of the moraine deposits and assess their behavior on water submergence, bulk samples were collected from the reservoir area for laboratory testing. Due to the heterogeneous nature of moraines, samples were collected from four locations in the reservoir area. In this context, a visual reconnaissance was performed to select the sampling areas. The samples were collected from a depth of 1.0 m to 2.0 m from the existing ground level at different elevations of the proposed reservoir slopes.

Grain size analysis and modified Proctor tests were first carried out on the collected samples. The grain size analysis showed that the samples were generally either sandy gravels or gravely sand, with fines varying from 40% to less than 10%. The samples were all generally well-graded. The modified Proctor compaction tests on these samples showed that the maximum dry density ranged between 1.99 Mg/m³ and 2.17 Mg/m³, and the optimum moisture content ranged from 8.5% to 9.2%.

The bulk samples were compacted for the laboratory direct shear tests to obtain the density and moisture content as computed in the field. The samples were tested for both dry and saturated conditions. The normal stress versus shear stress plots for dry and saturated tests are shown in Figure 4. The test results show that the failure envelope is slightly curved, as is generally obtained for over-consolidated soils. The results indicate that the angle of internal friction for the moraines varies between 30 and 42 degrees, with the majority of the results showing a range of 35 to 40 degrees. It is also observed that the angle of internal friction is not affected to a great degree by saturation. The cohesion intercept value ranges from zero to about 175 KPa, the lower values being associated with saturated tests, indicating a loss of cohesion on saturation. This diminishing or even disappearing cohesion has to be considered in the stability assessment as it could be the main factor for reducing the stability of the rim slopes upon impounding. The pronounced cohesion of the moraine material most likely originates from the chemical cementation of the grain surfaces. This cementation is most likely a chemical cohesion caused by coatings of calcite, silicate, or a similar bonding agent. The bonding most likely developed through the precipitation of salts during evaporation when the material was deposited. Without a doubt, the bonding dissolves to a large extent under saturation, as clearly demonstrated by the direct shear test results.

Figure 4.  Direct shear test results; dry state (left), saturated state (right).

DownLoad: Full-Size Img PowerPoint

The results of the direct shear tests were thoroughly analyzed and compared with characteristics of similar materials using available literature and past data. The shear strength parameters of various materials delineated in the geological cross sections were then selected using engineering judgment for performing the stability analysis and are provided in Table 1.

Table 1.  Engineering parameters for stability analysis.

Material Type	Material Composition	Unit Weight, γ (kN/m3)	C (Dry) (kN/m2)	C (Sat.) (kN/m2)	Φ (Degrees)
Glacial Till	Gravels ≈ 20% Sand ≈ 30% Fines ≈ 50%	21	150	0	35
ABGM-I	Gravels ≈ 60% Sand ≈ 25% Fines ≈ 15%	22	30	0	42
ABGM-II		22	30	0	42
RBGM		21	30	0	42
Glaciofluvial	Gravels ≈ 40% Sand ≈ 50% Fines ≈ 10%	21	30	0	40
Sand	Medium to Fine Negligible Fines	20	0	0	34
Clayey Silt	Fines ≈ 80%	18	150	0	30

---

 | Show Table

DownLoad: CSV

4.   Results and discussion

4.1.   Deterministic slope stability analysis

Slope stability analysis on the selected sections was performed using deterministic and probabilistic slope stability analysis techniques. First, deterministic slope stability analysis was performed under various stages of reservoir operation, and the SLOPE/W module of the limit equilibrium software GeoStudio was used. The Morgenstern-Price approach was deployed as it satisfies both force and moment equilibrium. Geological cross-sections developed at different locations of reservoir stretches during the project feasibility studies were used to define the stratigraphy of slopes for stability analysis. The analysis was carried out separately under static and earthquake scenarios. Earthquake loading was applied using pseudo-static and dynamic analysis approaches, and the results were later compared.

4.1.1.   Slope stability analysis under static loading

It is expected that the stability of slopes consisting of unconsolidated materials (moraines) will be greatly reduced during the first impounding and under further reservoir operation. As indicated by the direct shear test results, there is a significant loss of cohesion on saturation in these moraines. Considering this, the cohesion value for the saturated condition is adopted as zero for all soil layers. This strength attenuation, along with the increased pore water pressures, would greatly contribute to slope instability. Once impounding starts, it will have the following three detrimental effects on the stability of the slopes:

i. Disappearance of the cohesion at the toe of the slope that became saturated.

ii. Reducing the weight of the slope toe by setting it under uplift and reducing the support of the slope by the friction.

iii. Local mechanical erosion along the shoreline by waves.

Impacts (i) and (ii) will take full effect long before reaching the full supply level. This means that slopes that are not stable under the full supply level will fail during the earlier stages of impounding. Even the slopes with a factor of safety (FOS) of little more than 1.0 for a full supply level would slide during the impounding procedure. Hence, the impounding was done in three steps for each slope to raise the water level from the present state condition to the full supply level (1160 masl), and stability analysis was performed at each step. This procedure was performed for all the selected cross-sections. Images of software modeling and analysis of one sample cross-section (Section 10) have been presented in Figure 5. It should be noted that the assumed water levels and slide geometries are indicative only to illustrate the interaction of impounding and flattening of the slope by successive collapsing of smaller slide masses.

Figure 5.  Analysis results for the impounding of the reservoir—Section 10: (a) initial water level, (b) water impounding to an elevation of 1105 m, (c) new slope formation at 1105 m, (d) impounding to 1140 m, (e) new slope formation at 1140 m, (f) impounding to 1160 m, (g) new slope formation at 1160 m.

DownLoad: Full-Size Img PowerPoint

Figure 5a shows the result of stability analysis at the initial reservoir level, i.e., an elevation of 1070 masl. Here, the FOS is slightly greater than one. Hence, the slope is considered to be just stable. Figure 5b shows the analysis results after first impounding to an elevation of 1105 masl. The FOS falls just below one; hence, a slide is predicted. It is assumed that the sliding will occur along the failure surface indicated by the software. Figure 5c shows the new slope geometry that will consequently be maintained at this reservoir level as the slope angle flattens from 63° to 45°. Figure 5d shows the analysis result after the second impounding on an elevation of 1140 masl, where the FOS is below unity, and the slope is predicted to slide again. The new flatter slope geometry is shown in Figure 5e. Figure 5f shows the third impounding, elevating the reservoir to a full supply level of 1160 masl and reducing the FOS to 0.936. The rim slope will slide again, and Figure 5g shows the possible final geometry with the slope settled at 38° after the third slide.

The simulation of the impoundment process shows that the stability of the slope depends on how close its inclination approaches the average friction angle. The same pattern was observed in the analysis results of other cross-sections. Therefore, on average, the slopes steeper than 35° will generally slide, and their cross-section will be significantly flattened after the initial impoundment. High slopes with an overall inclination between 25° and 35° will at least deteriorate locally in the steeper sections. Some of these slopes may also slide globally depending on slope geometry, geologic features, and water level position. Slopes flatter than 25° on average will not be prone to any morphological changes due to impoundment apart from rim erosion and local collapses of overstep sections.

Due to the unique composition of reservoir slopes, there is a strong probability of slope instability in moraine deposits under a drawdown event. As described earlier, the strength attenuation and increased pore water pressures would result in slope failures during impounding. This effect would be further amplified during drawdown as the stabilizing effect of water on the slope face would also be lost. The load case rapid drawdown may occur upon a certain period of reservoir operation. Some of the slopes will have slid during the first impoundment, as shown by the analysis results, and therefore, the geometry attained after the first impounding was used for these slopes to perform the drawdown case. The first step in performing the rapid drawdown analysis was estimating the rapid drawdown rate. Figure 6 shows the methodology adopted for performing the rapid drawdown load case.

Figure 6.  Flow chart for rapid drawdown analysis.

DownLoad: Full-Size Img PowerPoint

For calculation of the rapid drawdown rate, a hypothetical scenario was assumed to decipher the maximum possible rate that the reservoir could be lowered from the full supply level (1160 masl) to a minimum operating level (1060 masl) if all the gates of the spillway and all outlet structures in the dam body were opened. This situation may involve a possible leakage in the dam body or in the case of an incoming emergency flood. For this, a cumulative discharge rating curve for all the outlet structures in the dam body was first developed. Using this cumulative discharge rating curve and the elevation-capacity curve for the Basha Reservoir, the rapid drawdown rate was estimated to be 16.67 m/day.

The SEEP/W module of GeoStudio was used to perform a realistic simulation of the drawdown of the reservoir at the calculated drawdown rate. For this, a transient seepage analysis was performed, and the actual position of the phreatic line for each slope section was determined as the water level was drawn from the full supply level (FSL) to the minimum operation level (MOL). An incremental time stepping sequence was used, and 6 time steps were deployed, each corresponding to 86, 400 seconds, i.e., 1 day. Slope stability analysis was carried out after lowering the water level for each day and using the phreatic lines thus obtained. The sequential drawdown analysis is shown in Figure 7. A graphical representation of combined results for the three load cases, i.e., present state, impounding to FSL, and rapid drawdown, are shown in Figure 8. The sections/cases with less than one FOS are considered failing.

Figure 7.  Analysis results for rapid drawdown—Section 10: (a) full supply level—El. 1160 m, (b) after day 1—reservoir El. 1143 m, (c) after day 2—reservoir El. 1127 m, (d) after day 3—reservoir El. 1110 m, (e) after day 4—reservoir El. 1093 m, (f) after day 5—reservoir El. 1077 m, (g) after day 6—reservoir El. 1060 m.

DownLoad: Full-Size Img PowerPoint

Figure 8.  Slope stability results under static loading.

DownLoad: Full-Size Img PowerPoint

From the results, it can be observed that three slopes, i.e., Sections 7, 10, and 12, fail on the first filling of the reservoir. This situation can be attributed primarily to the loss of cohesion on saturation. Five slopes, i.e., Sections 4, 5, 7, 10, and 12, fail in the case of the reservoir's rapid drawdown. Cross-section 6 has the highest FOS values, which is also the flattest section with a slope angle of 21°. Meanwhile, cross-section 12 has the lowest FOS values, which is also the steepest section with a slope angle of 58°.

Once the stability analysis for the rapid drawdown load case was performed, an attempt was made to estimate the safe drawdown rate of the reservoir without causing instability in the moraine slopes. For this, the back analysis method was deployed. Stability analysis was performed on all of the sections using the drawdown rates starting from 6.67 m/day and gradually decreasing to achieve a reasonable FOS for each slope. After several iterations, a safe drawdown rate of 2 m/day was estimated for the Basha Reservoir. This implies that if, during the operation of the dam, the reservoir has to be lowered rapidly due to an emergency, it should not be lowered at a rate greater than 2 m/day. Otherwise, it could trigger multiple landslides in the reservoir slopes.

4.1.2.   Slope stability analysis under earthquake loading

The moraine slopes are very close to the main fault; hence, these slopes are anticipated to be subjected to severe earthquakes in the future. For performing the analysis, it is logical to assume that the moraine slopes will be subjected to an earthquake once the first impounding of the reservoir has been completed. Hence, the load case earthquake was applied to the new geometry obtained after the first impounding. Two cases were analyzed to get a good picture of the possible hazards associated with the earthquake event. The first analysis was performed considering the slopes to be submerged at the minimum operating level, and the second analysis was performed at the full supply level, i.e., the least critical scenario and the worst possible scenario has been studied.

Earthquake loading was applied using both pseudo-static and dynamic analysis approaches. In the pseudo-static limit equilibrium analysis approach, a lateral force coefficient, called seismic coefficient K, is used to apply earthquake loading. The seismic coefficient of K = 2/3 of the maximum design earthquake (MDE) was used for the pseudo-static approach. Since the MDE for the dam site is 0.37 g, a seismic coefficient K of 0.25 g was applied.

The dynamic approach includes the application of earthquake loading as a real-time shaking force in the form of time history. Design time history was developed as follows:

• A site-specific design response spectrum was developed (Figure 9).

Figure 9.  Design response spectrum for the dam site.

DownLoad: Full-Size Img PowerPoint

• A sample time history from past earthquake data for similar site conditions was collected.

• Design time history was developed from the design response spectrum and sample time history using the spectral matching technique in the SeismoMatch software (Figure 10).

Figure 10.  Design time history for the application of earthquake loading.

DownLoad: Full-Size Img PowerPoint

Compared to Newmark's method, the QUAKE/W stress method was used to perform dynamic analysis. Newmark's method can be unreliable for situations where there is a large buildup of excess pore-water pressures or where there is more than about a 15% degradation in the shear strength during the earthquake shaking ^[27]. At first, the initial stress conditions in the reservoir slopes prior to the application of earthquake loading were developed using the SIGMA/W module of GeoStudio. Dynamic loading was then applied on the slopes using QUAKE/W. The applied dynamic loading deforms the structural mesh and generates new slope stress conditions. Figure 11 shows the slope modeled for the computation of in situ stresses and the application of dynamic loading.

Figure 11.  Dynamic analysis: (a) in situ stress calculation in SIGMA/W, (b) application of dynamic loading in QUAKE/W.

DownLoad: Full-Size Img PowerPoint

Slope stability analysis was then carried out in SLOPE/W on slopes under these newly induced stresses, and their FOS was calculated. Figures 12 and 13 compare FOS values achieved from pseudo-static and dynamic analysis approaches. The results show that 7 out of 15 slopes failed under earthquake loading when the reservoir was at the minimum operating level (1060 masl), while 11 out of 15 slopes failed when the reservoir was at the full supply level (1160 masl). Hence, in the case of a seismic event of considerable magnitude, large volumes of moraine slopes sliding into the reservoir can be expected. The comparison shows that the FOS achieved with dynamic loading is higher than with pseudo-static loading. Therefore, it is reasonable to suggest that the pseudo-static method is somewhat conservative compared to the dynamic method.

Figure 12.  Comparison of pseudo-static and dynamic analysis (at MOL).

DownLoad: Full-Size Img PowerPoint

Figure 13.  Comparison of pseudo-static and dynamic analysis (at FSL).

DownLoad: Full-Size Img PowerPoint

4.2.   Probabilistic slope stability analysis

The study of landslide stability on mountain slopes becomes challenging when the sliding materials are heterogeneous ^[28]. Due to the very nature of its deposition, the moraines are very heterogeneous, i.e., they contain a wide variety of materials ranging from very fine silts to large cobbles. Hence, assigning a single value for each material's angle of internal friction or cohesion would not be reasonable. Also, the moraines are distributed throughout the 105-kilometer-long reservoir stretch, while the direct shear test results were carried out on samples from only four different locations. Hence, this research considered the possibility of input parameter variability.

Probabilistic slope stability analysis using Monte Carlo simulation (MCS) was performed on cases where the FOS obtained from the deterministic approach was between 1.00 and 1.20, i.e., on slopes that are "just safe". To calculate the probability of failure, the FOS values below one were counted, and then a percentage of this counted number, combined with the total number of Monte Carlo trials converged, was taken. The probability of failure can be defined as the probability (in terms of percentage) of obtaining an FOS value smaller than unity. Based on the range of variability assigned to each shear strength parameter, 2000 Monte Carlo trials were deemed adequate and performed for each case.

Various probability density functions are available in GeoStudio, including normal, lognormal, uniform, triangular, and generalized spline functions. This study used a triangular probability density function for the input geotechnical parameters, as it best represents the expected distribution of parameter variability for moraine deposits. A triangular probability density function can be specified with three points: the min and max values together with the apex value. The actual values adopted for the cohesion and friction angle in the deterministic analysis were used as the max/apex values, and the smallest values expected for each material, based on the shear test results and literature, were used as the minimum values. The max and min values adopted for each material are provided in Table 2. Figure 14 shows a representative triangular probability density function for the cohesion and friction angle used for glacial till. Here, the value of cohesion used for the deterministic analysis approach was 150 kPa. However, according to the probability function assigned, the value of cohesion will vary from 150 kPa to 100 kPa. According to the function assigned, there is a 40% occurrence probability of 150 kPa, whereas for 100 kPa, it is just greater than zero percent. The software performs multiple analyses using all of these values according to the assigned distribution function and gives the results in the probability of failure.

Table 2.  Max and min shear strength parameters for probabilistic analysis.

Material Type	C - Dry (kN/m2)		ϕ (Degrees)
	Max Value	Min Value	Max Value	Min Value
Glacial Till	150	100	35	30
ABGM-I	30	0	42	37
ABGM-II	30	0	42	37
RBGM	30	0	42	37
Glaciofluvial	30	0	40	35
Sand	0	0	34	29
Clayey Silt	150	100	30	25

---

 | Show Table

DownLoad: CSV

Figure 14.  Probability density functions for (a) cohesion, C, and (b) friction angle, ϕ.

DownLoad: Full-Size Img PowerPoint

Figure 15 shows a typical probabilistic slope stability analysis result (Section 10). The FOS from the deterministic approach for this slope was greater than 1. However, the probabilistic analysis shows that if there is a slight deviation in selected soil parameters, the probability of failure for this slope is very high, i.e., 76.4 %. Hence, this slope cannot be deemed safe.

Figure 15.  Probabilistic slope stability analysis.

DownLoad: Full-Size Img PowerPoint

The same is true for many load cases for Sections 3, 7, and 10, where the FOS is just greater than one, and hence, these slopes cannot be deemed safe. Figure 16 shows the representation of these results in graphical form, which indicates that, in general, the load cases having an FOS between 1.00 and 1.06 have a probability of failure greater than 50%, and these slopes are considered unsafe.

Figure 16.  FOS vs. the probability of failure.

DownLoad: Full-Size Img PowerPoint

4.3.   Summary of analysis results

Using the results from both deterministic and probabilistic analysis, it is concluded that four (4) out of fifteen (15) slopes are unsafe on the first impounding of the reservoir, six (6) out of fifteen (15) slopes are prone to failure in the case of a rapid drawdown of the reservoir at a rate of 16.67 m/day. In comparison, eleven (11) slopes are deemed to fail in an earthquake event with a peak ground acceleration (PGA) of 0.25 g with the reservoir at the full supply level. Slopes 6, 8, 13, and 15 are safe and stable under all loading conditions.

The results show that the slopes that have an overall inclination steeper than 35° are generally susceptible to failure on the first filling of the Basha Reservoir. For a high slope, its failure will occur in several steps of partial slides until a slope with sufficient global stability has been established after impoundment. The overall inclination of such a slope after impounding should be approximately 30° as it is related to the friction angle of the loose material accumulating on the bottom of the slope. Submerged slopes and slopes with an impounded toe having an overall inclination of some 25° to 35° might be flattened to a certain extent by the collapsing of local steep or weak sections. Rim erosion may take place when prevailing water levels are elevated. Submerged slopes and slopes with an impounded toe with an overall inclination of less than 25° should remain globally stable under impoundment.

In comparison, slopes steeper than 28° are susceptible to failure under an earthquake of peak ground acceleration (PGA) of 0.25 g. The results show various depths and inclinations of the slide masses; however, most of them show global failures of the entire slope when earthquake loading is performed. Much of the slide mass will accumulate at the bottom part of the slope, leaving the total slope with an FOS higher than 1.0. The calculations have been performed for the MDE. However, the factors of safety far below 1.0 for some slopes indicate that earthquakes of lower magnitude and higher probability also have significant potential for morphological changes. It is pertinent to mention that the site observations of the moraines, which are presently subjected to inundation due to periodic high flows in the river, indicate that typical circular failures do not generally occur in moraines on inundation. Rather, failure takes the form of sloughing and collapse. Therefore, the failure surfaces obtained through the slope stability analysis described above can be considered the limit for a first-time failure, regardless of the type of failure mechanism.

Based on the stability analysis results, landslide susceptibility maps were generated for static and earthquake cases separately and are attached as an appendix. The allowable drawdown rate calculated for the Basha Reservoir is 2 m/day and should not be exceeded at any time during the dam operation. The research also indicates that the pseudo-static slope stability analysis method is somewhat conservative compared to dynamic analysis, i.e., it gives smaller FOS values.

The results of the slope stability analysis show that considerable failures can be expected in the reservoir slopes of the Diamer Basha Dam. Considering the massive quantities of potentially unstable materials, landslides could have severe adverse impacts on the safety of the dam and associated critical structures. Each of these potential impacts was evaluated in detail.

(1) Abrupt falling off of big chunks of slope mass can create high waves, which can overtop the dam. The heights of waves generated by the two largest landslides (Sections 3 and 7) were calculated using an empirical approach given by ^[29]. Table 3 shows the wave height generated due to slope mass falling into the reservoir. The results show that the heights of the resulting waves arriving at the dam would be much less than the available freeboard (10 m). Therefore, the dynamic waves generated do not pose a serious hazard to the dam.

Table 3.  Wave height generated due to slope mass falling.

Section No.	Distance from Dam Body (km)	Maximum Wave Height (m)
3	50	2.9
7	59	4.7

---

 | Show Table

DownLoad: CSV

(2) Similarly, the sediment load in the reservoir would increase when the moraine slopes collapse. An attempt was made to estimate the total volume and load that will collapse into the reservoir using the results of the slope stability analysis. The quantities shown in Table 4 were calculated:

Table 4.  Sediment load in the reservoir due to slope failure.

Sr. No.		Case	Total Volume (m3)	Sediment Load (Tons)
1	Under the first filling of the reservoir		46 million (46 × 106)	57 million
2	First filling of reservoir + seismic event		152 million (152 × 106)	188 million

---

 | Show Table

DownLoad: CSV

The reservoir sedimentation studies indicated that the Indus River would carry 200 million tons of sediment annually into the reservoir. Therefore, the sedimentation due to the collapse of moraines on the first filling is almost negligible. At the same time, the total sediment in the case of a seismic event would contribute less than one year of the river sediment.

(3) The total number of moraines that will collapse and deposit in the reservoir in the case of an earthquake is estimated to be about 152 million m³. Nevertheless, most of the slide mass originates below the full supply level, which means the elevation-area capacity curve will be modified. However, the overall storage capacity of the reservoir will not be reduced alarmingly.

(4) The reservoir width varies throughout the length, with a maximum width of about 2600 m and a minimum reservoir width of around 500 m. The critical assessment revealed that the reservoir valley width at the location of major unstable moraine deposits is sufficient. The failure of massive moraine deposits should not cause any major blocking in the reservoir. However, at the reservoir tail, there is a section where the valley is very narrow (less than 500 m), and there is also a potential for slope failures, which may locally block the reservoir. However, this blockage will be temporary, as floods will create high velocity at the narrow locations and will ultimately erode the landslide material. However, this local damming may cause back-flooding in the reservoir. Back-flooding is usually slow and typically presents little life hazard, but there could be substantial property damage in the upstream areas.

5.   Conclusions

The main objective of this work was to perform a slope stability assessment of the reservoir slopes of the Diamer Basha Dam and assess the potential impacts of these slope failures. Fifteen critical cross-sections were selected, and stability analysis was carried out under various stages of reservoir operation for both static and earthquake scenarios. Earthquake loading was applied using both pseudo-static and dynamic analysis approaches. A probabilistic analysis was then performed, and input parameters (C, ϕ) were assigned to a probability distribution function. Landslide susceptibility maps were generated using the results of these analyses. The findings of the study are as follows:

• The results of direct shear tests indicate a significant loss of cohesion on saturation. The cohesion of the moraines most likely originates from the chemical cementation of the grain surfaces by coatings of calcite, silicate, or a similar bonding agent. This bonding dissolves under saturation, resulting in a loss of cohesion. It is also observed that the angle of internal friction is not greatly affected by saturation.

• As explained in Section 4.1.1, the progressive rise of the water level during impounding, setting the toe of the slopes under uplift, and diminishing cohesion are the decisive load cases. Their effects were studied in detail for all the selected cross-sections and were described in detail for cross-section 10.

• The analysis results indicate that large masses of moraine slopes will collapse into the reservoir. About 25% of the total slopes analyzed are prone to failure during the dam's first impounding. However, the intensity of the critical slopes increases from 25% to 40% during the rapid drawdown condition at the rate of 16.67 m/day.

• In the event of an earthquake at peak ground acceleration, moraine slopes pose a serious threat to the reservoir, where 11 out of 15 slopes are at risk of failure. A slope angle greater than 35° is more susceptible to failure at the initial reservoir filling event. In comparison, slopes with angles greater than 28° are critical in the event of an earthquake.

• In terms of FOS, analyses indicate that the pseudo-static slope stability analysis method shows conservative results compared to dynamic analysis.

• Slope stability analysis shows that considerable failures in the Diamer Basha Dam reservoir slopes can be expected due to the massive quantities of potentially unstable materials. Hence, landslides occurring during any catastrophic event may have severe adverse impacts on the safety of the dam and associated critical structures. Dedicated instrumentation is necessary for regular monitoring of critical slopes.

• The study can be used as a reference for a detailed investigation of any dam where reservoir stability could be an issue. It can also help develop a better understanding of failure mechanisms that cause stability issues in unconsolidated glacial/moraine deposits. For further research, the hazards posed by the generation of dynamic waves and blocking of the reservoir may be studied using analytical methods or model studies to get a more realistic picture of the potential dangers associated with moraine instability.

Use of AI tools declaration

The authors declare that they have not used artificial intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors gratefully acknowledge the technical support and assistance from the Civil Engineering Department at the University of Engineering and Technology Lahore, and the Water and Power Development Authority of Pakistan (WAPDA).

Data availability statement

The authors declare that all utilized data is included in the manuscript.

Conflict of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.