Numerical Runout Modeling Analysis of the Loess Landslide at Yining, Xinjiang, China

Abstract

The Panjinbulake loess landslide is located in the western part of the Loess Plateau, in Yining County, Xinjiang, China. It is characterized by its long runout and rapid speed. Based on a field geological survey and laboratory test data, we used the DAN-W dynamic numerical simulation software (Dynamic Analysis Of Landslides, Release 10, O. Hungr Geotechnical Research Inc., West Vancouver, BC, Canada) and multiple sets of rheological models to simulate the whole process of landslide movement. The best rheological groups of the features of the loess landslide process were obtained by applying the Voellmy rheological model in the debris flow area and applying the Frictional rheological model in the sliding source area and accumulation area. We calculated motion features indicating that the landslide movement duration was 22 s, the maximum movement speed was 20.5 m/s, and the average thickness of the accumulation body reached 5.5 m. The total accumulation volume, the initial slide volume and the long runout distance were consistent with the actual situation. In addition, the potential secondary disaster was evaluated. The results show that the DAN-W software and related model parameters can accurately simulate and predict the dynamic hazardous effects of high-speed and long runout landslides. Together, these predictions could help local authorities make the best hazard reduction measures and to promote local development.

1. Introduction

Landslides are among the most destructive geological disasters with features of rapid speed, long runout distance, and entrainment effect [1,2,3]. Catastrophic landslide events are often triggered by heavy rainfall, earthquake, and engineering activities [4,5,6]. According to the spatial characteristics and trajectory of the sliding body, the entire process of landslide movement is mainly divided into three stages: The starting stage at the slide source area, the propagating stage, and the deposition stage [7,8,9,10]. In addition, entrainment, base liquefaction, and air cushioning occur during the landslide movement [11,12,13]. To reduce the landslide hazard loss, risk assessment is often requested [14,15,16]. The landslide runout analysis is a very effective method to assess a landslide hazard [17]. Landslide runout analysis involves two aspects: The simulation of previous landslides and the prediction of potential landslides [18]. Runout analysis could be used to design remedial engineering measures, such as barricades and berms [17]. The maximum runout distance, propagation velocity, and the deposit thickness, and provision of the basis for the design of remedial engineering measures, are obtained by landslide runout analysis.

Landslide runout analysis methods mainly include empirical-statistical methods and numerical models [19,20]. Empirical methods establish the geometric relationship between the landslide volume, height difference, and angle of reach (i.e., the angle of the line between the highest point of the rear edge and the farthest point of the sliding distance) to predict the sliding distance. The empirical methods could not precisely predict the runout distance in the different complex geological environments, including entrainment, friction resistance, and impaction. Compared with the empirical method, the numerical models could give more information about the dynamic features of the sliding mass under different geological environments, such as the scraping depth, the thickness of the accumulation area, velocity, and the scope of dangerous area. The numerical simulation methods include the discrete element methods and the continuum methods. The discrete element method is based on Newton’s second law and is used to analyze the interaction of particles constituting a landslide. It is suitable for landslides with debris flow patterns, such as the MatDEM (i.e., Fast GPU Matrix computing of Discrete Element Method, Nanjing University) and PFC (i.e., Particle Flow Code, Itasca) [21,22]. The continuum method, based on the momentum and mass equations incorporating the earth pressure theory, simulates the motion characteristics of the slip mass to obtain the velocity, position, and thickness of the slip mass [23,24]. The continuum methods have been successfully used to simulate previous hazards and predict potential hazards, such as debris flow, landslides, landslide bam, and avalanches [25,26,27]. The methods evolved into models such as the GeoFlow-SPH [28], LS-RAPID [29], Flow-2D [30], Kinematic model [31] and DAN model [23]. Based on the fluid continuity equation and motion equation, Hungr proposed the landslide dynamics model software DAN-W, which regards the sliding body as an equivalent fluid and can accurately calculate sliding motion characteristics [23]. These models provide a good method for the risk assessment of the geologic hazard.

The loess geological hazards frequently occurred in Tajik and Kazakhstan Tian shan area, which have also become a focus [32]. The loess has quite a widespread distribution in the world and it occupies approximately 10% of the total global land area. China is the country with most widely distributed loess area in the world. Loess is mainly distributed in the northwestern part of China, on the Loess Plateau, which covers an area of nearly 630,000 km², accounting for 4.4% of China’s land area [33,34] (Figure 1a). Due to its special geological structure, loess has high water sensitivity (i.e., loess undergoes a structural collapse when wetted) and is prone to geological hazards. The main types of hazards are loess landslides, such as the Heifangtai Landslide Group and Jingyang Landslide Group [35,36,37,38,39]. Among them, high-speed and long runout loess landslides have caused considerable losses in terms of human lives and property and have become an important research topic. Due to the porosity, weak cementation and water sensitivity of loess as the water content increases, the shear strength of loess declines sharply, and the loess structure is destroyed [40,41,42,43]. In addition, the strength loss in loess might also be a chemo-mechanical problem that involves volume and stress changes of the finest component due to changes in pore water salinity. Consequently, a loess slope loses stability and slides. Furthermore, the pore water pressure rises during the sliding process, and the phenomenon of motion liquefaction occurs, which readily forms a high-speed and long runout landslide [44,45,46,47].

In the study, we studied the characteristics of the Panjinbulake loess landslide through a field geological survey and aerial image analysis using drones. We used the landslide dynamics model DAN-W and multiple sets of rheological models to calculate the dynamic characteristics of this landslide. By using the different rheological models to simulate the different stages of the loess landslide (i.e., triggering in the sliding source area, propagation in the debris flow channel area, and deposition in the accumulation area), the best rheological model groups and parameters were obtained to improve the accuracy to analyze the loess dynamic characteristic. The potential secondary failure of the landslide was evaluated. This study could offer a basis to predict the potential landslide runout distance and define the hazard area, make necessary measures to prevent landslide induced damages (e.g., engineering measures, landslide early warning systems, and emergency response), and to favour local development.

2. Site Overview

The Panjinbulake loess landslide is located in the Karayagaq Township, Yining County, Xinjiang. The coordinates of the central point of the landslide are 81°30′32″ E and 44°11′48″ N, 53 km from downtown Yining City, and 43 km from Sailimu Lake (Figure 1b).

Among the landslide geological hazards in the Piliqinghe basin, most are high speed and long runout loess landslides, which pose a huge hazard to local agricultural and livestock production. Among these landslides, the Kezileisai landslide buried 402 cattle and 5 sheep, and the direct economic losses reached 490,600 Yuan (¥). There were some similar landslides in the Piliqinghe Basin (Figure 2b). The Panjinbulake loess landslide is a typical high speed and long runout loess landslide that occurred recently in the region. Through analysis of the induced factors and dynamic effects of the landslide, this study provides a reference to the landslide dynamic hazards in the loess area.

3. Post-Failure Behavior and Landslide Influential Factors

3.1. Post-Failure Behavior

The triggering of landslides was mainly from snow infiltration, which turn into long runout and rapid landslides, constitutes a typical disaster model in the loess area. The Panjinbulake landslide belongs to this typical geologic hazard model. Following the instability failure of the landslide, the front loess main body slipped due to the river flushing action. Pore-water pressure increased and soil saturation during the sliding process because of the snow infiltration. The landslide was transformed into debris flow, showing a flow state. Then, the landslide volume increased gradually with entrainment effects and flowed into the Pliqinghe Gully. Finally, the landslide struck the opposite side of the mountain and stopped.

3.2. Landslide Influential Factors

According to the field investigation, the occurrence of the Panjinbulake landslide was caused by a combination of factors, including geological structure and formation lithology, topography and hydrogeological condition.

3.2.1. Geological Structure and Formation Lithology

The landslide research area belongs to the western part of the Yili Valley. The terrain is generally high in the north and low in the south. It is gradually inclined from the northeast to the southwest. The elevation of the area is between 620 and 3700 m. It is a block-like eroded and uplifted mountain, covered with gravel and loess layers, showing a low mountain grassland landscape. The landslide area is located on the southwestern side of the West Tianshan Youdi trough fold belt in the southwestern Tianshan fold of the Tianshan-Xing’an trough fold area. It belongs to the junction of the Boroconu Mountain Complex Anticline and Yili Block and is located 2.7 km south of the Nalati deep fault zone (Figure 2a). The rocks in the study area are mainly Ishikirik group tuff, tuff lava, and gray-green coarse sandstone in the Carboniferous system and basalt in the Dahala Junshan Formation of the Carboniferous system. The surface layer is the Quaternary Holocene loess, with well-developed joints (Figure 2b). The structure was relatively loose, and the wormhole, large void structure, belonging to low-plastic silt. The silt (0.075–0.005 mm) content of the loess in the Piliqinghe area is high, reaching 69.8–86.0%; the fine sand (0.25–0.075 mm) content is 3.7–18.0%; and the clay (<0.005 mm) content is 10.3–12.2%. As a result, the loess expands and collapses after encountering water, and it is prone to motion liquefaction under certain static or dynamic water conditions, which provide good source conditions for landslides.

3.2.2. Topography

The Panjinbulake loess landslide was located on the south bank of the river and had long been subjected to the lateral erosion of the river, resulting in good conditions for the landslide front to be in the air. The hillslope was steep, with a slope of 40° (Figure 3a). Corresponding tensile stress condition occurred near the top of the mountain, and the cracks at the trailing edge of the landslide were gradually enlarged (Figure 3b). According to the Google Earth remote sensing image map from 18 May 2013, the front edge of the landslide had slipped. The sliding volume was about 9000 m³ (Figure 3a), and several tensile cracks appeared on the trailing edge. Consequently, a steep ridge of up to 1 m was formed (Figure 3c,d), which provided good topographic conditions for loess deformation and stress relief.

3.2.3. Hydrogeological Condition

The landslide area belongs to the temperate continental semi-arid climate. The average annual precipitation (for the period of 2011–2016) is 330.6 mm. Precipitation is highest from March to July, during which the monthly rainfall exceeds 30 mm, accounting for 52.5% of the annual rainfall (Figure 2c). Snowfall mainly occurs from October to the following March. The snowfall thickness can reach 94 mm per month. The fissure water inside the slope is frozen, and the vertical joints and cracks become enlarged due to the frost heaving action. From mid-March, the temperature rises, the snow that covered the surface begins to melt, infiltrating the cracks and joints and forming a certain transient water pressure and transient saturation zone in the surface layer of the slope. This results in a decrease in the anti-sliding force of the slope, thus inducing landslides. The rising water level of the river also causes the hydraulic gradient inside the slope to drop. Ice and snow meltwater can also be stored in the mountain for a long time and can continue to increase the slope sliding force and accelerate the formation of the potential slip surface. In addition, fissures are relatively developed and accumulated in the bedrock, which is exposed in the landslide. There is a large amount of ice and snow meltwater, which readily forms a “pipeline” channel that is in contact with the surface of the Quaternary aeolian loess and is discharged outward in the form of a spring. The flow volume of a spring was measured to be 5 × 10⁻⁵ m³/s. The mineralization of water is less than 1.0 g/L. Based on the soil test, the saturation of the soil is 85.7–91.2%, which shows that it is a very wet sliding body. The natural moisture content of the soil was 17.2–20.4%, the plastic limit was 15.8%, and the liquid limit was 27.6%. According to the measured data of groundwater level, the groundwater in the hill is shallow and buried in the range of 0–15 m. These provide good hydrological conditions for inducing landslides and also provide good groundwater conditions for the rapid conversion of loess landslides into high speed and long runout sliding landslides.

Qualitative analysis of the disaster-formation mechanism of the Panjinbulake landslide was conducted based on a field geological survey and remote sensing satellite images, but this was far from adequate for geological disaster prevention and control. Instead, quantitative analysis methods are required to fully investigate landslide movement and predict the secondary disaster, which can be explored using the landslide dynamic analysis software DAN-W (see Section 5).

4. Basic Characteristics of Landslides and Hazard Zoning

The landslide had a long-strip shape (Figure 4a). The slope before the landslide was close to 40°, and the main slip direction was N 69° E. The elevation of the trailing edge of the landslide was about 1280 m, the elevation of the landslide shearing edge was about 1190 m, and the horizontal distance reached 375 m (Figure 4b). Due to rainfall and snow infiltration, pore water pressure and soil saturation increased during the sliding process. The landslide was transformed into debris flow, showing a flow state. Based on the information obtained from the unmanned aerial vehicle data and remote sensing, geological field surveys, the landslide can be divided into the sliding source area, debris flow area, and accumulation area (Figure 5).

4.1. Sliding Source Area (Area I)

The sliding body was mainly Quaternary loess. There were many cracks in the upper part of the trailing edge, the crack width ranged from 23 to 54 cm, and there was an unstable body in the upper part of the trailing edge. Its volume reached 1.2 × 10⁴ m³. The average width of the sliding source area was 73 m. The area was about 1.2 × 10⁴ m², the average thickness of the sliding body was 4–7 m, and the volume was about 5.0 × 10⁴ m³.

4.2. Debris Flow Area (Area II)

The debris flow area presented a long and narrow shape, and bedrock appeared in the upper part of the debris flow area. The bedrock surface had signs of scratch, and there was an exposed water head. There was also a water head on the western side of the sliding source area. This might suggest that there was sufficient groundwater in the area that promoted the sliding of the landslide. In this area, the volume of the sliding body increased due to the entrainment. The area was 6.8 × 10³ m², the average thickness of the sliding body was 3–6 m, and the volume was about 25.0 × 10⁴ m³.

4.3. Accumulation Area (Area III)

The accumulation area had a fan shape, with a length of up to 110 m along the sliding direction and a maximum width of 100 m in the vertical sliding direction. The lithology of the accumulation area is dominated by Quaternary loess, which contains moderately weathered coarse sandstone scraped off the opposite mountain. The landslide struck the opposite side of the mountain and accumulated in the Piliqinghe Gully. The area was 1.46 × 10⁴ m², the thickness of the sliding body was 3–6 m, and the volume was about 6 × 10⁴ m³.

5. Dynamic Analysis

5.1. Theoretical Basis

DAN-W is numerical simulation software developed by Hungr to simulate the whole process of landslide movement and to study the dynamics of landslides [23]. The 3D numerical model was set up according to the two-dimensional simulation conditions provided by the calculation profile in Figure 3b. Based on the aerial views, the path widths of landslide were confirmed. In DAN-W, the Lagrangian analytical solution of the Saint–Venant equation is mainly used to treat the sliding body with the rheological features that are formed by a combination of several blocks with certain materials (Figure 6). In the curve coordinates, the corresponding physical equations and equilibrium equations are established for each block (Figure 6), as in Equations (1)–(7) [23].

(1)$\mathrm{F}={\mathsf{\gamma }\mathrm{H}}_{\mathrm{i}}{\mathrm{B}}_{\mathrm{i}}\mathrm{ds}\sin \mathsf{\alpha }+\mathrm{P}-\mathrm{T}$

Here, $\mathrm{F}$ is the sliding force (N); $\mathsf{\gamma }$ is the unit weight (KN/m³); $\mathrm{H}$ is the block height (m); $\mathrm{B}$ is the block width (m); $\mathrm{ds}$ is the nominal length of the block(m); $\mathsf{\alpha }$ is the slope foot (°); $\mathrm{P}$ is the internal tangential pressure (N); and $\mathrm{T}$ is the base resistance (N); i is the block index.

(2)$V_i=v_i^{\prime }+\frac{\mathrm{g}\left(\mathrm{F}\Delta \mathrm{t}-\mathrm{M}\right)}{{\mathsf{\gamma }\mathrm{H}}_{\mathrm{i}}{\mathrm{B}}_{\mathrm{i}}\mathrm{ds}}$

Here, $V$ is the new speed when sliding body movement. The new velocity at the end of a time step is obtained from the old velocity, $v$ (m/s); $\mathrm{g}$ is the gravitational acceleration(m/s²); $\Delta \mathrm{t}$ is the time step interval(s); $\mathrm{M}$ is momentum flux; and the other parameters are the same as in Equation (1).

(3)${\mathrm{h}}_{\mathrm{j}}=\frac{2{\mathrm{v}}_{\mathrm{j}}}{\left({\mathrm{S}}_{\mathrm{i}+1}-{\mathrm{S}}_{\mathrm{i}}\right)\left({\mathrm{B}}_{\mathrm{i}+1}+{\mathrm{B}}_{\mathrm{i}}\right)}$

Here, h is the average depth of the slip mass; $\mathrm{j}$ is block boundary index; i is block index; $\mathrm{S}$ is the curve displacement (m); and the other parameters are the same as in Equation (1).

(4)$\mathrm{V}={\mathrm{V}}_{\mathrm{R}}+{\displaystyle \sum }{\mathrm{V}}_{\mathrm{point}}+{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{Y}}_{\mathrm{i}}{\mathrm{L}}_{\mathrm{i}}$

Here, $\mathrm{V}$ is the entire volume of the loess landslide deposits(m³); ${\mathrm{V}}_{\mathrm{R}}$ is the volume of the initial landslide (m³); ${\mathrm{V}}_{\mathrm{point}}$ is the volume of the unstable body (m³); $\mathrm{Y}$ is the yield rate; $\mathrm{L}$ is the length of the i block; and $\mathrm{i}$ is the block index.

Momentum and mass during the entrainment of the path material could influence landslide kinematics. To describe the entrainment process quantitatively, an entrainment ratio (ER) could be offered to calculate the increase of the landslide volume for a specific entrainment zone in the DAN model [49].

(5)$\mathrm{ER}=\frac{{\mathrm{V}}_{\mathrm{Entrained}}}{{\mathrm{V}}_{\mathrm{Fragmented}}}=\frac{{\mathrm{V}}_{\mathrm{E}}}{{\mathrm{V}}_{\mathrm{R}}\left(1+{\mathrm{F}}_{\mathrm{F}}\right)}$

where ${\mathrm{V}}_{\mathrm{E}}$ (i.e., ${\mathrm{V}}_{\mathrm{Entrained}}$) is the volume of the entrained path material (m³$)$; ${\mathrm{V}}_{\mathrm{Fragmented}}$ is the volume of the fragmented material in the sliding source area(m³$)$. ${\mathrm{V}}_{\mathrm{R}}$ is the volume of the initial loess landslide (m³); and ${\mathrm{F}}_{\mathrm{F}}$ is the fractional amount of volume expansion due to fragmentation (0.25). The entire volume of the loess landslide deposits is equal to ${\mathrm{V}}_{\mathrm{R}}\left(1+{\mathrm{F}}_{\mathrm{F}}\right)+{\mathrm{V}}_{\mathrm{E}}$ [40]. In this study, ${\mathrm{V}}_{\mathrm{R}}$ equals 5.0 × 10⁴ m³ and ${\mathrm{V}}_{\mathrm{E}}$ equals 25.0 × 10⁴ m³. The length of the entrainment area was approximately 240 m. To simulate the phenomenon of entrainment, an ER equal to 4.0 was used in the DAN model of the loess landslide. According to the pore-water pressure increased and soil saturation during the sliding process because of the snow infiltration. The landslide was transformed into debris flow, showing a flow state, so the scraping volume was huge.

The movement speed of the sliding body and the thickness of the landslide accumulation body are calculated using Equations (1)–(5). In addition, the amount of resistance encountered during the movement of the sliding body is determined by different types of rheological models. In the DAN-W software, the resistance is mainly controlled by different base rheological models. DAN-W provides a range of rheological models. According to the existing research results and the trial-error method [50,51], the Voellmy model (V) and the Frictional model (F) are more suitable for landslide dynamic hazard research. The Frictional model is mainly used for landslides when the particle sizes of the residual body are large. The Frictional model is also used for mountains with open hillside cracks where the turbulent flow is not developed. The Voellmy model is suitable for the simulation of a landslide with fractured particles where there is a visible liquidized layer in the sliding mass. From Equation, it is evident that the rheological model is proportional to the velocity of the sliding body, so it could simulate the energy damage of the turbulent flow. This was caused by the liquefied material that has high moisture content, including the loose soil covering the flow path and a spring appearing in the path. This opinion has been accepted by Geotechnical Engineering Office (GEO) of Hong Kong [14,52].

Voellmy model: The expression of base resistance is as follows

(6)$\mathsf{\tau }=\mathrm{f}\mathsf{\sigma }+\mathsf{\gamma }\frac{{\nu }^2}{\xi },$

where f is the friction coefficient of the sliding body, σ is total stress perpendicular to the direction of the sliding path, γ is the material unit weight, ν is the moving speed of the sliding body, ξ is the turbulence coefficient, and τ is the resistance at the bottom of the sliding body. The constant friction coefficient (f) is a parameter that should be determined using the Voellmy model. The friction coefficient was modified by the pore pressure and could reach much smaller values when the path material shows wet features.

Frictional model: Assume that the flow of the sliding body is controlled by the effective normal stress acting on each block. The expression of resistance τ is as follows

(7)$\mathsf{\tau }=\mathsf{\sigma }(1-{\mathsf{\gamma }}_{\mathsf{\mu }})\tan \mathsf{\phi }$

where ${\mathsf{\gamma }}_{\mathsf{\mu }}$ is the pore pressure coefficient (specifically, the ratio of pore pressure to total stress); $\mathsf{\phi }$ is the internal friction angle; σ is the total stress perpendicular to the direction of the sliding path; and τ is the resistance at the bottom of the sliding body.

5.2. Model Selection

In the DAN-W software, the accuracy of the calculation result depends on three important factors: Sliding body motion trajectory, rheological model, and parameter selection. First, based on multi-period remote sensing images, aerial imagery of drones, and field geological surveys, topographic lines before and after landslides were determined (Figure 6). Second, the Panjinbulake loess landslide was divided into the sliding source area, debris flow area, and accumulation area (Figure 5). According to the hazard characteristics of different regions, it is critical to select suitable rheological models for different regions. Since the sliding source area started from the shearing exit, there are signs of scratch on the exposed area of the bedrock. According to the existing research results [23,49,50,51], the Frictional model(F) was suitable for the sliding source area. The Frictional model and the Voellmy model(V) were used in the debris flow area and the accumulation area, respectively. According to the landslide path sequence, four sets of the rheological model combinations, Frictional–Frictional–Frictional, Frictional–Frictional–Voellmy, Frictional–Voellmy–Voellmy, and Frictional–Voellmy–Frictioanl, were used to simulate the dynamic hazard effects of the Panjinbulake loess landslide, so as to select the most suitable rheological model combination to simulate the movement process of the loess landslide (see Table 1). The simulation results calculated from a combination of the four rheological models were compared to the features of the actual loess landslide. According to the results, we could find that the Frictional–Voellmy–Frictional model fit the above characteristics. From Equation, the Voellmy model is proportional to the velocity of the sliding body, and it could simulate the huge scraping force which removed the surface loess soil. Thus, the Voellmy rheology model had a better fit than the Frictional model in the debris flow area. Finally, we intended to simulate the movement characteristics of the Panjinbulake loess landslide using the Frictional–Voellmy–Frictional rheological model.

5.3. Parameter Selection

The dynamics of high speed and long runout landslides have been studied by researchers. The analysis of landslide dynamics depends, to a large extent, on the choice of parameters and the knowledge level of the author [23]. In this paper, the simulation parameters of the Panjinbulake loess landslide were mainly obtained by the field survey data and the existing research results [23,49,50,51]. For the Voellmy model, the main parameters were ξ = 400 m/s² (the software provides a range of 200–500 m/s²) and f = 0.05. When rain and snow melt water infiltrates, the groundwater level rises, and pore water pressure rises. The sliding body is close to the flow state, so the friction coefficient decreases gradually. The influence of underground groundwater on sliding body motion is realized by changing the friction coefficient in DAN-W (see Section 7.1). For the frictional model, the dynamic friction angle φ_b was set to 19° according the literature [23,49,50,51]. Due to the infiltration of ice and snow meltwater, the excess pore water pressure increases, and Ru was set to 0.7. Finally, according to the indoor geotechnical test and the engineering analogy method, the typical strength testing index φi was set to 20° and the unit weight ($\mathsf{\gamma }$) was set to 18 KN/m³. As shown in Table 2, based on the trial and error method and the existing research results, these rheological model combinations and parameters were used to simulate the dynamic hazard effects of the Panjinbulake loess landslide.

6. Results and Analysis

6.1. Speed Analysis

Using the DAN-W software and the Frictional–Voellmy–Frictional model, the total time of the Panjinbulake loess landslide movement was 22 s. It is assumed that the speed of the landslide was 0 m/s when starting in the sliding source area. The sliding body started from the shear exit (i.e., the leading edge of landslide). Due to the steep slope of the landslide, the speed was very fast, reaching 5 m/s very quickly (Figure 7). Under the action of gravitational potential energy, the speed of the sliding body increased rapidly and reached 21.5 m/s at X = 200 m. This process lasted 7 s, and the average acceleration reached 2.3 m/s². The sliding body moved from X = 200 m to X = 320 m in debris flow area. This area has exposed spring water (Figure 5), which provides suitable hydrogeological conditions for high speed movement of the sliding body. Thus, the sliding body was always in a high speed state. This phase lasted 8 s. The average speed reached 20 m/s and the sliding body entered the river terrace. Due to the friction of sandy gravel and gradual decrease of the slope, the moving speed of the sliding body dropped sharply, and the speed relative to the opposite mountain was 17.5 m/s. Due to the blocking of the opposite mountain, the sliding body stopped at a horizontal distance of X = 370 m, and the average acceleration in the deceleration phase reached −3.6 m/s², which lasted for 7 s.

6.2. Thickness Analysis of the Accumulation Body

As shown in Figure 8, when the landslide sheared and started at X = 140 m, the initial volume of the landslide was 5.62 × 10⁴ m³, and the average thickness of the sliding body in the sliding source area was 4.5–5 m. After sliding for 7 s, it reached X = 200 m. Due to the scraping effect, the landslide volume reached 13.3 × 10⁴ m³, and the thickness of the accumulation body reached 5 m. At 10 s, the sliding body moved forward to the slope toe, X = 250 m, and the volume of the sliding body reached 23.4 × 10⁴ m³. At 15.82 s, the sliding body moved to the opposite side of the slope, at X = 350 m, and the volume of the sliding body reached 27.7 × 10⁴ m³. Finally, at 22 s, the sliding body stopped moving, and the final total volume reached 32 × 10⁴ m³. The maximum thickness of the accumulation body at X = 255 m reached 6 m, and the average thickness of the accumulation body reached 0.2–1 m in the sliding source area, 1–2 m in the debris flow area, and 5–6 m in the accumulation area. The simulated results were less than the actual measured results because the dynamic model stretched a smooth two dimensional plane into three dimensions (Figure 6).

6.3. Typical Point Velocity Analysis and Accumulation Body Thickness Analysis

Combined with the actual situation of the field geological survey, the four points with horizontal distances X = 140, 200, 250, and 350 m were selected as the typical points for analysis. The calculation results based on the DAN-W software and F–V–F rheological model are as follows.

(1) The point X = 140 m is located at the landslide shear outlet (Figure 4b). Within 0–10 s, the speed at this point increased linearly from 0 to 17.5 m/s (Figure 9a), which indicates that the acceleration of the sliding body during the starting process increased, reaching 1.75 m/s². Also, typical negative terrain was present, and the lower part of the raised bedrock was exposed. This made it so that, after 12 s, the residual sliding body moved to the point where it was hindered by negative terrain with trough shape, the speed gradually attenuated, and the final thickness of the accumulation body at this point was 0.5 m (Figure 9b).

(2) The point X = 200 m is located in the debris flow area (Figure 4b). After the landslide slid for 7 s, it reached this point. Within 7–17.5 s, the speed at this location was always relatively faster, and the average speed reached 17.5 m/s (Figure 9c). In this area, the sliding body had a rapid speed and exposure to spring water, which provides conditions for the occurrence of a high speed and long runout landslide. At 8 s, the thickness of the accumulation body reached 4.53 m. After 17.5 s, the speed at this point gradually decreased, and the final thickness of the accumulation body reached 0.5 m (Figure 9d).

(3) The point X = 250 m is located at the foot of the slope (Figure 4b). After the sliding body slid for 9 s, it reached the foot of the slope. Due to the steep slope, the peak speed reached 19 m/s (Figure 9e), and the thickness of the accumulation body reached 5 m (Figure 9f). Scattering phenomenon was appeared at the foot of the slope where in the interval t = 20–24 s, and there was an increase of debris flow depth reaching a depth higher to the front passage. In addition, this was also a turning point for the speed of landslide movement. From this point on, the speed of the sliding body began to decrease due to the sudden slowing of the slope and the friction of the sliding body against the sandy gravel of the river terrace.

(4) The point X = 350 m is located at the foot of the opposite slope (Figure 4b). After the sliding body crossed the river terrace, it reached this point at 16 s; then, the speed dropped to 11 m/s (Figure 9g), and the thickness of the accumulation body reached 1.5 m (Figure 9h). From this point on, the sliding body began to hit the opposite side of the mountain. Because of this and even though a considerable amount of energy had been consumed, the sliding body continued to climb 20 m before stopping, indicating that a high-speed and long runout landslide has incredible energy and is severely devastating.

7. Discussion

7.1. Sensitivity Analysis

The DAN-W dynamic model was preferably used to predict the dynamic characteristics of the landslide accumulation area in the Loess Plateau. However, based on a field survey and laboratory tests, some parameters, e.g., the friction coefficient, were obtained by trial and error and the existing research results [23,49,50,51]. There are many factors affecting the friction coefficient, including the rate of water content, terrain, and ground temperature, with the water content the primary factor. Therefore, it is difficult to provide the friction coefficient in a more efficient way. When using the DAN-W dynamic model to predict the risk assessment, the variability range of the parameters should be considered. In this paper, the influence of friction on the landslide is simulated using multiple sets of working conditions. According to the numerical values recommended by the kinematic model software and existing research results, the friction coefficient is divided into three groups. The first group has a low friction coefficient (i.e., 0.05, 0.1, and 0.2) which means high water content, the second group a moderate friction coefficient (i.e., 0.2, 0.3, and 0.4) which means moderate water content, and the third group a high friction coefficient (i.e., 0.5, 0.55, and 0.6) which means low water content. The calculation results are as follows (Figure 10).

According to the numerical simulation results, it can be seen from Figure 10 that, if different friction coefficients are selected, the calculated moving distances are significantly different. For example, when a low friction coefficient (0.05–0.2) is selected, the moving distance of the sliding body is 350–375 m. When the moderate friction coefficient (0.2–0.4) is selected, the moving distance of the sliding body is 320–350 m. When the high friction coefficient (0.4–0.6) is selected, the moving distance of the sliding body is 280–300 m. According to the current geological survey in the field and the aerial view of the drone, the moving distance of the landslide reached a maximum of 366–375 m, which is in line with the calculation results based on the low friction coefficient. The low friction coefficient appears to be due to the saturation degree of its material, and there was sufficient groundwater in the slope area that promoted the sliding of the landslide.

7.2. Empirical-Statistical Model

The empirical formula, mainly based on simple geometric relations of landslides (Figure 11), is simple and effective in the prediction of long runout distance. Figure 11 reveals the geometric relationship between the landslide’s apparent friction angle (i.e., the angle between the trailing edge of the landslide and the farthest point of landslide movement), height difference, and motion distance [19]. Based on the concept of the apparent friction angle, Scheidegger proposed an empirical formula called the sled model to calculate the velocity of the sliding body [19]. The specific formula is as follows

(8)$V=\sqrt{2\mathrm{g}\left(\mathrm{H}-\mathrm{t}\times \mathrm{L}\right)}$

where $V$ is the velocity of the estimating point (m/s); $\mathrm{g}$ is the acceleration of gravity (m/s²); and $\mathrm{t}$ is the rake ratio between the highest point of the rear edge and the estimating point of the sliding distance (dimensionless); $\mathrm{H}$ is the height difference from the highest point of the rear edge and the calculated velocity point of the sliding distance(m); $\mathrm{L}$ is the horizontal distance between the landslide trailing edge and the calculated velocity point (m).

According to the calculation of the sled model (Figure 12), the maximum speed of the sliding body was 28 m/s, which occurred near the horizontal distance of 250 m. Similar to the calculation results of the DAN model, the sled model showed that the slide velocity increased sharply after the landslide occurred, and the speed decreased significantly when the slide moved to the front of the road and eventually struck the opposite side of the mountain. The maximum speed obtained by the sled model is far greater than that of the DAN-W dynamic model. Instead of taking the dynamic characteristic of the landslide into account, such as erosion and entrainment, the sled model only gives a preliminary description of the process of landslide movement variation. It could be seen that the calculation results of the DAN model are more accurate. While the concept of the sled model and the apparent friction angle tends to be conservative for landslide hazard prediction, they still comprise a qualitative and effective way to predict disasters.

The landslide under study is located in Piliqinghe Basin, located in the western part of the Loess Plateau and is part of the “Belt and Road” area. The location has many potential loess landslides, all of which pose a threat to agricultural production. It is especially important to research the prediction of potential landslide disasters, which is of great benefit to disaster prevention and devising mitigation measures. During the period 2017–2018, our team carried out field geological survey work in the area and counted 12 loess landslides that occurred. At the same time, the team measured the basic parameters of the landslide, and calculated the apparent friction angle of each landslide (Table 3, Figure 11). The statistical results indicate that the apparent friction angle of the loess landslide in this area is approximately 25°. Based on the concept of the apparent friction angle, the farthest running distance of the landslide (i.e., L_max) can thus be calculated by the formula which is as follows

(9)${\mathrm{L}}_{\max }=\frac{{\mathrm{H}}_{\max }}{\tan 25°}=2.15{\mathrm{H}}_{\max }$

where ${\mathrm{H}}_{\max }$ is the height difference from the highest point of the rear edge and the farthest point of the sliding distance (m).

The empirical formula and the DAN model can be used to predict and analyze the moving distance of potential landslides from qualitative and quantitative aspects separately.

7.3. Evaluation of Landslide Residual Risk

Based on the field survey and the aerial images from unmanned aerial vehicles, multiple tensile cracks appeared in the upper part of the trailing edge of the landslide and formed an independent unstable block. In the case of rainfall and ice and of snow melt infiltration, the unstable block would be extremely easy to slide, which poses a threat to agricultural production and road operation. To avoid secondary harm, in this paper, we use the DAN model to predict and analyze the movement trend of the unstable body. According to the results of the field geological survey, the unstable block area ranges from 3 to 6 m, with a volume of nearly 1.2 × 10⁴ m³.

The selected models and parameters are the same as in the previous landslide dynamic hazard analysis. Figure 13 shows the running velocity of the unstable body, the thickness of the deposit, and the predicted disaster threat zone. Since the slope of the unstable body surface is 30°, which is relatively gentle, most of the sliding body is deposited in the slide-source area after the unstable body slides. The average thickness of the sliding body in the sliding source area reaches 5 m, the longest distance reaches 300 m, and the maximum moving speed reaches 18.5 m/s. According to the sled model, the farthest distance of the unstable body motion is 362 m. Combined with the calculation results of the DAN-W model and sled model, it can be concluded that, if the unstable landslide body starts, it will be pose a threaten to road operation. Owing to its rapid movement speed, the landside mass also threatens the safe production of grazing herds of animals. Setting up a warning sign around the landslide to warn herders to locate grazing far from the area is recommended. The local should set engineering measures (such as garbion) around the road to ensure traffic safety.

8. Conclusions

Based on the geological survey in the field, multi-period historical remote sensing images and aerial images of the drone, combined with the geological conditions of the study area, we analyzed the inducing factors and runout process of the Panjinbulake loess landslide and predicted the secondary disaster. Furthermore, the DAN-W dynamic model and a set of combined basal rheological models (Frictional–Voellmy–Frictional models) can suitably simulate the dynamic hazard effects of the Panjinbulake loess landslide. We analyzed the influence of the landslide movement speed, typical point velocity, accumulation body thickness, and friction coefficient. The simulation results showed that the duration of the Panjinbulake loess landslide was 22 s, the maximum speed was 20.5 m/s, and the maximum thickness of the accumulation body was 5.5 m, which is in line with the actual situation based on the field investigation. The basal rheological model combination and parameters obtained through trial and error can be used to simulate and predict the long runout distance of loess landslides, and it is necessary in strengthening the early identification and prevention of loess landslide hazards using multi-precision observation technology and numerical techniques.

Author Contributions

L.Y. analyses, writing and dealing with data; Y.W. analyses the field geological phenomena; W.W. offers the technology of DAN-W; S.Z. analyses the field geological phenomena.

Funding

This study was supported by the National Key Research and Development Program of China [No.2018YFC1505404], China Geological Survey (DD20190647, DD20179609, DD20190637).

Acknowledgments

The authors are grateful to O.Hungr for supplying a copy of the DAN-W software. We thank MDPI (www.mdpi.com) for its linguistic assistance during the preparation of this manuscript. We also thank Gang Liu, Chinese Academy of Sciences, Feihang Qu, Northwest University and Jie Luo, University of Science and Technology of China, for their kind assistance. Finally, we thank editors and reviewers for their thoughtful review and valuable comments to the manuscript.

Conflicts of Interest

The authors declare no conflicts of interest.

Figures and Tables

Figure 1. Location map of the Panjinbulake loess landslide (a is modified from [33]; b is modified from [48]). (a) The distribution of loess in China. (b) The distribution of loess in the study area.

Figure 2. (a) The geological structure map in Yining Country. (b) The loess landslide distribution in the Piliqinghe Basin. (c) The monthly average precipitation of Yining County over time (2011–2016). (d) The monthly average temperature–snow depth map of Yining County over time (2011–2016).

Figure 3. The multi-temporal sensing images of the landslide.

Figure 4. (a) Overview of the Panjinbulake landslide. (b) The engineering geological section.

Figure 5. Three-dimensional Digital Elevation Model diagram of the landslide and typical pictures.

Figure 6. The 3D numerical DAN-W model of the landslide. The parameters are the same as in Equations (1)–(7) (It is modified from [23]).

Figure 7. Variation of speed versus sliding range.

Figure 8. Variation of the thickness of the accumulation body at different times.

Figure 9. Typical point velocity and variation trend of the accumulation body.

Figure 10. Variation trend of the moving speed versus the sliding rage of the sliding body under the action of different friction coefficients.

Figure 11. Sketch of the Empirical-statistical model (It is modified from [19]).

Figure 12. Speed contrast diagram of two models of the Panjinbulake landslide.

Figure 13. (a) The runout distance in the DAN-W. (b) Variation of the thickness of the unstable slope in the DAN-W. (c) Variation of the velocity of the unstable slope in the DAN-W. (d)The sled model for the two different reach angles. (e) Maximum extent of the unstable slope runout at different methods.

Table 1

Hydrodynamic model of the Panjinbulake landslide.

Model	Sliding Source Zone	Debris Flow Zone	Accumulation Zone
FFF	Frictional	Frictional	Frictional
FFV	Frictional	Frictional	Voellmy
FVV	Frictional	Voellmy	Voellmy
FVF	Frictional	Voellmy	Frictional

Table 2

Parameters of the Frictional–Voellmy model used for the Panjinbulake landslide.

Model	Unit Weight, $\mathsf{\gamma }$ (KN/m3)	Internal Frictional Angle, ${\mathsf{\phi }}_{\mathbf{i}}$ (°)	Friction Angle, ${\mathsf{\phi }}_{\mathbf{b}}$ (°)	Pore Pressure Coefficient, ${\mathsf{\gamma }}_{\mathsf{\mu }}$	Friction Coefficient, f	Turbulivity, $\mathsf{\xi }$ (m/s2)
Frictional	18	20	19	0.7	-	-
Voellmy	18	20	-	-	0.05	400

Table 3

Basic geometry of the loess landslides in the Piliqinghe basin.

Number	Loess Landslide	V (m3)	Lmax (m)	Hmax (m)	Hmax/Lmax	φ (°)
1	KS1#	5160	130	60	0.46	24.77
2	KS2#	108,000	137	83	0.605	31.21
3	AX1#	109,000	130	85	0.65	33.02
4	AX2#	10,000	75	37	0.49	26.10
5	KZ1#	12,000	335	145	0.43	23.27
6	KZ2#	17500	341	131	0.384	21.01
7	KZ3#	30,000	175	40	0.23	12.95
8	KZ4#	10,000	47	25	0.53	27.92
9	PL1#	66,000	175	85	0.49	26.10
10	PL2#	7000	74	40	0.54	28.36
11	PL3#	3440	136	55	0.41	22.30
12	Panjinbulake	300,000	375	160	0.43	23.26
Average					0.46	25