A Dynamic Strain-Rate-Dependent Contact Model and Its Application in Hongshiyan Landslide

An earthquake-induced landslide, mainly affected by seismic movement, is a frequent large-scale geological hazard in hydraulic engineering. This paper proposed a rate-dependent strain-softened micromechanical contact model and implemented it in discrete element method code, namely, PFC. Using the PFC-FLAC coupling scheme, the Hongshiyan earthquake landslide is analyzed as a case study. The influence of the strain rate, damping, and topographic effect is discussed. The results indicate that the rate-dependent micromechanical model can give a reasonable seismic-induced failure process compared with the in situ situation and provide a numerical technique for earthquake-induced landslide analysis and rockfall hazard prediction.


Introduction
Earthquake-induced rock landslide is a common type of geological disaster that causes great loss to human life and property. It is of great significance for the disaster prevention and mitigation of a slope to analyze the process of an earthquakeinduced landslide so as to reveal the formation mechanism of a rockslide and thus provide basis for the implementation of slope stability measures.
Current research on landslides is aimed at defining the degree of difficulty of predicting a landslide based on macroscopic factors of slope stability, as well as the intensity of earthquake and rainfall, so as to gain advanced warning. However, this method is only an approximation and is often inaccurate for specific slopes [1,2]. The dynamic response and deformation damage of a fractured rock mass under seismic loading are mainly controlled by structural planes, and its stability analysis is inaccurate due to the difficulty in determining mechanical parameters and the huge number of structural planes. Most of the current research methods are geomechanical analyses supplemented by engineering experience [3,4].
In recent years, numerical simulations using block discrete element, particle discrete element, discontinuous deformation analysis (DDA), and other discontinuous methods for feedback analysis of the seismic landslide process and mechanism have attracted wide attention [5][6][7]. The discrete element method (DEM), which fully takes into account the discontinuous and nonlinear dynamics characteristics of rocks, is widely used to solve rock engineering problems such as block caving and blasting, and to simulate the dynamic response of a slope under seismic wave; however, it is difficult to determine its parameters and is not suitable for determining the critical sliding surface of rock mass [8,9]. Therefore, the combined finite-discrete element method (FDEM) is often used in actual simulations, which can simulate the progressive failure process of materials from continuous to discontinuous states [10][11][12][13].
Once an earthquake-induced landslide begins, the motion behavior of the landslide body is closely linked to many factors such as the material heterogeneity, fracture propagation, strain rate, and topographic slope, in addition to seismic forces and gravity [14][15][16]. There are massive collisions and friction during the movement and accumulation process of the rock mass, which is characterized by strain softening [17,18]. In largescale landslides, motion parameters are strongly influenced by topographic factors, with concave regions having potential landslide hazards and gently undulating regions equally likely to produce large and frequent landslides when struck by an intense earthquake [19,20].
During seismic movement and accumulation, the strain rate of the rock changes dramatically with the movement rate, which in turn leads to differences in energy dissipation and crack propagation modes, and ultimately, postearthquake slope failure and accumulation patterns [21,22]. In seismic landslide analysis, it is important to consider the effect of dynamic strain rate changes on the fragment characteristics of the rock mass.
In this study, a dynamic soft-bond model considering the strain rate effect is developed to simulate the rock mass based on the Hongshiyan earthquake-induced landslide in Yunnan Province, China. According to the field investigation data, a continuous-discontinuous numerical model is established, through which a reasonable combination of damping values for the Hongshiyan landslide is discussed, and the start-up time, damage form, movement, and accumulation process of the landslide under the action of seismic waves are simulated and analyzed.

Dynamic Soft-Bond Model considering Strain Rate Effect
Particle flow code (PFC), one of the discrete element methods, has been widely used to analyze large deformations and failure of rock and soil mass in recent years. PFC allows the displacement and rotation of discrete particles and is not limited by the amount of deformation. This can be very helpful in dealing with the discontinuous mechanical problems and simulating the cracking, separation, and movement of medium such as slope rock mass [23,24]. Considering the limitation of calculation efficiency, the continuous-discontinuous coupling numerical simulation method is often used in the calculation and analysis of large-scale geotechnical engineering. The potential failure area is simulated by the PFC model, and the area without a large deformation and failure, such as bedrock, is simulated by the finite difference model, the Fast Lagrangian Analysis of Continua (FLAC). Coupling PFC and FLAC models can meet the requirements of both deformation and calculation efficiency [25].
2.1. Soft-Bond Constitutive Model. The linear parallel bond model is the most common contact model for simulating the failure process of a rock slope by using the particle flow method [26]. This model can transmit forces and moments between particles, and better reflect the mechanical response of cementation inside the rock mineral. However, due to the limitation of calculation scale, the size of particles used to simulate the failure of a slope can be up to several meters to tens of meters. Within this size range, there are a large number of internal fractures making the rock mass softened during the process of failure, which cannot be reflected by the linear parallel bond model. The bonded behavior of the soft-bond model proposed by Ma and Huang [27] is similar to that of the linear parallel bond model. However, contrary to the linear parallel bond model, the bond is not removed upon failure. Instead, it may enter into a softening regime until the bond stress reaches a threshold value at which the bond is removed. The slope and tensile breakage strength during softening are determined by the softening factor ζ and softening tensile strength factor γ, respectively. Therefore, the soft-bond model can be used to simulate the failure behavior of a large-sized granular rock mass such as large-scale landslides.
The force-displacement law for the soft-bond model updates the contact force and moment: where F is a linear force, F d is a dashpot force, and M is a linear moment.
As shown in Figure 1, the linear force is resolved into a normal and shear force, and the linear moment is resolved into a twisting and bending moment: where F n > 0 is tension. The update mode considering the where Δδ n and Δδ s are, respectively, the relative normaldisplacement and shear-displacement increments; Δθ b and Δθ t are, respectively, the relative bend-rotation and twistrotation increments; k n ð_ εÞ and k s ð_ εÞ are, respectively, the normal and shear stiffness; A, I and J are, respectively, the area, the moment of inertia, and the polar moment of inertia of the cross section. The cross-sectional properties of the soft-bond contact are updated as follows: Tensile failure

Geofluids
where R is the radius of the cross section, λ is the radius multiplier.
There is a high correlation between Young's modulus (E), Poisson's ratio (ν), effective modulus (E * ð_ εÞ), and the normal-to-shear stiffness ratio (k * ð_ εÞ) of the soft-bond model. E increases with the increase of E * ð_ εÞ, and ν increases to its limit with the increase of k * ð_ εÞ. With the relationship between normal stiffness k n ð_ εÞ, shear stiffness k s ð_ εÞ, and E * ð_ εÞ, k * ð_ εÞis given as follows: where R ð1Þ and R ð2Þ are, respectively, the radii of the two par-ticles in contact. The update mode of the maximum normal stress σð_ εÞ and shear stress τð_ εÞ considering the strain rate at the bond periphery is as follows: where β is the moment-contribution factor (0 ≤ β ≤ 1). By default, the bonded behavior of the soft-bond model is essentially that of the linear parallel bond. When ζ ≠ 0 and γ ≠ 1:0, a softening behavior in tension can occur. As shown in Figure 2, if the maximum normal stress at the bond periphery exceeds the bond tensile strength (σ > σ c ), then the bond enters in the softening regime, and the maximum stress is given by   where ζ is the softening factor, l c and l are, respectively, the critical and current bond elongation, which are set as follows:

Geofluids
where δl and jδθ b j are, respectively, the measure of the bond elongation and accumulated bending since softening started. If the current tensile stress is greater than the maximum stress (σ ≥ σ * ), then it is projected back on the softening envelope: Different from the linear parallel bond model, the tensile stress is reduced with continued unloading past the peak tensile strength, and softening continues until the tensile stress during unloading is less than the product of the softening tensile strength factor and peak tensile strength (σ ≤ σ c γ), after which the bond is broken in tension and F and M are reset to 0 (Figure 2(a)).

Geofluids
If the bond has not failed in tension, then shear failure is assessed. If the shear strength limit is exceeded ( τ > τ c , where τ c = c − σ tan ϕ), then the bond is broken in shear, and the frictional strength properties are reset as follows: where the critical bending and twisting moments are as follows: where λ b and λ t are, respectively, the bending and twisting friction multiplier (Figure 2(b)).

Strain Rate Effect.
The dynamic characteristics of rock are directly related to the strain rate, which can significantly affect the ultimate strength, fracture mode, fragmentation degree, and constitutive model of rock [16,28]. Related studies show that there is a logarithmic relationship between rock energy and strain rate. As the strain rate increases, both the strength of the rock and the energy absorbed per volume of rock increase, whose differences lead to the difference in mesoscopic crack propagation modes, and result in the difference in roughness of the fracture surface, and finally the difference in failure patterns [29,30]. When the failure process of a rock slope is simulated by the PFC model, the influence of the strain rate effect on the mechanical parameters of a rock mass should also be taken into account. The microscopic parameters of the discrete element model are closely related to the macroscopic mechanical characteristics, and there is a quantitative correlation between them [31]. It is necessary to analyze the relationship between normal stiffness (k n ), shear stiffness (k s ), tensile strength (σ c ), cohesion (c), and strain rate (_ ε) during simulation. Strain rate in PFC is defined as the ratio between the relative velocity of two contact particles and the centroid  30 40 Seismic wave (NS) Seismic wave (EW) where _ ε is the strain rate, v ! i and v ! j are the velocities of two particles in contact with one another, L is the centroid distance between two contact particles. Figure 3(a) shows the PFC3D numerical model, where the sample is a cylinder with a radius of 50 m and a height of 280 m, and the walls are applied around and up and down to fix it. There are 18773 particles with a radius of 1.5-3.0 m and a porosity of 0.35, whose microscopic parameters are set according to those of the slope model (Table 1). Firstly,  9 Geofluids the relationship between macroscopic parameters and strain rate is obtained through uniaxial and triaxial compression tests. Then, the relationship between micro-and macroparameters is obtained through a correlation study. Finally, the formula of microscopic parameters changing with the strain rate is obtained by combining the two [32]: where E * is the effective modulus (GPa); σ c and c are, respectively, the tensile strength (MPa) and cohesion (MPa); _ ε and _ ε 0 are, respectively, the dynamic and static strain rate (_ ε = 10 −5 when _ ε < 10 −5 ; _ ε 0 = 10 −5 ). Figures 3(b) and 3(c) show the stress-strain curves under different strain rates and confining pressures. The specimen exhibits adequate ductility after failure, which is in good agreement with the laboratory test results.  Figure 4(a), the surface layer of the left bank became loose and slid towards the riverbed, while the rock mass of the right bank toppled and collapsed at a high speed under the action of seismic waves, resulting in a large-scale landslide which rapidly accumulated towards the riverbed and blocked the river.

Analysis of the Hongshiyan Seismic Landslide
After the landslide induced by the Ludian earthquake, the topography of the right bank has changed significantly. Figures 4(b) and 4(c) shows the postearthquake topography and three main landslide runout areas of the right bank. The upper part of the slope is a huge remnant scarp with a height of 350 m and a width of approximately 800 m. The lower part is an original cliff with a height of 120 m, and the two parts are separated by a tilting platform dipping towards the river. Because of the influence of topography, the landslide is divided into three runout areas. In the damming body area, there are a large number of huge rocks, with the particle size of 5-10 m.
Considering that only the toe part of the slope on the left bank, with ancient landslide deposits, collapsed and slid under the action of seismic waves [33], and the river had no influence on the landslide process except for the final accumulation area, this paper only selects the rock mass on the right bank of Niulan River for the analysis of earthquake-induced rockfall and rockslide.
3.1.2. Model Building and Boundary Conditions. As shown in Figure 5(b), the three-dimensional model of the Hongshiyan slope is established by the continuous-discontinuous coupling method of FLAC3D and PFC3D. The source area consists of 96970 discrete element particles, with a radius of 1.5-3.0 m based on the actual situation. As described in Section  Table 1 shows the parameters of PFC particles and FLAC units. Figure 5(c) shows the velocity curves of the Hongshiyan seismic waves in three directions of EW, NS, and UD, which are provided by the China Earthquake Administration. The dynamic response of the Hongshiyan rock slope is analyzed by applying the X, Y, and Z three-different-direction seismic dynamic load and viscous boundary.

Evolution Process of Microscopic Contact Parameters of
Rock Mass. During the process of seismic loading, the deformation and strength of the rock contact, whose microscopic parameters are updated according to equation (13), change dynamically with the strain rate of each contact traversed every time step. For illustration, contacts with ID = 10000, 50000, and 200000 are selected randomly, with the strain rate and tensile strength curves of each shown in Figure 6. Figure 6 shows that the strain rates of the three contacts vary mainly between 10 -5 and 10 -2 , all with different patterns, indicating that the method can implement different temporal and spatial variations of the rock mass under seismic action, which can reflect the dynamic properties of the rock mass.
3.3. Earthquake-Induced Landslide Process. As shown in Figure 7(a), the Hongshiyan rock slope is characterized by    12 Geofluids the antidip stratified structure [34]. Under the action of the earthquake, the rock mass first presents an obvious progressive and continuous failure phenomenon, followed by body collapse along the main structural plane. The rock landslide undergoes a process from the key block failure to overall failure. As the blocks on steep rocks have different movable conditions and tightness degrees, it is inevitable that some blocks, which are subject to gravity, seismic forces, friction, and other factors, are first displaced outward and become the key blocks for block deformation and failure [35], as shown in Figure 7(b). Figure 4(b), the landslide is mainly divided into three runout areas, Area I, Area II, and Area III.  Figure 8(a). Figure 8(b) shows the vertical stress-time curves of measuring points A, B, and C. The vertical stress at measuring point A has several peaks, mainly 10.5 MPa at 5.7 s and 9.2 MPa at 8.9 s, reaching more than twice the initial stress, which is sufficient to cause the Area I landslide; the vertical stress at measuring point B has a huge peak between 5.4 s and 9.0 s, during which Area II absorbs a large amount of seismic wave energy and the rock collapses; the vertical stress-time curve at measuring point C is relatively gentle, with no obvious peaks, which shows that the landslide in Area III primarily originates from the failure of Area I and Area II when the stress change and the stable state is broken. Figure 9 shows the vertical velocity-time curves of the measuring points in different areas within 40 s of the start of the earthquake. Due to the high slope, the landslide has a large initial potential energy. Under the action of the earthquake, Area I where monitoring particle No. 1 is located is the first to be destroyed and slide. The potential energy is converted into kinetic energy through sliding, rolling, or free falling, so that the rock mass has a high velocity during its movement. The vertical velocity in Area I reaches a peak of 28.4 m/s at 14.1 s, followed by a collision with the slope and a significant dissipation of kinetic energy, and then enters a low-speed motion, with four major collisions during the entire landslide (Figure 9(a)). Area II, where monitoring particle No. 2 is located, slides after Area I slides for about 2 s. After a brief phase of sliding along the tilting platform to the edge of the cliff, it enters the falling phase, where the particles are subject to gravity and damping force only, and the velocity rises sharply, reaching a peak vertical velocity of 62.6 m/s when it falls to the bottom of the slope in 19.9 s. After a collision with the bottom, the kinetic energy is almost completely lost and the blocks rapidly stabilize (Figure 9(b)). Area III, where monitoring particle No. 3 is located, slides after the slide in Area I and Area II. During its movement, there are massive collisions and friction with the slope, and the velocity is fluctuating without any obvious peak, of which the maximum was only 15.0 m/s, and it does not fall to the bottom until about 91 s (Figure 9(c)).

Landslide Analysis in Different Areas. As shown in
Several typical landslide moments are shown in Figure 10. Figure 10(a) is the moment when Area I starts to slide, Figure 10(b) is the moment when Area I first collides with the slope, Figure 10(c) is the moment when Area II collides with the bottom, and Figure 10(d) is the moment when Area III slides along the slope for a certain distance. From  Figures 10(a)-10(c), the rock mass (red dash line in Figure 10(b)) where monitoring particle No.1 is located first slides due to seismic force and gravity, followed by a break in the overall stability of the slope, which leads to the collapse of the entire loose blocky structure and the eventual largescale landslide. The particle velocity in Area I reaches its peak and then collides with the slope uplift, after which it enters a continuous collision deceleration phase. From Figures 10(c)-10(d), the particles in Area II have reached the bottom of the slope when their velocity reaches its peak, after which they collide and lose all their energy, with little displacement. When Area I and Area II are largely stabilized, the particles in Area III have just slid to a position approximately 100 m below the source area.   14 Geofluids

Topographic Amplification Effect.
In landslide hazards, highly susceptible slopes represent the high-potential failure areas triggered by earthquakes. It was found that topographic irregularities can greatly influence the amplitude and frequency of ground motion, which is systematically amplified on convex terrain (e.g., hills and ridges) and attenuated on concave terrain (e.g., canyons and foothills), with complex amplification and attenuation patterns on mountain slopes leading to significant differential motions [36,37]. The Hongshiyan landslide is a super-large-scale  landslide, and the topography changes dramatically during the landslide, so the topographic effects need to be analyzed in relation to the actual landslide process. Figure 11 shows the displacement change of the sliding bed at different stages. During the seismic wave propagation stage, a significant seismic amplification effect occurs at the cliff. As the landslide proceeds, the cliff mass becomes more prominent and its curvature increases, followed with further enhancement of the seismic amplification effect and rapid increase of the displacement. When the majority of the sliding body falls to the bottom of the slope, the curvature of the cliff mass decreases instead, while the seismic amplification effect at the slope top gradually increases, resulting in obvious displacement. Figure 12(a) shows the displacement cloud map in the Y -direction 2 min after the start of the earthquake (the coordinate orientation is shown in Figure 4(b)). It can be seen that under the action of the seismic waves, the rock mass on the left side of the cliff (yellow dash line) uplifts by more than 2 m, and the effective free face narrows to form a spoon shape, which blocks the movement of Area I slides there and causes its velocity to drop from the peak [38]; the upper edge of the cliff (blue dash line) uplifts by approximately 2-4 m, blocking the movement of Area III. The postearthquake uplifts on the Hongshiyan slope due to the topographic amplification effect have the locking effect on the overall landslide. After the landslide, the sliding body composed of block stones accumulates on the slope surface and valley, partially occluded in the source area, which is proved by the postearthquake field investigation photographs (Figure 12(b)).   Mechanical energy is divided into two categories, namely, body energy and contact energy, and energy is accumulated in each cycle for balls, clumps, walls, and contacts. Mechanical body energy is the energy due to body force and dissi-pated by local damping and kinetic energy during the landslide; mechanical contact energy is defined by the contact models and dissipated by frictional sliding and viscous damping.

Discussions
4.1.1. Local Damping. Local damping acts on each ball or clump, which is applied with a damping force proportional to the magnitude of the unbalanced force, and the damping motion equation is [40] as follows:

Geofluids
where € x ðiÞ , _ ω ði−3Þ , F ðiÞ , and F d ðiÞ are, respectively, the mean acceleration, rotational acceleration, unbalanced force, and damping force of the particle; m and I are, respectively, the mass and rotational inertia of the particle; i is the degree of freedom of the particle, with i = 1, 2, 3 representing translation and i = 4, 5, 6 representing rotation. F d ðiÞ is given by where α is the local damping constant, sign ðxÞ is the sign function, and v ðiÞ is the generalized velocity as follows: Three different sets of local damping constants (α = 0:15, 0.3, and 0.5) are applied to the slope model under seismic conditions, and a continuous-discontinuous coupled numerical analysis is performed. The mean velocity-time and mean displacement-time curves for all particles are shown in Figure 13. Considering that 100 s after the earthquake starts, the velocity of particles tends to 0 and the displacement enters a gentle stage, the displacement cloud maps of the landslide for 100 s under different local damping constants are selected for comparative analysis ( Figure 14). Combining Figures 13 and 14, it can be seen that as the local damping increases, the velocity and displacement of particles show a downward trend and the amount of particles occluded in the source area gradually increases. By comparison, when the local damping constant α = 0:3, the location and scale of the model occluded area (Figure 15(b)) tallies well with that of the actual occluded area (Figure 15(a)).
where C ðiÞ is the viscous damping constant and _ δ ðiÞ is the relative normal (shear) velocity at contact; n represents normal, s represents shear. The relationship between the damping constant and critical damping constant is as follows: where m is the average mass of particles in contact, k ðiÞ and β ðiÞ are, respectively, the normal (shear) stiffness and critical viscous damping ratio at contact.

Geofluids
Nonzero local damping can be used for quasistatic deformation simulation. However, during a seismic landslide, the strain rate is high and there is a large number of collisions between particles; it is inappropriate to consider only local damping, and further comparative analysis should be carried out with different normal (shear) critical viscous damping ratios based on the soft-bond contact model.
Three different sets of critical viscous damping ratios (β n , β s = 0:3, 0.5, and 0.7) are applied to the slope model under seismic conditions, and a continuous-discontinuous coupled numerical analysis is performed. The mean velocity-time and mean displacement-time curves for all particles are shown in Figure 16, and the displacement cloud maps of the landslide for 100 s under different critical viscous damping ratios are shown in Figure 17. From Figure 16, the difference in the velocity and displacement of particles under different critical viscous damping ratios is not significant, and the peak velocity increases slightly with the increase of critical viscous damping ratio. Considering the huge kinetic energy of the landslide, which makes the sliding body eventually accumulate to form a damming body, the critical viscous damping ratio β n , β s = 0:5~0:7 is more appropriate.

Influence of Strain Rate Effect.
Many studies show that seismic cyclic loading usually induces the expansion of preexisting fractures in a rock mass, followed by the formation of isolated blocks, and finally the landslide under the action of gravity. Once the landslide starts, the movement of the sliding body is mainly affected by a variety of factors such as damping, block degree, rock disintegration, and topographic slope [41,42]. If constant material strength parameters are used and strain rate effect is not taken into account, the calculated results of disintegration and fragmentation of the landslide body, final runout distance, and deposition topography are very much inconsistent with that of field observation [18]. Figure 18 gives the vertical velocity-time and fracture number-time curves of the landslide with and without considering the strain rate effect, which shows that when the strain rate effect is considered, the peak velocity of particles decreases slightly and the number of fractures decreases by 5.5% compared to that without considering the strain rate effect. The strain rate effect has a certain influence on the propagation modes of rock cracks, resulting in differences in the degree of rock fragmentation and accumulation patterns. The simulation results are generally consistent with that of the field investigation, indicating that the simulation results using the dynamic soft-bond model can well reflect the movement and runout behavior of the Hongshiyan landslide. The final pattern of the landslide is shown in Figure 19.
Zhang et al. [43] used both the sphere-simplified model and the 3D reconstruction model to analyze the effect of block impact on the bedrock, and the results show that the 3D reconstruction model is beneficial in considering the combination effect of rockfall velocity and rotation during the collision. In this paper, the spherical model is used for the landslide body for simplification, and the influence of the microscopic morphology of rockfall on the motion state during impact needs further study.

Conclusions
Earthquake-induced slope failure is a dynamic and progressive cumulative failure process, in which the deformation and strength of the rock mass evolve with the strain rate. The dynamic soft-bond contact model which can reflect the strain softening and strain rate strengthening effect is established by using the continuous-discontinuous coupling method, and numerical simulations are carried out to analyze the Hongshiyan earthquake-induced landslide. The primary conclusions are summarized as follows: (1) A dynamic soft-bond contact model considering the influence of strain rate was developed based on the soft-bond model, which can meet both the strain rate effect and strain softening caused by the fracture network. The fragmentation, velocity, and runout distance of the landslide body are reduced by less than 10% when the strain rate effect is considered compared to when it is not considered.
(2) During the earthquake landslide, the damping parameters have a significant influence on the movement and accumulation process of the landslide. The calculated results show that the landslide runout area with low local damping (0.3) and high viscous damping ratio (0.7) is very much consistent with that of the field investigation (3) The Hongshiyan landslide begins from the occurrence and development of fractures under seismic action, then to the formation of a structural plane separating the rock mass, and ends with a rockfall and a rockslide under seismic forces and gravity. The landslide body experiences continuous collision with the uplifted rock mass, during which the energy is continuously lost, and it finally accumulates at a certain range of the valley, with some occluded in the source area

Data Availability
Data is available on request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.