Numerical simulation of rainfall-induced rock mass collapse and debris flow

Geologic disaster of rainfall-induced rock collapse and debris flow is an essential part of research in the field of geotechnical engineering. Based on the discrete element method, the article derived the constitutive model of particles bond-damage fracture-move. A two-dimensional slope model of rock mass collapse was established to simulate slide of rock mass collapse and dynamic evolution process of debris flow. The change of mechanics parameters of six monitoring points were tracked and analyzed. The results showed that, nonlinear motion of microscopic particles was obvious during the process of bond-damage fracture-move. Based on the geological conditions of rock mass collapse and debris flow on Greenland in Denmark, a three-dimensional slope model was established. The results showed displacements and velocities of measuring points both displayed linear relationship with slope when slope increased from 1.0 to 2.0; while non-linear capacity was strong under the slope with a degree above 2.0. The research will be the foundation for nonlinear movement of debris flow and this kind of disaster induced by different factors.


INTRODUCTION
Slope sliding and debris flow is a typical geological disaster whose forming mechanism is influenced by many factors.As movements of geological disasters induced by rainfall are complex, some technical problems have not been broken through on dynamic problems.Researches on dynamics of debris flow are currently hot and difficult problems for academia.
In order to explore the physical mechanism of debris flow, scholars at home and abroad have done much research on theoretical exploration and practical application (Iverson, 1997;Kane and Ponten, 2012;Neal et al., 2013).As a kind of special fluid, debris flow possesses of tremendously complex features.Scholars mainly establish dynamic model of debris flow basing on the theory of SH and their researches on processes of debris flow always focus on numerical simulation (Tang, 1994;Yu, 1996;Ma et al., 2008).The United States Geological Survey (USGS) set up the institutions on debris flow and their statistical model has been applied.Many scholars in the domestic such as Du (1987), Wu (1990), Kang (2004) and so on, have made pioneering research on formation mechanism of debris flow in China *Corresponding author.E-mail: jikunzhao_2006@163.com,Tel: +8613675132073.Fax: +862558606573.Author(s) agree that this article remain permanently open access under the terms of the Creative Commons Attribution License 4.0 International License and accumulated valuable observation data.
Gravity is the main driving force of rainfall inducing rock collapse and debris flow.The dynamic mechanism of liquid is very complicated that contains many forms of motion, such as collision, friction, sliding and rotation between particles.An important approach to studying its movement is numerical calculation.Discrete element method has been utilized in many researches for slope stability (Federicoa and Michaele, 2007;Lei and Wang, 2006).It also has a big superiority in solving the formation and evolution process of debris flow and has solved the difficult problem on movement evolution of discontinuous medium.Sitharam and Vinod (2010) attained the degree of impact of other factors on dynamic shear modulus and damping ratio of particle materials through Discrete Element Method (DEM).Mohsin et al. (2013) established a one-dimensional model by using Hydrologic Engineering Center River Analysis System (HEC-RAS) and studied debris flow in view of dam break.Wang and Fulvio (2011) analyzed the influence of factors on rock collapse by DEM.The parametric inversion analysis on the rolling motion was completed by Jiang et al. (2008) through DEM.Li et al. (2009) analyzed stability of rock slope by DEM.Tang and Xu (2007) studied problems on geotechnical multi-scale by coupling the discrete element and shell finite element.Zhou et al. (2010) realized the two dimensional discrete-continuous coupling analysis for problems on geotechnical engineering deformation by Particle Flow Code (PFC) and Fast Lagrangian Analysis of Continua (FLAC) Socket I/O interface; numerical examples showed the effectiveness of the method on coupling analysis.Fan et al. (2010) simulated dynamic process of entrainment along the way, they found that the entrainment action could increase the motion volume of debris flow and the debris flow would be more destructive and harmful.Based on the Bingham theological theory, Zhang et al. (2012) built a model of debris flow, simulated the gully bed evolution and carried out the bed erosion-deposition processes.Nicholas et al. (2014) analyzed seven debris flows initiated in proglacial gullies.The high slopes and elevations of debris flow-producing proglacial areas reflect positive slope-elevation trends for the Mount Rainier volcano.
Rock collapse is a key cause secondary disaster.Such as sediment deposition, debris flow erosion, impact and stacking problems.Erosion and the impact of disasters is one of the main causes of the disaster along the way.Sediment deposition is one of the main causes of the local lake.Because of the above reasons, the paper studied rock collapse and sediment deposition.Research on effects of microscopic mechanic parameters on the macroscopic characteristics in process of rock mass collapse and debris flow is less.Based on the discrete element method, the article derived the constitutive model of particles connect-damage fracture-move.A slope model was established to simulate slide of rock mass collapse and dynamic process of debris flow.The article also analyzed mechanism characteristics of rock collapse and evolution of debris flow under different slopes.

THE CONSTITUTIVE MODEL OF PARTICLES BOND-DAMAGE FRACTURE-MOVE
There are two kinds of connected models on discrete element.One is contact connection model which is equal to the rigid ring joint at contact points; there is a pair of springs with a constant normal stiffness and the tangential stiffness which can only transmit force.It cannot transmit bending moment; otherwise it will break when the force is greater than the connection strength.The other is parallel connection model which is equal to the flexible rubber connection on the contact surface and has rigidity connection itself.It will fracture when the force is greater than the connection strength.However, it cannot describe the particle damage stage.The paper deduced the constitutive model of damage of the particle.Microscopic parameter of contact connection model consists of normal connection strength and tangential connection strength.Zhao et al. (2013) have derived constitutive model of particles bond-damage fracture before.The paper would derive the equation of movement after the particles being damaged.
Assume that A and B are two adjacent particles of rock mass and the connection damage mechanical behavior is described by normal connection strength; tangential connection strength; normal connection stiffness; tangential connection stiffness and connection radius.The unit normal vector i n of particle -particle interface is represented as: x are the center vector of particles A and B, respectively; center distance d between particles is: ) R are the radius of particles A and B, respectively.The contact point of two particles is determined by the following expression: ) Total contact force is decomposed into tangential component and normal component along the contact surface: In the formula, F and 3 M are initialized to zero when the connection is formed later, elastic force and torque increment in the contact caused by the displacement increment and rotational increment will be stacked in the current value.Within one time step t ∆ , the elastic force increment is: In the formula, ) are instantaneous speed and instantaneous rotational speed of particles j Φ respectively.The j Φ can be determined by the following formula: Within one time step t Δ , the elastic torque increment is: In the formula, A is contacting connection area; J is polar moment of inertia of contact surface; I is moment of inertia of contact surface relative to the contact point along the s i θ ∆ -axis direction.Parameters above can be determined by the following formula: Through the force and torque between particles A and B, calculating the maximum tensile stress and maximum shear stress of connection: In the above equation, when the inequality is equal, namely, the maximum tensile stress max When the maximum tensile stress and maximum shear stress is smaller than normal strength and shear strength, the connecting force and torque among particles will be applied to the connecting particles A and B respectively.The expression is: In the formula, When shear stress is stronger than yield stress, the particle begins to flow and stress would reduce.Its translational motion and rolling movement mainly depends on the force and bending moment.The velocity and direction of the particle are controlled by the integral calculation from Newton's law.The resultant force and translational motion satisfy the following equation: where, m is quality, i x  is accelerated velocity, and i g is acceleration vector of volume force.i F is resultant force on the particle that can be expressed as the following: where, Resultant moment and rotation movement are satisfied in the following relation: where, I is inertia moment; ω  is vector of angular acceleration.The acceleration of the particle at time t can be obtained by central difference method: Translational velocity and angular velocity at time ( ) 2 / Δ + t t can be achieved by substituting Equation (4) into Equations ( 1) and (3), respectively: The center position of the particle at time ( ) 2 / Δ + t t can be obtained by Equation ( 21): The parameters such as translational velocity ( ) and rotational velocity ( ) and moment ( ) on particles at next moment can be achieved by loop iteration process of the method.
Based on constitutive model of particles, the author wrote a c++ program to establish user-defined constitutive model.Dynamic link library (DLL) file was output and called in the main particle discrete element program.

2D simulations
In order to study the geological disasters of the rock collapse and debris flow induced by rainfall, a 2D rock   5.The 2D numerical model is 110 m long and 30 m high.According to the initial void ratio with 0.1, particle generates randomly in the slope area and its radius changes within 4~6 cm.The contact mechanics parameters of microscopic particles have greater influence on the macroscopic mechanical properties.The mechanics parameters shown in Table 1 have been analyzed by author (Zhao et al., 2013), wherein the relative change rate of each index after rain refers to literature (Li et al., 2009).In the initial state, the 2D model was stable.After fully rainfall, the parameters of rock were weakened.The initial damage data of stiffness was 0.1407.This value was substituted into Equation ( 14), the rainfall stiffness value were calculated.The change rate of intensity was 0.2113.Then the data of t ξ was substituted into Equation ( 12), the rainfall intensity could be calculated.The overall strength of rock mass was induced by initial damage weakening, which resulted in the late rock collapse and debris flow.
The weakening of parameters in the damage process of rock mass follows the constitutive equation.Rock collapse model was shown in Figure 1a.Six different monitoring points on the slope surface and internal points were set to monitor and track motion of particles, the stress and deformation, respectively.There were three sections to monitor the total discharge of debris flow.Frist section is at the toe.Every five meters were set two sections in its base.
The top of the model slope is seen as the origin of coordinates, the vertical downward is negative and the horizontal position to the right is positive.Due to rainfall infiltration, index of shear strength of rock mass significantly became lower.The instability of rock mass would be damaged under the action of gravity, and change of vertical position was negative.
Figure 2 presents curves of vertical location versus time with monitoring points.There are many differences among them with rock collapse-damage fracture-debris flow.The processes of rock collapse and debris flow were shown in Figure 3. Taking the movement of Point 1 as an example, before turning inflection Point 1, the bonding force between the particle and adjacent particles were strong.With the process of collapse, the particle collided with other particles and the bedrock, and then damage and fractures where generated at inflection Point 1. Then the evolutionary phase came.The results showed that not all particles of rock mass were damaged and fractured; those unbroken particles mainly depended on their weight and moved along with others.Often this part of rock mass is large, and would have great impact on the downstream.Figure 3(a)-(e) present three stages of slope failure and evolution: (1) slope with overall stability; (2) initiation and perforation of sliding surface of slope; (3) evolution of debris flow.When time is before 0.8 h and the calculation step is 2.7 × 10 4 , slope is at the stable stage.Loose blocks which attached on the slope, dropped to the slope toe under the action of weight, while some particles dropped a little farther away from the slope toe (Figure 3 (a)).When time reached 1.5 h, the slope was at the early stage of sliding and connection between particles turned from damaging to fracturing gradually.Micro sliding surface first appeared at Point 3, with large angular velocity and speed horizontally along the slope surface which weakened integral strength of the slope.When time reached 3.5 h, debris flow moved with a good liquidity.
Due to the particles at monitoring Point 3 was located outside of the rock cliffs and the lower was impending, they fell early under the influence of weight (Figure 4).When time is 1.5 h, the monitoring point presents the large angular velocity and horizontal velocity.The characteristics of nonlinear motion of slope slide and debris flow is obvious.The main reason is that at early collapse of rock mass, the particles which were connected by the coherence force first collapsed and then collided with adjacent rock in the whereabouts particles, and then the path changed.According to different rotation directions of particles in the process of movement, the clockwise rotation is positive, counterclockwise rotation is negative.In the early stage of debris flow, nonlinear motion was also obvious.However, the motion tended to be stable in later evolution of debris flow.This is because the broken grain went into the hole among large particles, the overall porosity reduced, and then the overall viscosity increased.
The comparison was shown in Figure 5 between the simulation and the experiment.From Figure 5(a), the three section total discharge simulation results agree with the experimental results.Before 2 h, the three section numerical simulation results agreed with experimental results.Then total discharge and impact pressure results were slightly larger than the experimental values in Section 1. Debris flow total discharge results slightly smaller than the experimental values in Section 2. Although, there were some fluctuations, but the overall were regularity.Figure 5(b) presents that in Section 1, 1.9 h occurred at the peak impact pressure 2915KN, then the intensity values produce small fluctuations.It was suggested that the time rate of change of flow rate reaches a maximum debris flow.In the other two sections, the peak impact pressure was clearly lagging behind, but overall were slight fluctuations up and down in the experimental.

3D simulation studies
In order to study the geological disasters of the rock mass collapse and debris flow induced by rainfall, and compared with debris flow of Greenland in Denmark (Costard et al., 2002), four kinds of 3D models were established.They are 60 m wide and 150 m high.The slopes of bedrock are 1.0, 1.5, 2.0 and 2.5, respectively (Figure 6).According to the initial void ratio with 0.1, particle generated randomly in the slope area and its radius changes within 2~8 cm.The mechanics parameters were also shown in Table 1.In the initial state, the 3D model was stable.The particles within the boundary two XY displacements is constrained in the Z direction, and is allowed to occur displacement  movement in the XY plane.After fully rainfall, the parameters of rock were weakened.The initial damage data of stiffness was 0.1407.A weakening of the parameters at this time was 0.2113.The overall strength of rock mass was induced by initial damage weakening, which resulted in the late rock collapse and debris flow.Two monitoring points in the slope were set to monitor and track displacement and velocity.The processes of rock mass collapse and debris flow were shown in Figure 7.There are also three stages: (1) slope with overall sliding; (2) stage of breaking and rolling; (3) evolution of debris flow.By comparison, the deposition form of numerical simulation agrees with Greenland's form.Figure 7 presents that integrity of the slope is strong at the beginning of the slope slide and the whole rolling of the slope body occurred on the surface of the bedrock.Figure 8 shows the curves of angular velocity of monitoring points versus time.The angular velocity in clockwise is positive and in counterclockwise is negative.In early slide, the angular velocity is about 0.2 rad/s.Figure 9 presents that within less than 4 h, number of crushing particle has nothing to do with slope because the aggregate particles are damaged by particle collision force in the prophase of evolution.However, as time went on, energy for decline was released, so the number of crushing particles and slope increased the linear after 4.2 h.

Number of crushing particles
The rock rolled along colliding downstream, the slope body itself would be broken and generated about 81 particles, which is shown on the curve at 4 h in Figure 9.The whole strength of the slope body decreased as well as the cohesive strength.The curve of work of gravitational force in 4~8 h (Figure 10) shows that partial block separated and increment of gravity power was obvious.Later in the process of debris flow, it continued to flow along the bedrock, glutinousness and viscosity increased; the angular velocity of rock mass is close to zero (Figure 8).The curves of work of friction that produced by the viscous force are flat (Figure 10) in this phase.Then, debris flow reached the bottom of the bedrock with the advancement of upstream particles layered erosion accumulation was generated by front-end particles (Figure 11).
Figure 11 presents the contrast between contours of accumulation body by the numerical simulation and Greenland when the slope is 2.5.The numerical 0 and 10 m contours in front of accumulation area were slightly smaller than that of Greenland.The main reason is that in simulation, the established rock particles are sphere particles, the porosity of accumulated particles are large; while the density of actual rock mass is great after medium-sized particles filled the voids among large particles in the process of debris flow resulting in flat deposits.Taken together the contours of accumulation body by simulation fit in geological disasters of Greenland, so the rationality of the model has been verified.
Figure 12 presents the increase of slope, displacement and the speed of monitoring points along the slope present nonlinear trend for increasing.When the slope increases from 1.0 to 2.0, the curves of displacement increase by linear ratio.When the slope increases from 2.0 to 2.5, displacement and velocity increase significantly.Under the same slope, nonlinear motion of displacement and the speed of monitoring points are obvious.The reason is that the collision among particles and crushing changed their trajectory and orientation which would affect the displacement and velocity components.The above explain that the nonlinear characteristic of microscopic particles is closely combined with macroscopic changing process of rock mass.

DISCUSSION
After a rainfall, there is rock collapse and the occurrence and debris flow evolution of geological disasters.Research on this problem has great superiority by discrete element method.This method can be used to study micro body particle movement.Real time tracking effect is good.However, there are still problems and needs improvement.In the debris flow motion process, further research needs associated with slurry.In Figure 5, the two-dimensional simulation results and experimental results of node factors exist differences in slurry.
Numerical simulation of the discrete element method can achieve large-size 3D slope model at the same time need to generate a lot of particles.However, the calculation takes more computing time is currently drawbacks.How to achieve ultra-performance parallel computing is one of the ways to expand its field of application.

Conclusions
Based on the discrete element method, the article derived the constitutive model of particles connectdamage fracture-move and simulated rock mass collapse and dynamic process of debris flow.The article analyzed (1) The process of rock mass collapse and debris flow is complex kinetic processes, and presents progressive dynamic evolution process from continuum to discontinuous.The experiments of the model have practical significance on engineering to describe the nonlinear motion kinetic characteristics of the whole process of slope, (2) There are three stages of rock mass collapse and debris flow with different kinetic characteristics.In the evolution of debris flow, slope on cliffs is easy to slide, with great energy and obvious nonlinear movement, (3) Rainfall is one of the main factors that induce geological disasters.The numerical results agree well with the experimental results.Tracking the evolution of the dynamic parameters of rock particles can compensate experiment deficiencies.
increment and tangential relative displacement increment respectively between two contacting particles within one time step

σψU
is equal to normal strength c σ .The maximum shear stress max τ is equal to shear strength c τ .Under the rainfall condition at time t , the value of is the original value of c σ and c τ .tξis the relative change rate after rainfall.This data is weakened due to rainfall parameters.When the maximum tensile stress and maximum shear stress exceed the normal strength and shear strength, will happen.Normal damage and shear damage variable are defined as: are normal displacement and shear displacement of the particles respectively at t moment.The damage variable is 0 when the relative displacement reaches threshold displacement, and its value will become 1 when the relative displacement reaches n U .The normal stiffness n k and tangential stiffness s k at damage stage are calculated by the following formula: the torque of particles A and B.

F
and s is scalar tangential component of force; i n and i t are unit normal component and unit tangential component, respectively.

Figure 1 .
Figure 1.The model of rock collapse (a) Two dimensional model (b) Plan of the area (1:25000).

Figure 2 .
Figure 2. Curves of vertical location versus time with different monitoring points.

Figure 3 .Figure 4 .Figure 5 .Figure 6 .
Figure 3. 2D simulation of rock collapse and debris flow.(a) Early rock collapse (b) Stage of rock collapse (c) Early evolution of debris flow (d) Middle evolution of debris flow (e) Stable accumulation of debris flow.

Figure 7 .
Figure 7. 3D simulation of rock collapse and debris flows (a) Stage of slope slide (b) Stage of debris flow (c) Accumulation of numerical simulation (d) Accumulation of Greenland.

Figure 8 .
Figure 8. Curves of angular velocity of monitoring points versus time.

Figure 9 .
Figure 9. Curves of number of crushing particles versus time under different slope.

Figure 10 .
Figure 10.Curves of work of gravitation and viscous friction versus time, respectively.

Figure 11 .
Figure 11.Contours of accumulation body of numerical simulation and of Greenland.

Figure 12 .
Figure 12.Curves of dynamic parameter of monitoring points versus time (a) Displacement of monitoring points (b) Velocity of monitoring points.