Dynamic Process and Mechanism of the Catastrophic Taihongcun Landslide Triggered by the 2008 Wenchuan Earthquake Based on Field Investigations and Discrete Element Method Simulations

The Taihongcun landslide, which was a remarkable geological disaster triggered by the 2008 Wenchuan earthquake, had a volume of about 2 × 106 m3 and killed about 23 people. Through detailed field investigations, basic information of topography, geological structure and stratigraphy for the landslide were acquired and key kinetic characteristics of the landslide were identified. On the basis of filed investigations, 2D numerical models with discrete element method (DEM) were established to simulate the kinematics and failure process of the landslide. To ensure the validity of the dynamic calculations, the free-field boundary condition was developed and introduced intro the DEM models. According to filed investigations and DEM simulations, the dynamic processes of the Taihongcun landslide can be divided into four phases: fragmentation, projection motion, scraping, and granular debris flow and accumulation. In addition, the parameter analysis showed that the particle bond strength had a significant influence on the runout distance and landslide debris morphology. Finally, the possible mechanism of the Taihongcun landslide was determined: a rock mass of poor quality provided the lithological basis for this landslide formation; a joint set J1 in the back scarp and a weak interlayer of carbonaceous slate and shale between the upper sliding mass and the bedrock formed the rupture boundaries of the upper source area; a strong seismic ground motion was the external excitation that triggered the destructive landslide event; additionally, hypermobility was caused by the high elevation and topographical conditions of the landslide.

The Taihongcun landslide, which was a remarkable geological disaster triggered by the 2008 Wenchuan earthquake, had a volume of about 2 × 106 m 3 and killed about 23 people. Through detailed field investigations, basic information of topography, geological structure and stratigraphy for the landslide were acquired and key kinetic characteristics of the landslide were identified. On the basis of filed investigations, 2D numerical models with discrete element method (DEM) were established to simulate the kinematics and failure process of the landslide. To ensure the validity of the dynamic calculations, the free-field boundary condition was developed and introduced intro the DEM models. According to filed investigations and DEM simulations, the dynamic processes of the Taihongcun landslide can be divided into four phases: fragmentation, projection motion, scraping, and granular debris flow and accumulation. In addition, the parameter analysis showed that the particle bond strength had a significant influence on the runout distance and landslide debris morphology. Finally, the possible mechanism of the Taihongcun landslide was determined: a rock mass of poor quality provided the lithological basis for this landslide formation; a joint set J 1 in the back scarp and a weak interlayer of carbonaceous slate and shale between the upper sliding mass and the bedrock formed the rupture boundaries of the upper source area; a strong seismic ground motion was the external excitation that triggered the destructive landslide event; additionally, hypermobility was caused by the high elevation and topographical conditions of the landslide.

INTRODUCTION
On May 12, 2008, a magnitude Mw 7.9 earthquake occurred in the Longmenshan Mountains tectonic zone of Sichuan Province, China. The earthquake's epicenter was located at Yingxiu Town, Wenchuan County (31. 0°N and 103.4°E), and the hypocenter was at a depth of approximately 19 km. Surface displacements produced by the fault rupture reached several meters, propagating from the epicenter for approximately 240 km along the mountain range (Fan et al., 2018). According to the statistical data provided by the Chinese government, this earthquake claimed the lives of more than 87,000 people and injured approximately 459,000 others (Qi et al., 2010). More than 60,000 rock avalanches and landslides were triggered by the 2008 Wenchuan earthquake (Gorum et al., 2011;Wang et al., 2014), and several hundred landslide dams were formed.
Landslides triggered by strong earthquakes have always been a popular topic of discussion. The dynamics and failure mechanisms of earthquake-induced landslides are critical to elucidate for improving hazard prevention and mitigation measures. Field investigations indicate that the kinematic processes of earthquake-induced landslides are characterized by high speeds, long runouts, and high mobilities Yin et al., 2009;Chigira et al., 2010;Dai et al., 2011;Gorum et al., 2011;Tang et al., 2011;Huang et al., 2012;Zhang and Yin, 2013;Zhang et al., 2016). Considerable kinetic energy, high toe to rupture surface elevations, elastic energy release due to grain fragmentation, and high excessive pore water pressure generated by undrained loading were underscored as the principal controlling factors for the long travel distance, high speed, and high mobility Chigira et al., 2010;Zhang and Yin, 2013;Wang et al., 2014). It was concluded that the occurrence of an earthquakeinduced landslide is controlled mainly by seismic, terrain, and geological factors (Huang et al., 2012). Strong ground motion generated by the earthquake and topographic amplification of seismic waves due to features of the terrain were emphasized as the main triggers of rock avalanches and landslides (Havenith et al., 2003;Murphy, 2006;Bourdeau and Havenith, 2008;Yin et al., 2009;Huang et al., 2012;Zhou et al., 2013).
Field investigation and numerical simulation are two fundamental techniques for analyzing the dynamics and revealing the possible mechanism of individual landslides induced by strong earthquakes. Through field investigation, detailed information about a single landslide can be obtained, such as its geological and topographical features (Has and Nozaki 2014;Wang et al., 2014;Zhang et al., 2016;Zhuang et al., 2018). Furthermore, field investigation can provide evidence of the dynamic processes for landslides (Zhang et al., 2011, Zhang et al., 2016Xu et al., 2013). Numerical simulation is an important approach for modelling the dynamic process of rock avalanches and exploring the formation mechanism of earthquake-triggered landslides. At present, numerical simulation methods for the seismic response analysis of rock or soil slopes are mainly divided into two categories. The first category is composed of traditional continuum methods, including the finite element method (FEM) (Charatpangoon et al., 2014;Che et al., 2016) and finite difference method (FDM) (Bouckovalas and Papadimitriou, 2005;Bourdeau and Havenith, 2008;Zhou et al., 2013). The continuous simulation methods are based on the assumption of small deformation and displacement; therefore, using this type of approach, it is difficult to describe large displacements along discontinuities and rotations of blocks in a jointed rock mass under seismic loading. The second category consists of discontinuum methods such as discontinuous deformation analysis (DDA) (Wu and Chen, 2011;Zhang et al., 2015;Huang et al., 2016;Fu et al., 2019) and the discrete element method (DEM). With an explicit scheme-based force balance solution, the DEM exhibits a high computational efficiency in the seismic analysis of largescale rock mass engineering. In the DEM, blocks or particles resulting from different discretizations of the domain of interest can lead to different approaches. The most popular block-based DEMs are the commercialized UDEC and 3DEC for solving twoand three-dimensional (2D and 3D) problems; these programs have been extensively used in seismic response analyses for rock engineering (Kveldsvik et al., 2009;Cui et al., 2016;Song et al., 2018). The particle-based DEM method introduced by Cundall and Strack (1979) simulates the mechanical behavior of a material by representing it as an assemblage of discs (2D) and spheres (3D). The DEM was originally used to investigate granular materials such as soils/sands and was eventually implemented in a commercial code, i.e., Particle Flow Code (PFC) (Itasca Consulting Group Inc., 2008). To solve problems related to a rock mass, a bonded particle model (BPM), in which intact rocks are represented as assemblages of cemented rigid particles, was developed (Potyondy and Cundall, 2004), and a smooth joint contact model (SJCM) was proposed to model the joints (Ivars et al., 2011).
As previously documented, many earthquake-triggered landslides are characterized by flow-like movement, and the sliding material disintegrates into fragments during the highspeed movement, ultimately forming a granular deposit. The particle-based DEM method simulates the geomaterial as assemblies of rigid discs (in 2D) or spheres (in 3D) bonded at their contacts and allows finite displacements, rotations and separations of the discrete particles. The failure process in the synthetic material can be represented by the breakage of bonds and propagation and coalescence of these microcracks. Therefore, the particle-based DEM is more appropriate than DDA or the continuum methods for simulating the dynamic failure process and kinematic characteristics of landslides triggered by earthquakes. Huang et al. (2012), Tang et al. (2013), Zhou et al. (2013), Yuan et al. (2014), and Deng et al. (2017) developed 2D DEM models to simulate the dynamic process of landslides triggered by earthquakes and explore their seismic failure mechanism. In their simulations, the sliding surface and bedrock were represented by wall elements, and the seismic ground motion histories were applied to the wall boundaries. One drawback of this approach is that the wall boundaries cannot absorb the scattering waves from the upper sliding body, and the radiation of waves from the sliding mass into the bedrock cannot be modeled. Therefore, nonreflecting boundary conditions instead of wall boundaries are required to simulate the failure process of earthquake-induced landslides.
In this study, detailed geological and geomorphological surveys were conducted for the catastrophic Taihongcun landslide triggered by the 2008 Wenchuan earthquake; this landslide is located in Chenjiaba town, Beichuan County, Sichuan Province, China. From the detailed field investigation, the main characteristics of the Taihongcun landslide were identified. Subsequently, 2D DEM models were established to simulate the kinematic and failure processes of the Taihongcun landslide through the PFC2D code. To ensure the validity of the dynamic calculations, the free-field boundary condition was developed and introduced into the DEM models. The effects of the vertical and horizontal seismic forces on the dynamic failure process of this landslide were considered. Finally, based on field investigations and numerical simulations, we discussed the possible mechanism contributing to this landslide event.

GEOMORPHOLOGICAL AND GEOLOGICAL SETTINGS
The catastrophic Taihongcun landslide occurred on the eastfacing slope of Shuijingyan Mountain (31°58′25.4″N, 104°35′55.6″E, see Figures 1, 2), which is located northeast of Taihong Village, Chenjiaba Town, Beichuan County in Sichuan Province, China. The landslide, triggered by the 2008 Wenchuan earthquake, killed approximately 23 people, destroyed one village, and blocked the Duba River. The Taihongcun landslide lies on the right bank of the Duba River, originating from the north-south trending ridge of the Shuijingyan mountain ( Figure 2B). Due to the longterm river erosion, the west bank is concave, causing the slope near the bank to be unstable ( Figure 2A). The geomorphology of the Taihongcun coseismic landslide is chair-shaped and can be divided into two parts by the existence of a notable cliff on the edge of the upper sliding bed ( Figure 3A). The upper source area is between 960 and 1,300 m in elevation, with an estimated volume of 1.5 million m 3 , and the original slope was between 20°and 30°(   Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 710031 6 thickness of 40 m; the deposits are gray, with moderately weathered slate fragments present at the surface. The lower cliff induced by the secondary landslide is between 780 and 960 m in elevation and 400 m in width, with an average slope of approximately 50°. Field investigations show that the Yingxiu-Beicuan fault is located at the toe of the landslide (Figure 3). The floodplain of the Duba River is below the lower cliff, and the depth from the bottom of the river to the bank is 35-41 m. The main sliding direction along line I-I′ in Figure 3B is 110°. Figure 3C displays the geological longitudinal section through the main slding direction I-I′ for the Taihongcun landslide. According to field investigations, the Taihongcun slope consists of the Cambrian Qiujiahe Formation (Єq), Cambrian Youfang Formation (Єy) and Silurian Hanjiadian Formation (Sh). The Cambrian Qiujiahe Formation is mainly composed of carbonaceous and siliceous slate ( Figure 5D); the Cambrian Youfang Formation is composed of extensively jointed sandstone ( Figure 5E); and the Silurian Hanjiadian Formation is composed of thin layers of argillaceous slate. The Cambrian Qiujiahe Formation outcrop, in which the slate is strongly weathered, is observed at the top of the back scarp, with a thickness of approximately 96 m. A 0.4 m wide tension crack is found at the top of the back scarp ( Figure 4A). The lower part of the back scarp consists of the Cambrian Youfang Formation, in which the sandstone is intercalated with limestone ( Figure 4B), with a thickness of approximately 80 m. Karst erosion has developed along the preexisting fissures in the limestone, weakening the stability of the back scarp. The bedding planes of the strata around the back scarp of the slope gently dip away from the slope surface at an angle of approximately 10°. Two dominant and steep joint sets (J 1 , J 2 ) are widely developed in the back scarp ( Figure 4B). Set J1 dips away from the back scarp surface at an angle of approximately 80°and controls the surface deformation of the back scarp.
The displaced material covering the sliding bed of the upper source area is mainly composed of gray slate blocks and fragments, as observed on the surface. The bedrock of the upper source is the extensively jointed sandstone of the Cambrian Youfang Formation ( Figure 5E). Several gullies incised into the deposit by rainfall runoff enabled the internal structure of the deposit and sliding surface to be examined and sampled. As shown in Figures 5A-C, according to observations from a gully, a weak interlayer of carbonaceous slate and shale exists between the landslide deposits and the bedrock, with a dip angle ranging from 17°to 23°. The thickness of this layer is approximately 0.2-1.5 m. An obvious sliding surface can be seen between the displaced material and carbonaceous shale ( Figures  5A-C). Therefore, the carbonaceous shale is the strata on top of which the sliding surface occurred. The lower cliff was determined to have formed after the Wenchuan earthquake by comparing Figures 2A,B; this area was subsequently covered by Quaternary colluvial and eluvial deposits with an average thickness of 70 m. The basement rock of the cliff consists of the carbonaceous and siliceous slate of the Cambrian Qiujiahe Formation and sandstone of the Cambrian Youfang Formation.  Figure 6A. (C) Compact landslide deposits composed of angular fragments on the location shown as the red circle A in Figure 6A.
August 2021 | Volume 9 | Article 710031 Based on the field observations, the strata in the shallow part of the cliff are highly weathered. Due to the long-term strong tectonic activities of the Yingxiu-Beicuan fault, which lies across the foot of the cliff, the rock mass of the entire slope is fractured and is vulnerable to weathering. Thus, the strength of the slope rock mass is low. The Yingxiu-Beicuan fault is the interface between the Cambrian strata (in the hanging wall) and Silurian strata (in the footwall) ( Figure 3).

KINEMATIC CHARACTERISTICS OF TAIHONGCUN LANDSLIDE
As shown in Figure 3, the study area can be divided into three regions: the upper source area, the secondary landslide area, and the debris deposit area. The elevation difference between the toe and the crest of the slope in the upper source area is 260 m, and the corresponding horizontal distance is 460 m. The original slope of the upper source area failed due to the substantial seismic forces. After the landslide initiation in this area, the source mass detached from the parent rock mass and slid along the sliding surface at an increasing speed. Hightemperature carbonization observed near the sliding surface indicates that a high-speed shearing effect between the sliding body and the sliding surface was produced during the highspeed movement of the displaced mass ( Figures 5A-C). The initial velocity was obtained when the upper sliding body slid in front of the toe of the sliding surface and was projected into the air. After the failed rock mass was projected into the air from the toe of the upper sliding surface, this mass would exhibit a parabolic airborne course, and the gravitational potential energy would be converted into kinetic energy in the falling stage. The secondary landslide area is characterized by the cliff below the upper source area. A large elevation difference of 200 m between the toe of the sliding surface and the cliff base allowed the airborne rock mass from the upper source area to acquire considerable momentum and thus to scrape the preexisting deposit overlaying the cliff. The loose deposit material was entrained by the rock debris due to the impact and pushed downwards with it.
In the deposit area, the granular debris flow reached the valley bottom below the 180 m high cliff, blocked the Duba River, and climbed the opposite bank of the river ( Figure 6A). The horizontal travel distance was approximately 1,200 m, and the vertical fall height was approximately 550 m. The thickness of landslide deposits is affected by terrain. The deposits in the gullies and rivers are thickest, and the average thickness is approximately 50 m. Specifically, the deposits are the thickest at the landslide dam, with a thickness of up to approximately 95 m. The thickness of the debris flow deposits along the river course is generally approximately 30 m, while the thickness at the edge of the accumulation area is only 2-5 m. Field investigations show that the deposits are composed mostly of highly fragmented sandstone and shale; the grain size is in the range of sand, gravel, and cobble, with some large boulder-sized blocks ( Figure 6C). The main grain size is gravel and cobble, and the majority of the deposits are subangular to angular. As shown in Figure 6B, in the front of the deposits on the opposite bank, wheat stalks were pushed down and arranged in radial lines by the powerful air blast induced by the displaced material, indicating the high-speed movement of the granular debris flow.

DISCRETE ELEMENT METHOD MODELING
The dynamic processes of landslides under seismic loading conditions are very complicated. The DEM is well suited for simulating the cracks and subsequent large displacements of catastrophic landslides or rock avalanches and is therefore a useful tool for analyzing their complex kinematics and failure mechanisms. In this section, based on the two-dimensional code PFC2D, DEM models were used to investigate the kinematical process and failure mechanism of the Taihongcun landslide triggered by the 2008 Wenchuan earthquake.

Free-Field Boundary Condition for Discrete Element Method
In a numerical analysis of the seismic response of rock slopes, earthquake waves such as plane body waves are usually applied at the bottom of the model and propagate upward. The boundary conditions at the sides of the model should account for the freefield motion that would exist in the absence of a structure, such as a slope or a tunnel. The free-field boundary algorithm proposed by Zienkiewicz et al. (1989) is selected to "enforce" the free-field motion in such a way that the boundaries retain their nonreflecting properties (i.e., outward waves originating from the structure are properly absorbed). In this approach, the lateral boundaries of the main grid domain are coupled to the free-field domain by viscous dashpots to simulate an absorbing boundary, and the unbalanced forces from the free-field grid are applied to the main grid boundary. The free-field grid is solved in parallel with the main grid. For a lateral boundary of the main grid domain with its normal in the direction of the x-axis, the free-field boundary condition for 2D problems can be expressed as follows: where σ x and τ xy are the normal and shear components of the traction applied to the lateral boundary, respectively; _ u m x and _ u m y are the x and y components of the main grid point velocity at the lateral boundary; _ u ff x and _ u ff y are the x and y components of the velocity of the grid point in the free-field domain; and _ u ff x and _ u ff x are the tractions contributed from the free-field grid in the x and y directions.
Eqs 1, 2 are defined based on the continuum condition. To apply the free-field boundary condition in the DEM models, Eqs 1, 2 need to be discretized. For regular packing geometries, it may be possible to develop analytical expressions of the equivalent discrete free-field boundary condition based on Eqs 1, 2. For a square array of disks (each particle has four neighbors) of unit Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 710031 9 thickness in 2D, the particles in the assemblies have a common diameter D b . The regular square packing is subjected to an unconfined compression test. The axial stress, k n,SJ , is related to the axial force, k s,SJ , as follows: Substituting Eq. 3 into Eqs 1, 2, the free-field boundary condition for the DEM in 2D is obtained by where f x and f y are the normal and shear components of the force applied to the particles at the lateral boundaries of the main grid domain; f ff x and f ff y are the forces of a particle in the free-field domain in the x and y directions; and β 2p and β 3s are the nondimensional calibration factors for optimum wave absorption. Figure 7A shows the 2D DEM model of the Taihongcun landslide, constructed based on the geological longitudinal section shown in Figure 3C. The sliding mass in the upper source area can be detailed in Figure 7B, and the secondary landslide area in the model is shown in Figure 7C. To capture the main features of the large rock avalanche, the following assumptions were made:

Numerical Model
1. The sliding mass in the source area is cut by two sets of discontinuities, including joint set J 1 and the bedding planes. 2. Given the radius difference between the bedrock and the sliding mass, an 8 m thick transition layer is produced between the bedrock and sliding mass (see the enlarged view in the box in Figure 7B). 3. The rupture surface, including the back scarp and the sliding surface, is preset along the interface between the transition layer and sliding mass. 4. The bedrock below the sliding surface is considered homogeneous. 5. To avoid a large computational burden, the bedrock under the riverbed and east bank are simulated by wall elements. The seismic ground motion of the wall elements is in accordance with the movements of disk A (see Figure 7D).
The sliding mass in the upper source area is composed of slate and sandstone (see Figure 7B). In PFC2D, a joint is identified as a contact between particles lying on opposite sides of the joint plane. Joint plane contacts follow a contact model with properties that are different from those assigned to particle contacts elsewhere in the rock mass.
The most commonly used method for the generation of rock joints is to set the bond strengths of the joint plane contacts to low values (or zero for the lowest shear strength) (Kulatilake et al., 2001;Wang et al., 2003;Huang et al., 2015). The shortcomings of this method are that the surface of a generated joint has an inherent microscale roughness that is due to the circular shape and the nonuniform size and distribution of the particles (Bahaaddini et al., 2015). Therefore, it is difficult for the particles to slide along the discontinuities, especially along the sliding surface. To solve this problem, the smooth joint contact model (SJCM) was developed and introduced into PFC to simulate joints on a mesoscopic scale (Ivars et al., 2011). The SJCM simulates the behavior of an interface regardless of the local particle contact orientations along the interface. The behavior of a frictional or bonded joint can be modeled by assigning the SJCM to all contacts between particles that lie on opposite sides of a joint. With this approach, the shear behavior and mechanisms of rock joints under direct shear tests were studied in detail (Bahaaddini et al., 2013(Bahaaddini et al., , 2015Gutiérrez-Ch et al., 2018). In this study, the J 1 joint set, bedding planes, sliding surface and back scarp are modeled with the SJCM ( Figure 7B). Thus, the sliding mass in the upper source area can slide or separate along the sliding surface and back scarp.
The dynamic boundary conditions and monitoring disks for the DEM model are shown in Figure 7D. To absorb outward propagating waves and avoid plane waves propagating upward from being distorted at the lateral boundaries of the model, the free-field boundary conditions described in Free-Field Boundary Condition for Discrete Element Method are applied to the model. To ensure that the free-field domain has the same elastic behavior as the main grid domain, the microscale properties of the freefield domains are the same as those of the bedrock. The free-field domain is composed of four soil columns. The soil columns FFD-I and FFD-II are located on the left of the main grid domain, and the soil columns FFD-III and FFD-IV are located on the right. Note that only a seismic P-wave is applied on the bottom boundaries of soil columns FFD-I and FFD-IV and that only a seismic S-wave is applied on the bottom boundaries of soil columns FFD-II and FFD-III, while seismic P-and S-waves are applied to the bottom boundary of the main grid domain. When only an S-wave (or a P-wave) is required as input for the model, the free-field domain is merely composed of soil columns FFD-II and FFD-III (or FFD-I and FFD-IV). Viscous boundary conditions are applied to the bottom boundary of the main grid and free-field domains.

Model Parameters and Seismic Input
PFC2D simulates the intact material as an assembly of circular particles bonded at their contacts with the contact bond model or parallel-bond model and is thus referred to as a bonded particle model (BPM) (Potyondy and Cundall, 2004). The macroscopic mechanical behavior of the synthetic material depends on the microscopic parameters of the particles and bonds between particles. However, the microscopic parameters cannot, in general, be related directly to a set of relevant material properties (Itasca Consulting Group Inc., 2008). In general, the microscopic parameters can be obtained by means of a calibration process in which a series of numerical tests are performed to reproduce the desired macroscopic properties measured in laboratory experiments. In this study, numerical uniaxial compressive strength tests are carried out on BPMs to calibrate microscopic parameters against the uniaxial compressive strength UCS, elastic modulus E and Poisson's ratio v. The derived microscopic parameters of the BPMs for the landslide are listed in Table 1. A comparison of the macroscopic properties between the BPMs and laboratory experiments is conducted. The results from the BPMs are very close to those from the laboratory experiments. Numerical direct shear tests under constant normal stress are employed to determine the strength and stiffness parameters of the SJCM interfaces. The SJCM is applied to the contacts between particles that lie on opposite sides of the shear plane. To simulate the direct shear tests, the lower block is kept stationary during the tests, and a horizontal velocity V, which is low enough to ensure that the sample remains in quasistatic equilibrium, is applied to upper walls. The normal stress is applied vertically to the upper block, and this normal stress is kept constant during the tests by using a servo-control mechanism (Itasca Consulting Group Inc., 2008). The shear and normal displacements are obtained by tracing the horizontal and vertical displacements of the upper block, respectively, and the shear stress is calculated by dividing the average force reaction on the upper left and lower right walls by the joint length. The calibrated SJCM parameters for the sliding surface, bedding planes and joint set J 1 are listed in Table 2.
Because no strong motion station is present at the landslide site, records of the Wenchuan earthquake waves from the strong  (Figure 7) are constructed based on a geological longitudinal section ( Figure 3C), it is necessary to project the horizontal acceleration records in the east-west and north-south directions into the main sliding direction of the landslide using the following equation: where a H is the acceleration in the main sliding direction, a E−W is the acceleration in the east-west direction, a N−S is the acceleration in the north-south direction, and 110°is the azimuth of the main sliding direction. The ultimate acceleration histories after filtering and baseline correction are shown in Figure 8. The velocity time-history records are transformed into force time-history records and applied to the bottom boundaries of the main grid and free-field domains.
To effectively dissipate the kinetic energy of the assembly system, numerical damping must be used in the simulations. Local damping and viscous damping are available in DEM models to dissipate kinetic energy. Local damping applies a damping force, with a magnitude proportional to the unbalanced force, to each particle, while viscous damping adds normal and shear dashpots at each contact. Local damping is not applicable for landslide simulation (Deng et al., 2017). Therefore, viscous damping is used in the simulations. The parameters of viscous damping can have a significant impact on the results, but they cannot be explicitly related to any physical mechanism. The parameters of viscous damping in DEM models can also influence the grain motion and reflect energy dissipation during collision. The damping parameters are correlated with the coefficients of restitution, which can be measured in the laboratory or field. The results of field restitution coefficient tests conducted by Giani (1992) are used to define the damping parameters. The values of the normal and shear damping coefficients are thus set to 0.4 and 0.2, respectively, (Deng et al., 2017).

NUMERICAL SIMULATION RESULTS AND ANALYSIS
The dynamic process of the Taihongcun rock avalanche can be classified into four stages: upper sliding mass fragmentation, projection of the upper displaced mass, scraping of the upper sliding mass toward the Quaternary colluvial and eluvial deposit, and granular debris flow and accumulation. Each of the stages is discussed in detail in Upper Sliding Mass Fragmentation-Granular Debris Flow and Accumulation. The numerical parameter of internal bonding strength is the most crucial factor influencing landslide transport and deposition in the DEM model . Therefore, the influence of the parallel-bond strength on the runout distance and deposition patterns is investigated in Influence of the Bond Strength on the Runout Distance and Deposition Patterns. Upper Sliding Mass Fragmentation Figure 9 shows the progressive fragmentation of the upper sliding body under the combined horizontal and vertical seismic loadings; the cluster plots indicate fractured blocks (Itasca Consulting Group Inc., 2008). In the static gravity equilibrium situation ( Figure 9A), damage can be observed in the rock mass near the scarp and sliding surface and on the top part of the slope. Cracks or fissures develop along the bedding planes and joint set J 1 in the damaged regions, and the rock mass is broken into small blocks. The entire model can be brought to equilibrium under gravity so that the slope is stable before the earthquake. After the seismic force is applied to the base of the model, the slope is not stable (Figures 9E-H). In the initial stage of the Taihongcun rock avalanche, a rock slide occurred, and an obvious settlement is observed near the top of the slope when the shear stress at the sliding surface reaches the shear strength ( Figure 9B). The development of cracks or fissures in the rock mass, induced by the shaking and sliding, is observed ( Figures 9E-H). The high-speed shearing effect between the sliding body and the sliding bed resulted in fragmentation of the rock blocks into small particles ( Figures 9B,C). The jointed rock blocks in the slope begin falling and fragment into smaller sizes as a result of the internal forces exceeding the bonding strength ( Figures 9E-H). Under strong seismic loadings, tensile cracks develop in the rear part of the sliding body, forming a tensile stress concentration zone ( Figures 9E-D). Then, the cracks widened laterally and propagated vertically to the intersection between the scarp and the sliding surface, generating horizontal and vertical displacements. The large elevation difference between the slope and the slip plane facilitated the rockslide, triggered by the violent quaking.

Projectile Motion of the Upper Displaced Mass
Under seismic and gravity forces, the upper sliding mass achieved a high initial speed of approximately 40 m/s over a short period as it slid away from the toe of the sliding surface ( Figure 10A).
There is sufficient open space in the movement direction of the landslide, allowing the upper sliding mass to run out from the top of the 180 m cliff at a high speed, along a parabolic path ( Figures  9E-G, Figure 10A). After the upper sliding mass was projected into the air, the speed increased rapidly during the vertical drop from the toe of the sliding surface to the valley bottom. The high initial velocity, high elevation of the toe of the sliding surface and sufficient open space in the movement direction induced projectile motion. Scraping of the Upper Sliding Mass Figure 10B shows the scraping of the upper sliding mass into the Quaternary colluvial and eluvial deposit on the lower cliff. As the gravitational potential energy would be converted into kinetic energy in the falling stage, with a substantial kinetic energy, the thrown upper sliding mass collided with and intensely scraped off the Quaternary deposit, completely entraining and removing the deposits. The bedrock below the Quaternary deposits was completely exposed after the event, and the lower cliff formed. The lower part of the landslide was triggered into a secondary failure due to the earthquake activity and the effects from the upper part of the landslide. Figure 11 shows the velocity histories of the monitoring particles at different locations in the upper sliding mass and the lower Quaternary deposit (shown in Figure 7D). As shown in Figure 11A, the maximum velocity of the monitoring particles is approximately 60 m/s, which was recorded at point 5, and the peaks of the recorded velocities decrease from the front of the upper sliding mass to the rear. Particles in the rear of the upper sliding mass (particles 1, 2, and 3) start moving slowly, and their movement continues for a short time. However, particles in the front of the upper sliding mass (particles 4 and 5) start moving quickly, and their movement continues for a longer period. The peaks in the recorded velocities in the front of the upper sliding mass (particles 4 and 5) appear later than those in the rear of the upper sliding mass (particles 1, 2, and 3). The seismic response of the upper sliding mass is most intensive at t 6-50 s, and the direct cause of the landslide is the earthquake. As illustrated in Figure 11B, monitoring particle 6 in the top of the lower Quaternary deposit was displaced by the thrown upper sliding mass at approximately 13 s. The particles in the top of the lower Quaternary deposit were displaced earlier and farther than those in the bottom. Therefore, the velocity peaks of the monitoring particles decrease from the top of the lower Quaternary deposit to the bottom. The speeds of the monitoring particles in the lower Quaternary deposit peaked at approximately 30 s, then gradually decreased, and reached zero at approximately 43 s. The upper sliding mass was disintegrated into fragments due to the strong vibration, intensive collision, and violent scraping action. Then, the fragmented geomaterial transformed into granular debris with a rapid, flow-like motion. The lower Quaternary deposit was entrained due to the rock debris collision on the slope surface. The topography after sliding is shown in Figure 12A, which indicates that there are two deposition zones of the Taihongcun landslide: the first is in the upper source area, where fragmented rock debris accumulates on the sliding surface; the second is in the valley of the Duba River. In the first deposition zone (Figure 12B), the average thickness of the deposit is approximately 40 m, which is in accordance with the field observations, and the deposit can be divided into two layers: an upper layer of slate and a lower layer of sandstone. In the second deposition zone, the granular debris flow reached the valley bottom below the 180 m high cliff, blocked the Duba River, and climbed over the opposite bank of the river ( Figure 12C). The second deposition zone has a maximum thickness of approximately 35 m on the opposite bank, and the greatest height of the landslide dam in the Duba River valley is approximately 95 m. The landslide deposit in the second deposition zone is mainly composed of sandstone and Quaternary deposits. The Quaternary deposits are covered by sandstone debris. The distance between the summit of the mountain before the earthquake and the edge of the deposit on the opposite bank is approximately 1,200 m in the horizontal direction, and the vertical height from the summit of the mountain to the bottom of the Duba River is approximately 600 m.

Influence of the Bond Strength on the Runout Distance and Deposition Patterns
Some researchers noted that the basal friction coefficient and the internal bonding strength were the most crucial factors influencing landslide transport and deposition in PFC2D models (Tang et al., 2009;Deng et al., 2017;Zhang et al., 2017). The influence of the basal friction coefficient on the runout distance and deposition patterns can be understood easily. However, the influence of the internal bonding strength on the runout distance and deposition patterns is Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 710031 disputed and poorly understood. A series of numerical simulations that were performed to investigate the runout distance and deposition area of the Taihongcun landslide in relation to internal bond strength are presented in this section. Three cases for the parallel-bond normal strength of the particles are used in the numerical simulations: Case 1. The parallel-bond normal strength for sandstone is 4 MPa, and the parallel-bond normal strength for slate is 1.6 MPa.
Case 2. The parallel-bond normal strength for sandstone is 12.5 MPa, and the parallel-bond normal strength for slate is 5 MPa.
Case 3. The parallel-bond normal strength for sandstone is 25 MPa, and the parallel-bond normal strength for slate is 10 MPa.
Figures 13A-C shows the final topography of the first deposition zone for the three cases, and Figures 13D-F shows the final topography of the second deposition zone for the three cases. As depicted in Figure 13, the final topography of the deposits in the simulations indicate that the internal bonding significantly influences the resultant runout distance and debris flow avalanche deposit. Figure 13 show that Case 1 results in the minimum thickness of the first deposition zone, the maximum thickness of the second deposition zone, and the longest runout distance, while Case 3 results in the maximum thickness of the first deposition zone, the minimum thickness of the second deposition zone and the shortest runout distance. Moreover, the deposits for Case 1 are the densest of the three cases. The series of simulation results indicates that a lower bond strength enables the upper sliding mass to more easily fragment (Figure 13A), while a high bond strength will cause the main body mass to behave as a rigid block, thus hindering fragmentation ( Figure 13C). With respect to the back analysis of the runout distance and debris accumulation topography, Case 2 provides the best simulation results for the Taihongcun landslide dynamics process.

FORMATION MECHANISM ANALYSIS
Based on field investigations and discrete element simulations, we conclude that the occurrence of the Taihongcun landslide is the result of a combination of seismic, terrain, and geological factors. The formation mechanism of the Taihongcun landslide will be discussed in this section.

Weak Rock in the Slope
Rock mass quality is the basis for the stability of rock slopes. The strength and deformation modulus of a rock mass with a good quality are high, which is beneficial to the slope stability, while that with a poor quality are low, which challenges slope stability. The rock mass in the study area was of poor quality due to the shearing and brecciation associated with past fault movement, rendering the rock slopes prone to failure and contributing to the high degree of disintegration of the failed rock mass. The strength of the Cambrian sandstone and slate in the upper source area is low, and the rocks are vulnerable to weathering, leading to fracturing of the upper sliding mass. The weak rock mass of the slope is the lithological basis for this landslide formation.

Unfavorable Geological Structure
Based on the field investigations, there exists a group of joints J 1 , which steeply dip outward, in the residual slope of the upper source area, and the back scarp developed along this joint set. As observed in the outcrop ( Figure 5), a weak interlayer composed of carbonaceous slate and shale is present. Under strong ground motion induced by earthquakes, the slope is easily shears along the weak interlayer, and a large-displacement slide can occur. The joint set J 1 and the weak interlayer formed the rupture surfaces of the upper source area.

Strong Seismic Ground Motion
The seismic intensity in the study area is more than X (Chinese seismic intensity) (Huang and Li 2009), and the seismogenic Yingxiu-Beicuan fault passes through the foot of the valley slope. The original PGA in the area is greater than 600 Gal, a seismic intensity greater than X . Furthermore, the main source area is located in the hanging wall of the fault, and the hanging wall effect in the near-fault ground motion Huang and Li, 2009) may be obvious in the study area. The strong seismic ground motion is an external excitation that triggers the destructive landslide event.

High Elevation and Topographical Conditions
The vertical distance from the toe of the upper sliding surface to the valley bottom is approximately 300 m. Therefore, the upper sliding mass can gain substantial kinetic energy, transformed from the gravity potential energy, in the falling stage. There is enough transport space downslope to facilitate the projection of the upper sliding mass into the air. Three faces of the Shuijingyan mountain, where the landslide originated, are free, and the terrain is steep, high and convex. The topographic amplification effect on seismic waves is strong at or near ridge crests in this kind of topography, which favors slope failure. The high elevation and topographical conditions provided the basis for hypermobility.

CONCLUSION
This paper focuses on the dynamic process and formation mechanism of the catastrophic Taihongcun landslide triggered by the 2008 Wenchuan earthquake. Through detailed field surveys, basic information of the topography, geological structure and stratigraphy for the landslide was acquired, and key kinetic characteristics of the landslide were identified. Based on field investigations, 2D DEM models were established to simulate the kinematics and failure process of the landslide.
To ensure the validity of the dynamic calculations, the freefield boundary condition was developed and introduced into the DEM models. The modeled scenario had a reasonably good fit to the actual topography, adding credibility to the modeling results. Based on the field investigations and numerical simulations, the dynamic processes of the Taihongcun landslide was divided into four phases: upper sliding mass fragmentation, projection of the upper displaced mass, scraping of the upper sliding mass toward the Quaternary colluvial and eluvial deposit, and granular debris flow and accumulation. The progressive fragmentation and projectile motion of the upper sliding body under the combined horizontal and vertical seismic loadings were effectively simulated. The mechanism of the Quaternary deposit on the lower cliff scraped by the upper sliding mass was revealed through these simulations. In addition, a parameter analysis showed that the internal bond strength of the particles had a significant influence on the runout distance and landslide debris morphology.
According to field investigations and numerical simulations, a possible mechanism of the Taihongcun landslide was determined: the rock mass in the study area was of poor quality and was vulnerable to weathering, which provided the lithological basis for this landslide formation; the joint set J 1 , which steeply dips outward in the residual slope, and the weak interlayer, which was comprised of carbonaceous slate and shale, formed the rupture boundary of the upper source area; the strong seismic ground motion was an external excitation that triggered the destructive landslide event; additionally, the high elevation and topographical conditions provided the bases for hypermobility.
Notably, the seismogenic Yingxiu-Beicuan fault that ruptured during the Wenchuan earthquake lay across the foot of the valley slope. The role of oblique-thrust faulting in the initiation of earthquake-induced landslides needs to be quantitatively assessed. However, it is difficult to account for oblique-thrust faulting in our simulations at present. Therefore, further research is needed to consider the effects of oblique thrusting on the initiation of near-fault landslides in such simulations.

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.

ETHICS STATEMENT
Written informed consent was obtained from the relevant individuals for the publication of any potentially identifiable images or data included in this article.