Evaluating pre- and post-failure behavior of a large waste dump by combining LEM and DEM approach: a case study from a hydropower station, China

1 Introduction

The civil engineering sector remains highly active in China and continues to generate substantial volumes of construction waste annually. A vast number of waste dumps already exist and are expected to expand rapidly to accommodate increasing waste materials from various sub-sectors, including mining (Wang et al., 2022; Liu et al., 2019), hydropower projects (Sun et al., 2020), highway construction (Zhang et al., 2025), and municipal solid waste (Zhan et al., 2018). A prominent example is the upcoming hundreds of pumped-storage hydropower projects over the next 10–20 years, a critical initiative for achieving China’s carbon neutrality goals by 2030 and 2060. Each of these huge projects will produce millions of cubic meters of excavated rock and soil waste, necessitating permanent disposal in designated valley-fill waste dumps (Li et al., 2024). However, these waste dumps often have relaxed design standards compared to main structures, often resulting in poorly graded materials, suboptimal geotechnical properties, and insufficient compaction. Under extreme hydrological or seismic conditions—particularly in cases of drainage system failure—such deficiencies may trigger large-scale slope instability (Sun et al., 2020). Valley-type spoil landslides can transform into high-mobility debris flows, causing ecological damage, soil erosion, and construction delays. Extremely, downstream river blockages may form catastrophic barrier lakes. Thus, rigorous pre-failure stability assessment coupled with post-failure long-runout landslide simulation plays an imperative role for effective risk mitigation.

For pre-failure instability assessment, the limit equilibrium method (e.g., Pan et al., 2021; Hu et al., 2022) or finite element-based strength reduction method (e.g., Feng et al., 2015) are still the most used techniques for obtaining the factor of safety (FS) at present. For example, Wang et al. (2021) employed the limit equilibrium method (LEM) to back-analyze the failure mechanism of a waste rock dump from a mining in Miyi County, China, which occurred due to the combination of inadequate compaction and the succumbed underlying stratum under heavy surcharge. Sun et al. (2020) evaluated slope stability at a slag disposal site using the Swedish arc method, examining factors such as site scale, slope gradient, location characteristics, and environmental conditions. Zhang et al. (2022) employed the finite element strength reduction method to investigate the effect of varying interstitial fillers on the stability of soil-rock mixture slopes in slag dump sites. The routine use of these two types of methods is often limited to a two-dimensional (2D) cross-section that is “subjectively selected” by the analyst. Thus, the resulting FS can only reflect the safety margin under any pre-selected conditions. However, waste disposal gullies typically exhibit pronounced three-dimensional (3D) topographic features, which provide multidirectional lateral confinement and enhance the overall stability of the waste dump. Thus, it becomes infeasible to select a 2D cross-section that conforms to the actual potential sliding mechanism. For instance, an arbitrarily chosen 2D profile cannot simultaneously follow a curvilinear and undulating 3D terrain while aligning with the probable sliding direction, inevitably resulting in a distorted FS.

The gully-type waste dump concerned in this study demands 3D slope stability analysis more critically compared to other scenarios, which could be achieved using 3D limit equilibrium or strength reduction methods. The latter, however, introduces significantly higher computational complexity and is more susceptible to non-convergence issues. In contrast, 3D LEMs retain much of the simplicity of their 2D counterparts while yielding more realistic FS. Nevertheless, achieving an ideal—or at least a viable trade-off—method to harmonize the simplicity of limit force-balance principles with the intricate mechanics of inter-slice forces within a rigorous theoretical framework remains an unresolved challenge, despite decades of research pursuit (Anagnosti, 1969; Lam and Fredlund, 1993; Zhu and Lee, 2002). While numerous studies have demonstrated their potential applications in various engineering contexts (e.g., Kalatehjari et al., 2015; Qi et al., 2021), these techniques have yet to transition from theoretical frameworks to widespread applications due to limited commercial software support. Particularly, there is limited evidence in the literature demonstrating the necessity of 3D limit equilibrium analysis for hydropower waste dump stability–a critical knowledge gap that this study aims to address.

Post-failure runout simulations of waste dumps have garnered significant research attention due to their critical role in quantitatively assessing risk consequences (Wang et al., 2022; Wang et al., 2021; Gao et al., 2019; Li et al., 2024). For instance, Liu et al. (2020) employed the Tsunami-Squares method to simulate the rapid flow dynamics of tailings sand in Chengui Town, China. Their study generated hazard maps by analyzing flow paths and deposition thickness, enabling a comprehensive evaluation of economic losses, fatalities, and environmental contamination for different hazard zones. Further advancing numerical techniques, Wang et al. (2022) developed an enhanced Smoothed Particle Hydrodynamics (SPH) model to investigate the long-runout landslide mechanisms in an iron mine dump, with particular emphasis on improving prediction accuracy in failure magnitude and impact extent through a modified kernel function. Gao et al. (2019) applied SPH to simulate a catastrophic landfill landslide with a specific focus on the propagation and fluidization processes during long-distance travel. These studies collectively contribute to a deeper understanding of flow dynamics and risk assessment. However, most existing simulation have focused on the post-failure behavior of waste dumps from either mining activities (e.g., Wang et al., 2021; Cho and Song, 2014; Liu et al., 2019) or municipal solid waste (e.g., Gao et al., 2019; Zhu et al., 2018) with continuum-based modeling frameworks. This modelling preference is well justified, as wastes from these two sources typically consist of fine particles that exhibit fluid-like behavior. Such behavior arises from mechanisms such as the development of excess pore-water pressure (Zhan et al., 2018) or air-cushion effects (Wang H. B. et al., 2020; Wang S. et al., 2020), as implied by Martín-Moreno et al. (2018). The implemented continuum-based constitutive models, e.g., the modified Cam-Clay (Wang et al., 2022) and the frictional-Bingham (Gao et al., 2019), have proven effective in capturing those phenomena.

Nevertheless, these continuum-based approaches are ill-suited for simulating rock fragment-dominated debris from hydropower projects, as they do not explicitly account for interparticle interactions—such as sliding, separation, and collisions—that govern the mobility of coarser granular flows. In such cases, the localized fluidization of fine particles between larger rock fragments becomes negligible. For these types of granular assemblies, the Discrete Element Method (DEM) demonstrates its unique advantages (Spreafico et al., 2017; Zhu et al., 2022), as it is specifically designed to model discrete particle interactions. However, the DEM application to the study of hydropower-related waste dumps remains scarce in the literature. Recent work by Li et al. (2024) demonstrated the feasibility of DEM in this context, which also indicated the challenges in dynamically calibrating mesoscale parameters (e.g., interparticle friction coefficients) to reduce the simulation uncertainties, according to specific conditions and environments. Thus, further DEM-based investigations are required to establish its applicability for hydropower-related waste dumps.

The aim of this study is twofold. First, to explore the effects of pronounced three-dimensional (3D) topographic features on the stability of a gully-type waste dump, we develop a state-of-the-art 3D LEM model for analyzing both global and localized stability under various potential failure mechanisms. This self-coded 3D LEM was recently proposed to overcome critical limitations inherent in most existing frameworks, including oversimplified assumptions regarding inter-column forces and incomplete satisfaction of force/moment equilibrium conditions. Subsequently, we construct a 3D DEM framework to investigate the constraints from meandering valley on the runout behavior of failed waste dump materials. To address the model calibration challenge, a two-stage calibration strategy is proposed: The first stage was used to identify secondary DEM parameters with limited influence on runout behaviour and to establish a reasonable model configuration, whereas the second stage focused on the microscopic friction coefficient μ because of its dominant influence on interparticle sliding and post-failure mobility. A key highlight here is the qualitative cross-model comparison between the pre-failure 3D LEM and post-failure DEM simulations at the engineering scale, for ultimately refining the representative parameter. This project-tailored parameter determination ensures more reliable and targeted stability assessment and runout prediction, particularly given the enhanced rigor of 3D methodologies. We demonstrate the modelling strategies of both the pre-failure 3D LEM and post-failure DEM simulations in evaluating the potential risk of a waste dump in a pumped storage power station project in Yongjia County, Zhejiang Province, China. Results provide cross-model insights into the importance of targeted reinforcement measures, as well as critical data for future hazard assessments.

2 Project overview and site conditions

2.1 Project general information

The Yongjia Pumped Storage Power Station is one of more than hundreds key project under China’s “Medium- and Long-Term Development Plan for Pumped Storage Energy (2021–2035).” Located approximately 320 km from Hangzhou, Zhejiang Province, this project comprises upper and lower reservoirs with a 577 m hydraulic head (Figure 1a). The two reservoirs feature 60 m- and 81 m-high concrete-faced rockfill dams with storage capacities of 5.91 and 6.15 million m³, respectively. Major earthwork quantities for the main structures include 5.96 million m³ of earth\rock excavation, 2.34 million m³ of earth\rockfill and 1.22 million m³ of underground excavation. Excluding dam fill material, the surplus soil and rock fragments are >3,000,000 m³ and will be permanently stored in two designated waste dumps (Figure 1a), close to the upper and lower reservoirs, respectively.

This study analyzes the waste dump adjacent to the lower reservoir (1,740,000 m³ in total), which presents three primary hazards: (1) slope instability, (2) soil erosion, and (3) environmental contamination. The significant volume of accumulated waste threatens critical infrastructure at the valley mouth, including the stakeholder’s camp, steel pipe plant, construction facilities, and Yangwan Village (Figures 1a,c,d). The most severe risk involves potential waste mass landslides that can be triggered by extreme precipitation events through groundwater level elevation, which could catastrophically impact downstream areas and obstruct the 15-m-wide Zhang’ao Creek, significantly exacerbating flood risks. Without proper proactive mitigations integrated into the design, these hazards may lead to severe ecological degradation and long-term environmental pollution.

2.2 Topography, geology, and waste dump

The studied waste dump is situated downstream of the dam’s left bank within Zhukeng Gully (Figure 1), oriented roughly North-South and is approximately 450 m from Yangwan Village. The terrain features slopes of 25°–45°, with gentler slopes along ridges, saddles, and gully mouths, while the left bank exhibits steeper cliffs. The overburden layers (1–3 m thick) consist of residual-colluvial gravelly clay and slope-derived boulder-silt mixtures. The gully bottom contains alluvial pebble-cobble deposits (2–5 m thick) underlain by rhyolitic welded tuff. The weathered zone (2–8 m thick) is discontinuous with exposed bedrock in the gully.

The waste materials are characterized from three primary sources. (1) excavated earth predominantly composed of gravelly silty clay (thickness: 0.5–4.5 m) that is sourced from the dam’s right bank. (2) Quarried rocks consist of welded tuff extracted from fully to strongly weathered layers (0.5–4.0 m thick) and (3) moderately weathered zones (8–10 m thick). The waste dump occupies an area with elevations ranging from 50 to 125 m, extending 560 m along the gully with an average longitudinal gradient of 134‰. The total designed storage capacity is approximately 2.02 million m³. The engineered slope features a crest elevation of 120.00 m with a 1:2.0 (H: V) front slope ratio, 3-m-wide berms are incorporated at 20-m vertical intervals. During disposal, the spoil surface is graded at 3%–5% from the centerline toward both sides to ensure effective drainage. Stability is enhanced through layered placement with compaction, particularly in areas designated for temporary surface facilities. A reinforced concrete retaining wall is constructed at the toe.

2.3 Meteorological and seismic conditions

The site has a subtropical monsoon climate with distinct seasonal precipitation. Typical rainfall intensities for precipitation events with different return periods at the site are illustrated in Figure 2. The mean annual rainfall reaches 1,567.2 mm, predominantly (80%) concentrated between March and September. The site is located within the seismically active South China Coastal Zone, where 10 potential seismic sources have been identified (3 with M max = 6.0 and 7 with M max = 5.5). Destructive seismic impacts primarily come from moderate to strong earthquakes nearby or large distant earthquakes. The site has a basic seismic intensity of VI, as per Chinese seismic intensity scale.

3 Numerical method and model

3.1 Instability initiation computation

3.1.1 2D and 3D limit equilibrium method

A reliable stability assessment of waste dump slopes depends on three key factors: (1) simple yet efficient calculation methods; (2) comprehensive evaluation scenarios, and (3) appropriate geotechnical parameters. Given that the two-dimensional (2D) LEM remains the prevailing methodology, this study first applies the well-established Morgenstern-Price method to determine the FSs of the waste dump under two conditions: (1) perennial groundwater levels and (2) 50-year return period rainfall events, as it has capability of accommodating slip surfaces of arbitrary geometry. The method accommodates slip surfaces of arbitrary geometry, and its theoretical basis and engineering reliability have been extensively verified through benchmark solutions and practical applications (Morgenstern and Price, 1965).

It is important to emphasize that conventional 2D methods fail to account for critical 3D effects inherent in slope stability problems, typically yielding overly conservative results. The waste dumps are situated within a topographically complex mountain gully along the left bank of Zhang’ao Creek (Figure 3a: 3D perspective view, and Figure 3b: plan view). The site’s geomorphology is characterized by significant surface irregularities, with a distinctive boundary configuration analogous to a tri-digit palmate structure. It features three converging sub-valleys that intersect at the gully’s central axis and is flanked by two moderately inclined ridges at the downstream terminus. Such complex 3D geomorphic features necessitate 3D method to properly account for the spatial distribution of stability-influencing factors.

However, current 3D LEMs exhibit several theoretical and practical limitations: (1) incomplete consideration of inter-column forces, (2) partial satisfaction of equilibrium conditions, and (3) computationally intensive. Thus, commercially available 3D platforms have yet achieved broad adoption. To address these constraints, this investigation employs a novel self-developed 3D methodology (Qi et al., 2021). The methodology extends the fundamentals of 2D Morgenstern-Price Method into 3D space by preserving its rigorous treatment of inter-slice forces (Figure 4). This methodology offers three significant advantages: (1) incorporation of all inter-column normal and shear forces (Figure 4a), (2) strict adherence to complete static force/momentum equilibriums, and (3) demonstrated validation against both numerical benchmarks and field observations (Qi et al., 2021). This approach provides closed-form solution via a high-efficient triple-nested iterative algorithm. It is also versatile to enable quantitative assessment of hydrogeological effects through arbitrary groundwater table specification and seismic loading conditions. The method has been also validated against several field cases, including the Chengdu Municipal Landfill facility (Qi et al., 2021) and the catastrophic Shuicheng landslide in Guizhou Province (Figure 4b; Zhou et al., 2021). These previous validations provide the methodological basis for its application to the present valley-type waste dump. It should be noted that these previous benchmark and field-case comparisons validate the adopted 3D LEM methodology rather than the site-specific predictions of the present waste dump.

The phreatic surface used in the LEM analysis was obtained from a saturated-unsaturated seepage analysis governed by mass conservation and Darcy’s law, expressed in vector form as ∂ θ / ∂ t = − ∇ · q and q = − K θ ∇ h + z, where θ is volumetric water content, q is the water-flux vector, K θ is hydraulic conductivity, h is pressure head, and z is elevation head. The nonlinear hydraulic behavior under unsaturated conditions was described using the Fredlund and Xing (1994) soil–water characteristic curve (SWCC):wherein which θ h is the volumetric water content at matric suction head, h θ s is the saturated volumetric water content, and a f, b f, and c f are fitting parameters controlling the shape of the SWCC.

θ v = C h θ s ln e 1 + h a f b f c f

C h = 1 − ln ⁡ 1 + h h r ln ⁡ 1 + 10 6 h r

For the waste material, the fitting parameters were estimated from the grain-size distribution following the procedure of Houston et al. (2006), with c f = 5.02 and h r = 100. The corresponding unsaturated hydraulic conductivity was represented by the relative hydraulic-conductivity functionwhere K is the hydraulic conductivity under unsaturated conditions and K s is the saturated hydraulic conductivity. The parameters used for the seepage analysis in PLAXIS 3D are summarized in Table 1.

K r h = K h K s

TABLE 1. Hydraulic parameters used in the saturated–unsaturated seepage analysis.

θ r and θ s denote the residual and saturated volumetric water contents, respectively; a f and b f are Fredlund–Xing fitting parameters; and K s is the saturated hydraulic conductivity. For the waste material, c f = 5.02 and h r = 100 were estimated from the grain-size distribution.

For the rainfall scenario, the 50-year return-period rainfall shown in Figure 2 was applied to the exposed ground and dump surfaces as a time-dependent normal flux boundary. The perennial groundwater field was used as the initial condition, while the remaining external boundaries were assigned according to their hydraulic roles as drainage at model laterals or no-flow boundaries at model bottom; the zero-pressure-head contour obtained from the seepage calculation was then adopted as the phreatic surface in LEM analyses. The pore-water pressure acting on the slip surface was calculated from the water head, defined as the vertical elevation difference between the phreatic surface and the corresponding point on the slip surface.

3.1.2 Column-discretized computational model

3D slope stability analysis critically depends on an accurate computational mode, high-resolution topographic data of the waste dump area were first acquired using unmanned aerial vehicle (UAV) and processed in ArcGIS to generate a digital elevation model. To optimize computational efficiency, we delineated this area by a rectangular domain encompassing merely the waste dump (Figures 5a,b). Finally, the raster-based approach in ArcGIS was adopted to discretize the model and each 3D column is defined by a 4 m × 4 m horizontal projection. Note that the model does not consider elevation-specific geological parameters, since the critical slip surfaces primarily traversed the waste dump mass or along the pre-existing ground surface, rather than the underlying strata as demonstrated in strength parameters discussed later (Table 2). The discretized domain comprises a 131 × 119 grid (Figure 5c). Preliminary refinement trials showed that further reducing the column size produced only negligible changes in the global FS and the principal spatial distribution of the local FS, while substantially increasing the computational cost. Therefore, the adopted discretization is considered sufficient for the present engineering-scale analysis; minor uncertainties associated with discretization are not expected to affect the representativeness of the calculated FS, while the column-scale local FS spatial distribution remains effective in characterizing three-dimensional variations in sliding-mass thickness, basal inclination, and lateral confinement. The area exhibits pronounced 3D terrain characteristics, featuring three converging valleys in the mountainous hinterland. Near the outlet, two prominent ridges create a distinct morphological constriction, resulting in an S-shaped thalweg.

TABLE 2. Shear strength parameters.

Regarding this, three primary potential sliding directions were identified (Figure 5b) and five representative slip surfaces were selected for analyses (Figures 5d–h):

• Sliding Direction I incorporates two 3D slip surfaces (A and B), where the sliding mass is unimpeded by topographic constraints. Two cross-sections 1-1 and 2-2 were selected for comparative 2D analysis.
• Sliding Direction II involves Slip Surface C, characterized by relatively shallow failure depths. It Corresponds 2D analysis on Cross-section 3-3.
• Sliding Direction III includes Slip Surfaces D and E, situated above a left-bank topographic protrusion. This failure mechanism involves shallow sliding with significant right-bank geometric confinement. Corresponding 2D cross-sections are donated by 4-4 and 5-5.

Additionally, a basal failure mode was considered along the pre-existing ground surface (Slip Surface F, Figure 5c). All the 2D cross-sections are shown in Figure 3b.

The geological survey indicated that the stratigraphy of the waste dump comprises (from top to bottom): gravelly silty clay (thickness: ∼3 m), strongly weathered rock mass (∼5 m), moderately weathered rock mass (∼35 m), and slightly weathered rock mass. Table 2 summarizes the triaxial testing-derived strengths. The waste dump has a relatively lower strength than other geological strata. Yet, compared to conventional soil-rock mixed slopes—or even well-graded embankment dam fills of analogous composition—the waste dump material exhibits higher strength variability due to its inferior lithology, suboptimal gradation, and inadequate compaction. Triaxial tests indicate that its shear strength is marginally lower than that of heavily compacted rockfill but remains within documented ranges for similar waste materials (Dong, 2006; Yan, 2003; Zhai, 2009; Tan and Deng, 2009). Note that the concrete retaining wall interfaces with the natural cohesive foundation. Following the Code SL379-2017, the base-soil interface parameters are selected as 0.9 φ and 0.2 c, where φ and c denote the natural soil’s strength.

3.2 Post-failure simulation

3.2.1 Discrete method and constitutive model

As discussed in the Introduction section, most existing studies adopt the smoothed particle hydrodynamics (SPH) method to assess the post-failure flow of waste dumps. This continuum-based approach, incorporating the frictional Bingham model (Gao et al., 2019), is suitable for simulating the flow-like motion of loose mining or landfill spoils of very fine particles. However, the waste dump studied herein consists of loosely packed rock fragments with weak interparticle bonds (Figure 1b), falling under the category of discrete media. A DEM is required to capture the large discontinuous deformations, such as particle sliding, overturning, and separation during motion. Here, we employ the DEM-based particle flow code, i.e., PFC3D, to simulate the full failure progression along this gully, encompassing stages from initial instability and particle disintegration through flow propagation, energy dissipation via particle collisions, to final deposition. Through systematic investigation of how natural gully morphology influences flow paths and quantitative assessment of terrain effects on deposition patterns, this simulation is dedicated to offering a rigorous mechanistic evaluation of runout hazards.

The computational cycle in PFC3D sequentially alternates between the force-displacement law and Newton’s second law to determine interparticle contact forces and the resultant translational/rotational motion of particles, respectively. It offers diverse modeling strategies to represent particles of arbitrary geometries, as well as various constitutive laws governing particle-particle interactions. However, given the immense scale of the waste deposit, complex particle morphologies and advanced constitutive models would demand prohibitive computational costs—an impractical trade-off for their marginal gains in simulation accuracy. Thus, spherical particles are adopted for simulating the dynamic response of such large-scale waste deposits under complex topographic loading paths.

The linear bond model that extends the basic linear contact formulation by introducing a bond mechanism is implemented. This mechanism can capture cohesive effects arising from weak interparticle bonding, such as those attributable to fine-grained cohesive soils filling the interstitial voids among rock particles. It combines a linear elastic contact model with a parallel bonding mechanism to capture both normal/tangential forces and interparticle cohesion. The normal and tangential contact forces are governed by: where, k n and k t are normal and tangential stiffness, respectively, and determined by the Young’s modulus E and particle radii. δ n and δ t are normal and tangential overlap displacements, respectively. μ is friction coefficient. The tangential displacement follows a stick-slip behavior constrained by the Coulomb friction limit. The parallel bond introduces a finite-sized bonding area between particles R, and the bond stresses are computed as where I = π R 4 4 and J = π R 4 2 are moments of inertia. The bond will break if σ > σ c or τ > τ c.

F n = k n · δ n

F t = min k t · δ t , μ F n

σ = F n b π R 2 + M t b R I

τ = F t b π R 2 + M n b R J

3.2.2 DEM of the waste dump and its calibration

The DEM model integrates two key components: the high-resolution terrain model and the granular waste representation. The complex high-resolution terrain model was constructed in Rhinoceros by interpolating elevation data from UAV and subsequently converted into an STL-format triangulated surface. The imported mesh, consisting of 25,948 facets, was transformed into rigid wall elements. The domain was large enough to encompass the dump area and its potential sliding path while ensuring sufficient resolution to accurately capture terrain effects with the linear particle-wall contact model. Given that the waste material mainly composed of coarse rock fragments, the bulk waste mass was simulated as an assembly of 125,000 spherical particles with radii 0.3–0.6 m, a bulk density of 1950 kg/m³ and a porosity of 0.3. Preliminary simulations suggested that further particle refinement (e.g., explicitly modeling small soil particles and/or increasing particle number) would provide limited improvements in accuracy, indicating an appropriate trade-off between model fidelity and computational efficiency for this large-scale waste dump of around 1,740,000 m³.

The DEM’s constitutive model requires calibration of microscopic parameters that determine the macroscopic kinematic process of granular system. Their values depend inherently on particle geometry and physical properties, yet their experimental determination remains challenging, and rigorous micro-macro theoretical relationships are lacking. Virtual triaxial simulations are commonly adopted to calibrate these parameters by matching stress-strain curves between simulated and experimental data, this approach has several limitations for our purpose: (1) it requires intensive trial calculations or sophisticated mathematical optimization (Kim et al., 2021); (2) inherent particle assembly’s heterogeneity means the calibration would result in non-unique solutions rather than a single optimal parameter set; and (3) virtual triaxial simulations that characterizes pre-failure deformation and strength properties under quasi-static conditions cannot capture the post-failure dynamics. Given the fundamentally critical role of microscopic parameters and the strikingly limited DEM simulations specifically addressing dump sites in hydropower projects, this study systematically calibrates DEM parameters through an integrated approach. The methodology starts with parametric sensitivity simulations informed by existing simulations on analogous geotechnical scenarios (e.g., general landslides, mining and hydropower waste dumps). The stability-analysis results were then used for a qualitative cross-model consistency assessment of the DEM simulations.

Existing numerical simulations of both natural soil landslides (Zhang et al., 2012; Zhou et al., 2019; Wu et al., 2020) and artificial waste dumps (Li et al., 2024; Lai, 2019) have established the following parameter ranges: E = 10⁸–10⁹ Pa, η = 1.0–2.5 for linear spring; σ _*n* = 10⁷ and τ = 10⁷ Pa for parallel bonds, and μ = 0.1–0.65 for frictions. It should be emphasized that μ is the microscopic friction coefficient assigned to particle contacts in the DEM, whereas φ is the macroscopic friction angle describing the bulk shear strength of the waste material in the LEM; the two parameters act at different scales and cannot be directly equated. Note that k _*n* and k _*t* are derived from particle radii (R), E and η, rather than being directly specified, to ensure more realistic collision responses. The viscous damping ratio was explicitly set to 0.3 to account for energy dissipation during particle interactions. This viscous damping coefficient was determined based on its well-established relationship with the restitution coefficient (Ai et al., 2022), which governs the inelastic behavior of contacting particles. Based on these ranges, preliminary simulations were conducted to assess the relative influence of key parameters on critical output metrics, including sliding scale, velocity, and influence range. It is shown that both E and η = k n / k t exhibit negligible effects on post-failure behavior, particularly on runout distance and deposition morphology (Figures 6a–e). This is because these microscopic parameters predominantly govern particle deformation during collisions rather than the macroscopic dynamics. Notably, while E does not significantly alter outcomes, it, however, impacts the computational efficiency; thus, E = 10⁸ Pa was adopted subsequently to balance accuracy and efficiency. Further, the remaining parameters exhibited moderate effects except for μ.

Comparison between Figures 6f,g and Figure 6a illustrates the distinct post-failure behaviours by varying μ. For μ = 0.65 and 0.375, only superficial particle mobilization occurred within the dump area, demonstrating a fairly stable system under the control of retaining structure. However, a reduction to μ = 0.1 triggered basal failure initiation, characterized by rapid propagation and substantial runout distances. The particle cascades extended through the gully before depositing on the foothill plain. This aligns with the DEM investigations of Li et al. (2024) regarding friction-controlled sliding mechanisms of a hydropower waste dump. This can be attributed to the fact that μ is the dominant parameter governing post-failure dynamics through its influence on inter-particle slippage behavior. Thus, the potential hydrological effects in the deposition zone as previously introduced warrant particular attention. Intensive rainfall in this area (Figure 2), as the major triggering factor for waste landslides, may increase saturation or reduce effective stress within the waste. The lubrication effects induced by water films can significantly weaken interparticle interactions, which will be further declined by the development of excess pore water pressure. These two factors complementarily indicate clearly lower μ in field conditions (An et al., 2021; Li et al., 2012; Wei et al., 2022). Therefore, refined simulations under low-μ range are required to better characterize the actual failure processes. Subsequently, we will employ the calibrated DEM model with computational fidelity and efficiency to focus on effect of low-μ on failure propagation mechanisms.

4 Computational results

4.1 Stability computation results

For 2D stability analysis, the five representative cross-sections were strategically selected considering both pre- and post-deposition topographical characteristics. For each section, the entry-exit search algorithm in GeoSlope was employed to determine the most critical slip surface. Three predetermined failure mechanisms were considered to guide the automatic search algorithm: (i) entry from the crest and exit from the slope surface; (ii) both entry and exit from the slope surface; and (iii) entry from the crest and exit from the bottom of the concrete wall. These three mechanisms were selected to represent the principal entry–exit configurations permitted by the geometry of the present waste dump, including shallow slope failure, deeper failure extending from the crest, and toe-controlled failure associated with the retaining wall. They are intended as representative search scenarios for the present site rather than an exhaustive description of all possible failure mechanisms. Table 3 presents the computed 2D FS under (1) perennial groundwater levels and (2) 50-year return period rainfall events. The phreatic surfaces were derived from coupled unsaturated-saturated finite element seepage analyses by subjecting the slope to the corresponding rainfall scenarios (Figure 2), the modelling details are omitted here for brevity. Figures 7a–e illustrates the representative computed results for five cross-sections.

TABLE 3. Computed FS for 2D and 3D conditions.

Pre-/post-“/” values show FS, for perennial groundwater and 50-year rainfall cases.

All 2D cross-sections maintain FS > 1.0 for all potential failure mechanisms with the FS values under the perennial groundwater conditions systematically greater than those under a 50-year return period. Comparatively, Failure Mechanism iii represents the most danger mechanism for all five 2D sections, and the most critical case emerges at Section 2-2 with the minimum FS of 1.239. The FS variations among cross-sections are closely associated with pre- and post-deposition topographic conditions. Section 1-1 has moderate inclinations in both the native terrain and emplaced waste material. Section 3-3 features near-horizontal native terrain adjacent to the toe. Sections 4-4 and 5-5 display convex-upward basal geometries near to the toe. Conversely, Section 2-2 presents an inclined original terrain and a substantially steeper overlying deposit inclination, which constitutes the principal instability factor. Parametric studies on phreatic surface variations reveal that Failure Mechanism iii demonstrates pronounced sensitivity, in particular, Section 2-2 shows the most FS reduction. A conservative trial calculation, assuming complete saturation of waste mass due to prolonged high-intensity precipitation concurrent with drainage system malfunction, results in a near-critical safety factor (FS = 1.06), leaving minimal stability margin.

The most 3D critical slip surfaces are automatically located through a simulated annealing algorithm in each local area (Figure 5). Under both hydrogeological conditions specified above, all computed 3D FSs surpass their corresponding 2D values in the same region. This observation is unequivocal and can be attributed to the enhanced sliding resistance from the 3D terrain features along both flanks of the primary sliding direction, which superimposes onto the above 2D cross-sectional topography effect—A finding consistent with well-documented evidence (Qi et al., 2021).

Our 3D LEM can not only compute the global FS but also use the spatial distribution of the column-scale local FS indicator to illustrate the stability variations associated with topographic changes in sliding-mass thickness, basal inclination, and lateral confinement. Figures 8a,b shows the spatial FS distribution across discretized columns for Slip Surface B (correlating with 2D Section 2-2 with the minimum FS) under both hydrogeological scenarios, which reveals a distinct trough-shaped profile, with minimized values at the central portion and progressively higher values toward both left and right extremities. This pattern directly reflects the 3D geometry of failure mass with progressively thinner thickness and gentler slipping inclinations towards both lateral boundaries (see Figures 5a,b,e). For the potential slip surface along the foundation-waste interface, Figures 8c,d shows that the majority of columns on the left-side demonstrate exceptionally high FSs, contrasted by moderately elevated values for a limited of columns at the right terminus. This asymmetry arises from the constraining effect of bilateral protruding ridges flanking the gully on waste mass mobilization (Figures 5a,b). Although this geometric configuration produces a seemingly favorable FS of 1.9, closer examination reveals that the right-side ridge aligns geometrically parallel to the gully outlet orientation, leaving an unobstructed kinematic pathway for right-side waste migration toward the lower-left discharge area and accounting for the diminished FS for this particular slip mechanism of slip surface B. If this adverse topography coincides with the extreme scenario of fully saturated wastes, the 3D FS drops sharply to 1.15. Despite the nominal safety margin, potential for localized slope instability persists under these compounding effects, considering also the significant inherent variability in the waste strengths due to possible inadequate compaction.

4.2 DEM computation results

Previous trial simulations have indicated that the post-collapse behavior of the waste body is likely to be controlled by a relatively low μ. Following Li et al. (2024), who suggested that μ is unlikely to be <0.1 considering the natural geomorphological characteristics of the gully, we systematically varied its value from 0.1 to 0.2 to examine the variations in the runout dynamics (Figures 9–11), including velocity evolution (Figure 12a), runout distance (Figure 12b), and final deposition morphology. The three values represent a low-friction condition with extensive mobilization (μ = 0.10), an intermediate condition (μ = 0.15), and a transition condition between long-runout failure and localized movement (μ = 0.20), respectively.

At μ = 0.1, a two-stage failure developed along the meandering gully, comprising an initial rapid toe failure and subsequent delayed upslope mobilization (Figures 9a–h show temporal evolution from 10 s to 3,000 s). The initial stage featured rapid acceleration of the mobilized mass, initiating from the toe and attaining a peak velocity of 15.31 m/s in 10 s (Figure 12a). Subsequently, the mass decelerated at gentler terrain (10–40 s) before undergoing renewed acceleration along steeper gully segments (40–60 s). The leading runout particles from the primary failure eventually decelerated upon reaching the alluvial plain at the gully mouth, where deposition occurred gradually (60–480 s). The secondary stage involved the delayed mobilization of upslope material. This subsequent failure propagated over the previously deposited primary mass (480–1,200 s) before reaching complete stabilization at 3,000 s). Although the secondary failure substantially increased depositional thickness along the gully and expanded the distal inundation area, it did not significantly alter the maximum runout distance (680 m) achieved by the leading particles within 120 s (Figures 12e–h).

At μ = 0.15, the failure initiated similarly at the slope toe, with kinematic characteristics analogous to that at μ = 0.1 (Figures 10a–h show temporal evolution from 10 s to 1800 s). After overtopping the retaining wall, the failed mass exhibited three distinct mobility phases: initial acceleration (0–10 s), transitional deceleration across moderate gradients (10–30 s), and secondary acceleration along steeper gully segments (30–120 s), see Figure 12a. Subsequent deposition on the alluvial plain at the gully mouth (120–540 s) was characterized by velocity fluctuations with net deceleration, ultimately resulting in complete stabilization within 600 s. Comparatively, there is a significant reduction in mobilized volume without secondary slope failure. Terminal deposition patterns consisted of only several isolated, poorly sorted accumulations across the fan surface (Figures 10f,g,h). Runout prediction revealed an earlier attainment (120 s) of the reduced maximum travel distance (667 m, Figure 12b).

Further increasing μ to 0.2 significantly reduced the mobilized volume at the toe (Figures 11a–f show temporal evolution from 2 s to 1800 s). The entire sliding process completes in 120 s with three characteristic velocity phases: rapid acceleration (0–8 s), sharp deceleration (8–30 s), brief recovery (30–40 s), before progressive stabilization in mid-slope colluvial gullies with moderate gradients (Figure 12a). The resultant maximum runout distance decreased substantially to just 80 m (Figure 12b).

Physically, increasing μ enhances basal frictional resistance, which reduces the net driving force available for acceleration and increases the rate of frictional energy dissipation along the travel path. Consequently, both the peak velocity and the final runout decrease with increasing μ. However, there are only negligible differences in runout distance (Δ < 3%) between the μ = 0.1 and μ = 0.15 cases. This is because the initially high-velocity toe-originated masses decouple from the main sliding body behind. The kinematics of these leading-edge masses are primarily governed by local damping coefficients and topographic constraints rather than intergranular friction. At μ = 0.2, rapid stabilization occurs with runout and impact zone confined to proximal slopes, effectively eliminating potential threats to downstream infrastructures on the piedmont plain.

5 Discussion

Previous DEM simulations have underscored the critical importance of representative μ selection in reproducing the runout and deposition behavior of waste materials: mid-slope confinement (μ ≥ 0.2), partial gully mouth reach (μ = 0.15), and extensive downstream impacts (μ = 0.1). Attempts have been made to correlate this micro-scale friction coefficient with macro-scale shear strength, which, however, exhibits limited generalizability beyond test conditions. Laboratory-derived cross-scale relationships typically cannot be directly transferred to field-scale scenarios.

Given this context dependence and the fact that these two parameters respectively control pre-failure instability and post-failure runout behaviors at the same field scale, it would be more reasonable to examine the cross-scale macro-micro relationships directly within this actual engineering context. Here, we estimate the appropriate range of μ for this particular problem through comparative analysis between pre-failure stability LEM calculations and post-failure dynamic DEM simulations. Figures 13a–f illustrates the DEM-simulated incipient failure dynamics of waste materials for different values of μ (=0.10–0.65). At μ = 0.2, the incipient failure (denoted by red coloration with velocities largely exceeding ambient particles) is confined at left lateral margin of the valley outlet (Figure 13d). This localized failure pattern shows remarkable correspondence with both 2D and 3D limit equilibrium computations where the Cross-section 2-2 and Slip Surface B possess most critical safety margins. Any other external disturbances would initiate this specific local failure mode (FS ≈ 1.0) prior to other potential failure mechanisms. The remarkable DEM–LEM consistency qualitatively supports μ = 0.2 as a representative value for the present model and engineering conditions, rather than as a uniquely calibrated field value. This is also broadly consistent with the macro-mesoscopic parameter relationship proposed by Liu (2019): where I is an intermediate variable. Substituting φ = 32° yields μ = 0.21. If μ falls below 0.2, slope failure could extend over a much wider region (Figures 13e,f), which might represent more extreme conditions such as the combination of sustained intense rainfall, drainage system failure, and even seismic events. In addition, insufficient compaction could somehow reduce the shear strength, if to an extent that representative μ drops to approximately 0.1. Both the scale of potential sliding and its consequences would escalate, posing a significant threat to adjacent farmland and residential areas.

μ p = − 2 2 + 2 I 2 + 2 I

I = 1 + tan ⁡ φ 2 + tan ⁡ φ 2

Although the current design provides an overall stability margin, localized risks may persist under extreme conditions. The numerical results suggest that risk mitigation should primarily focus on improving surface and subsurface drainage, limiting rainfall infiltration, and strengthening monitoring of the locally critical areas identified by the LEM and DEM analyses. Erosion control and ecological restoration may also be considered as supplementary measures. These recommendations are directional rather than site-specific design prescriptions, and their detailed configurations and effectiveness should be evaluated through dedicated engineering design and further numerical or field validation.

6 Conclusion

This study integrates three-dimensional limit equilibrium analysis and discrete element modeling to comprehensively evaluate the pre-failure stability and post-failure mobility characteristics of a specific valley-type waste dump. The two approaches characterize the stability condition of the waste body and the subsequent runout process, respectively, while their results are further correlated based on the spatial correspondence of critical failure regions.

Under the prescribed site geometry, material properties, and hydrological conditions, the calculated factor of safety (FS) values for all five analyzed regions exceeded 1.0. The two-dimensional limit equilibrium analyses identified failure mechanism III as the most critical scenario. Under the 50-year return period rainfall condition, the minimum FS occurred at cross-section 2-2, with a value of 1.239. The three-dimensional limit equilibrium analysis also identified the corresponding slip surface B as the most critical region, with an FS of 1.698, and all 3D FS values were higher than their corresponding 2D results. Under the fully saturated condition represented by sustained heavy rainfall combined with drainage system failure, the minimum FS values from the 2D and 3D analyses decreased to 1.06 and 1.15, respectively. The relatively higher FS values obtained from the 3D analysis are likely attributable to the multidirectional lateral confinement provided by the meandering gully and surrounding ridges. This indicates that a single 2D cross-section may not adequately capture the three-dimensional confinement effects induced by complex topographic conditions.

The DEM simulations demonstrated that, with friction coefficients (μ) of 0.10, 0.15, and 0.20, the maximum runout distances of the waste body were 680, 667, and 80 m, respectively, while the corresponding stabilization times were approximately 3,000, 600, and 120 s. These results indicate that μ has a significant influence on the mobility range and duration of post-failure movement under the simulated conditions. At μ = 0.10, the waste body exhibited a two-stage failure process, in which secondary failure mainly increased the mobilized volume and deposition extent but did not significantly extend the maximum runout distance. At μ = 0.20, the mobility range and potential impact area were confined to the proximal slope area. Furthermore, the critical region identified by the limit equilibrium analysis showed a good spatial correspondence with the localized initial mobilization zone obtained from DEM simulations at μ = 0.20, highlighting the complementary capabilities of the two approaches in identifying potential failure-prone areas.

The results indicate that the studied waste dump maintains a certain stability margin under the conventional conditions considered in this study. However, under extreme hydrological conditions involving sustained heavy rainfall and drainage system failure, the FS decreases substantially, and localized areas may approach a critical state. Meanwhile, DEM simulations suggest that reduced μ values may result in longer runout distances and larger impact zones. Therefore, maintaining the long-term effectiveness of drainage systems is essential, and downstream potential impact areas should be incorporated into risk management and emergency planning based on the predicted runout pathways under low-μ conditions. It should be noted that the conclusions of this study are influenced by the specific site geometry, material properties, constitutive models, and loading scenarios considered herein. Therefore, the quantitative results should not be directly generalized to other waste dumps, particularly in cases where fine-scale hydro-mechanical coupling effects and fluid-particle interactions must be explicitly accounted for. Nevertheless, the proposed LEM–DEM integrated framework provides a methodological approach for simultaneously evaluating pre-failure stability and post-failure mobility, and its applicability to other valley-type waste dumps require further validation under different engineering conditions.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Hans-Balder Havenith, University of Liège, Belgium

Reviewed by: Xiaochao Li, China Institute of Water Resources and Hydropower Research, China

Ahmed Chetti, University Center of El-Bayadh, Algeria