Energy and structural evolution process of high-altitude and long-runout landslides induced by a strong earthquake

Restoring the dynamic evolution and energy transfer mechanisms of ancient, earthquake-induced long-runout landslides remains challenging. This study investigates the fluidization and momentum transfer mechanisms of the Mogangling landslide. We reconstructed the high-precision pre-failure topography using a contour-continuity restoration method to accurately define the initial energy state. A three-dimensional discrete element code was then employed to simulate the landslide kinematics over a duration of 87 s. During this process, the sliding mass achieved a maximum velocity of 27.5 m s⁻¹. To elucidate internal momentum transfer, we quantitatively identified effective block collisions by applying an optimal prominence threshold of 30 % to velocity-distance curves. Furthermore, the Alpha Shape algorithm was utilized to extract the continuous structural evolution of the landslide mass. The simulation demonstrates that while the global mechanical energy of the system conforms to the law of conservation, energy at the individual block scale is non-conserved. Intense internal collisions facilitate a pushing effect, transferring kinetic energy from the rear to the frontal mass, thereby sustaining hyper-mobility and prolonging the runout distance. Crucially, this study quantifies the critical thresholds marking the transition from solid-phase sliding to granular-phase flow. The results indicate that macroscopic fluidization and maximum kinetic energy variation occur when the Volume Swelling (VS) rate reaches 29.46 % and the Area Growth (AG) rate reaches 319.59 %. These findings provide a quantitative link between microscopic energy dissipation and macroscopic structural fragmentation, offering vital dynamic parameters for hazard-chain modeling.

1 Introduction

Landslide is a common and widespread geological disaster, which causes great loss of life and property. In recent years, the chain disasters of high-level rock landslides in high-intensity earthquake areas have shown a trend of frequent occurrence (Carpignano et al., 2009; Fan et al., 2019), such as donghekou landslide debris flow dammed lake induced by the 2008 Wenchuan earthquake (Sun et al., 2011), Hongshiyan landslide induced by the 2014 Ludian earthquake (Hu et al., 2017). Strong earthquakes in high-altitude canyon regions often trigger catastrophic long-runout landslides. These geological hazards are characterized by high-speed ejection, dynamic fragmentation, and hyper-mobility, posing severe threats to downstream communities and critical infrastructure like hydropower stations. The 1786 Kangding earthquake (M = 7.75), which triggered the massive Mogangling landslide and subsequently blocked the Dadu River, serves as a classic example of such a hazard chain. Understanding the dynamic evolution and energy mechanisms of such paleo-landslides is crucial for future risk assessment.

At present, in terms of surface terrain restoration, the terrain restoration before a landslide is less. For landslides that have not occurred too far away, the key information of landslide characteristics is often obtained through photos, field surveys, optical remote sensing images and UAV impact, so as to further build geomorphic maps through GIS (Geographic Information System) for subsequent analysis (Zhang et al., 2024; Song at al.2017). For most numerical simulations involving terrain restoration, a two-dimensional simulation is more often used. The simplified model is adopted to simulate the overall process of landslide by inputting external stress (Du et al., 2021; Cheng and Yu-Sheng, 2014; Li et al., 2024; Bernardie et al., 2021; Dahal et al., 2022). This paper uses the method of restoring the continuity of contour lines to establish the terrain data before a landslide, which can provide a new idea of terrain restoration before a landslide for the same type of research.

Existing studies have extensively investigated landslide kinematics using field investigations, flume experiments, and numerical simulations (Huang et al., 2012; Intrieri et al., 2017). Continuum mechanics-based methods (Zhu et al., 2010; Zhang et al., 2020a, b) have successfully reproduced runout distances and deposition areas. However, a critical scientific gap remains: the quantitative link between microscopic energy evolution and macroscopic structural changes is not yet fully understood. While it is widely accepted that fragmentation induces fluidization (dynamic disintegration), traditional methods struggle to quantify the real-time evolution of the landslide's structure (e.g., volume expansion and specific surface area growth) during the sliding process. Furthermore, for ancient landslides, the pre-failure topography is often missing. Most back-analyses rely on simplified slope assumptions or current topography, which fails to accurately represent the initial gravitational potential energy, leading to significant uncertainties in dynamic modeling.

To address these challenges, this study proposes an integrated methodological framework combining topographic reconstruction, discrete element modeling, and morphological quantification. First, unlike simplified geometric assumptions, we apply the contour-continuity restoration method to reconstruct the high-precision pre-failure topography of the Mogangling landslide, ensuring the accuracy of the initial energy state. Second, we utilize the 3D Discrete Element Code (3DEC) to simulate the discontinuous motion, allowing for the extraction of energy evolution data at both the global and individual block scales. Third, and most notably, we introduce the Alpha Shape algorithm to quantitatively monitor the structural evolution of the landslide mass. By calculating the time-varying Volume Swelling (VS) and Area Growth (AG) rates, we aim to establish a quantitative relationship between energy dissipation and structural fragmentation.

This study focuses on the Mogangling landslide to achieve the following objectives: (1) to reconstruct the dynamic fragmentation process driven by seismic loading; (2) to reveal the mechanism of momentum transfer via effective collisions; and (3) to quantify the critical thresholds of structural indices (VS and AG) that mark the transition from solid-phase sliding to granular-phase flow. The findings provide new theoretical insights into the fluidization mechanism of rock avalanches and offer dynamic source parameters for hazard-chain modeling in the Dadu River basin.

2 Study Area and Geological Backgrounds

As shown in Fig. 1, the Mogangling landslide is located at the upstream of Caihong bridge, Jinguang village, Detuo Township, Luding County, Sichuan Province, and the right bank of Dadu River is close to the mouth of Moxi River, a tributary. The coordinates of the landslide center are: $\mathrm{29}\mathit{°}{\mathrm{37}}^{\prime }{\mathrm{30.726}}^{\prime \prime }$ north latitude, $\mathrm{102}\mathit{°}{\mathrm{09}}^{\prime }{\mathrm{41.190}}^{\prime \prime }$ east longitude. It connects Luding County to the north, Shimian County to the south, Moxi town, and Kangding County to the West. The middle reaches of the Dadu River are located at the intersection of “Y”-shaped faults composed of the Xianshuihe fault, Longmenshan fault, and Anninghe fault. It is an area of high seismic intensity, and nearly 20 large and super large earthquake landslides are densely distributed here (from the Detuo Guza river section). The lithology of the landslide area is mainly gray white, strongly to weakly weathered Jinning plagioclase granite (${\mathit{\gamma }}_{{\mathrm{O}}_{\mathrm{2}}}$). Quartz diorite (O₂) and a small amount of basic magmatic rocks are distributed on the western slope of Mogangling. The sand slate of the Triassic Baiguowan formation (T_3bg) is exposed on the left bank of the Dadu River. The quaternary system mainly includes landslide accumulation ($Q_{\mathrm{4}}^{\text{del}}$), colluvial Deluvial deposit ($Q_{\mathrm{4}}^{\text{col}+\text{dl}}$) and alluvial deposit ($Q_{\mathrm{4}}^{\text{al}}$), of which colluvial Deluvial deposit is located at the rear edge of the landslide, the mountain side of the left bank of Dadu River and the surrounding steep slope, and alluvial deposit is located at the river valley terrace. Historical earthquakes have caused extremely heavy losses to southwest China, especially the series of losses caused by earthquakes, landslides, river blocking, and a breach super flood disaster chain. The Mogangling landslide induced by the Kangding earthquake (M_s 7 3/4) in 1786 burst the Dadu River after 9 d. The sudden flood made the downstream area a vast ocean, and tens of thousands of people drowned.

Figure 1 Position of the study area: (a) the location of the Mogangling landslide; (b) remote sensing image of Mogangling; (c) engineering geological map and (d) engineering geological profile.

In the construction of major projects in Western China, the effect of seismic force on slope instability has been underestimated to a certain extent. For example, the stability investigation of the reservoir slope of large-scale power stations in Western China mostly focuses on the stability of existing landslides and the evaluation of the surging waves after instability, and the analysis of the instability probability of potentially unstable high and steep slopes under the action of seismic force and the evaluation and analysis of dam break or dam overturning caused by abnormal bedrock landslide such as surging waves are relatively few. The Dadu River high-intensity area is a major hydropower project in China. There are a large number of paleoseismic landslides in this area. When the development background and laws of these major earthquake landslides have not been fully studied in detail, the construction of high-risk reservoirs and dams has great potential risks. Therefore, taking the Mogangling landslide as an analysis case, exploring the dynamic propagation characteristics can provide an evaluation basis for the potential landslide in the high-intensity area of the Dadu River, which is of great significance to the development of Western China.

3 Methodologies

In order to describe the energy transfer process of earthquake-induced landslide, the representative Mogangling landslide is selected as a case study. As shown in Fig. 2, the research method involves three parts: (1) surface terrain restoration: the terrain before the landslide is deduced, and the landslide model is established based on the existing terrain of the Mogangling landslide by restoring the continuity of contour lines. It includes strata, geological structure, boundary conditions, seismic force input and mechanical parameters; (2) Analysis of Landslide Evolution Process: carry out discrete element numerical simulation of the landslide model, infer the landslide evolution process according to the changes of motion parameters and stress in the numerical simulation process, and add effective collision extraction to analyze the motion law; (3) Energy evolution analysis: the potential energy, kinetic energy and dissipation energy of the landslide are extracted by monitoring the speed and displacement of the whole landslide and some individuals, and their laws are observed. At the same time, through the comparison of effective collision and motion trajectory, the influence of effective collision frequency of different parts of landslide mass on energy transfer is analyzed; (4) Extraction of volume and area of landslide accumulation area: the envelope of landslide mass is generated by alpha shape, so the volume expansion rate and surface area growth rate of landslide mass in the process of landslide can be calculated, and the evolution law can be obtained by combining the energy change.

3.1 Geo-mechanical model building

In order to study the energy transfer in the process of landslide, numerical simulation technology is used in this study. Before model construction, it is necessary to establish a geometric model of the ground surface to obtain the sliding bed and sliding mass. According to the DEM data obtained by UAV scanning, the point coordinates are extracted, and the landslide area is delimited for precision screening to reduce the point density. Convert the geographical coordinate system to the projection coordinate system to facilitate the export of the filtered point cloud. After reducing the noise caused by vegetation and repairing the cavity, the complete topographic surface data of the Mogangling area after the landslide were obtained. The obtained terrain data is corrected by the normal vector to generate the existing model as the sliding bed. Because the landslide occurred too long ago, only the Tiezhuang Temple tablet unearthed in the 1980s has records of earthquake landslide, river dam break, and water discharge. Therefore, this study uses the method of restoring contour lines to simulate the terrain before the landslide. The scanned terrain data is output as contour lines, and then the contour lines corresponding to the sliding bed and accumulation area are modified to be straight. The restored contour lines are transformed into point clouds for subsequent processing, and the topographic map before the landslide is obtained after noise reduction, cavity repair, normal correction, and other processing.

Figure 2 Flow chart of the proposed method of this study.

As shown in Fig. 3a, the geometric model of the terrain surface is established according to DEM. The whole landslide occurrence area can be divided into four parts: the rock wall area, colluvial slope area, slope debris flow area, and main accumulation area. There is a reverse fault along the right bank of the Dadu River, namely the Dadu River fault (Detuo fault). According to the reconstructed topographical surface geometry, the Mogangling landslide is a typical armchair-shaped terrain, with a main sliding direction of 75°, distributed along the NEE-SWW direction, about 450 m long, about 1000 m wide along the river, and a planar area of 0.45 km². The average thickness of the landslide mass is about 100 m, and the volume is about 2.4×10⁷ m³. The elevation of the slope toe is 1120 m (about 5 m higher than the river surface), and 1120–1330 m is the front edge of the landslide accumulation mass, with a slope of about 50°; the elevation of landslide platform is between 1330 and 1380 m, and the slope is about 5–12°; the elevation of the rear wall of the landslide is between 1380 and 1890 m, of which 1380–1600 m are covered by colluvial deposits, with a slope of about 35°, and 1600–1890 m are landslide walls, with a slope of about 57°; the elevation of the back wall ridge is about 2000 m. As shown in Fig. 3b, the contour line is modified along the extension direction of the mountains on the left and right sides to repair the sliding bed generated by the landslide into the terrain of the mountains. The modified contour is converted into a point cloud and repaired. As shown in Fig. 3c, make a section along $a-a^{\prime }$ in Fig. 1c. It can be seen that although compared with the real terrain, the restored terrain is smooth and the surface does not fluctuate, the mountains become continuous, which can be regarded as the terrain before the landslide.

Figure 3 Landslide restoration diagram: (a) 3D views of the study area; (b) comparison map of surface topography before and after restoration, and (c) terrain restoration comparison map.

Due to the need to simulate the large deformation and displacement of rock mass in the process of landslide, the three-dimensional discrete element program (3DEC), one of the most widely used numerical methods in rock mechanics, is used to simulate the motion process of the Mogangling landslide. According to the extracted terrain data after the landslide, the three-dimensional modeling is carried out, and the landslide mass is constructed according to the restored terrain surface.

As shown in Fig. 4, the established model includes the landslide mass and base block. The length of the landslide mass in the sliding direction is 1413 m, and the vertical height is about 871 m; the bottom of the landslide mass is 797 m away from the Dadu River, and the vertical height is about 491 m. Because the occurrence distribution data of the structural plane of the landslide mass before the landslide can not be obtained, the impact of the structural plane on the landslide process has to be ignored in this study. The landslide mass is evenly divided into 8.7 m×8.7 m blocks on the xoy surface. The height is determined by the vertical height from the surface of the landslide mass to the sliding bed. A uniform structural plane is set inside. According to the geological engineering drawing, different strata are generated and marked with different colors. The whole model contains 25 345 blocks, 2788 m long, 3341 m wide, 1121 m high, and the volume of the landslide is 2.4×10⁷ m³.

Note that, since the Mogangling landslide was triggered by the 1786 Kangding earthquake, the considerable time elapsed since the event and the lack of instrumental ground-motion records for this historical earthquake make it impossible to directly obtain the seismic input parameters from the original incident. This study adopts the acceleration records measured at the Luding station during the Wenchuan earthquake as the input ground motion for the dynamic numerical simulation analysis of the Mogangling landslide. The justification for this substitution is threefold: (1) Source Mechanism: both the 1786 Kangding earthquake (Xianshuihe Fault) and the 2008 Wenchuan earthquake (Longmenshan Fault) occurred within the convergent tectonic system of the eastern Tibetan Plateau. They share similar crustal stress environments within the Y-shaped fault zone, resulting in comparable high-energy release characteristics. (2) Propagation Path and Site Effects: although the epicentral distances differ, the recording station (Luding) and the Mogangling landslide are situated in the same Dadu River canyon corridor. Seismic waves propagating through this specific topography undergo similar scattering and topographic amplification effects. Using the Luding station record inherently preserves these critical local valley-site response characteristics that synthetic waves or records from other regions would miss. (3) Damage Potential: previous studies indicate that the damage distribution and intensity of the Kangding earthquake are comparable to the effects of the Wenchuan earthquake in this region, making the selected record a representative proxy for the seismic loading experienced by the Mogangling slope.

Figure 4 Discrete numerical model with fine topography: (a) Stratigraphic model; (b) Structural surface model.

In the course of a discrete-element simulation-particularly when modeling nominally semi-infinite bodies such as the use of fixed or elastic boundaries in a static analysis will inevitably reflect outbound seismic waves into the model during a subsequent dynamic analysis. Consequently, when employing 3DEC for seismic-dynamic modeling, it is essential to impose non-reflecting viscous boundaries in conjunction with free-field boundaries. Specifically, viscous boundaries are applied at the model base where the incident seismic waves enter, while free-field boundaries constrain the lateral faces that are not subjected to seismic input. Dynamic loading must be introduced only after the model has attained static equilibrium; i.e., the dynamic analysis is predicated on a prior static analysis that has reached convergence. For the static phase, boundary conditions consist of vertical-direction velocity restraints at the base and horizontal-direction velocity restraints along the lateral faces.

In the numerical simulations, the blocks were assigned a Mohr–Coulomb constitutive model. The physical and mechanical properties of both the rock mass and the discontinuities critically influence the deformation-failure mechanisms of the rock slope and the evolution of key response variables such as displacement, velocity, and acceleration. In this study, the selection of physical-mechanical parameters was based primarily on laboratory testing of site samples, empirical analogy from similar engineering projects in the Dadu River basin (Zhao et al., 2021), and parameter inversion based on the landslide's final deposition. The slope is underlain by plagioclase granite, which is subdivided into strongly weathered and moderately weathered zones. The discontinuities are represented by both fully developed potential rupture planes and joint sets, all of which are governed by the Mohr–Coulomb failure criterion. In the numerical model, the parameters for the rock mass and discontinuities used in the analyses are those listed in Table 1.

Table 1 Physical and mechanical parameters of the rock mass and discontinuities (parameters derived from laboratory tests and back-analysis based on Wu et al., 2013).

3.2 Effective collision identification

During the landslide, numerous block collisions occurred, altering the kinematic behavior of the rock mass. Some impacts produced substantial acceleration of individual blocks, while others caused deceleration due to energy dissipation. To elucidate the anomalous dilation of the sliding mass, this study examines those collisions that markedly increase block velocities-hereafter termed “effective collisions” (Ge et al., 2021).

Effective collisions of rock blocks are identified by their characteristic velocity changes. Figure 5a presents a typical velocity-travel-distance curve for an individual block. A local velocity peak (denoted by red or purple markers) corresponds to an acceleration phase, but only those pronounced peaks (highlighted in red) yield a substantial velocity increase and, consequently, a markedly longer travel distance. Hence, effective collisions can be detected by pinpointing the significant local maxima on the velocity-distance profile.

Figure 5 Effective collisions identification: (a) definition of effective collisions; (b–f) threshold=10 %, 20 %, 30 %, 40 %, and 50 %.

Significant local maxima are defined as points whose velocity markedly exceeds that of their neighbors. To pinpoint these prominent peaks on the velocity-distance curve, a suitable prominence threshold must be chosen-one that quantifies each peak's standout height, width, and position relative to adjacent peaks. This threshold critically influences the detection of effective collisions. As illustrated in Fig. 5b–f, low thresholds (10 % and 20 %) spuriously flag minor undulations – often caused by numerical fluctuations in contact forces rather than physical impacts – as peaks, leading to an overestimation. Conversely, high thresholds (40 % and 50 %) fail to detect significant momentum-transfer events. Consequently, a 30 % threshold was adopted as the optimal value to distinguish effective physical collisions from numerical noise. At a threshold of 30 %, however, the identification of effective collisions is deemed acceptable. Under this setting, 439 effective collisions were detected from the travel-distance versus velocity profile.

3.3 Energy evolution

During a landslide, dynamic energy propagation involves three forms of energy: kinetic energy (KE), potential energy (PE), and dissipated energy (DE), with the sum of KE and PE defined as mechanical energy (ME). For an individual rock block within the sliding mass (i.e., a single block in the numerical model), at time t_i, its mechanical energy is given by:

$$\begin{array}{}\text{(1)}& \text{ME}\left(t_i\right)=\text{KE}\left(t_i\right)+\text{PE}\left(t_i\right)\end{array}$$

where, at time t_i, PE(t_i), and DE(t_i) can be expressed as:

$$\begin{array}{}\text{(2)}& {\displaystyle }\text{KE}\left(t_i\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}m\left[v\right(t_i){]}^{\mathrm{2}}\text{(3)}& {\displaystyle }\text{PE}\left(t_i\right)=mgh\left(t_i\right)\text{(4)}& {\displaystyle }\text{DE}\left(t_i\right)=\text{ME}\left(t_{\mathrm{0}}\right)-\text{ME}\left(t_i\right)\end{array}$$

In these expressions, m is the mass of the monitored block (kg); v(t_i) is the block's velocity at time t_i (m s⁻¹); h(t_i) is the block's elevation at time t_i (m); ME(t₀) is the mechanical energy of the block in its initial state (J). Since the initial velocity is 0, ME(t₀) reduces to the initial potential energy. Because this study focuses exclusively on energy transfer during the landslide, the potential-energy calculation uses the model's initial height as the reference datum rather than true sea level; only the change in potential energy due to relative displacement is considered. For the landslide system as a whole, the riverbed of the Dadu River is selected as the datum plane for potential energy calculation. For specific monitoring blocks, the minimum detected vertical distance is used as the datum.

As shown in Fig. 6, ten monitoring points were installed in the lower, lower-middle, upper-middle, and upper portions of the slide mass to capture representative kinematics and deposition characteristics of the landslide blocks. In addition, the maximum unbalanced force was tracked to verify that the model reached post-failure stability. The properties of each monitoring point are listed in Table 2.

Figure 6 Location of 10 monitoring points.

Table 2 The basic information of blocks where the 10 monitoring points are located.

For the global analysis, the overall displacement, velocity, and stress variations of the entire slide mass were recorded and interpreted in conjunction with the stages of the failure process. For the local analysis, displacement, velocity, and stress histories were extracted for each of the ten monitored blocks in all relevant directions and used, together with the governing equations, to compute the evolution of their energies. Following these energy-evolution calculations and kinematic observations, the sequence of failure and run-out at Mogangling was inferred by integrating the numerical results with field-investigation findings.

3.4 Structural extraction of sliding body

To examine the relationship between the dispersal degree of the landslide deposition zone and the evolution of energy during the sliding process, we quantify dispersion by extracting the volume and surface area of the deposition zone. An alpha-shape algorithm is applied to the slide mass to generate an enclosing envelope, from which its volume and area are computed. To facilitate a direct comparison between the changes in volume and area and the corresponding energy evolution, two dimensionless indices are introduced: volume swelling (VS) and area growth (AG):

$$\begin{array}{}\text{(5)}& {\displaystyle }\text{VS}={\displaystyle \frac{V_t-V_{\mathrm{0}}}{V_{\mathrm{0}}}}\times \mathrm{100}\mathit{\%}\text{(6)}& {\displaystyle }\text{AG}={\displaystyle \frac{S_t-S_{\mathrm{0}}}{S_{\mathrm{0}}}}\times \mathrm{100}\mathit{\%}\end{array}$$

In these expressions, V₀ denotes the initial volume of the landslide mass (m³); V_t denotes the volume of the enclosing envelope at time t (m³); S₀ denotes the initial surface area of the landslide mass (m²); S_t denotes the surface area of the enclosing envelope at time t (m²).

Considering that the choice of the alpha parameter (α) in the Alpha Shape algorithm directly affects the accuracy of the envelope fitting as well as the computed volume and surface area, this study uniformly sets the α value to 1. This approach ensures a tight fit of the envelope model while minimizing geometric deviations introduced by parameter tuning.

Building on this foundation, to further explore the relationship between rock mass fragmentation and energy transfer during the landslide process, the evolution of the landslide system's kinetic energy (KE), potential energy (PE), and dissipated energy (DE) concerning the dimensionless indices volume swelling (VS) and area growth (AG) is visualized using line plots. Specifically, the horizontal axis represents either VS or AG, while the vertical axis shows the instantaneous values of the three energy components (KE, PE, DE). Each energy curve is distinguished by a unique marker and line style to ensure clarity and facilitate visual comparison. By connecting discrete data points at successive time steps, a continuous evolution trajectory is formed, allowing dynamic trends in energy conversion and dissipation to be observed across different levels of mass dispersion. This visualization approach not only intuitively illustrates the influence of rock fragmentation and dispersal on energy distribution but also provides a quantitative basis for understanding energy transfer mechanisms in the dynamics of landslides.

4 Results

4.1 Motion and deposition of the Mogangling landslide

The numerical simulation successfully reproduced the entire process of the Mogangling landslide, which lasted for approximately 87 s. Figure 7 illustrates the kinematic evolution through velocity and stress distribution cloud maps.

• a.

  Initial State (0 s). As shown in Fig. 7a and h, the slope is initially stable with zero velocity before the seismic input.

• b.

  Vibration and Tension-cracking Stage (0–3 s). As shown in Fig. 7b and i, the maximum velocity reaches 6.1 m s⁻¹. During this phase, tensile stress concentration is observed along the steeply dipping structural planes at the rear of the slope. Cracks begin to develop along the pre-existing boundaries, although the main sliding mass has not yet detached.

• c.

  High-speed Initiation Stage (3–7 s). As shown in Fig. 7c and j, the landslide mass accelerates rapidly, with the maximum velocity reaching 8.2 m s⁻¹. Continuous tensile cracks form along the rear and lateral margins, effectively separating the landslide mass from the stable bedrock. The stress monitoring indicates significant shear failure developing along the sliding zone.

• d.

  Collapse and Ejection Stage (7–9 s). As shown in Fig. 7d and k, the maximum velocity peaks at 21.2 m s⁻¹. The landslide mass disintegrates into fragmented blocks. The frontal part of the slide collides with the lower bedrock slope, while the middle and rear portions collapse downslope, losing their original structural coherence.

• e.

  Impact and Deposition Stage (9–15 s). As shown in Fig. 7e and l, the maximum velocity reaches 27.5 m s⁻¹. The entire sliding mass impacts the valley floor. The simulation shows that the landslide material spreads laterally upon reaching the Dadu River valley.

• f.

  Accumulation Stage (15–87 s). As shown in Fig. 7f and n, the velocity gradually decreases. By t=87 s, the velocity approaches 0 m s⁻¹. The landslide material completely blocks the river channel, forming a dam-like accumulation.

Figure 7 Simulation results of Mogangling landslide: (a–g) cloud map of velocity when t=0, 3, 7, 9, 12, 15 and 87 s; (h–n) cloud map of stress when t=0, 3, 7, 9, 12, 15 and 87 s.

4.2 Energy evolution of the landslide

4.2.1 Energy evolution for the whole debris mass

As shown in Fig. 8, the energy evolution of the entire sliding mass can be divided into distinct phases. Initially (0–3 s), the potential energy (PE) dominates with minimal kinetic energy (KE). From 3–12 s, the landslide initiates rapidly, with KE increasing steadily to its peak. Around 13 s, corresponding to the impact with the riverbed, the KE drops sharply. Subsequently, the DE (Dissipated Energy) continuously increases while PE decreases until the mass comes to rest.

Figure 8 Evolution of energy during the landslide.

4.2.2 Energy evolution for individual rock blocks

The energy curves for individual blocks (Fig. 9) exhibit localized kinematic characteristics that differ significantly from the smoothed global average. Following the data correction, the mechanical energy of discrete blocks demonstrates realistic fluctuations driven by terrain variations and internal collisions. The localized dissipated energy curves display a characteristic staircase accumulation pattern. Rather than experiencing unphysical decreases, the dissipated energy surges rapidly at specific moments. These sudden leaps in energy dissipation coincide precisely with the effective collisions during the movement, indicating that the kinetic energy of individual blocks is violently consumed through localized impacts, fragmentation, and intense basal friction.

Figure 9 Energy analysis for 10 individual blocks in Mogangling landslide: (a–j) energy evolution of point 1, 2, 3, 4, 5, 6, 7, 8, 9, and 10.

Figure 10a shows the trajectories of the 10 monitored blocks, while Fig. 10b presents the frequency of effective collisions experienced by blocks at different locations during the landslide process.Combining this with Fig. 10b, it is clear that as the migration distance increases during the landslide, the number of collisions also rises continuously, reaching a peak once the blocks enter the river channel. At this stage, the frequency at which individual blocks receive energy transfer from surrounding blocks also increases.

Figure 10 (a) Movement trails of 10 monitoring points and (b) the plot of travel distance vs. collision frequency.

4.3 Structural evolution as a function of energy changes

Figure 11 illustrates the Alpha Shape envelopes of the deposit mass at 10 representative time points. As the simulation progresses, both the envelope volume and surface area exhibit a continuous increasing trend.

Figure 11 Comparison of the overall envelope fitting of the landslide at 10 different times and the shape of the landslide mass in numerical simulation: (a–h) when time=0, 3, 7, 9, 12, 15, 32 and 87 s.

Figure 12 presents the kinetic energy (KE), potential energy (PE), and dissipated energy (DE) as functions of the dimensionless indices: Volume Swelling (VS) and Area Growth (AG). The data indicates a critical threshold behavior, signifying a fundamental transition in the mechanical state of the landslide system. When the variation amplitude of KE reaches its maximum, the Volume Swelling rate (VS) reaches approximately 29.46 %, and the Area Growth rate (AG) reaches about 319.59 %. This phenomenon mirrors the mechanical evolution of granular flows described by Takahashi (1978), where internal shear and dispersive stresses facilitate the transition from a locked static assembly to a mobilized fluid-like state. Simultaneously, the Dissipated Energy (DE) exhibits a sharp increase corresponding to these geometric expansion thresholds. This surge in energy dissipation reflects the complex solid-fluid interactions and the shift in energy partitioning that characterize the physics of debris flows (Iverson, 1997).

Figure 12 Overall energy variation of the sliding mass: (a) as a function of volume; (b) as a function of area.

5 Discussion

5.1 Dynamic evolution and energy transfer mechanism

The numerical simulation reveals that the Mogangling landslide underwent a complex dynamic process involving seismic triggering, high-speed ejection, and collision-induced disintegration. The initiation phase was governed by the dynamic response of the slope to seismic loading. As observed in the early stages of the simulation, the amplification of seismic waves at higher elevations, combined with the “back-slope effect,” significantly increased the Peak Ground Acceleration (PGA) at the slope crest. This dynamic stress concentration triggered intense tensile cracking along the rear structural planes (Fig. 7), a phenomenon consistent with the failure mechanism of the Donghekou landslide triggered by the 2008 Wenchuan earthquake (Sun et al., 2011; Wang et al., 2014).

A critical finding of this study is the mechanism of long-runout propagation driven by energy transfer. While the mechanical energy of the global landslide system conforms to the law of conservation (Fig. 8), the energy evolution at the individual block scale is non-conserved and highly heterogeneous. Our monitoring data (Fig. 9) and effective collision analysis (Fig. 10) demonstrate that blocks in the rear portion of the slide mass experienced a sharp increase in dissipated energy upon impacting the frontal blocks. Crucially, this collision process facilitated a significant transfer of momentum and kinetic energy from the rear to the front. This “pushing effect” effectively energized the frontal rock mass, enabling it to maintain high velocities and achieve a runout distance that exceeds frictional predictions. This observation quantitatively supports the theory of momentum transfer in rock avalanches (Heim, 1932; Davies, 1982), explaining why the frontal debris of the Mogangling landslide could travel extensively into the valley.

5.2 Structural fragmentation and solid-to-fluid phase transition

The coupling between energy dissipation and structural evolution offers a novel perspective on the fluidization mechanism of rock avalanches. Traditional studies often describe disintegration qualitatively; however, our application of the Alpha Shape algorithm allows for a quantitative definition of the transition from a solid-phase slide to a fluid-phase flow.

As shown in Fig. 12, there is a distinct synchronization between the peak kinetic energy and the geometric expansion of the landslide mass. Specifically, when the Volume Swelling rate (VS) surpasses a critical threshold of approximately 29.46 % and the Area Growth rate (AG) exceeds 319.59 %, the landslide exhibits behavior characteristic of a granular fluid.

The fundamental significance of the observed VS and AG indices lies in their ability to provide a quantitative definition for the evolutionary transition of the landslide from a coherent solid mass to a dispersed granular flow. This physical transformation is schematically illustrated in Fig. 13. As depicted in Fig. 16a, the initial rock mass behaves as a large, coherent block. In this phase, the calculated VS and AG values remain low (below the critical thresholds), indicating minimal internal fragmentation and volume expansion. The motion is primarily governed by basal sliding friction with high resistance due to block interlocking. The identified thresholds (VS>29.46 % and AG>319.59 %) mark the critical tipping point of macroscopic disintegration and dilatancy. Crossing these values signifies that the internal structure has been sufficiently disrupted to overcome interlocking. Once these thresholds are surpassed, the landslide evolves into the state shown in Fig. 13b, characterized by a highly fragmented, dispersed granular mass. The consistently high VS and AG values quantify this state, where significant dilation reduces inter-particle contacts, allowing collisional rheology to dominate and enabling long-runout propagation.

Figure 13 Schematic illustration of the solid-to-granular phase transition: (a) the coherent solid phase; (b) the dispersed granular-fluidized phase.

5.3 Uncertainties and limitations of the numerical modeling

While the discrete element simulation successfully reproduces the macroscopic behavior of the Mogangling landslide, several uncertainties and limitations associated with the modeling choices must be acknowledged.First, regarding the pre-failure topography, the contour-restoration method inevitably introduces a smoothing effect, potentially filtering out micro-geomorphic features that contribute to basal friction. Nevertheless, considering the massive volume of the landslide (2.4×10⁷ m³), the global kinematics and runout path are controlled primarily by the valley-scale topography rather than local surface roughness.Second, the use of the seismic record from the Luding station (2008 Wenchuan earthquake) as a proxy for the 1786 Kangding event introduces uncertainty regarding the exact spectral content. However, this substitution is justified by the tectonic similarity (both events occurred within the Xianshuihe-Longmenshan fault system) and, crucially, the preservation of the Dadu River canyon's specific valley-site effects in the recording, which are vital for accurate slope response analysis. Third, the discretization strategy in 3DEC involved simplified structural planes and a relatively uniform block size distribution due to computational constraints. This simplification restricts the model's ability to simulate the continuous generation of fine-grained gouge and the complex interaction of varying block sizes. Furthermore, as a distinct element code, 3DEC simulates debris-flow-like behavior through the interaction of dry blocks rather than fluid dynamics equations such as involving pore water pressure. Therefore, the fluidization discussed in this study refers to the mechanical regime transition from a coherent solid mass to a granular flow state driven by collisions, rather than a saturated slurry flow. Despite these constraints, the extracted energy evolution trends and the identification of effective collisions remain robust indicators of the landslide's dynamic disintegration process.

5.4 Implications for hazard-chain evolution

Unlike gravitational landslides, earthquake-induced landslides are driven by strong external dynamic loading. This study utilizes the dual innovations of discrete element energy extraction and structural morphological analysis to quantitatively characterize how this loading alters the landslide's motion state. (1) Quantitative Characterization of Motion State: the strong seismic force is not merely a trigger but the driver of the energy-structure evolution. As captured by our energy monitoring, the seismic input is rapidly converted into the kinetic energy of the rock mass. This energy surge drives intense internal collisions, leading to the structural disintegration quantified by our Alpha Shape indices. The crossing of the VS and AG thresholds (as shown in Fig. 13) marks the definitive transition from a friction-dominated solid state to a collision-dominated fluid state. (2) Implications for Hazard Scope and Prediction: this structure energy coupling mechanism has profound implications for predicting the disaster scope. The structural fragmentation directly facilitates the momentum transfer mechanism, which sustains the landslide's hyper-mobility. Consequently, the runout distance and impact area are significantly amplified compared to intact rock slides. For future hazard forecasting, this suggests that prediction models must move beyond simple topography-based frictional laws. They must account for the degree of structural fragmentation and energy dissipation efficiency induced by the earthquake intensity. Understanding this phase transition mechanism is key to accurately delineating the potential danger zones of catastrophic long-runout landslides.

The simulation results confirm that the Mogangling landslide completely blocked the Dadu River, forming a landslide dam with a deposition geometry characterized by a significant upstream protrusion (Fig. 7). Although the post-damming hydraulic processes were not modeled in this study, the kinematic data implies severe hazard-chain consequences. The rapid deposition (velocity dropping to 0–87 s) suggests an instantaneous river blockage, which would lead to a rapid rise in the upstream water level and subsequent backflow flooding. Historical records confirm that the breach of this dam 9 d later caused a catastrophic flood downstream. For the ongoing hydropower development in the Dadu River basin, these findings underscore the critical need to consider the “strong earthquake–landslide–impulse waves–damming–outburst flood” hazard-chain. Future risk assessments should integrate the dynamic parameters extracted in this study into hydraulic models to predict potential dam-breach floods and optimize the siting of dams and reservoirs in high-intensity seismic zones.

6 Conclusions

This study investigated the dynamic evolution and energy transfer mechanisms of the Mogangling high-altitude long-runout landslide triggered by the 1786 Kangding earthquake. By integrating pre-landslide topographic reconstruction based on contour continuity, 3D discrete element numerical simulation and Alpha Shape structural quantification, we successfully reproduced the landslide process. This work establishes a quantitative framework linking microscopic block collisions to macroscopic morphological changes, providing a theoretical basis for analyzing energy conversion in complex rock avalanches.

Based on the numerical results and theoretical analysis, the main findings are summarized as follows:

1. The simulation results reveal distinct behaviors between the global and slide mass and individual blocks. For the landslide mass as a whole, the conversion among potential, kinetic, and dissipated energy follows the law of energy conservation; specifically, the dissipated energy increases continuously, and the sum of the three energy components remains equal to the initial potential energy. However, for individual blocks, energy is not conserved. This is due to the intense interaction between blocks, where collisions and compression facilitate energy transfer. This phenomenon is quantitatively reflected by the temporary fluctuations observed in the dissipated energy curves of specific blocks at different stages of the movement.

2. The landslide process involves significant geometric evolution. The collision and friction between blocks during energy transfer lead to rock mass fracturing and disintegration. The generated debris fills the interstitial voids and coats the surface of the deposit, resulting in a sudden expansion of the overall envelope volume and a sharp increase in the specific surface area. These structural changes are strongly correlated with the peaks in energy dissipation, indicating that mechanical work is primarily consumed by the fragmentation and disintegration of the rock mass.

3. Based on the analysis of effective collisions, a pushing effect mechanism was identified. During the collapse and sliding stages, the rear rock mass experiences fewer collisions but retains higher potential energy, which is converted into kinetic energy to impact the frontal blocks. This collision leads to a sharp increase in the dissipated energy of the rear blocks and a transfer of kinetic energy to the front rock mass. This energy transfer mechanism effectively pushes the frontal material forward, enabling the landslide to achieve a longer runout distance.

4. Despite the insights gained, this study has certain limitations. First, the numerical simulation model involves simplifications of the complex geological structures and micro-parameters, which may introduce uncertainties compared to the actual heterogeneous rock mass. Second, the current model primarily focuses on dry granular flow and does not fully account for the coupling effects of groundwater or fluid-solid interactions such as impulse waves, which are critical in specific landslide scenarios. Future research will focus on refining the parameter calibration process and incorporating fluid-mechanical coupling to enhance simulation accuracy.

5. A preliminary restoration method for the ancient landslide topography was established based on contour continuity. For the source area, the trajectory equations of the stable contour lines outside the landslide boundaries were extracted to reconstruct the missing terrain. For the deposit area, where reference features are scarce, the restoration was guided by the geomorphological characteristics of surrounding mountain areas within the same stratigraphic unit. This approach provides a reliable geometric basis for subsequent numerical simulations.

Appendix A: Nomenclature

Table A1 Nomenclature of key parameters used in this study.