Critical area identification and dynamic process simulation for landslide hazard chain formation in the upstream Jinsha River

The upper reaches of the Jinsha River, with their complex terrain and active tectonic movements, are vulnerable to landslide-induced hazard chain events, which endanger the safety of residents and infrastructure construction. Based on the analysis of the development background of the hazard chain in the upstream area of the Jinsha River, five factors, including the lithology, distance to faults, distance to rivers, peak ground acceleration, and slope degree, were selected to identify the critical landslide-prone areas. Principal component and grey correlation analyses were then conducted to determine the contributions of these different factors. Based on ArcGIS, the study zone was categorized into five classes of landslide susceptibility: very high, high, moderate, low, and very low. The identification of the critical target areas for landslide hazard chain formation showed satisfactory accuracy. The very high- and high-susceptibility areas are concentrated along the Jinsha River. The dynamic process of a typical landslide in a very high-susceptibility area was numerically simulated using OpenLISEM. The high-precision Baige landslide data of the study area were used to calibrate the practicality of the input mass parameters, including cohesion, internal friction angle, D50, and D90. The movement and accumulation processes of a typical landslide were then numerically simulated with the verified data. The entire landslide accumulation covers an area of 0.45 k m 2 , with a length of 1,600 m and a width of 270 m. Thus, the OpenLISEM model, which combines mass, topography, and landcover parameters, is feasible for the numerical simulation of landslide dynamic processes. The prediction of the dynamic processes and accumulation morphology of landslides can provide a reference for the formation processes and mechanisms of the landslide-induced hazard chain in the upper Jinsha River.


Introduction
Due to the compressed active fault zone and severe erosion of the valley shear, the Jinsha River basin exhibits unique fold mountains and deep valley geomorphology Zhu et al., 2021). With the complicated lithology and steep slope degree, the upper reaches of the Jinsha River are vulnerable to landslide-dam breach-flood hazard chain events, e.g., the Baige , Sela (Zhu et al., 2021), Temi , and Guili (Xu et al., 2022) landslides, which threaten the safety of people and property upstream and downstream (Wei and Siming, 2020;. Therefore, the identification of critical areas and simulation of the dynamic processes of landslide hazard chain formation is critical in the upstream Jinsha River basin. Field surveys are exceedingly difficult to perform in this region because the landslide disaster site is in a remote alpine canyon area with a harsh environment and limited access. Thus, high-precision ground observation technologies, including remote sensing, interferometric synthetic aperture radar (InSAR), and light laser detection and ranging (LiDAR), are used to identify landslide hazards (Lu et al., 2019). Moreover, Landsat, ALOS, Quick Bird, and other satellites with high-resolution images have been widely used for landslide identification; however, they only identify landslide hazards with clear deformation (McDonald and Grubbs, 1975;Sato and Harp, 2009;Youssef et al., 2009;García-Davalillo et al., 2014). InSAR with centimeter-high precision and all-weather surface observation technology can perform large-scale landslide identification and monitoring to build a landslide hazard inventory (Woods et al., 2020;Urgilez Vinueza et al., 2022;Zhang et al., 2022). Because of the satellite lateral imaging mode, the terrain may cause geometric distortion of SAR images, including perspective shortening, overlap, and shadow (Jie et al., 2018). LiDAR can not only directly obtain three-dimensional terrain coordinates to provide high-accuracy topographic images but also remove the effective influence of vegetation to obtain true ground elevation data (Gorsevski et al., 2016;Abdulwahid and Pradhan, 2017;Xu et al., 2019). However, LiDAR is not suitable for large-scale disaster identification due to its harsh operating conditions and high costs (Lu et al., 2019). Thus, while high-accuracy remote sensing technology has become an important means to obtain landslide hazard information, the reliability of remote sensing data interpretation and the accessibility and costs of high-precision image data remain issues in landslide risk assessment. In addition to satellite sensing technology, scholars have proposed statistical models for landslide susceptibility assessment, with remarkable results (Song et al., 2012;Xu et al., 2013;Khan et al., 2019;Liu D et al., 2021;Yong et al., 2022). Models such as the information quantity model (Yang et al., 2018), the weight-of-evidence model (Wang et al., 2016), the logistic regression model (Xing et al., 2021), the analytical hierarchic process (Yoshimatsu and Abe, 2006), and the principal component analysis (Chang et al., 2014) have demonstrated excellent accuracy. Based on the rich results in the selection of evaluation factors, the weight of the factor determination, and the construction of the evaluation model, the present study applied mathematical statistics to identify the key target areas of landslide hazards to not only objectively analyze the assessment of regional landslide disaster susceptibility but also the relationships and impact degrees of slope failure factors.
Many numerical simulation tools for the simulation of landslide dynamic processes are effective for the analysis of landslide hazards. Pastor et al. (2021) and Ouyang et al. (2019) used the depth-integrated continuum method for the dynamic simulation of debris and landslides. Based on the discrete element method An et al. (2021) established an adapted Hertz-Mindlin contact model between particles and the ground surface to accurately simulate the landslide dynamic process. Zhang et al. (2012) simulated the entire process of failure and instability of the Jiweishan high-speed remote landslide in Chongqing via PFC 3D software. Smoothed particle hydrodynamics (SPH) based on Lagrangian particle-based meshless methods have been widely employed in applications in geotechnical engineering (Peng et al., 2019;Zhu et al., 2020). The depth-integrated continuum and discrete element methods consider the landslide body as a fluid and can efficiently simulate the landslide motion and accumulation processes. However, they cannot easily simulate the initial failure mechanism and triggering factors such as the erosion and volume expansion processes due to the disaggregation and fragmentation of the rock mass (Wen- Jie et al., 2021). The landslide movement process is affected by many factors, including mass composition, trigger mechanism, and vegetation environment. The above numerical simulation methods do not consider the interaction between the material source and the environment during the movement process. Aiming to precisely assess the landslide dynamic process, the present study used OpenLISEM to perform the numerical simulations. OpenLISEM divides the research area into several grids of equal size; inputs corresponding terrain, vegetation, and material source parameters into different grids; couples the distributed basin hydrological and two-phase flow models; compares the slope, ground roughness, and other parameters of adjacent grids; and simulates disaster processes such as landslides, debris, and flash floods (OpenLISEM manual 2017, https://lisemmodel.com).
The complete geohazard chain includes potential hazards, primary hazards, secondary hazards (series), and hazard-bearing bodies (Peng and Jian, 2021). This study identified potential geatators by hazard susceptibility assessment and researched the dynamic process of the primary hazard by numerical simulation. Based on remote sensing and field investigation data, the weights of the landslide factors were determined via principal component analysis and the grey correlation degree method. A landslide susceptibility evaluation model was then constructed to identify high-susceptibility landslide areas. Based on the vulnerability assessment, the key target areas of landslide hazards were numerically simulated in OpenLISEM to obtain the movement process and accumulation form of landslides, and to provide a reference for a landslide-dam breach-flood disaster chain prevention in the upstream area of the Jinsha River.

Study area
The study area is situated in the upper Jinsha River valley, at the intersection of Tibet and Sichuan, and covers an area of 60,352 km 2 . The geographic coordinates are 97.33°E-99.94°E and 28.68°N-32.73°N (Figure 1). The terrain is characterized by "V"shaped gullies in the high-mountain region, which is severely affected by river erosion. The elevation of the study area ranges from 2092 to 6,088 m. Based on the genesis, the geomorphic types in the study area can be divided into five types: erosion accumulation, structural erosion, structural denudation, structural dissolution, and glacier geomorphology (Bai et al., 2014). The main fault strikes have a northwest-southeast orientation. The lithology of the study area is
The study area has highland climates with an average temperature of -4.9°C-7.8°C. Furthermore, due to the high terrain and strong solar radiation, the area is cold and dry with 53-67% relative humidity, an annual average rainfall of 387.0-657.6 m, and an average of 100.3-169.8 days of precipitation. The rainy season is from May to October, accounting for more than 90% of the annual rainfall. Snow is the main form of precipitation in the region, and the daily rainfall in most areas is 50 mm. Additionally, the vegetation is mainly grassland and swamp meadows (Wu, 2007).

Database and method
3.1 Susceptibility assessment 3.1.1 Factor sources Landslide susceptibility mapping is essential for identifying areas with high landslide risk. In this study, the landslide inventory with a total of 635 hazard points was extracted from the Resource and Environment Science and Data Center (RESDC) of the Chinese Academy of Science. Based on previous research, the specific characteristics of the study region, including lithology, distance to the fault, distance to the river, and peak ground acceleration (PGA), were selected to assess the landslide susceptibility (Table 1). The lithology and fault were obtained from the RESDC of the Chinese Academy of Science, while the river data were obtained from the Open Street Map. The fault and river buffer distances were calculated using the surface analysis tools in ArcGIS. The slope data were extracted from the digital elevation model (DEM) downloaded from the advanced land observation satellite (ALOS) with a resolution of 12.5 m × 12.5 m. The influence degree of seismic intensity was reflected by PGA extracted using the seismic ground motion parameter zonation map of China (GB18306 2015).

Factor classification
The geographic information system (GIS) was used to show the terrain features of the lithology and other factors. The hazard inventory, combined with the five evaluation factors, was digitized and stored in the GIS. The classification of the five factors in this study is shown in Figure 2. The lithology composition in the study area is sophisticated. Herein, it is divided into five categories: soft rock (I), half-soft rock (II), halfhard rock (III), hard rock (IV), and others (V), which include loose deposits, complex structural surfaces, and unknown rock. Landslides are mainly located in soft rock composed of siltstone, sandstone, etc. because these rocks are vulnerable to external forces ( Figure 2A). The faults are densely distributed in the study area. This study selected the Jinsha-Honghe fault as the main research target. The rock near the fault zone is relatively broken and has poor stability. The landslide activity intensity decreases with increasing distance from the fault ( Figure 2B). The surface stream deepens in the valley in the study area through erosion and transportation, producing favorable conditions for hazards. More than 50% of the hazards are located within 4 km of the river ( Figure 2C). PGA was used to represent the magnitude of the earthquake intensity and was divided into <0.15, 0.2, and 0.3 g. According to the seismic intensity map in the spatial distribution, most landslides are mainly concentrated in areas of high earthquake intensity ( Figure 2D). The slope gradient was obtained using the 3D-analysis slope tool in ArcGIS. The landslides are concentrated in regions ranging from 10°t o 40°( Figure 2E).

Factor weight calculations
According to the landslide intensities for the selected factors, a score was given to the different intervals of each evaluation factor ( Table 2). The discriminant coefficient matrix of the landslide inventory in the study area and factors is shown in Table 3.
The contribution degrees of the five factors were obtained via principal component analysis (Chang and Tang et al., 2014). The discriminant coefficient matrix (Table 3) was imported into SPSS software for principal component analysis to obtain the contributions of the factors. The first four factors, i.e., lithology, faults, rivers, and PGA, contributed to >80%. Lithology was the most influential factor and was, therefore, selected to determine the contribution of the other factors (Table 4).
The grey correlation analysis (Wei et al., 1998) mainly determines the correlation between the dominant factor and the other four factors. When the variable was transformed into dimensionless data, the outcome was a coefficient matrix of factors and initialized data. From this, the absolute D-value Eq. 1was calculated: where i=1,2. . . . . .m; k=1,2,. . . . . .n; x i (k) =coefficient matrix of factor lithology; and Δ i (k) =absolute value. The extreme of the matrix (Δ min and Δ max ) was obtained via Eq. 2 The correlation coefficient of the evaluation factors was calculated using Eq. 3:   where ξ i (k) is the correlation coefficient and ρ is the distinguishing coefficient, which takes 0.5 as the value. The size of the value can control the data transformation and significant differences in the correlation coefficients. The greater the correlation between other factors and the dominant factor, the greater the impact on geological disasters. The weights of the evaluation factors are shown in Table 5.

Model generation
Five-factor raster maps with 12.5 m resolution were constructed using GIS. The values of the factor raster were calculated using the Raster Computing Tool. The final model was established using these factors.
where i = 1, 2. . .. . .5, R is the value of landslide susceptibility evaluation, x i (k) is the value of the landslide in Table 2, and ω i is the weight of the landslide in Table 5.

Landslide numerical simulation 3.2.1 Theory
To better simulate the dynamic process of critical landslide hazards, a discrete numerical modeling method was applied using the open-source software OpenLISEM, which requires the subdivision of both space and time into a discrete set of locations (https://lisemmodel.com). Originally, OpenLISEM was a physically based numerical model designed to simulate event-based runoff, flooding, and erosion on a catchment scale. By combining solid and water runout flow equations, OpenLISEM includes a series of dynamic hydrological processes, such as precipitation, interception, surface flow, splash detachment, erosion, and sediment transportation. Many revisions and additions have subsequently been incorporated into the OpenLISEM application. OpenLISEM was further developed as a multi-hazard model, including groundwater flow, slope stability, slope failure, mass movements, deposition, entrainment, and earthquake effects (Pudasaini, 2012;Bout et al., 2018;Scaringi et al., 2018). Furthermore, the model incorporates the iterative slope failure method based on a modified infinite slope mode. The conventional infinite slope model predefines the bottom of the soil layer as the potential slip surface, while the iterative method iteratively searches the potential slip surface. The equation for determining the factor of safety (FOS) is

FOS
c + c′ + γ − mγ ω z + mγ ω z cos 2 β tan φ′ γ − mγ ω z sin β cos β where c and c′ (kpa) are the effective soil cohesion and root cohesion, respectively; γ and γ ω (kg/m 3 ) are the soil and water densities, respectively; m is the effective saturation level of the soil; z (m) is the soil depth; β (°) is the slope angle; and φ′ (°) is the effective internal friction.
The two-phase runout flow within OpenLISEM is a combination of water and solid dynamics (Pudasaini, 2012). Using it, landslides, water flow, and debris flow can be simulated, including their interactions. The full momentum source terms for both the fluid and solid phases are as follows: where S s and S f (m/s 2 ) are the momentum source terms for the solid and fluid phases, respectively; α s and α f are the volume fractions for the solid and fluid phases, respectively; P b (kg/ms 2 ) is the pressure at the base surface; b (m) is the basal surface of the flow; N R is the Reynolds number; N RA is the quasi-Reynolds number; C DG is the drag coefficient; γ is the density ratio between the fluid and solid phase; χ (m/s) is the vertical shearing of fluid velocity; ε is the aspect ratio of the model; and ξ (1/m) is the vertical distribution of α s .

Data input and calibration
This study ignores the interception model. Therefore, the input data of the OpenLISEM model can be divided into three categories: landcover, mass, and topography parameters.
The landcover parameters include land-use type, vegetation cover (vegc), vegetation height (ch; m), and leaf area index of the plant cover in a grid cell (lai; m 2 /m 2 ). Land-use type and ch were obtained from Wang et al. (2012). The upper reaches of the Jinsha River are mainly grassland and farmland. Based on ArcGIS software, vegc was obtained from the linear range of NDVI (normalized difference vegetation index) from Landsat remote images (http:// Frontiers in Earth Science frontiersin.org www.gscloud.cn/search). Finally, lai was derived as follows (Choudhury, 1987;Choudhury et al., 1994): Topography parameters such as the slope, sine of slope gradient in the direction of the flow (Gradient), random roughness (RR), local surface drainage direction network (LDD), and main catchment outlet corresponding to LDD (Outlet) can be derived from the DEM using the PCRaster program. Manning's index, another topography parameter, was obtained using the OpenLISEM Manual (2017), according to land use (Wang, 2018).
The mass parameters numerically used in the simulation include mass depth (mm), initial moisture (-), cohesion (kpa), internal friction angle (radians), porosity (-), density (kg/m 3 ), and D50 and D90 (cm). The mass depth is calculated based on the empirical formula reported by Tang et al., 2012: where T is the average soil depth (m) and s l is the landslide area (m 2 ).
It is difficult to directly obtain accurate values for the other mass parameters due to the limitations of field testing technology. On October 10 and November 3, 2018, two large landslides near Baige village occurred in the same location on the right bank of the Jinsha River ( Figure 5). Many scholars performed field investigations, and research on the slope failure mechanisms and parameters showed inversion in the Baige landslide, which provided a valuable opportunity to adjust parameters (Ouyang et al., 2019;Xu et al., 2021;Zhou et al., 2020;Zhao et al., 2020;Zhang et al., 2020;Zhou et al., 2022;Liu X et al., 2021;Sun, 2021;Wang et al., 2019).
The first Baige landslide was successfully simulated using the OpenLISEM program. The results were consistent with the modeling performed by Ouyang et al., 2019 (Figure 3). The appropriate parameters were obtained by numerically simulating the first Baige landslide event, as shown in Table 6.

Critical landslide identification
To identify the critical areas for landslide hazard chains, factors such as lithology, fault, river, PGA, and slope were selected for the landslide susceptibility analyses. The five landslide factors were weighed via principal component and grey correlation analyses. Table 5 shows that lithology has the greatest impact on the landslide hazard chain, with a weight value of 0.33, while PGA has the least impact, with a weight value of 0.07. Fault, river, and slope have weights of 0.2, 0.27, and 0.13, respectively. Consequently, the R-value reflecting the landslide susceptibility ranges from 1.33 to 5, which is classified into five groups (Table 7). Very high and high susceptibility areas account for 45.42% of the total area. However, the area accounts for 75.28% of the landslides. To verify the reliability and applicability of the model for the vulnerability evaluation of landslide hazard points in the study area, ROC (receiver operating characteristic) curve and AUC (area under the curve) values were selected for testing. The ROC curve is an effective method to evaluate the performance of the classification algorithm; that is, the relationship between the simulated and sampled values. The horizontal axis is the cumulative value of the false positive rate (FPR), which indicates the proportion of susceptible areas, and the vertical axis is the true positive rate (TPR), which  indicates the cumulative value of the proportion of disaster points. The AUC value represents the area between the ROC curve and the abscissa axis and, with values in the range of [0.5,1], it indicates the good fitting effect of the model on the processed data. Figure 4 shows that the AUC value of this model is 0.70. These results demonstrated the high evaluation accuracy of the susceptibility assessment model. The landslide susceptibility distribution map in the Jinsha river upstream is shown in Figure 5. The very low and low susceptibility areas are mainly located in the northeast and southwest of the study area, mainly covering the Dege and Mangkang counties. The moderate susceptibility area is very dispersed and primarily located in central Batang county and the western area of Dege county. The very high and high susceptibility areas with the largest floor areas are located near the main stream and tributaries of the Jinsha River.
Based on the susceptibility assessment results, the Litang landslide area with very high susceptibility was selected as a critical hazard for the numerical prediction simulation (Figure 6). The landslide shape was obtained from Cui et al. (2020).

Dynamic numerical simulation
The Litang landslide, with an area of 1.32 km 2 , is located near the Baige landslide. Based on Eq. 7, the mass average soil depth is 15.2 m. Based on the GIS platform, the mass, topography, and landcover parameters were converted into a "map" file by QGIS and PCRaster. These map data were then input into OpenLISEM for numerical simulation.
According to the simulation results, Figure 7, which is plotted at 20s intervals, displays the entire dynamic process of the landslide from startup to relative stability, which occurs in approximately 100 s. At the initial stage of the landslide dynamic process, the rock mass encounters stability failure and begins to slide along the bedrock surface. From time t = 0-20 s, the landslide body slides downstream in the northwest direction with an average velocity of about 12 m/s ( Figure 7A). ROC curve.

FIGURE 5
Landslide susceptibility map in the study area.

Frontiers in Earth Science frontiersin.org
From t = 20-40 s, the landslide movement is uniform with an average velocity of 15 m/s ( Figure 7B). The mass at the front edge of the landslide reaches the Jinsha River at t = 40 s. The main acceleration stage refers to the period from t = 40-80 s, and the peak value of the velocity reaches approximately 59 m/s ( Figure 7C). During this period, the rock mass is divided into two parts,i.e., upstream and downstream, due to the influence of the watershed in the central section of the landslide. The deceleration stage refers to the period between t = 80-100 s Figure 7D). During this stage, the average velocity of the rock mass near the upstream part approaches zero, while the average velocity of the rock mass near the downstream part is <10 m/s.

FIGURE 6
Critical landslide location in the study area. Simultaneously, some of the mass materials reach the Jinsha River, and the upstream part of the mass begins to form a barrier dam. When t = 100 s, the accumulation stage is complete, and the landslide nearly stops. A large volume of the mass is deposited in the valley within this period ( Figure 7E). As shown in Figure 8, the entire landslide dam covers an area of 0.45 km 2 , with a length of 1,600 m and a width of 270 m. Due to the terrain, the final accumulation of the landslide includes two barrier dams with an average thickness of 17 m, which are thin on the sides and thick in the middle. The maximum accumulation thickness is higher on the upstream side relative to the downstream side of the dam (upstream, 55.32 m; downstream, 44.33 m).

Discussion and conclusion
This study used susceptibility assessment to identify an area with high susceptibility for a landslide, which was numerically simulated using an OpenLISEM model, to provide a reference for geological hazard prevention in the upper reaches of the Jinsha River.
For susceptibility assessment, the landslide inventory and five hazard factors, including lithology, fault, river, PGA, and slope data, were digitalized and categorized in the ArcGIS model. Based on the landslide intensity in the classification of the five factors, the discriminant coefficient matrix of landslides in the study area was constructed. Principal component and grey correlation analyses were performed to calculate the weights of the factors, which indicated that lithology had the largest impact on landslides. Furthermore, an assessment model was established. The susceptibility results indicated that the high-susceptibility zone accounts for 45.42% of the total area but comprises 75.28% of the landslide numbers, which are located near the main stream and tributaries of the Jinsha River. Therefore, the susceptibility map not only can be used as a basic tool for critical landslide identification but also helps in land use planning. A lack of accurate factor measurements may affect the precision of the factor data. Despite the limitations mentioned above, the ROC curve results showed that the proposed method has the potential for risk reduction in the study area. Moreover, we performed susceptibility analysis for one model. Quantitative and qualitative models are increasingly applied to research on landslide susceptibility, with continued improvements. Landslide hazard data provide important information for susceptibility assessment. However, due to technical limitations, additional research is needed on hazard risk evaluation in areas lacking landslide hazard data.
According to the susceptibility assessment results, the Litang landslide from the high-susceptibility area was selected as a typical hazard for simulation by OpenLISEM. Using the landslide that occurred in Baige on October 10, 2018, as an example, the mass parameters were calibrated to be inconsistent with the movement of the Baige landslide. Then, the mass, topography, and landcover parameters were input into OpenLISEM. The landslide dynamic process simulation was completed in approximately 100 s and comprised four stages: initial start-up, acceleration, deceleration, and accumulation. In the acceleration period, the landslide body was divided into two parts due to the terrain. The entire deposition covered an area of 0.45 km 2 , with a maximum thickness of 55.32 m. The numerical simulation analysis of the Baige and Litang landslides showed that OpenLISEM can be applied to the research and analysis

Frontiers in Earth Science
frontiersin.org of landslide movement. Using the PCRaster platform, OpenLISEM has good compatibility with the GIS platform, which can provide a reference for further analysis of secondary disaster evolution research (based on the GIS platform). OpenLISEM model simulation requires a large amount of accurate data, and some areas cannot be used for simulation studies of landslide dynamic processes without good simulation data.

Data availability statement
The original contributions presented in the study are included in the article/supplementary material. Further inquiries can be directed to the corresponding author.