Characteristics and Dynamic Process Modeling of the Rainfall-Induced Landslide on August 21, 2020 in Hanyuan County, China

On August 21, 2020, a landslide occurred in Zhonghai village of Hanyuan County, Sichuan Province, China. The landslide is triggered by two successive rounds of heavy rainfall. This landslide can be clearly divided into an initial landslide and a main landslide. The main landslide is activated by the intense impact and overloading of the initial landslide. The depth-integrated continuum method is adopted to simulate the dynamic process. Due to the complicated failure process, it is found that there is no proper unified parameters in Coulomb model which could well reproduce the two successive landslides. It implies that the dynamic process of landslides is highly associated with the characteristics of sliding bodies. Here, an implement of variable frictional coefficients for different parts is proposed and the parameters are calibrated. It is demonstrated that results from numerical modeling match well with the field investigation. The complicated landslides in two different stage can both be efficiently revealed by depth-integrated continuum modeling.


INTRODUCTION
Rainfall especially continuous rainfall and heavy rainfall is the most common cause of landslides. In the past several decades, rainfall-induced landslides have occurred frequently all over the world (Ayalew 1999;Lin and Jeng 2000;Chen and Lee 2004;Guzzetti et al., 2004;Huang 2009;Yin 2011). Recently, there are several catastrophe rainfall-induced landslides in China, which have caused devastating disasters and is a huge challenge to residents and project (Xu et al., 2012;Yin et al., 2016;Yin et al., 2017;Ouyang et al., 2018;Fan et al., 2020;Xue et al., 2021). For instance, the rainfall triggered landslide on June 28, 2010 in Guanling County destroyed two villages and killed 99 people (Xing et al., 2014). With frequent extreme weather (Wei et al., 2020;Yang and Huang 2021), it is essential to understand the failure mechanism and dynamic process under heavy rainfall.
Landslide runout analysis is mainly implemented in three ways: physical experiments, empirical-statistical methods and numerical models. Physical experiments are widely used to study on the runout mechanism and the deposition characteristics of landslides (Kokelaar et al., 2014), but it is usually limited by the size effect. Empirical-statistical methods mainly focus on geometrical correlations: 1) the correlation between the angle of reach and the landslide volume (Scheidegger 1973;Corominas 1996;Hunter and Fell 2003); 2) the correlation between deposition area of landslide and the landslide volume (Iverson et al., 1998;Griswold and Iverson 2008); 3) the correlation between the travel distance and other geometrical parameters (Finlay et al., 1999). By the empirical formula based on statistical data, the runout distance and deposition area can be predicted simply and effectively. In addition, empiricalstatistical methods are useful for risk assessment mapping in large regional scales (Griswold and Iverson 2008;McDougall 2017). For all that, empirical-statistical methods cannot provide dynamic characteristics such as speed and deposition depth of landslides at a specific moment.
With the improvement of computer operation efficiency, various numerical models have been used to analysis the runout and dynamic characteristics of landslides. Common numerical models for landslide include: 1) discrete element method (Tang et al., 2009;Lin and Lin 2015); 2) the smoothed particle hydrodynamics (Pastor et al., 2009;Pastor et al., 2014); 3) discontinuous deformation analysis (Beyabanaki et al., 2016;Zhang et al., 2016); 4) the depth-integrated continuum method (Savage and Hutter 1989). Compared with other methods, the depth-integrated continuum method has a advantages of higher computational efficiency and less parameters (Iverson and Ouyang 2015), so it is widely used to simulate runout of landslides (Hungr 1995;Pirulli 2005;McDougall 2006;Beguería et al., 2009;Iverson and George 2014). Due to the complexity and diversity of landslides, it is difficult to determine motion parameters of different landslides and predict the runout of landslides accurately. Therefore, back analysis of the landslide parameters is very meaningful for understanding and predicting the runout of similar landslides. As a general rule, previous studies on numerical simulation of landslides mostly focused on the motion parameters inversion of single landslide and the parameters are used to predict the landslides at the same location. However, motion parameters of landslide at the same area could be different and even different zones of a landslide may have different motion parameters. Therefore, this study focuses on the motion parameters and processes of different landslides parts.
At 3:50 AM on August 21, 2020, a landslide (located at E 102°41′44″ and N 29°20′35″) occurred at Zhonghai village in Hanyuan County of Sichuan Province, China. The landslide killed nine fatalities, destroyed eight houses and caused traffic disruption of provincial road S435. During the heavy rainfall, deformation signs were found on the slope, and the landslide occurred after the rainfall. In this paper, the remote sensing method based on unmanned aerial vehicle (UAV) is used to obtain geographic reference digital surface model (DSM) and digital orthophoto model (DOM) and the failure mechanism of the landslide is analyzed. The terrain data is processed based on the pre-and post-disaster terrain and SLBL (Jaboyedoff et al., 2004;Jaboyedoff et al., 2009) is used to determine the sliding plane. Finally, the depth-integrated continuum method is used to inverse the dynamic process of the initial and main landslides.

REGIONAL GEOLOGICAL AND GEOMORPHOLOGY SETTING
The landslide occurred in Hanyuan County, which is 285.7 km away from Chengdu and 147 km away from Ya'an City ( Figure 1A). It is located in the transitional area between the Sichuan Basin and the Qinghai-Tibetan Plateau and its landform is mainly mountainous ( Figure 1B). The Dadu River runs from west to east and the Liusha River runs from north to south through Hanyuan County. The center of the study area is river valley area and surrounded by high mountains. The highest peak reaches an altitude of 4021 m, the lowest altitudes is 550 m, and the maximum relative elevation difference is 3471 m in the study area. The region is prone to landslides due to high mountains, steep slopes and incision of rivers.
Due to frequent tectonic deformation in history, the soil mass is loose and the integrity of rock mass is poor in Hanyuan County. The overlying soil of the landslide area is cultivated soil and sandy soil. The bedrock is strongly weathered sandmud interbeds, which have bad physical and mechanical properties. In past decade, multiple earthquakes such as the Lushan Ms 7.0 earthquake in 2013 ( Figure 1C) occurred in nearby faults, which could further cause degradation of geological properties of rock and soil mass. Generally, The features of rock-mass body in the landslide area are relatively bad.
Hanyuan County have a subtropical monsoon-influenced humid climate with four distinct seasons. Annual mean temperature in the study area is 17.7°C (the maximum temperature is 40.9°C and the minimum temperature is −29°C). The study area has an mean annual rainfall of 741.8 mm, which has apparent seasonal mal-distribution. The rainfall in summer and autumn (from July to August) occupies 80% of the annual rainfall, while the rainfall in winter and spring (from November to April of the next year) occupies only 12% of the annual rainfall. There were two rounds of heavy rainfall and an accumulated rainfall 167.5 mm FIGURE 2 | The rainfall in 1 h intervals and cumulative rainfall before the slide. The data is provided by Anle Town rainfall station which is 2.4 km away from the landslide.
Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 720872 ( Figure 2) before the landslide occurred. With the infiltration of rainwater, the weight of rock and soil in the slope increased and the mechanical properties decreased, resulting in slope failure. Hence, rainfall is deemed as the direct triggering factor of the landslide.

GENERAL DESCRIPTION OF THE ZHONGHAI LANDSLIDE
On August 26, 2020, field investigations on Zhonghai landslides in detail were carried out. A total of 803 photographs were taken by UAV (unmanned aerial vehicle) covering an area of 0.3 km 2 . DSM and DOM with an accuracy of 0.055 m were established. The DOM of the landslide is shown in Figure 3. This landslide is relatively complicated and owns a different stage. It is not feasible to describe the landslide in a single process. Based on the field investigations, comparison from pre-and post-disaster, sliding plane analysis by SLBL (Jaboyedoff et al., 2004;Jaboyedoff et al., 2009), the landslide is partitioned into the initial landslide, the main landslide and two unstable mass. Zhonghai landslide has a runout distance of 590 m, the highest elevation of 1013 m and the lowest elevation of 840 m with a total area of 7.8 × 10 4 m 2 . The slide direction of the initial landslide was 16°NE. The initial landslide has the maximum elevation difference of 53 m, a length of 120 m and a maximum width of 134 m, covering an area of 1.1 × 10 4 m 2 . The initial landslide volume is about 8.5 × 10 4 m 3 and there are scarps about 3-5 m high at the upper edge of the initial landslide. Two unstable mass and the main landslide is located below the initial landslide. The unstable mass Ⅰ with an area of 1.27 × 10 4 m 2 is located on the left side of the main  The topography of the original slope provided a favorable terrain condition for the initiation of the landslide. The original slope was high in the south and low in the north with irregular plane shape. The slope topography was summarized as "steep-gentle-steep". The upper part of the slope is under good free face condition with a slope gradient of about 35°and the middle slope is gentle with a slope angle of 10°-20°as well as the lower slope is steep with a slope angle of 20°-25°. This topographic feature made the upper slope become unstable firstly and caused the chained reaction.
The overlying soil of Zhonghai landslide is cultivated soil and sandy soil, with loose structure, large gap and strong permeability which is beneficial to rainwater infiltration. The rock under the soil layer is weathered sand-mud interbeds. Under the long-term moistening of infiltration water, the mudstone interlayer is obviously softened and argillized, so its resistant shear strength is greatly reduced. The rock and soil in the landslide area lay the material foundation for the deformation and failure of the slope.
The landslide occurred in rainy season and experienced longterm rainfall before the landslide. According to the data of rainfall stations near the landslide (Figure 2), there are the two rounds of heavy rainfall and the most intense hourly rainfall is more than 30 mm/h. The continuous rainfall lead to increase the weight of the slope and decline the mechanical properties of rock and soil, especially weak mudstone interlayer, which is the direct triggering factor of the landslide.
At 3:50 am, on August 26, 2020, after two rounds of heavy rainfall, the steep slope in the upper area became unstable and formed the initial landslide. Under the influence of the impact and loading of the initial landslide, the slope in the middle area had intense deformation, resulting in two unstable mass and triggering the main landslide. The landslide mass in the source zone was relatively deep and its movement is constrained. It imposed strong impact and loading on the bottom slope. The landslide mass in the fluidization zone moved fast so it destroyed the house and buried the provincial road.

Model Description
To analyze the depth and velocity distribution of landslide in the process of sliding, the depth integral continuous statistical method is adopted. The governing equation is as follows where ρ means the averaged density, h is flow height, u and v are the averaged velocities in the x and y directions, respectively. (τ zx ) b and (τ zy ) b are the basal shear stress components in the x and y directions, respectively. k active and k passive present the lateral earth pressure in the active elongation and passive compression states, respectively. The Coulomb frictional model is adopted and it can be written as follows: where c and δ are the cohesion and basal friction angle of the mass, respectively. The above-metioned equations are solved by MacCormark-TVD method and it has been implemented in Massflow software (Ouyang et al., 2013). The Massflow has been successful applied different kinds of earth-surfaced flows, such as landslides (Ouyang et al., 2017;Wei et al., 2020), debris flows (Ouyang et al., 2015a;Horton et al., 2019) and floods (Ouyang et al., 2015b).

Analysis of the Numerical Results
The pre-disaster geomorphology and building distribution were obtained from the historical Google image. The pre-disaster terrain is obtained from Sichuan Bureau of Surveying, Mapping and Geoinformation with a resolution of 5 m. The post-disaster terrain is build based on UAV images with resolution of 0.055 m. In the computational procedure, all terrain data were resampled to 1 m grid. The ratio of mobility (H/L) of the initial landslide and the main landslide is determined to be 0.44 and 0.26. First, the ratio of motion are used as the friction coefficient in numerical simulation and the numerical results don't match with the actual condition very well. Then, the parameters are adjusted for trial calculation. In addition, it is difficult to obtain good results by only using one parameter in the simulation of the main landslide. Finally, by parameter calibrations, the basal friction coefficient is set to be 0.5 for the initial landslide; the basal friction coefficient 0.31 (in the source zone) and 0.25 (in the fluidization zone) is adopted for the main landslide.
Deposition thickness contour map of the initial landslide at t 0, 2, 4, 6, 8 s are shown in the Figure 4. At t 2 s, the landslide mass moved down fast along the direction of 15.5°NE. At t 4 s, the front of the landslide reached the mild slope. After 4 s, the landslide movement was blocked by smooth terrain and the landslide mass began to deposit. After 6 s, most of the mass was deposited in the mild slope and sliding surface. The final deposition area of the numerical simulation is compared with data acquired by DSM and field measurements. The result shows that the simulated deposition zone is mostly located in the actual boundary indicated by red line and the part without deposition matches well with the exposed scarp at the upper edge of the landslide. It is proved that the numerical and measured results are generally in agreement.
Flow speed contours of the initial landslide at t 0, 2, 4, 6, 8 s are shown in the Figure 5. At t 2 s, the speed of landslide front reached a speed 7 m/s with effect of geopotential. At t 4 s, the front of the landslide was blocked, but the landslide mass on two sides still had a slow speed. After 4 s, the landslide mass began to FIGURE 6 | The maximum thickness and speed of the initial landslide during the sliding process.
Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 720872 deposit and the speed of mass on both sides decreased continuously. At t 8 s, the most mass became static except for some on the steep scarp. The maximum deposition thickness and maximum flow speed contours of the initial landslide are shown in the Figure 6. The speed almost increased and then decreased from top to bottom. The slope below the initial landslide suffered from the violent impact and compression of the initial landslide, resulting in two unstable mass and triggering the main landslide.
Deposition thickness contours of the main landslide at t 0, 8, 16, 24, 32 and 40 s are shown in the Figure 7. At t 8 s, the front of the landslide destroyed houses on both sides of the provincial road. At t 16 s, the front edge of the landslide moved to the bank and destroyed two houses near the bank. After 16 s, the landslide no longer moved forward and began to deposit mass affected by relatively flat terrain. At t 40 s, the most landslide mass is completely deposited. The final deposition area of the numerical simulation is compared with  The maximum deposition thickness and maximum speed of the main landslide are shown in the Figure 9. The speed almost decreased and then increased from top to bottom. The maximum speed was near the free surface on the upper edge of landslide. The high kinetic energy of the landslide destroyed houses in the fluidizaition zone.
A comparison between simulation result and terrain aquired from UAV along the main sliding profile A-A' and the transverse profiles is shown in Figure 10. For the main profile ( Figure 10A), the computed and measured results are generally in agreement, except for the location of provincial road and houses. The landslide mass on the road has been cleared and the destroyed buildings have increased the deposition depth, leading to an difference between the simulated results and the real depth. For the several transverse profiles ( Figures 10B-D), the computed results are generally consistent with the measured profile. In order to further quantify the accuracy of the simulation, the average error between the computed results and the real depth is calculated to be 2.2 m. It is proved that the deposition depth of simulation result matches with the actual depth very well and the depth-integrated continuum method can invert the dynamic process of the Zhonghai landslide.
In order to justify the parameters, numerical results for the initial landslide with μ 0.45 and 0.55 are presented in Figures  11A,B. It is seen that the mass moves downstream further with μ 0.45 and the runout distance is shorter with μ 0.55. Obviously, the friction coefficient of 0.5 is the most reasonable. In addition, numerical results for the main landslide with μ 0.27 and 0.31 are presented in Figures 11C,D. Although the front edge of the landslide is well matched with the real result relatively well, there is little mass in the upper part of the landslide in Figure 11C. Correspondingly, the deposition depth in the source zone of the main landslide and measured result are generally in agreement, but the landslide runout distance is shorter in Figure 11D. Hence, it is necessary to use different friction coefficients for different zones. On the other hand, pore pressure commonly increases with the decrease of slope height, so the lower sliding mass has greater mobility. Friction coefficient of 0.5, 0.31 and 0.25 is used from top to bottom as a result of the increase of pore water pressure. Therefore, from the aspect, the distribution of friction coefficient is reasonable.

CONCLUSION AND DISCUSSION
The landslide occurred at Zhonghai village in Hanyuan County is comprehensively analyzed by field investigation, aerial image recognition and numerical modeling. The failure mechanism and movement characteristics are unique and worthy to be noticed. Some findings are summarized and discussed as follows.
1. For Zhonghai landslide, the two successive rounds of heavy rainfall was recognized as the main cause of its occurrence. Meanwhile, the steep terrain, loose soil and crashed rock mass with poor mechanical properties were its vulnerable conditions. Loose cultivated soil and sandy soil is beneficial to rainwater infiltration, which increases the weight of the slope. The weak mudstone intercalation is softened and argillized, under the longterm moistening of infiltration water, so its resistant shear strength is greatly reduced. The "steep-gentle-steep" terrain promoted the initiation of the landslide and lead a series of chain reactions. The initial landslide in the high elevation can trigger a serial of corresponding consequences. Thus, for risk evaluation of this kind of landslides, the chained reaction need attention.
2. The rainfall triggered landslides are a complicated process. Even if the landslides are located at the same area, the dynamic processes are apparently different. For the initial landslide, the speed almost increased and then decreased from top to bottom. Correspondingly, the back edge of the main landslide with good free surface condition and the fluidization zone with higher mobility had high speed as well as the middle part moved slowly. It brings a huge challenge for numerical modeling. Here, we adopted a scheme of variable friction coefficients to reveal the internal mechanism. Friction coefficient of 0.5, 0.31 and 0.25 is used from top to bottom and it obviously has the correlation with the pore pressure which commonly increases with the decrease of slope height. In the future, a unified model coupled with infiltration and pore pressure evolution might be a creditable attempt.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.