Earthquake Research in China  2020, Vol. 34 Issue (1): 96-109     DOI: 10.19743/j.cnki.0891-4176.202001005
A Discrete Element Analysis of the Sliding Friction Heat in High-Speed and Long-Runout Landslides
ZHU Yao1, ZHU Chengguang1, LIU Chun1,2, LIU Hui1, DENG Shang2
1. School of Earth Sciences and Engineering, Nanjing University, Nanjing 210046, China;
2. Nanjing University High-tech Institute at Suzhou, Suzhou 215123, Jiangsu, China
Abstract: Earthquake is one of the main causes of high-speed and long-runout landslides. Generally, the heat generated in the sliding zone is significant in such devastating landslides. In this study, we establish a two dimensional slope model which includes 0.2 million elements to simulate the development of high speed and long-runout landslides using the discrete element software MatDEM. The model not only suggests that heat is produced by friction and fracturing, but also simulates the process of tension generation in cracks and the generation of a high heat zone near the sliding region. Besides, the heat field graph indicates a banded high heat belt that is related to the location of the thickest sliding body. The logarithms of the total calorific value and the highest value in the heat zone during the sliding process are linearly related to the logarithm of the landslide height.
Key words: MatDEM     Discrete element method     Energy     Heat     Landslide

INTRODUCTION

Earthquake is one of the main causes of high-speed and long-runout landslides (Guo Changbao, 2018). These landslides attract attention of geologists because of the risk of secondary disasters (Tang Huiming et al., 2015; Chen Lujun et al., 2008). During the rapid sliding process, the high-speed and long-runout landslides are associated with shear heating, thermal pressurization, and thermal decomposition of carbonates at the slide shear zone. Under the influence of this mechanism, the generated heat may lead to the rise of pore pressure, which may further result in the reduction of frictional resistance to sliding (Goren L. et al., 2010). As a result, the sliding distance of some landslides may exceed the theoretical value due to the reduced sliding friction, with the speed reaching as high as 100 m/s. Examples are the Italian Vajont landslide and a series of landslides triggered by Wenchuan earthquake (Breien H. et al., 2008). High-speed and long-runout landslides accompanied by fluid in the rock mass can cause long-distance displacement of materials, which may run for several kilometers to flood cities. Therefore, it is of great significance to study friction heat generated by high-speed and long-runout landslides.

Investigating the characteristics of earthquake-induced landslides are beneficial for monitoring the landslide process as well as preventing future disasters. The destructiveness of landslides limits researchers' ability to verify these events with field-monitoring data. In addition to some equivalent bump tests (Wang Y. F. et al., 2017), most tests focus on torsional shear apparatus with the use of only small sample volumes. Therefore, scholars start using computer simulation methods (Huang Runqiu et al., 2008), such as element flow code, dynamic analysis, and computational fluid dynamics, to analyze earthquake-induced landslide mechanisms. In particular, some researchers combine the response surface methodology with computer simulation. For example, Fei (Yin Yueping et al., 2012) performs Monte Carlo simulation using ground-penetrating radar responses to calculate the system failure probability of slopes. Li (Hu Mingjian et al., 2015) reviews various applications of response surface methods in different slope reliability problems. Several mechanisms are proposed to interpret response surface, such as the landslide fluidization theory, gas lubrication theory, and friction weakening theory (Ge Yunfeng et al., 2016).

The high-speed torsional shear test (De Blasio F. V., 2008) is widely used to investigate the mechanisms of earthquake-induced landslides. In this study, the change in effective friction coefficient and temperature under high-speed shear is investigated. Wang Y. F. et al., (2017) obtain data on the temperature rise on the sliding surfaces of rock samples during high-speed shearing. Their results indicate that the temperature is mainly determined by the shearing rate. Dalguer L. A. et al., (2003) conduct shearing tests on rock samples to simulate energy distribution on inner planes of rock during sliding deformation. They also apply the DEM and used a thermal infrared imager to track heat distribution during shearing. The researchers find that the temperature on the shear surface increase acutely. Thus, a large amount of heat is generated during high-speed sliding in a short time due to high friction on the sliding surface. However, the heat cannot be dissipated in a similarly short period, and thus a sudden increase in temperature occurs around the sliding zone. Rocks and soils can melt at high temperatures (Zhang Lijiao et al., 2008); such phenomenon may reduce the strength of geomaterials and change the friction coefficients of rock and soil. High temperatures may also cause evaporation of pore water, resulting in formation of an air cushion at sliding interfaces (Habib P., 1975). Travelletti J. et al., (2014) monitor the variation of temperature and displacement during the process of an earthquake-induced landslide using a thermal infrared imager and a three-dimensional laser scanner; but the researchers only achieve limited success in monitoring internal zones of the slope. Rice J. R., (2006) observes that a considerable amount of heat is generated during faulting, which may reduce the shear strength of rock.

In this study, the sliding friction heat of high-speed and long-runout landslides are investigated by using the discrete element method (DEM). The model of a large-scale slope is built in the discrete element software, MatDEM. The dynamic landsliding process is simulated, and landslide slip characteristics are investigated based on the velocity field. The heat generated during the landsliding is quantified using a heat detection module. The influence of landslide height on the heat of the sliding zone is explored on a preliminary basis.

1 THE DISCRETE ELEMENT METHOD 1.1 The Contact Model

In the discrete element model, rock or soil is composed of a series of granular units that satisfy the Newtonian equation of motion. The units are affected by normal and shear spring forces (Jiang M. J. et al, 2003; Liu Chun et al., 2013, 2015). The system can be represented by the following expressions:

 ${F_{\rm{n}}}{\rm{ = }}{K_{\rm{n}}}{X_{\rm{n}}}$ (1)
 ${F_{\rm{s}}}{\rm{ = }}{K_{\rm{s}}}{X_{\rm{s}}}$ (2)

where Fn represents the normal force (N); Kn is normal stiffness of the bond (N/m); Xn is normal relative displacement of the spring (m); Fs is tangential force (N); Ks is shear modulus (N/m); and Xs represents shear relative displacement (m).

 ${F_{{\rm{nmax}}}}{\rm{ = }}{K_{\rm{n}}}{X_{\rm{b}}}$ (3)

When the normal force Fn exceeds the Fnmax in Eq.(3), a normal fracture is generated between two adjacent units. The maximum displacement between them is Xb. When the two elements are compressed again, the normal force between the elements is restored.

 ${F_{{\rm{smax}}}}{\rm{ = }}{F_{{\rm{s0}}}} - {\mu _{\rm{p}}}{F_{\rm{n}}}$ (4)
 $F{'_{{\rm{smax}}}} = - {\mu _{\rm{p}}}{F_{\rm{n}}}$ (5)

Where Fs0 is the interelement shear resistance (N), μp is the interelement coefficient of friction, and Fn is the normal force (compressive force is negative). The intact bond between two particles breaks when the external force exceeds the maximum shear force (Fsmax) determined in Eq.(5).

1.2 Calculation of Energy

The mechanical energy in the discrete element model is mainly composed of the kinetic energy, the elastic potential energy, and the gravitational potential energy. The element's elastic potential energy is defined by the stiffness and relative displacement of the spring between elements (Liu Chun et al., 2017):

 ${E_{\rm{e}}} = 0.5 \cdot {K_{\rm{n}}}{X_{\rm{n}}}^2 + 0.5 \cdot {K_{\rm{s}}}{X_{\rm{s}}}^2$ (6)

where Kn and Ks represent normal stiffness and shear stiffness, respectively, and Xn and Xs are their corresponding relative displacements. The gravitational potential energy (Eg) of the element is:

 ${E_g} = mgh$ (7)

where m is the element's mass (kg); g is the gravitational acceleration (g=-9.8 m/s2), and h is the height compared with a reference level. The kinetic energy of a cell can be expressed as follows:

 ${E_{\rm{k}}} = 0.5 \cdot m{v^2}$ (8)

Where m is the element's mass and v is the velocity (m/s).

In a discrete element numerical simulation, mechanical energy gradually transforms into heat of many forms: friction heatis generated when sliding friction occurs between the elements; fracturing heat is released when the connection between two elements breaks; and viscous heat is produced by the damping motion of the elements (Liu Chun et al., 2017). The total heat in this model is the sum of the viscous heat (Qv), fracturing heat (Qb), and frictional heat (Qf):

 $Q = {Q_{\rm{v}}} + {Q_{\rm{b}}} + {Q_{\rm{f}}}$ (9)
1.3 MatDEM Software

Based on the contact principle and energy relation presented, the high-performance discrete element software MatDEM is developed. By using an innovative matrix computing algorithm, the speed of discrete element computation is improved dramatically in the software, which may process millions of elements in one workstation. MatDEM exhibits strong pre- and post-processing capabilities and can automatically generate a variety of cloud figures and demos of the simulation results. Its primary technical advantages, compared with those commercial softwares, are the capacity of large-scale calculation, high computational efficiency, and the ability to calculate heat generation (Liu Chun et al., 2017). Furthermore, the software can be downloaded free from the website http://matdem.com, the source code of this study is included in the software.

The system can simulate soil deformation, fracture, destruction, and other behaviors using three-dimensional discrete element modeling. Besides, the software can be used to calculate the conversion of mechanical energy to heat and provides statistics regarding factors involved in heat generation. At present, the software is equipped with the functions of simulating formation of compaction bands (Liu Chun et al., 2015), three-dimensional landsliding (Scaringi Gianvito et al., 2018) and land subsidence (Liu C. et al., 2019), etc.

2 MODELING A SLOPE

A numerical model of slope dynamics is also formulated using MatDEM. The modeling process covers the establishment of a random accumulation model, formulation of a geometric model, and assignment of material parameters. In Fig. 1, a simplified model with less number of elements is used to introduce the modeling process. While in the following simulation process, a model consisting of 220 000 elements is used to simulate the process of sliding.

 Fig. 1 (a) Initial random-packing assemblage. (b) Additional gravity element deposition. (c) The coordinate figure in layers. (d) Cutting slope profile and giving layer characteristics to slope model
2.1 Random Accumulation Model

The initial model is constructed by simulating the natural deposition process of particles. As denoted in Fig. 1(a), a 220 m × 120 m rectangular box model is configured. The average particle size of the elements is 0.1 m. Elements are first distributed regularly within the region, and random initial velocities are assigned to the elements, which move to random position via numerical iterations. Then, with the addition of gravity, the elements fall down and accumulate in the box. Finally, a flexible band consisting of numerous particles is used to load the model. The particles are compacted repeatedly for the purpose of obtaining closely packed model elements. As a result, a two-dimensional packing model is built.

2.2 Geometry Modeling and Attaching Materials

The modeling of a slope in MatDEM is based on digital elevation models, which record the coordinates of a series of points of three-dimensional surfaces. Layers of a slope can be defined by these surfaces. Under the circumstances of two-dimensional models, these layer surfaces can also be defined by aline chart in the office EXCEL software. As shown in Fig. 1(c), the coordinate file of the line chart can be imported to MatDEM, and be further used to define layers of the slope. The geometric model of the slope is shown in Fig. 1(d), three layers are distinguished according to the four surfaces illustrated in Fig. 1(c), in which layer 2 is a weak layer with lower strength.

Two materials are defined in this model and the characteristic parameters are shown in Table 1. As shown in Fig. 1, the layer formation of A corresponds to material 1 and formation B corresponds to material 2. In detail, layer B is a weak layer corresponding to material 2 with lower strength. The element properties are determined by the soil properties and the conversion formulas (Liu Chun et al., 2017). In a discrete element-compacted packing model, there is an analytical solution between the mechanical parameters of the element and the mechanical properties of the overall model. Based on the conversion formula given by the linear elastic model (Liu Chun et al., 2017), the mechanical parameters of the discrete elements can be determined by the macro-mechanical properties of the rock and soil, and the specific parameters are shown in Table 2.

Table 1 Average mechanical properties of two soil layers of the model

Table 2 Mechanical properties of soil

When actual material properties are assigned to the elements, the element stiffness and the stress between elements may change, resulting in imbalance forces in the model. Therefore, the model requires a re-balancing operation before the numerical simulation. In the operation, the slope elements are glued and temporally assigned strong tensile strength to avoid failure. After hundreds to thousands of iterations, the imbalance forces between elements will be released.

3 NUMERICAL SIMULATION AND ANALYSIS OF EARTHQUAKE-INDUCEDLANDSLIDING 3.1 Numerical Simulation

According to the above modeling method, a two-dimensional slope model is created (Fig. 2(a)). The slope model has a height of 100 m and an average inclination angle of 60°, and contains approximately 212 000 elements. The radius of the elements ranges from 0.16 to 0.21 m, and the average radius is 0.18 m. The numerical simulation is conducted on a GPU workstation (Tesla GPU K20c). The slope slides occur along the weak surface because of gravity. The simulation of the sliding process takes 6.5 h, which corresponds to 12 s in actual situation. Through numerical simulation, the horizontal and vertical displacement, vertical velocity, shear stress, energy value, and heat distribution are obtained.

 Fig. 2 (a)-(c) Distribution of velocity field at 1, 6, and 12 s. (d)-(f) State of element bond at 1, 6, and 12 s
3.2 Sliding Characteristics

Figs. 2(a)-2(c) denote the horizontal velocity of the landslide element at 1, 6, and 12 s, respectively. The different colors in the figure indicate different velocity values. Fig. 2(a) represents the initial status of the landslide at 1 s. At the foot and waist, the maximum speed is observed (4 m/s). At 6 s (Fig. 2(b)), the horizontal velocity of the upper field increases, and the sliding bed becomes almost stationary. The landmass slides along the weak layer are at high speed, with the velocity at the waist and lower part of the slope reaching a maximum value of 15 m/s. At 12 s (Fig. 2(c)), the landslide speed gradually slows down such that the velocity at the trailing edge near the slip zone is 0.5-0.8 m/s. The sliding bed then stops.

The boundary of the weak layer (B layer) is outlined in Figs. 2(a)-2(c). As sliding proceeds, the weak layer of the leading edge becomes thicker and deforms due to the shift in landslide movement from downward sliding to horizontal and upward extrusion. The resistance force is increasingly mobilized, and the soil layer is deformed by compression.

Figs. 2(d)-2(e) illustrate the states of connection elements of the landslide at 1, 6, and 12 s. The black line delineates the boundary of the landslide. The green line in the figure represents the complete connection among the particles (cementation), while the blank area stands for the disconnected links. Fig. 2(d) presents the situation at 1 s at the waist and the foot of the landslide surface, in which numerous faults are evident. The soil disaggregates into numerous falling particles. Simultaneously, fracture zones are apparent in the weak layer. An upper weak layer is subsequently formed with the sliding body. Bonding between the elements of the sliding body gradually diminishes during sliding.

At 6 s (Fig. 2(e)), almost all the bonds are broken in soil in the following areas: foot of the slope, back edge of the slope, around the slip zone, and around the slope surface. Only a fraction of the middle region retains complete connection, and the connection decreases as sliding progresses. At 12 s (Fig. 2(f)), the entire earthquake-induced landsliding mass is essentially a loose accumulation of soil and there are some fractures along the slip zone (Figs. 2(b), (c)) aligned at approximately 30° to the surface (dotted line). These fractures are attributed to soil failure from the buildup of a substantial shear force on the surface of the slip zone. According to the Mohr-Coulomb failure criterion, maximum shearing occurs at approximately 30° alignment between the shear failure direction and the minimum principal stress direction. When the shear stress is large, tensile stresses and associated tensile fractures may appear in the direction of the minimum principal stress. The results of this simulation show that the investigated earthquake-induced landsliding on the slope is mostly due to the influence of gravitational forces on the soil mass. At the waist and foot of the slope, stress concentrations destroy the cohesion of the soil mass, acting as an early trigger of the landslide. With the progression of failure, a slip zone is formed. Finally, the back of the slope is destroyed, leading to the eventual collapse of the slope. Numerous particles fall into the foot of slope, corresponding to a typical gravity landslide (Parise M. et al., 2000; Xu Qiang et al., 2016).

3.3 Energy Transformation in Sliding Process

MatDEM provides data regarding variety types of energy, especially heat, in landslides. Fig. 3 shows the energy conversion curve of sliding. The discrete element model enables assessment of energy conservation. During the sliding process, gravitational potential energy gradually decreases and transforms into kinetic energy, which reaches a maximum value at 5 s and decreases gradually to a stable value after 10 s. Potential energy fluctuates during earthquake-induced landsliding due to the requirement of energy to establish a compensatory relationship of conservation between kinetic energy and elastic potential energy. Large vibrations promote soil structural adjustments that can disrupt cementation of grains. This further reduces soil density and enhances soil fluidity.

 Fig. 3 Energy conversion curves during sliding

A substantial amount of heat is generated during sliding. Fig. 4 shows the heat generation curves, including those for damping heat, frictional heat, and fracture heat. Frictional heat at the level of 5 × 108 J is generated by shear friction among soil particles, reflecting the same trend as total heat. In addition, gravitational energy converts to heat during sliding. The results show that the mechanical energy is mainly dissipated in frictional heat generation within 12 s. Frictional force is the dominant factor in heat generation on the sliding surface, as demonstrated by Liu Chun et al. (2017) on monitoring heat generation in a quartzite element stacking model.

 Fig. 4 Heat generation curves during sliding
3.4 Simulation of Heat Generation

Figs. 5(a)-5(d) shows the results of the application of the model for estimating heat distribution in a landslide at 1, 6, 9, and 12 s. The warm colors represent higher heat levels. In Fig. 5(a), heat is generated at the foot and waist of the slope at the initial stage of sliding. Heat induced by cracking and friction increases due to soil disaggregation, and sliding shear of the weak layer produces an obvious heat belt. At 6 s (Fig. 5 (b)), the sliding belt generates the most heat in a horizontal distance range of 40-80 m. As shown in Fig. 5(a), this region is in the low-to-middle part of the slip zone. During the sliding process, these elements in the zone slide the most and undergo maximum pressure, resulting in the highest heat generation and the formation of the corresponding heat belt. At 9 s (Fig. 5(c)), the heat belt moves forward, following the thick zone in the landslide. The heat belt slides from 40-80 m to 50-100 m horizontally. This movement corresponds to the trajectory of the maximum thickness of the sliding soil.

 Fig. 5 (a)-(d) Distribution of heat during sliding at 1, 6, 9, and 12 s

For the same displacement and friction conditions, greater amounts of pressure result in larger generation of heat. The landslide surface is somewhat circular such that the volume of soil in the middle of the sliding body is the largest, whereas the volume in the foot and head is small. Soil volume directly determines the pressure around the slip zone as well as the magnitude and distribution of heat. At the end of the slide (Fig. 5(d)), the peak heat value of the unit reaches approximately 1.9 × 105 J. This process continues for only 12 s. The produced heat can increase the temperature of 1 kg of clay soil (1406 J/(kg·℃)) to 152 ℃, which changes the properties of the soil (Zhang Fanyu et al., 2018). At this temperature, moisture can evaporate from the soil. Moreover, associated thermal pressure can also change the mechanical characteristics of the soil at a microscopic level (Wang Y. F. et al., 2017; Noda H. et al., 2011).

3.5 Heat Distribution Characteristics of Landslides at Different Heights

To explore the influence of landslide height on the heat value of the slide zone, the model is operated with inputs of constant values of dip angle and height. For the height values, 50, 100, 150, 200, 300 and 400 m are used. The number of elements in each modeling is approximately 200, 000 and six modeling scenarios are simulated. High heat zones are observed in the landslide model. The maximum heat values of the corresponding elements are 7 900, 1.9 × 105, 6.3 × 105, 4.8 × 106, 1.6 × 107 and 3.8 × 107 J. The total heat generated during the sliding process is 9 × 106, 4.8 × 108, 2.3 × 109, 6.4 × 109, 3.3 × 1010 and 7.6 × 1010 J. Fig. 6 indicates the relationship between heat value and land heights on the double logarithmic table. The total heat value of heat belt and the maximum heat value of elements are exponentially related to landslide height.

 Fig. 6 Relationship between the total and maximum values of heat in landslides of different heights

The empirical formula between heat value and landslide height is expressed as follows:

 ${Q_1} = 0.7719{H^{4.2916}}$ (10)
 ${Q_2} = 0.0008{H^{4.1427}}$ (11)

Where Q1 is the peak value of the heat element, Q2 is the total value of the heat belt, and H is the landslide height.

The numerical simulation shows that the location of the peaks of heat elements varies with landslide height. As denoted in Fig. 7, when landslide height is small, the peak of the heat element appears in the front of the heat belt. As the height increases, the heat peak moves to the band's rear. This is caused by the larger slip distance that is generated at the higher portions of the slope.

4 DISCUSSION

Numerical simulations indicate that the simulation accuracy is influenced by the element number. After testing, the results show that the landslide scale model has certain requirements for the number of model particles. When the average particle size is greater than 1/100 of the model length, the simulation results of the landslide will change significantly with the increase of the particles. By contrast, when the particle size is smaller than this value, the simulation results will not change significantly. Obviously the simulation results will be closer to reality as the number of particles increases.

As shown in Fig. 4, the heat generated during the entire sliding process is mainly the frictional heat. The rate of increase of frictional heat has a positive correlation with the sliding speed. It is shown that the long runout phase is important for the heat generating in the development process of a landslide. Furthermore, the simulation experiments show that the total heat value of heat belt and the maximum heat value of element are exponentially related to landslide height. It should also be noted that the highest heat value of elements realized in the heat zone is 1.9 × 105 J, the density of the cohesive soil is 1.7 × 103 kg/m3, the average particle size is 0.1 m, and the mass is about 0.88 kg. For an idealized 0.88 kg of cohesive soil with a specific heat of 1 406 J/(kg·℃), the soil temperature may increase to approximately 172℃. Such temperature is sufficient to vaporize the liquid water in the cohesive soil to reduce the friction of the slip zone.

However, the model still needs some improvement. Firstly, this model explores the process of heat generation in sliding belts, and points out that the water in the sliding belts may vaporize and reduce the sliding frictional force. This process is not reflected in the model. At present, a scheme is proposed to achieve this process by reducing the sliding friction of particles whose temperature reaches a certain limit, which requires later testing. The current maximum sliding speed of landslides can reach 25 m/s. After adding this change, this speed will be faster. Secondly, the model also needs to be extended to three-dimensional level and the constitutive model between particles should be changed to adapt to different landslides.

5 CONCLUSIONS

In this paper, a slope model is established using the MatDEM software. The sliding process and the heat generated during high-speed and long-runout landslides are simulated. The law of conservation of energy is satisfied by the discrete element model, and the heat distribution characteristics around the sliding zone are analyzed for landslide models at different heights.

(1) The heat is mainly composed of frictional heat that accumulates near the slip zone, and probably promotes the long-runout landslide.

(2) During landsliding, a banded high heat belt forms along the sliding area. This heat belt is related to the location of the maximum thickness of the sliding body. As the landslide height increases, heat increases exponentially.

(3) For a modeled slope height of 100 m, the highest heat value of elements realized in the heat zone is 1.9 × 105 J. For an idealized 0.88 kg of cohesive soil with a specific heat of 1 406 J/(kg·℃), the soil temperature may increase to approximately 172 ℃. This temperature is high enough to evaporate water and induce changes in the soil from thermal pressure.

(4) In the double logarithmic table, the peak value of the heat element and the total value of the heat belt have a strong linear relationship with landslide height, which influences the position of the core of the heat belt. In detail, the distance between core's position to the rear of the slip zone decrease with the rise of the landslide height.

REFERENCES
 Breien H., De Blasio F.V., Elverhøi A., Høeg KBreien H., De Blasio F.V., Elverhøi A., Høeg K. Erosion and morphology of a debris flow caused by a glacial lake outburst flood, Western Norway[J]. Landslides, 2008, 5(3): 271-280. Chen Lujun, Xing Aiguo, Chen Longzhu, Yin Yueping, Tong LiqiangChen Lujun, Xing Aiguo, Chen Longzhu, Yin Yueping, Tong Liqiang. Numerical analysis on the flying of highspeed and long runout landslide[J]. Hydrogeology & Engineering Geology, 2008, 35(5): 1-6 (in Chinese with English abstract). Dalguer L.A., Irikura K., Riera J.DDalguer L.A., Irikura K., Riera J.D. Simulation of tensile crack generation by three-dimensional dynamic shear rupture propagation during an earthquake[J]. Journal of Geophysical Research:Solid Earth, 2003, 108(B3): 2144. De Blasio F.VDe Blasio F.V. Production of frictional heat and hot vapour in a model of self-lubricating landslides[J]. Rock Mechanics and Rock Engineering, 2008, 41(1): 219-226. Ge Yunfeng, Tang Huiming, Wang Liangqing, Xiong Chengren, Zhang Shen, Wang DingjianGe Yunfeng, Tang Huiming, Wang Liangqing, Xiong Chengren, Zhang Shen, Wang Dingjian. Strain energy evolution of penetrative rock joints under shear loading[J]. Chinese Journal of Rock Mechanics and Engineering, 2016, 35(6): 1111-1121 (in Chinese with English abstract). Goren L., Aharonov E., Anders M.HGoren L., Aharonov E., Anders M.H. The long runout of the Heart Mountain landslide:Heating, pressurization, and carbonate decomposition[J]. Journal of Geophysical Research:Solid Earth, 2010, 115(B10): B10210. Guo ChangbaoGuo Changbao. Recognition and prevention of high-speed long-distance landslides[J]. Land Resources Science and Culture, 2018(1): 32-35 (in Chinese). Habib PHabib P. Production of gaseous pore pressure during rock slides[J]. Rock Mechanics, 1975, 7(4): 193-197. Hu Mingjian, Pan Huali, Zhu Changqi, Wang FawuHu Mingjian, Pan Huali, Zhu Changqi, Wang Fawu. High-speed ring shear tests to study the motion and acceleration processes of the Yingong Landslide[J]. Journal of Mountain Science, 2015, 12(6): 1534-1541. Huang Runqiu, Pei Xiangjun, Li TianbinHuang Runqiu, Pei Xiangjun, Li Tianbin. Basic characteristics and formation mechanism of the largest scale landslide at Daguangbao occurred during the Wenchuan earthquake[J]. Journal of Engineering Geology, 2008, 16(6): 730-741 (in Chinese with English abstract). Jiang M.J., Konrad J.M., Leroueil SJiang M.J., Konrad J.M., Leroueil S. An efficient technique for generating homogeneous specimens for DEM studies[J]. Computers and Geotechnics, 2003, 30(7): 579-597. DOI:10.1016/S0266-352X(03)00064-8 Liu C., Fan H. M., Zhu C G, S hiLiu C., Fan H. M., Zhu C G, S hi. Discrete element modeling and simulation of 3-Dimensional large-scale landslide-taking Xinmocun landslide as an example[J]. Journal of Engineering Geology, 2019, 27(6): 1362-1370 (in Chinese with English abstract). Liu Chun, Pollard D.D., Gu Kai, Shi BinLiu Chun, Pollard D.D., Gu Kai, Shi Bin. Mechanism of formation of wiggly compaction bands in porous sandstone:2. Numerical simulation using discrete element method[J]. Journal of Geophysical Research:Solid Earth, 2015, 120(12): 8153-8168. DOI:10.1002/2015JB012374 Liu Chun, Pollard D.D., Shi BinLiu Chun, Pollard D.D., Shi Bin. Analytical solutions and numerical tests of elastic and failure behaviors of close-packed lattice for brittle rocks and crystals[J]. Journal of Geophysical Research:Solid Earth, 2013, 118(1): 71-82. DOI:10.1029/2012JB009615 Liu Chun, Xu Qiang, Shi Bin, Deng Shang, Zhu HonghuLiu Chun, Xu Qiang, Shi Bin, Deng Shang, Zhu Honghu. Mechanical properties and energy conversion of 3D close-packed lattice model for brittle rocks[J]. Computers & Geosciences, 2017, 103: 12-20. Noda H., Kanagawa K., Hirose T, Inoue ANoda H., Kanagawa K., Hirose T, Inoue A. Frictional experiments of dolerite at intermediate slip rates with controlled temperature:rate weakening or temperature weakening?[J]. Journal of Geophysical Research:Solid Earth, 2011, 116(B7): B07306. Parise M., Jibson R.WParise M., Jibson R.W. A seismic landslide susceptibility rating of geologic units based on analysis of characteristics of landslides triggered by the 17 January, 1994 Northridge, California earthquake[J]. Engineering Geology, 2000, 58(3/4): 251-270. Rice J.RRice J.R. Heating and weakening of faults during earthquake slip[J]. Journal of Geophysical Research:Solid Earth, 2006, 111(B5): B05311. Scaringi Gianvito, Fan Xuanmei, Xu Qiang, et alScaringi Gianvito, Fan Xuanmei, Xu Qiang, et al. Some considerations on the use of numerical methods to simulate past landslides and possible new failures:the case of the recent Xinmo landslide (Sichuan, China)[J]. Landslides, 2018, 15(7): 1359-1375. Tang Huiming, Liu Xiao, Hu Xinli, Griffiths D.VTang Huiming, Liu Xiao, Hu Xinli, Griffiths D.V. Evaluation of landslide mechanisms characterized byhigh-speed mass ejection and long-run-out based on events following the Wenchuan earthquake[J]. Engineering Geology, 2015, 194: 12-24. DOI:10.1016/j.enggeo.2015.01.004 Travelletti J., Malet J.P., Delacourt CTravelletti J., Malet J.P., Delacourt C. Image-based correlation of Laser Scanning point cloud time series for landslide monitoring[J]. International Journal of Applied Earth Observation and Geoinformation, 2014, 32: 1-18. DOI:10.1016/j.jag.2014.03.022 Wang Y.F., Dong J.J., Cheng Q.GWang Y.F., Dong J.J., Cheng Q.G. Velocity-dependent frictional weakening of large rock avalanche basal facies:Implications for rock avalanche hypermobility?[J]. Journal of Geophysical Research:Solid Earth, 2017, 122(3): 1648-1676. Xu Qiang, Li Yanrong, Zhang Shuai, Dong XiujunXu Qiang, Li Yanrong, Zhang Shuai, Dong Xiujun. Classification of large-scale landslides induced by the 2008 Wenchuan earthquake, China[J]. Environmental Earth Sciences, 2016, 75(1): 22. Yin Yueping, Xing AiguoYin Yueping, Xing Aiguo. Aerodynamic modeling of the Yigong gigantic rock slide-debris avalanche, Tibet, China[J]. Bulletin of Engineering Geology and the Environment, 2012, 71(1): 149-160. DOI:10.1007/s10064-011-0348-9 Zhang Fanyu, Kong Ran, Peng JianbingZhang Fanyu, Kong Ran, Peng Jianbing. Effects of heating on compositional, structural, and physicochemical properties of loess under laboratory conditions[J]. Applied Clay Science, 2018, 152: 259-266. Zhang Lijiao, Ye Shulin, Liu QingchengZhang Lijiao, Ye Shulin, Liu Qingcheng. Research on the application of soil natural thermoluminescence survey in prospecting landslide boundary[J]. Uranium Geology, 2008, 24(5): 298-301 (in Chinese with English abstract).