A BP neural network-based micro particle parameters calibration and an energy criterion for the application of strength reduction method in MatDEM to evaluate 3D slope stability

<!-- metadata commented in wiki content


'''A BP neural network-based micro particle parameters calibration and an energy criterion for applying strength reduction method in MatDEM to evaluate 3D slope stability'''

Wei Jiang*<sup>1, 2</sup>, Yihong Tan<sup>1</sup>, Jinzhou Yan<sup>3</sup>, Ye Ouyang<sup>1</sup>, Zhaoyu Fu<sup>1</sup>, Qiang Feng<sup>1</sup>

1.Hubei Key Laboratory of Disaster Prevention and Mitigation, China Three Gorges University, Yichang, Hubei, 443002, China

2. Key Laboratory of Geological Hazards on Three Gorges Reservoir Area, Ministry of Education, Yichang, Hubei, 443002, China

3. Yichang Construction Investment & Development Co., Ltd., Yichang, Hubei, 443002, China

Correspondence: Wei Jiang ([mailto:jiangweilion@ctgu.edu.cn jiangweilion@ctgu.edu.cn])
-->
==Abstract==

To enhance the applicability of discrete element method in 3D slope stability analysis, a BP neural network-based micro parameter calibration method and an energy criterion are proposed by taking MatDEM as an example. Firstly, the relationship between the micro particle parameters and the shear strengths of particle aggregate are represented by using the BP neural network. And then the micro particle parameters are obtained for the given shear strengths by using a correction calibration. Next, the energy conversions are investigated for the stable and instable slope models in MatDEM. From a view of practical application, the abrupt in variation tendency and magnitude of the kinetic energy is selected for indicating the emergence of the limit equilibrium state of a slope. Finally, the effectiveness of the proposed improvements is testified by taking Baijiabao landslide as an example. Results verify that the calibration method established in this study is applicable to provide the micro particle parameters when the shear strength is constantly reduced, and the factor of safety determined by the kinetic energy criterion reflects the landslide stability at the global level.

'''Keywords''': 3D slope stability analysis, MatDEM, strength reduction method, micro parameters, failure criterion

==1. Introduction==

Slope stability analysis is a classic problem in geotechnical engineering. Three dimensional (3D) analyses has some significant advantages over two dimensional (2D) analyses in the slope stability analysis that the geometry structure, boundary conditions and spatial variations of a slope is fully took into account. The limit equilibrium method (LEM) mainly developed in the last century is still the most popular methodology among engineers. A considerable amount of research has developed 3D LEM from different point of views such as the method of columns, rigorous solutions [1] or utilizing lower bound theorems [2]. The progress in the development of 3D LEM theory till 2021 has been sumarized by Firincioglu [3]. To achieve an accurate result by 3D LEM, slope and slip surface geometry including geological structures must be precisely determined in the model, which is not always possible in natural instabilities with complex geology. 

With the advances of computers, numerical analysis methods have become affordable to engineers, such as the finite element method (FEM) [4-5], discrete element method (DEM) [6-7] and others [8-15]. These numerical methods fall into two categories, continuum-based and discontinuum-based numerical methods. The latter kind of methods, e.g. DEM, overcome the limitation of the assumption of macro-continuity in the former kind of methods, and are capable of simulating the generation and development process of the slip surface. Due to these advantages, DEM is becoming increasing popular for slope stability analysis in situation where the continuum-based numerical model is inapplicable [16-18].

When analyzing the slope stability using numerical methods, the strength reduction method (SRM) [19] is a commonly used strategy for determining the factor of safety (FOS) of a slope. The SRM was initially proposed for FEM analysis of slope stability [20]. SRM continuously reduces the shear strength of soil and rock until a slope reaches a limit equilibrium, and then the reduction factor is considered as the FOS of a slope. Theoretically, SRM can be apllied in DEM similarity to the application of SRM in FEM. But, two critical issues are expected to be addressed when analyzing the slope stability by applying SRM in DEM.

The first critical issue is how to reduce the shear strength of soil and rock in the DEM modeling. DEM simulates the mechanical properties of soil and rock from the micro point of view, and the relationship between the macro parameters of soil and rock and the micro parameters of particles is complicated. Theoretical and numerical approaches were used to reveal this relationship for soil and rock simulated by different models, e.g. the bonded particle model [21-22], the flat jointed model [23], the close-packed lattice model [24], and the coarse grained particle model [25]. But, a rigid theory is still not available to ensure the magnitude of the relationship between macro parameters and micro parameters [16]. The trial-and-error calibration remains the usual approach to obtain the micro parameters of particles corresponding to the macro parameters [26], even if it exhibits high computing costs. In slope stability analysis by using SRM in DEM, the shear strength is constantly reduced, which means that the micro parameters are necessary to be adjusted correspondingly. To avoid the cumbersome calibration on the micro parameters, the shear strength reduction is usually simplified as the reduction of several micro parameters [27-28] based on the regression analysis of the effect of micro parameter on macro property. The simplification may results in the FOS deviating from its original definition, in view of the complicated relationship between the shear strength and the micro parameters. 

With the rapid development of neural network technology, recently many researchers attempted to establish an intelligent method to represent the relationship between macro and micro parameters for DEM simulation. Albeit some sophisticated algorithms have been employed for accomplishing this object, e.g. the improved simulated annealing algorithm [29], a sequential quasi-Monte Carlo filter [30] and the non-dominated sorting genetic algorithm [31], the backpropagation (BP) neural network is still the most commonly used algorithm in the previous attempts [32-35] due to its simplicity and widely available software. However, BP neural network is possible to obtain a result with a low precision when the output layer has more neurons than the input layer, even if a considerable amount of samples are collected [35]. Therefore, BP neural network has potential to be used in the reduction of the shear strength in the DEM modeling, while a reliable way is to be developed.

The second critical issue is how to determine the FOS when applying SRM in DEM to analyze the slope stability. The target is to establish an appropriate criterion to indicate that a slope is in the limit equilibrium state. When analyzing the slope stability by using SRM in FEM, whether the slope is in the limit equilibrium state is usually evaluated by using three criteria: (1) penetration of the plastic zone in the slope, (2) convergence of the unbalance force, and (3) an abrupt change in the displacements of selected characteristic points. These criteria cannot be simply duplicated into the slope stability analysis by using SRM in DEM. In many DEM modelings, the particles are treated as rigid bodies, and thus there will be no plastic zone, which induces that it is impossible to employ the first criterion for the DEM analysis of the slope stability. In case some particles fall or a small region of the slope is deformed, the unbalance forces in the DEM modeling will not convergence, while actually the slope is still stable overall. Thus, it is unreasonable to judge whether the slope is unstable based on the second criterion [36]. Similarly, in views of the fact that the displacement characteristic of some selected points cannot reflect the overall behavior of the entire slope, it is improper to determine whether the entire slope is unstable based on the third criterion [37]. To avoid the limitations of the traditional criteria, the abrupt change in the vector graphic of global displacements [38], the variation coefficient of global displacements [36], and the slope failure extent [16] have been used as the failure criterion in the slope stability analysis by using SRM in DEM. These criteria are capable to represent the slope stability at the global level, but they are inapplicable in some cases. For example, the vector graphic of global displacements may be quite disordered for a complex slope, and it is difficult to evaluate significant changes. Thus, a reasonable criterion is expected to determine the FOS when analyzing the slope stability by using SRM in DEM.

The fast GPU matrix computing discrete element method (MatDEM) [39] is a newly developed code based on the DEM theory and MATLAB, which has been verified as an effective numerical simulations by many examples [40-43]. Since matrix operations and high-performance GPU calculations are used in the code, its computational efficiency is significantly improved. By taking MatDEM as an example, two improvements are proposed in this study to address the aforementioned two issues. Firstly, to precisely represent the reduction of the shear strengths of the geological material in the MatDEM modeling, the relationship between the macro shear strengths and the micro particle parameters is represented by an overdetermined BP neural network, and a correction calibration is developed to obtain the micro particle parameters corresponding to the original and reduced shear strengths. Then, to reasonably determine the FOS of a slope at the global level, the energy conversion and heat generation during the failure of a slope are fully investigated, and thus an energy criterion is developed for the application of SRM in MatDEM to evaluate the slope stability. Finally, the validity of the proposed strategies is testified by taking Baijiabao landslide situated in Zigui County, Hubei Province of China as a case study.

==2. Fundamental principles of the MatDEM==

===2.1 Contact model of particles===

The basic model of the MatDEM code is a series of 3D close-packed elastic particles as shown in [[#img-1|Figure 1]]a, and these particles interact through the spring forces ([[#img-1|Figure 1]]b).

{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
|
{|
|+
|-
|style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image1.png|250x250px]]
|-
|style="text-align: center;font-size: 75%;padding-bottom:10px;"|(a) 3D close-packed elastic particles
|-
|style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image2.png|474px]]
|-
|style="text-align: center;font-size: 75%;padding-bottom:10px;"|(b) Normal and tangential springs between two particles
|}
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 75%;"| '''Figure 1'''. The contact model of particles in MatDEM
|}


The normal force <math>F_n</math> and the normal relative displacement <math>X_n</math> between two adjacent particles are simulated by a normal spring:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>F_n=\left\{\begin{array}{c}
K_nX_n,\\
0,\\
K_nX_n,
\end{array}\mbox{ }\begin{array}{c}
X_n\leq X_b,\mbox{ intact }\\
X_n\geq 0,\mbox{ broken}\\
X_n<0,\mbox{ broken}
\end{array}\mbox{ }\begin{array}{c}
\qquad \mbox{(a)}\\
\qquad\mbox{(b)}\\
\qquad\mbox{(c)}
\end{array}\right.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

where <math>K_n</math> is the stiffness of the normal spring, and <math>X_b</math> is the breaking displacement. Initially, the particles are interconnected with their adjacent particles and subjected to tensile or compressive spring forces as Eq. (1a). In case of <math>X_n</math> between the particle pair exceeding <math>X_b</math>, the spring breaks and the tensile force ceases to exist between them (Eq. (1b)). However, the compressive force may act between them when they return to a compressive status (Eq. (1c)). It is noticed that the tensile force is positive, and the compressive force is negative in Eq. (1).

The shear force <math>F_s</math> and the tangential relative displacement <math>X_s</math> between two particles are simulated by a tangential spring:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>F_\mbox{s}=K_sX_s</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

where <math>K_s</math> is the stiffness of the tangential spring. The tangential spring also has a fracture criterion based on the Mohr-Coulomb yield criterion:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>F_{smax}=\left\{ \begin{array}{c}
F_{s0}-\mu_pF_n,\\
-\mu_pF_n,
\end{array} \mbox{ }\begin{array}{c}
\mbox{intact }\\
\mbox{ broken}
\end{array}\mbox{ }\begin{array}{c}
\qquad \mbox{(a)}\\
\qquad\mbox{(b)}
\end{array}\right.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

where <math>F_{smax}</math> is the inter-particle shear resistance, <math>F_{s0}</math> is the initial inter-particle shear resistance when no normal force exists, and <math>\mu_p</math> is the inter-particle coefficient of friction. When the tangential force <math>F_s</math> exceeds <math>F_{smax}</math> in Eq. (3a), the tangential spring will break, whereupon the shear force <math>F_s</math> is limited to be less than or equal to the inter-particle shear resistance <math>F_{smax}</math> of the broken spring. If the tangential spring is broken and the magnitude of external shear force exceeds the limit <math>F_{smax}</math> in Eq. (3b), two particles begin slipping, and the slipping friction between them is <math>F_{smax}</math> in Eq. (3b).

===2.2 Energy system of the MatDEM modeling===

The MatDEM modeling follows the law of energy conversation and Newton’s law of motion. The energies of a model fall into two categories, mechanical energy and heat. A brief introduction of their compositions is given as follows.

(1) The mechanical energy in the system consists of elastic potential energy, gravitational potential energy and kinetic energy. The elastic potential energy <math>E_e</math> is the sum of the strain energy of normal and tangential springs between particles:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>E_e=\frac{1}{2}K_nX_n^2+\frac{1}{2}K_sX_s^2</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}


The gravitational potential energy <math>E_g</math> of a particle is:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>E_g=mgh</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

where <math>m</math> is the particle mass, <math>g</math> denotes the gravity acceleration, and <math>h</math> represents the height above the reference level. The kinetic energy <math>E_k</math> of a particle is computed as:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>E_k=\frac{1}{2}mv^2</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}

in which <math>v</math> is the scalar of particle velocity.

(2) The heat in the system comes from the viscous damping, the spring breaking and the friction. Damping is used in the MatDEM to weaken the elastic wave energy of the model and dissipate the kinetic energy in the system. The damping force <math>\boldsymbol{F}_d</math> is given by

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{F}_d=-\eta \boldsymbol{v}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}

where <math>\eta</math> is the damping coefficient and <math>\boldsymbol{v}</math> is the vector of particle velocity. Because the time step of the simulation is very small, the particle velocity is assumed to be constant in a step. Viscous heat <math>Q_d</math> generated by viscous damping is calculated by

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>Q_d=-\boldsymbol{F}_d\cdot d\boldsymbol{x}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

in which <math>d\boldsymbol{x}</math> denotes the particle displacement in the current time step.

When an intact spring breaks, the spring force reduces, and the elastic potential energy of the inter-particle normal spring and/or tangential spring will dissipate into heat. Thus, breaking heat is equal to the reduction of the elastic potential energy <math>E_e</math>. The calculation of breaking heat is based on the status of the inter-particle normal force. If the inter-particle normal force is tensile, both the normal and tangential spring forces will reduce to zero, and thus the breaking heat <math>Q_b</math> is the sum of elastic potential energy of both normal and tangential springs:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>Q_b=\frac{1}{2}K_nX_n^2+\frac{1}{2}K_sX_s^2</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|}


If the normal force is compressive, the normal spring force will not be changed, and the tangential spring force reduced from <math>F_{smax}</math> in Eq. (3a) to <math>F_{smax}</math> in Eq. (3b). The breaking heat <math>Q_b</math> is the reduction of the elastic potential energy of the tangential spring, which can be calculated by

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>Q_b=\frac{\left({\left(F_{s0}-{\mu }_pF_n\right)}^2-{\left({\mu }_pF_n\right)}^2\right)}{2K_s}=</math><math>\frac{F_{s0}^2-2{\mu }_pF_nF_{s0}}{2K_s}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}


Friction heat <math>Q_f</math> generated during the sliding process is defined as the product of the average sliding friction and the effective sliding distance:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>Q_f=\vert 0.5\left(\boldsymbol{F}_{s1}+\boldsymbol{F}_{s2}\right)\cdot d\boldsymbol{S}\vert </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}

in which <math>\boldsymbol{F}_{s1}</math> and <math>\boldsymbol{F}_{s2}</math> are the sliding friction forces at respectively the beginning and the end of the current time step, and <math>d\boldsymbol{S}</math> denotes the effective sliding distance. More details about Eq. (11) can be found in [39].

The total energy of a MatDEM model is the sum of all mechanical energy and heat:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>E_{total}=E_e+E_g+E_k+Q</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

where <math>Q</math> represents the sum of viscous heat <math>Q_d</math>, breaking heat <math>Q_b</math> and friction heat <math>Q_f</math>:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>Q=Q_d+Q_b+Q_f</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|}


Based on the law of energy conversation, the total energy of an isolated system is constant. If the system deforms under external force, the increment of the total energy must be equal to the work done by the external force.

==3. A BP neural network-based micro particle parameters calibration==

===3.1 The structure of the BP neural network to predict the shear strength===

In this study, the BP neural network is used to establish the relationship between the shear strength and micro parameters for MatDEM modeling. BP neural network is a multilayer feed-forward neural network which consists of an input layer, implicit layer and output layer as shown in [[#img-2|Figure 2]] .

<div id='img-2'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
|
{|
|+
|-
|style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image16.png|354px]]
|}
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 75%;"| '''Figure 2'''. Typical structure of BP neural network
|}


When the value of neurons in the input layer is given, neurons in the implicit layer and the output layer are calculated according to the following forward propagation:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>u_i^l={\sigma }_i^l\left(\sum_{j=1}^{k_{l-1}}w_{ji}^lu_j^{l-1}+ b_i^l\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}

where <math>l</math> is the number of layers in the implicit layer, <math>u_i^l</math> denotes the <math>i</math>-th neuron in the <math>l</math>-th layer, <math>{\sigma }_i^l</math> is the activation function, <math>k_l</math> is the number of neurons in the <math>l</math>-th layer, <math>w_{ji}^l</math> is the weight and <math>b_i^l</math>  is the deviation. The training of BP neural network indicates constantly adjusting the weights and deviations to minimize the differences between the output value and the true value.

The structure of the BP neural network prediction model for the shear strength is plotted in [[#img-3|Figure 3]]. The micro parameters are placed in the input layer and the shear strengths are placed in the output layer. It is noticed that only the micro parameters in terms of the particle contact model is involved in the BP neural network. Although some other micro parameters, e.g. the particle radius and the particle density, and the particle packing pattern have been verified to influence the macro mechanical properties of geotechnical materials in the numerical model [31], their effect is ignored in this study. The reason is that the particle radius and the particle density will be constant when the numerical model of the slope has been established, and all particles are prescribed to be close-packed in MatDEM.

<div id='img-3'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
|
{|
|+
|-
|style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image18.png|396px]]
|}
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 75%;"| '''Figure 3'''. BP neural network prediction model for the shear strength
|}


Considering that an ultimate aim of this study is to acquire the micro parameters corresponding to the original and reduced shear strength, the established BP neural network seems unfavorable. The reasons are explained as follows. If placing the micro parameters in the output layer and the shear strengths in the input layer, the output layer have more neurons than the input layer, which forms an underdetermined mathematical problem. In this situation, the BP neural network is possible to result in a result of low precision. For example, to calibrate micro particle parameters for Particle Flow Code (PFC) modeling, a BP neural network was built by placing three macro geotechnical mechanical parameters in the input layer and four micro parameters in the output layer [35]. Although the BP neural network database consists of four hundred samples, results of some testing samples were still unsatisfactory. Thus, the micro parameters and the shear strengths are placed in the input layer and the output layer respectively, to promise the validity of the BP neural network. Base on the BP neural network predict model of the shear strength, a calibration of the micro particle parameters will be developed similar to the work in [34].

===3.2 Virtual 3D direct shear tests to measure the shear strength===

The direct shear test and the triaxial compression test are two commonly used tests to measure the shear strength of geotechnical materials in laboratory, and their virtual counterparts have been established in DEM simulation by researchers [30,44]. Under the same micro parameter conditions, the two virtual tests are possible to result in different values of the shear strength [44]. Considering that the shear strength of geotechnical materials in the latter case study is obtained by using direct shear test, virtual 3D direct shear test is used to measure the shear strength in this study.

As shows in [[#img-4|Figure 4]], the virtual 3D direct shear test model consists of two parts, the specimen and the shear box. The specimen is composed of 7636 spherical particles, whose radiuses vary from 0.83 mm to 1.20 mm. The radius of the specimen is 30 mm and the height is 20 mm. The shear box consists of 7374 walls. To simulate the shearing process, the imposed boundary conditions are the horizontal velocity of the lower box and the normal pressure on the top wall. By prescribing the horizontal displacement of the lower box in each time step to be 0.05 mm, the shear stress was recorded until the specimen was finally damaged.

<div id='img-4'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
|style="padding:10px;text-align:center;"| [[Image:Review_852894090255-image19.png|360px]]
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 75%;"| '''Figure 4'''. Virtual 3D direct shear test model for measuring the shear strength of particle aggregate
|}


Given the micro particle parameters as <math>K_n =5.0 \times 10^6</math> N/mm, <math>K_s =1.0\times 10^6</math> N/mm, <math>X_b =1.0\times 10^{-5}</math> m, <math>F_{s0}=1.0\times 10^6</math> Pa and <math>\mu_p=0.20</math>, four normal pressures, 100 kPa, 200kPa, 300kPa and 400 kPa, are applied on the top wall to perform the virtual 3D direct shear test. The shear stress−shear displacement curves under the four normal pressures are plotted in [[#img-5|Figure 5]]a. The trend of the four curves is consistent with each other. The increasing of the imposed normal pressures induces a higher peak value of the shear stress, and the shear displacement where the shear stress reaches its peak value increases slightly, which is similar to the result of the previous virtual direct shear test [44]. The Coulomb formula

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\tau =C + \sigma \tan\varphi </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|}

is employed to determine the shear strength of the particle aggregate as shown in [[#img-5|Figure 5]]b. Finally, the friction angle <math>\varphi</math> is 32.4°, and the cohesion <math>C</math> is 76 kPa.

<div id='img-5'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
|
{|
|-
| style="text-align:center;padding:10px;"| [[Image:Review_852894090255-image21.png|270px]]
|-
|style="text-align:center;font-size: 75%;padding-bottom:10px;"| (a) Shear stress−shear displacement curves
|-
|  align="center" padding:"10px"| [[Image:Review_852894090255-image22.png|246px]]
|-
|  style="text-align:center;font-size: 75%;padding-bottom:10px;"|(b) Determination of the shear strengths
|}
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 75%;"| '''Figure 5'''. Virtual 3D direct shear test model for the shear strength measurement
|}

===3.3 A correction calibration of micro particle parameters based on the BP neural network===

Based on the previous studies on the influence of the micro particle parameters in MatDEM on the mechanical properties of the particle aggregate [39], four levels are considered for <math>F_{s0}</math> and <math>\mu_p</math>, and three levels are considered for <math>K_n</math>, <math>K_s</math> and <math>X_b</math> when establishing the BP neural network database. [[#tab-1|Table 1]] lists the level values of micro particle parameters. By using the full orthogonal combination, 432 combinations of micro particle parameters are generated. And then, the virtual 3D direct shear test is executed to obtain their corresponding shear strengths. Finally, the BP neural network database is composed by 432 groups of micro particle parameters and the corresponding shear strengths.

<div class="center" style="font-size: 75%;">'''Table 1'''. The level values of micro particle parameters to establish the neural network database</div>

<div id='tab-1'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;width:auto;" 
|-style="text-align:center"
! <math>K_s</math> (N/m) !! <math>K_n</math> (N/m) !! <math>F_{s0}</math> (Pa) !! <math>\mu_p</math> !! <math>X_b</math> (m) 
|- style="text-align:center"
| <math>1.0\times 10^5</math>, <math>1.0\times 10^6</math>, <math>1.0\times 10^7</math>
| <math>5.0\times 10^5</math>, <math>5.0\times 10^6</math>, <math>5.0\times 10^7</math>
| <math>1.0\times 10^5</math>, <math>5.0\times 10^5</math>, <math>1.0\times 10^6</math>, <math>5.0\times 10^6</math>
|<math>0.10, 0.20, 0.30, 0.40</math>
|<math>1.0\times 10^{-6}</math>, <math>1.0\times 10^{-5}</math>, <math>1.0\times 10^{-4}</math>
|}


The training of the BP neural network actually successfully establishes the mapping relationship between the micro particle parameters in the input layer and the shear strengths in the output layer, which can be denoted as the following formula.

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\left(C,\tan\varphi \right)=f\left(K_n,K_s,X_b,F_{s0},{\mu }_p\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
|}


Assuming that the shear strengths (outputs) are continuous in the micro particle parameter space, the micro particle parameters (inputs) corresponding to the target shear strengths could be resolved similarly to the resolution of the nonlinear equations, because the gradient at any given input could be computed in the BP neural network. In the first place, an objective function is defined as

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>Error={\left(\frac{C-C^\ast }{C^\ast }\right)}^2+{\left(\frac{\tan\varphi -\tan{\varphi }^\ast }{\tan{\varphi }^\ast }\right)}^2</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
|}

where <math>C^\ast</math> and <math>{\varphi }^\ast</math> are the target cohesion and the target friction angle, respectively. And then, the gradient descent method is used to adjust the inputs until the objective function have a value smaller than the prescribed threshold value. For a simplicity expression, the input (<math>K_n</math>, <math>K_s</math>, <math>X_b</math>, <math>F_{s0}</math>, <math>\mu_p</math>) is denoted by <math>\boldsymbol{X}</math>, the output (<math>C</math>, <math>\tan\varphi</math>) is denoted by <math>\boldsymbol{Y}</math>, and thus Eq. (16) can be written as <math>\boldsymbol{Y}=f(\boldsymbol{X})</math>. The resolution process is explained as follows:

(1) Determine the initial exploration input <math>\boldsymbol{X}_0</math> by using the K-means clustering analysis [45]. K-means clustering is a vector quantization method originally from signal processing and can partition <math>n</math> observations into <math>k</math> clusters where each observation belongs to the cluster with the nearest mean. The samples in the neural network database are partitioned into <math>k</math> clusters by using K-means clustering analysis. K is takens as 15 in this study. Distances between the cluster centers and the target are computed by Eq. (17), and the nearest cluster center is selected as <math>\boldsymbol{X}_0</math>.

(2) Calculate the response outputs <math>\boldsymbol{Y}_i</math> for the exploration point of <math>\boldsymbol{X}_i</math> at step <math>i</math>, and then invoke the objective function. When the objective function has a value smaller than the threshold value, the resolution is terminated and <math>\boldsymbol{X}</math> takes the value of <math>\boldsymbol{X}_i</math>. If not, continue to perform step (3).

(3) Calculate the gradient of the objective function at <math>\boldsymbol{X}_i</math>, and update <math>\boldsymbol{X}</math> along the descent direction of the objective function. The formulation to update <math>\boldsymbol{X}</math> is

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\Delta \boldsymbol{X}=-\eta \frac{\partial \left(Error\right)}{\partial \boldsymbol{X}}\vert _{\boldsymbol{X}=\boldsymbol{X}_i},\boldsymbol{X}_{i+1}=</math><math>\boldsymbol{X}_i+\Delta \boldsymbol{X}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
|}

in which <math>\Delta \boldsymbol{X}</math> is the increment used to update <math>\boldsymbol{X}_i</math> and <math>\eta</math> is the step length. <math>\Delta \boldsymbol{X}</math> is obtained by multiplying the negative gradient of the objective function by <math>\eta</math>. On account of the nonlinearity of forward propagation rule, <math>\eta</math> is introduced for performing a step search approximation along the error gradient descent direction.

(4) Let <math>i=i+1</math> and return to step (2).

Forward propagation of the established BP neural network can be regarded as a substitution of the virtual 3D direct shear test. The samples in the neural network database are always finite, which indicates that the substitution is certain to be approximate. Therefore, the inputs obtained from the mathematical resolution may have a low accuracy. An additional step called as “verifying test” is used to determine whether the inputs are acceptable. The verifying test performs virtual 3D direct shear test by using the inputs obtained from the mathematical resolution, and calculates the error between the resulted shear strength and the target by using Eq. (17). If the objective function has a value smaller than the threshold value, the inputs are acceptable. Otherwise, a correction should be introduced for obtain a more reasonable input.

The failure of the verifying test indicates that the current BP neural network is still not accurate enough to obtain a fine mathematical resolution. A feasible strategy to improve the precision of BP neural network is adding some new samples to the database. Although the current BP neural network is not accurate enough, it provides a rough search direction for the resolution and the inputs resulted by the last round of the mathematical resolution must be a bit closer to the optimal solution than the initial inputs. Considering that, a new sample consisting of the resulted inputs in the last round and the corresponding shear strength is added to the database, and then the neural network is retrained. After the retraining of the BP neural network, the mapping relationship in Eq. (16) has been replaced by a new relationship. Then, the mathematical resolution is executed again by selecting the resulted inputs in the last round as the initial exploration input, and the verifying test is performed again. The “mathematical resolution−verifying test−neural network correction” flow is repeated till the objective function has a value smaller than the threshold value in verifying test. The entire flow of the correction method to calibrate the micro particle parameters is plotted in [[#img-6|Figure 6]] .

<div id='img-6'></div>

{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
|
{|
|+
|-
|style="text-align: center;padding-bottom:10px;"| [[Image:Review_852894090255-image26.png|600px]]
|}
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 75%;"| '''Figure 6'''. The entire flow of the correction method to calibrate the micro particle parameters
|}

===3.5 Validity of the proposed correction calibration===

To verify the validity of the proposed calibration method, four different target shear strengths are given, and the optimal micro particle parameters are searched by prescribing the threshold value as <math>2.0\times 10^{-4}</math> for Eq. (17). After three or four rounds of mathematical resolution, the resulted micro particle parameters passed the verifying test. [[#tab-2|Table 2]] lists the micro parameters provided by the mathematical resolution and the shear strengths provided by the verifying test in each round.

<div class="center" style="font-size: 75%;">'''Table 2'''.  Micro particle parameters and shear strengths in each round</div>

<div id='tab-2'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;width:auto;" 
|-style="text-align:center"
! Target shear strength !! Round !! <math>K_n</math> (N/m) !! <math>K_s</math> (N/m) !! <math>X_b</math> (m) !! <math>F_{s0}</math> (Pa) !! <math>\mu_p</math> !! <math>C</math> (kPa) !! <math>\varphi (^o)</math>
|-style="text-align:center"
|  rowspan='4' | <math>C^\ast=30.0</math> kPa <br> <math>{\varphi }^\ast =19.6^o</math>
|  1
|  <math>0.27\times 10^7</math>
|  <math>5.27\times 10^6</math>
|  <math>1.95\times 10^{-5}</math>
|  <math>1.25\times 10^{5}</math>
|  <math>0.320</math>
|  <math>21.7</math>
|  <math>30.2</math>
|-style="text-align:center"
|  <math>2</math>
|  <math>0.47\times 10^7</math>
|  <math>1.84\times 10^6</math>
|  <math>2.72\times 10^{-5}</math>
|  <math>3.67\times 10^{5}</math>
| <math> 0.220</math>
|  <math>28.4</math>
|  <math>22.8</math>
|-style="text-align:center"
|  <math>3</math>
|  <math>2.03\times 10^7</math>
|  <math>2.32\times 10^6</math>
|  <math>1.84\times 10^{-5}</math>
|  <math>6.42\times 10^{5}</math>
|  <math>0.143</math>
|  <math>29.5</math>
|  <math>18.9</math>
|-style="text-align:center"
|  <math>4</math>
|  <math>1.27\times 10^7</math>
|  <math>2.09\times 10^6</math>
|  <math>2.04\times 10^{-5}</math>
|  <math>5.95\times 10^{5}</math>
|  <math>0.141</math>
|  <math>30.2</math>
|  <math>19.6</math>
|-style="text-align:center"
| rowspan='3' |<math>C^\ast=28.5</math> kPa <br> <math>{\varphi }^\ast = 18.6^o</math>
|<math> 1 </math>
|<math> 2.15\times 10^7 </math>
|<math> 4.92\times 10^6 </math>
|<math>2.36\times 10^{-5}</math>
|<math>3.29\times 10^{5}</math>
|<math>0.114</math>
|<math>26.3</math>
|<math>16.7</math>
|-style="text-align:center"
|<math>2</math>
|<math>1.63\times 10^7 </math>
|<math>4.40\times 10^6 </math>
|<math>2.02\times 10^{-5}</math>
|<math>3.17\times 10^{5}</math>
|<math>0.127</math>
|<math>29.2</math>
|<math>17.9</math>
|-style="text-align:center"
|<math>3</math>
|<math>2.38\times 10^7 </math>
|<math>4.75\times 10^6 </math>
|<math>1.61\times 10^{-5}</math>
|<math>4.52\times 10^{5}</math>
|<math>0.129</math>
|<math>28.6</math>
|<math>18.5</math>
|-style="text-align:center"
| rowspan='4' |<math>C^\ast=45.0</math> kPa <br> <math>{\varphi }^\ast =28.0^o</math>
|<math>1</math>
|<math>4.03\times 10^7 </math>
|<math>3.82\times 10^6 </math>
|<math>5.25\times 10^{-5}</math>
|<math>2.07\times 10^{5}</math>
|<math>0.324</math>
|<math>48.9</math>
|<math>31.5</math>
|-style="text-align:center"
|<math>2</math>
|<math>2.94\times 10^7 </math>
|<math>4.43\times 10^6 </math>
|<math>3.85\times 10^{-5}</math>
|<math>2.18\times 10^{5}</math>
|<math>0.296</math>
|<math>45.7</math>
|<math>29.7</math>
|-style="text-align:center"
|<math>3</math>
|<math>3.16\times 10^7 </math>
|<math>5.26\times 10^6 </math>
|<math>4.19\times 10^{-5}</math>
|<math>1.95\times 10^{5}</math>
|<math>0.268</math>
|<math>45.3</math>
|<math>28.9</math>
|-style="text-align:center"
|<math>4</math>
|<math>3.25\times 10^7 </math>
|<math>4.47\times 10^6 </math>
|<math>4.37\times 10^{-5}</math>
|<math>2.23\times 10^{5}</math>
|<math>0.272</math>
|<math>44.8</math>
|<math>27.8</math>
|-style="text-align:center"
| rowspan='4' |<math>C^\ast =12.0 </math> kPa <br> <math>{\varphi }^\ast =33.0^o</math>
|<math>1</math>
|<math>3.27\times 10^7 </math>
|<math>4.13\times 10^6 </math>
|<math>1.74\times 10^{-5}</math>
|<math>2.95\times 10^{5}</math>
|<math>0.330</math>
|<math>18.5</math>
|<math>30.5</math>
|-style="text-align:center"
|<math>2</math>
|<math>3.65\times 10^7 </math>
|<math>3.29\times 10^6 </math>
|<math>1.86\times 10^{-5}</math>
|<math>3.20\times 10^{5}</math>
|<math>0.306</math>
|<math>16.1</math>
|<math>28.5</math>
|-style="text-align:center"
|<math>3</math>
|<math>2.45\times 10^7 </math>
|<math>2.94\times 10^6 </math>
|<math>1.66\times 10^{-5}</math>
|<math>2.41\times 10^{5}</math>
|<math>0.358</math>
|<math>14.3</math>
|<math>32.0</math>
|-style="text-align:center"
|<math>4</math>
|<math>2.36\times 10^7 </math>
|<math>3.15\times 10^6 </math>
|<math>1.57\times 10^{-5}</math>
|<math>2.46\times 10^{5}</math>
|<math>0.379</math>
|<math>11.9</math>
|<math>33.2</math>
|}


Results in [[#tab-2|Table 2]] shows that <math>Error</math> in Eq. (17) decreases whith the mathematical resolution round increasing, which implies that the target shear strengths are approximated gradually. Taking the case <math>C^\ast</math>=30.0 kPa, <math>{\varphi }^\ast=19.6^o</math> as an example, the relative errors of the resulted shear strengths are respectively 27.7% and 63.5% for the cohesion and the friction coefficient in the first round, by using the target shear strengths as the reference. In the second round, they decrease to 5.3% and 18.1% for the cohesion and the friction coefficient, respectively. Ultimately, they decrease to 0.7% and 0.2% in the fourth round for the cohesion and the friction coefficient, respectively. Thus, the proposed method has the ability to obtain the micro particle parameters that well reflects the shear strengths of particle assembly.

When performing the slope stability analysis by applying SRM in MatDEM, the shear strengths of the geological material in the MatDEM modeling are constantly adjusted. Once a new value of the shear strengths is given, a new combination of micro particle parameters is to be searched. At the moment, the established neural network database will provide a new initial exploration input <math>\boldsymbol{X}_0</math> by using the K-means clustering analysis, and then the proposed method is performed to obtain the micro particle parameters.

==4. An energy criterion to evaluate the slope stability==

As mentioned in Section 2.2, the MatDEM modeling follows the law of energy conversation and Newton’s law of motion. When slope instability happens, the displacements of soil and rock will induce the conversion between various kinds of energies. Consequently, the energy conversion of an instable slope in MatDEM must be different with that of a stable slope. In this section, a simple 2D slope is taken as an example to investigate the energy conversion difference between the stable and instable slope. And then, an energy criterion is established to evaluate the slope stability.

===4.1 Investigation on the energy conversion difference between the stable and instable slope===

As  shows in [[#img-7|Figure 7]], the slope is 7 m height and the ratio of slope is 7: 3. And, the slope is assumed to be homogeneous and composed by one kind of geotechnical material for simplicity. The MatDEM model of the slope contains 6305 particles. The radius of particles varies from 0.064 m to 0.096 m. The density is taken as 2.04&#x00d7;10<sup>3</sup> kg/m<sup>3</sup>, and the gravity acceleration <math>g= 9.8</math> m/s<sup>2</sup>. All particles at the left, right and bottom boundaries are fixed as the boundary conditions of the model.

<div id='img-7'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
|style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image27.png|300px]]
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 75%;"| '''Figure 7'''. A simple homogeneous slope 2D model
|}


In the numerical simulation of the slope model, three groups of micro particle parameters listed in [[#tab-3|Table 3]] are selected, which are denoted as Case 1, 2 and 3 respectively. Considering that the chief purpose here is to reveal the energy conversion difference between the stable and instable slope, the shear strengths are not measured for the three groups of micro particle parameters by using the virtual numerical test.

<div class="center" style="font-size: 75%;">'''Table 3'''. Micro particle parameters to be used in the investigation of energy conversion</div>

<div id='tab-3'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;width:auto;" 
|-style="text-align:center"
! Parameter !! <math>K_n</math> (N/m) !! <math>K_s</math> (N/m) !! <math>X_b</math> (m) !! <math>F_{s0}</math> (Pa) !! <math>\mu_p</math>
|-style="text-align:center"
|  Case 1
|  <math>7.80\times 10^6</math>
|  <math>3.10\times 10^6</math>
|  <math>1.10\times 10^{-4}</math>
|  <math>5.46\times 10^5</math>
|  <math>0.020</math>
|-style="text-align:center"
|  Case 2
|  <math>2.80\times 10^8</math>
|  <math>1.50\times 10^8</math>
|  <math>2.60\times 10^{-5}</math>
|  <math>4.90\times 10^5</math>
|  <math>0.120</math>
|-style="text-align:center"
|  Case 3
|  <math>9.60\times 10^7</math>
|  <math>2.80\times 10^7</math>
|  <math>9.30\times 10^{-6}</math>
|  <math>3.90\times 10^5</math>
|  <math>0.309</math>
|}


After 10 time steps of the computation, the slope is verified to be instable in Case 1 but stable in Case 2 and 3. It is noticed that one time step in MatDEM consists of a number of sub-time steps, and the value of a sub-time step is self-determined by the software. In the computation of the slope model, one time step is composed of 63 sub-time steps, and a sub-time step takes 0.0001<math>s</math>, 0.00015<math>s</math> and 0.0002<math>s</math> for Case 1, Case 2 and Case 3, respectively. The displacement cloud charts are plotted in [[#img-8|Figure 8]]a and [[#img-9|Figure 9]]a for Case 1 and Case 2, respectively. The displacement cloud chart of Case 3 is not provided because the displacement difference between Case 2 and Case 3 is not noticeable. [[#img-8|Figure 8]]b and [[#img-9|Figure 9]]b illustrate the variations of the total energy, elastic potential energy, kinetic energy, gravitational potential energy and heat with the time steps for Case 1 and Case 2, respectively. Some interesting phenomena are observed in the energy-time steps curves, and a further investigation is valuable for the establishment of an energy criterion to evaluate the slope stability.

<div id='img-8'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
|
{|
|+
|-
| style="align:center;padding:10px;"| [[Image:Review_852894090255-image28-c.png|324px]]
|-
|style="align:center;font-size: 75%;padding-bottom:10px;"|(a) Displacement cloud chart
|-
|style="align:center;padding:10px;"| [[Image:Review_852894090255-image29.png|312px]]
|-
|style="align:center;font-size: 75%;padding-bottom:10px;"|(b) Energy variations with time step
|}
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 75%;"| '''Figure 8'''. Displacements and the energy conversion for the slope in Case 1
|}


<div id='img-9'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
|
{|
|+
|-
|style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image30-c.png|330px]]
|-
|style="text-align: center;font-size: 75%;padding-bottom:10px;"|(a) Displacement cloud chart
|-
|style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image31.png|306px]]
|-
|style="text-align: center;font-size: 75%;padding-bottom:10px;"|(b) Energy variations with time step
|}
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 75%;"| '''Figure 9'''. Displacements and the energy conversion for the slope in Case 2
|}


Firstly, the conclusion that MatDEM follows the law of energy conversation is confirmed again by the energy-time steps curves. For both two cases, the total energies are almost constant during the computation, no matter whether the slope is stable or not. However, the total energy of Case 2 has a larger value than that of Case 1, which can be explained as follows. In the MatDEM modeling, the gravitational potential energy is evaluated by the relative vertical displacement of particles. It is initialized to be zero for a particle, and becomes negative when the particle falls while positive when the particle rises. The kinetic energy is zero for a particle unless the motion of the particle happens. The heat induced by the viscous damping, the spring breaking and the friction is zero in the initial state. Therefore, the total energy should be equal to the elastic potential energy initially. Since the elastic potential energy is computed as Eq. (4), <math>K_n</math> and <math>K_s</math> dominate the value of the elastic potential energy. Because the values of <math>K_n</math> and <math>K_s</math> in Case 2 is greater than that in Case 1, the total energy of Case 2 has a larger value than that of Case 1.

Secondly, the variations of the kinetic energy and the gravitational potential energy in [[#img-8|Figure 8]]b and [[#img-9|Figure 9]]b are apparently different. With the time step increasing, the kinetic energy increases but the gravitational potential energy decreases in [[#img-8|Figure 8]]b. But, the gravitational potential energy and the kinetic energy are nearly constant in [[#img-9|Figure 9]]b. [[#img-10|Figure 10]]  plots the variations of the kinetic energy and the gravitational potential energy with time steps for all three cases.

<div id='img-10'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
|
{|
|+
|-
|style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image32.png|306px]]
|-
|style="text-align: center;font-size: 75%;padding-bottom:10px;"|(a) Kinetic energy
|-
|style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image33.png|306px]]
|-
|style="text-align: center;font-size: 75%;padding-bottom:10px;"|(b) Gravitational potential energy
|}
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 75%;"| '''Figure 10'''. Variations of the kinetic energy and the gravitational potential energy with time steps for all three cases
|}


For Case 3, the kinetic energy is observed, and the gravitational potential energy is positive. That may be caused by the slope deformation. The kinetic energy and the gravitational potential energy caused by the deformation can be easily distinguished from that caused by the failure. Above all, the magnitude of the energy induced by the deformation is tiny compared to that induced by the failure. For Case 3, the kinetic energy and the gravitational potential energy have peak values of 2800 J and 4700 J, respectively. For Case 2, they have peak values of 116 J and 945 J, respectively. Next, when the two energies are caused by deformation, their variations show the convergent trend. But, when the two energies are induced by the failure, their variations show the divergent trend. Since the clear difference in the variation of the kinetic energy and the gravitational potential energy can be observed between a stable slope and an instable slope, they have the potential to be used in the evaluation of the slope stability.

Finally, the heat is investigated. The variations of heat are illegible in [[#img-8|Figure 8]]b and [[#img-9|Figure 9]]b because the magnitude of heat is too small. [[#img-11|Figure 11]]a plots the variations of the heat with time steps for all three cases. Obviously, the heats in Case 2 and Case 3 are much greater than that in Case 1. But it is incorrect to have a conclusion that the heat in a stable MatDEM slope model is more than that in an instable MatDEM slope model. The heat in the system is composed of the contributions of the viscous damping, the spring breaking and the friction force, in which the latter two sources are calculated based on the micro particle parameters. In [[#img-11|Figure 11]]a, the heat curve for Case 1 is situated in the lowest position, which maybe resulted by the micro particle parameters in [[#tab-3|Table 3]]. Thus, the magnitude of the heat is unreasonable to be chosen as a reference in the slope stability assessment. But, the shapes of the three heat curves are quite different. For a better demonstration, the heat curve is repainted for Case 1 in [[#img-11|Figure 11]]b. The heat curves of Case 2 and Case 3 are roughly convex, but the heat curve of Case 1 is concave. If the slope is stable, the position change of the particles induced by the deformation will stop rapidly, which promises the heat to be stable after several time steps. If the slope is sliding, the positions of the particles are still continuously changed, which results in the on-going increasing of the heat. Therefore, the shape of the heat curve may be effective in the evaluation of the slope stability.

<div id='img-11'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
|
{|
|+
|-
|style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image34.png|312px]]
|-
|style="text-align: center;font-size: 75%;padding-bottom:10px;"|(a) Three cases
|-
|style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image35.png|300px]]
|-
|style="text-align: center;font-size: 75%;padding-bottom:10px;"|(b) Case 1
|}
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 75%;"| '''Figure 11'''. Variations of the heat with time steps [68]
|}

===4.2 A criterion based on the kinetic energy of the MatDEM model===

As mentioned in the Section 4.1, the variation of the kinetic energy and the gravitational potential energy and the shape of the heat curve are verified to be different between a stable slope and an instable slope when performing the simulation by MatDEM. Their validities are compared in this section from a view of practical application.

When performing the slope stability analysis by using SRM in MatDEM, the micro parameters of the particles will be modified according to the reduction factor. Both the gravitational potential energy and the kinetic energy are independent of the micro particle parameters. However, the heat is calculated based on the micro particle parameters. When the energy variations resulted by using different micro particle parameters are collected together, it is more suitable to compare the results of the gravitational potential energy and the kinetic energy than that of the heat. Additionally, the number of the time steps has a great influence on the reliability of the conclusion if estimating the stable state of a slope by using the shape of the heat curves. As shows in [[#img-11|Figure 11]]a, the overall trend of the heat curve is convergent for Case 2, while a slight turning happens at the sixth time step. A large number of time steps seem to be feasible to overcome the effect of the local abnormality, but it will cause another issue. In fact, the heat curve must be convergent ultimately no matter whether the slope is stable or not. In case of an instable slope, the position of the particles will be relocated finally when the sliding is terminated, and thus the heat will be also stable. Therefore, an incorrect conclusion may be drawn based on the shape of the heat curve when the number of time steps is inappropriate. In short, the gravitational potential energy and the kinetic energy are more suitable to be compared than the heat when estimating the stable state of a slope.

There is an essential association between the gravitational potential energy and the kinetic energy. The motion of the particles results in the change of their positions, so the kinetic energy always emerges before the gravitational potential energy in MatDEM. This deduction can be verified by the variations of the kinetic energy and the gravitational potential energy for Case 3. As shown in [[#img-10|Figure 10]] , for Case 3 it is the fourth time step that the kinetic energy begins to decrease, while it is the seventh time step that the gravitational potential energy begins to decrease. This phenomenon indicates that the convergence of the kinetic energy should be earlier than that of the gravitational potential energy for a stable slope. Moreover, the variation of the gravitational potential energy is noticed to be small when the time step ranges from four to seven for Case 1. The reason is that the motion of the particles in the model maybe temporarily dominated by the horizontal motion during the failure, but only the vertical motion of the particles induces the variation of the gravitational potential energy. Thus, the gravitational potential energy may have a more complicated curve than the kinetic energy. Based on the above analysis, the kinetic energy is taken as the main reference to estimate the stable state of a slope.

When performing the slope stability analysis by using SRM in MatDEM, the variations of the kinetic energy are collected under different reduction factors. Then, the variation tendencies and the magnitude of the kinetic energies are compared. The abrupt in the variation tendencies and the magnitude is considered as a sign for the emergence of the critical point, and the safe factor of the slope is determined finally. At this point, a criterion has been proposed to evaluate the slope stability based on the kinetic energy.

==5. Case study: Baijiabao landslide==

===5.1 Geological background===

Baijiabao landslide is situated in Guizhou village of Zigui County, Hubei Province of China on the right bank of the Xiangxi River (30<sup>o</sup>58′59.9′′N, 110<sup>o</sup>45′33.4′′E). A photograph of the landslide that was taken from the opposite bank of the Xiangxi River is shown in [[#img-12|Figure 12]]. The elevation of the landslide ranges from 125 m to 265 m above MSL, and the slope of the landslide ranges from 10<sup>o</sup> to 20<sup>o</sup>. The upper boundary of the landslide is defined by the interface between the bedrock and the soil. The left and right boundaries are defined by two natural gullies that contain many surface cracks. The toe of the landslide varies in elevation between 125 m to 135 m. The primary sliding direction of Baijiabao landslide varies between 75<sup>o</sup> and 85<sup>o</sup> (SW-NE). This landslide is 550 m long and 400 m wide, covers an area of 20 ha, and has an estimated volume of nearly 1 million m<sup>3</sup>. The landslide mass is primarily composed of loose Quaternary deposits. Its sliding zone is mainly silty clay. And the underlying bedrock contains quartz sandstone and argillaceous siltstone of the Early Jurassic Xiangxi Formation, which dips into the hill at an angle of 30<sup>o</sup> to 40<sup>o</sup> ([[#img-13|Figure 13]]).

<div id='img-12'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
|style="padding:10px;"| [[Image:Review_852894090255-image36.png|600px]]
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 75%;"| '''Figure 12'''. A photograph of Baijiabao landslide, taken from the opposite bank of the Xiangxi River
|}


<div id='img-13'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
|
{|
|+
|-
| style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image37.png|414px]]
|-
|  style="text-align: center;font-size: 75%;"|(a) Engineering geological map    
|-
|  style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image38.png|414px]]
|-
|  style="text-align: center;font-size: 75%;padding-bottom:10px;"|(b) Geological cross-section
|}
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 75%;"| '''Figure 13'''. Geological background of Baijiabao landslide [46]
|}


Baijiabao landslide poses significant threats to public safety. The Zi-Xing road crosses the central part of the landslide mass as [[#img-13|Figure 13]]a shows. Before Three Gorges Reservoir (TGR) impoundment, 165 residents lived in the landslide area, while only 20 residents live there today. The deformation rate of the landslide was observed to increase significantly after TGR impoundment. Monitoring system installed in October 2006 consists of four GPS monitoring sites, which are noted as ZG323, ZG324, ZG325 and ZG326 in [[#img-13|Figure 13]]a. The detailed deformation history of Baijiabao landslide can be referred to the related literatures [46-48].

===5.2 Landslide model in the MatDEM modeling===

This subsection describes the establishment of the 3D MatDEM model for the stability analysis of Baijiabao landslide. Firstly, the geomorphy of the landslide is simulated by using ArcGIS software based on the contour line measured in 2006 ([[#img-14|Figure 14]]a). Then, a 3D model of Baijiabao landslide is built in MatDEM. The model is 564 m long, 497 m wide and 186 m height, covering the whole landslide mass ([[#img-14|Figure 14]]b). Next, geotechnical materials are allocated according to the engineering geologic investigation data. The sliding zone in the model ([[#img-14|Figure 14]]c) is positioned by using the interpolation of locations where the sliding zone is encountered in the drilling exploration. The landslide mass has a thickness varying between 20 m to 30 m in the front area, and has a thickness varying between 10 m to 40 m in the trailing area.

<div id='img-14'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
|
{|
|+
|-
|style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image39.png|198px]]
|-
|  style="text-align: center;font-size: 75%;"|(a) Geomorphy model
|-
|  style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image40.png|360px]]
|-
|  style="text-align: center;font-size: 75%;"|(b) 3D MatDEM model
|-
|  style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image41.png|318px]]
|-
|  style="text-align: center;font-size: 75%;padding-bottom:10px;"|(c) An enlargement to exhibit the sliding zone
|}
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 75%;"| '''Figure 14'''. 3D MatDEM model for the stability analysis of Baijiabao landslide
|}


Finally, the 3D MatDEM model of Baijiabao landslide is composed by 1,207,985 particles, in which the landslide mass and the sliding zone are represented by 733,104 removable particles, and the bedrock is represented by 474,881 fixed particles. The radius of particles varies between 0.55 m and 0.83 m. In [[#img-14|Figure 14]]c, the landslide mass, the sliding zone and the bedrock are presented as brown particles, green particles and blue particles, respectively.

===5.3 Mechanical parameters and micro particle parameters===

According to the report of engineering geological exploration, the unit weight and shear strengths are listed in [[#tab-4|Table 4]] for the geotechnical materials in the landslide mass, the sliding zone and the bedrock, respectively. All parameters are measured in the native state, and taken their statistical mean value. It is noticed that the stability of Baijiabao landslide is actually affected by the rainfall and the reservoir water fluctuation. In view that the main purpose of this study is to improve the ability of MatDEM in 3D slope stability evaluation, the rainfall and the reservoir water fluctuation are ignored temporarily.

<div class="center" style="font-size: 75%;">'''Table 4'''. Unit weight and shear strengths of geotechincal materials in different regions</div>

<div id='tab-4'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;width:auto;" 
|-style="text-align:center"
! Region !! Unit weight <math>\gamma</math> (kN/m<sup>3</sup>) !! <math>C</math> (kPa) !! <math>\varphi</math> (<sup>o</sup>)
|-style="text-align:center"
|  Landslide mass
|  18.0
|  17.1
|  20.2
|-style="text-align:center"
|  Sliding zone
|  20.5
|  24.6
|  15.2
|-style="text-align:center"
|  Bedrock
|  26.0
|  1,200
|  32.1
|}


Using the BP neural network-based calibration method established in Section 3, micro parameters of particles representing various geotechnical materials are obtained and listed in [[#tab-5|Table 5]]. Because the particles have a much larger radius than the particles used in Section 3, the size of the virtual 3D direct shear test model is adjusted. Previous studies [21,49] verified the effect of the particle size on the macro properties in DEM modeling, and declared that the ratio <math>L</math>/<math>R</math> is a critical variable to be carefully selected when performing virtual strength test. Here <math>L</math> is the size of the strength test model, and <math>R</math> denotes the average radius of particles. According to Su’s suggestion [49], <math>L</math>/<math>R</math> should be more than 200 to avoid the effect of the particle size on the macro properties of the particle assembly in DEM. Finally, the specimen to be used in virtual 3D direct shear test has a radius of 21 m and a height of 14 m.

<div class="center" style="font-size: 75%;">'''Table 5'''. Resulted micro particle parameters and corresponding shear strangths</div>

<div id='tab-1'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;width:auto;" 
|-style="text-align:center"
! Region !! <math>K_n</math> (N/m) !! <math>K_s</math> (N/m) !! <math>X_b</math> (m) !! <math>F_{s0}</math> (Pa) !! <math>\mu_p</math> !! <math>C</math> (kPa) !! <math>\varphi</math> (<sup>o</sup>)
|-style="text-align:center"
|  Landslide mass
|  <math>3.27\times 10^7 </math>
|  <math>5.27\times 10^6</math>
|  <math>1.95\times 10^{-5}</math>
|  <math>1.25\times 10^{5}</math>
|  <math>0.220</math>
|  <math>17.0</math>
|  <math>20.2</math>
|-style="text-align:center"
|  Sliding zone
|  <math>2.49\times 10^7 </math>
|  <math>1.84\times 10^6</math>
| <math> 2.72\times 10^{-5}</math>
|  <math>3.67\times 10^{5}</math>
|  <math>0.120</math>
|  <math>24.6</math>
|  <math>15.2</math>
|-style="text-align:center"
|  Bedrock
|  <math>4.27\times 10^7 </math>
|  <math>2.09\times 10^6</math>
|  <math>2.04\times 10^{-5}</math>
|  <math>5.95\times 10^{5}</math>
|  <math>0.441</math>
|  <math>37.5</math>
|  <math>32.0</math>
|}


With respect to the bedrock, the micro particle parameters in [[#tab-5|Table 5]] resulted in the shear strengths different to the actual shear strengths in [[#tab-4|Table 4]]. The reason is that it is hard to represent the three geotechnical materials by using micro parameters of the same order of magnitude. The cohesion of the bedrock is nearly 100 times that of the landslide mass and the sliding zone. Because all the particles representing the bedrock are fixed in the model, the difference between the cohesion of the bedrock in [[#tab-4|Table 4]] and that in [[#tab-5|Table 5]] will have little influence on the computation. Moreover, to reduce the workload, the micro particle parameters of the bedrock will not be changed when the stability of the landslide is evaluated by using SRM.

===5.4 Computational results===

As concluded in Section 4, the abrupt in the variation tendencies and the magnitude of the kinetic energy can be utilized to determine the FOS when estimating the slope stability by using SRM in MatDEM. Under different reduction factors, the variations of the kinetic energy with time steps are recorded for the landslide model and plotted in [[#img-15|Figure 15]]. The kinetic energy curve resulted by a greater reduction factor always has a higher situation in [[#img-15|Figure 15]]. But, if the reduction factor is less than 1.06, the difference in the magnitude of the kinetic energy between different reduction factors is relatively small at each time step. When the reduction factor is increasing from 1.06 to 1.07, the magnitude difference has a dramatically increasing. Meanwhile, the kinetic energy curves present a convergent trend when the reduction factor is less than 1.06, but present a divergent trend when the reduction factor arrives at 1.07. Thus, the FOS of Baijiabao landslide is considered as 1.06.

<div id='img-15'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
|style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image42.png|318px]]
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 75%;"| '''Figure 15'''. Variation of kinetic energy with time steps under different reduction factors
|}


By using the proposed BP neural network-based calibration method, the micro particle parameters are constantly adjusted during the computation. [[#tab-6|Table 6]] lists the micro particle parameters and the shear strengths under different reduction factors.

<div class="center" style="font-size: 75%;">'''Table 6'''. Micro particle parameters and corresponding shear strangths under different reduction factors</div>

<div id='tab-6'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;width:auto;" 
|-style="text-align:center"
! Region !! FOS !! <math>K_n</math> (N/m)!! <math>K_s</math> (N/m) !! <math>X_b</math> (m) !!<math>F_{s0}</math> (Pa) !! <math>\mu_p</math> !! <math>C</math> (kPa) !! <math>\varphi</math> (<sup>o</sup>)
|-style="text-align:center"
|  rowspan='7' |Landslide mass
|  <math>1.00</math>
|  <math>3.27\times 10^7 </math>
|  <math>5.27\times 10^6</math>
|  <math>1.95\times 10^{-5}</math>
|  <math>1.25\times 10^{5}</math>
|  <math>0.220</math>
|  <math>17.0</math>
|  <math>20.2</math>
|-style="text-align:center"
|  <math>1.10</math>
|  <math>3.35\times 10^7 </math>
|  <math>4.89\times 10^6</math>
|  <math>1.61\times 10^{-5}</math>
|  <math>1.03\times 10^{5}</math>
|  <math>0.175</math>
|  <math>15.5</math>
|  <math>18.5</math>
|-style="text-align:center"
|  <math>1.03</math>
|  <math>3.29\times 10^7 </math>
|  <math>5.13\times 10^6</math>
|  <math>1.86\times 10^{-5}</math>
|  <math>1.18\times 10^{5}</math>
|  <math>0.207</math>
|  <math>16.5</math>
|  <math>19.7</math>
|-style="text-align:center"
|  <math>1.04</math>
|  <math>3.18\times 10^7 </math>
|  <math>4.97\times 10^6</math>
|  <math>1.82\times 10^{-5}</math>
|  <math>1.16\times 10^{5}</math>
|  <math>0.202</math>
|  <math>16.3</math>
|  <math>19.5</math>
|-style="text-align:center"
|  <math>1.05</math>
|  <math>3.22\times 10^7 </math>
|  <math>5.05\times 10^6</math>
|  <math>1.77\times 10^{-5}</math>
|  <math>1.14\times 10^{5}</math>
|  <math>0.198</math>
|  <math>16.2</math>
|  <math>19.3</math>
|-style="text-align:center"
|  <math>1.06</math>
|  <math>3.53\times 10^7 </math>
|  <math>5.02\times 10^6</math>
|  <math>1.75\times 10^{-5}</math>
|  <math>1.12\times 10^{5}</math>
|  <math>0.194</math>
|  <math>16.0</math>
|  <math>19.1</math>
|-style="text-align:center"
|  <math>1.07</math>
|  <math>3.43\times 10^7 </math>
|  <math>4.95\times 10^6</math>
|  <math>1.71\times 10^{-5}</math>
|  <math>1.09\times 10^{5}</math>
|  <math>0.190</math>
|  <math>15.9</math>
|  <math>19.0</math>
|-style="text-align:center"
|  rowspan='7' | Sliding zone
|  <math>1.00</math>
|  <math>2.49\times 10^7 </math>
|  <math>1.84\times 10^6</math>
|  <math>2.72\times 10^{-5}</math>
|  <math>3.67\times 10^{5}</math>
|  <math>0.120</math>
|  <math>24.6</math>
|  <math>15.2</math>
|-style="text-align:center"
|  <math>1.10</math>
|  <math>2.45\times 10^7 </math>
|  <math>1.81\times 10^6</math>
|  <math>2.36\times 10^{-5}</math>
|  <math>3.05\times 10^{5}</math>
|  <math>0.098</math>
|  <math>22.4</math>
|  <math>13.9</math>
|-style="text-align:center"
|  <math>1.03</math>
|  <math>2.69\times 10^7 </math>
|  <math>1.86\times 10^6</math>
|  <math>2.61\times 10^{-5}</math>
|  <math>3.48\times 10^{5}</math>
|  <math>0.113</math>
|  <math>23.9</math>
|  <math>14.8</math>
|-style="text-align:center"
|  <math>1.04</math>
|  <math>2.38\times 10^7 </math>
|  <math>1.83\times 10^6</math>
|  <math>2.57\times 10^{-5}</math>
|  <math>3.42\times 10^{5}</math>
|  <math>0.111</math>
|  <math>23.7</math>
|  <math>14.6</math>
|-style="text-align:center"
|  <math>1.05</math>
|  <math>2.33\times 10^7 </math>
|  <math>1.96\times 10^6</math>
|  <math>2.54\times 10^{-5}</math>
|  <math>3.36\times 10^{5}</math>
|  <math>0.109</math>
| <math>23.4</math>
|  <math>14.5</math>
|-style="text-align:center"
|  <math>1.06</math>
|  <math>2.41\times 10^7 </math>
|  <math>1.92\times 10^6</math>
|  <math>2.50\times 10^{-5}</math>
|  <math>3.29\times 10^{5}</math>
|  <math>0.107</math>
|  <math>23.2</math>
|  <math>14.4</math>
|-style="text-align:center"
|  <math>1.07</math>
|  <math>2.46\times 10^7 </math>
|  <math>1.87\times 10^6</math>
|  <math>2.47\times 10^{-5}</math>
|  <math>3.23\times 10^{5}</math>
|  <math>0.104</math>
|  <math>23.0</math>
|  <math>14.2</math>
|}


[[#img-16|Figure 16]] illustrates the final displacement of Baijiabao landslide model when the reduction factor is taken as 1.06. Some particles in the front area of the model have the displacements of nearly 60 m, and most particles have the displacements less than 10 m. Thus, the FOS determined by the abrupt in the variation tendencies and the kinetic energy magnitude in fact reveals the landslide stability at the global level. The displacements when the landslide is close to its limit equilibrium state may be valuable in the landslide failure prediction. The field data [46] verified that at the monitoring period from January 2007 to December 2019, the ZG326 site ([[#img-13|Figure 13]]a) in the middle area of Baijiabao landslide recorded the max cumulative displacement of 1.8 m. Its displacements are still observed when the reservoir water level changes and the rainfall happen, while Baijiabao landslide hasn’t slipped entirely till now. Therefore, the surface movements before its failure are uncertain. In views that this study aims to establish a micro particle parameter calibration method and a criterion for the application of SRM in MatDEM to analysis 3D slope stability, the influence of the reservoir water and the rainfall on the stability of Baijiabao landslide is left out of consideration. Hence, the displacements in the limit equilibrium state obtained in this study are not suitable for the failure prediction for Baijiabao landslide.

<div id='img-16'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
|style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image43-c.png|324px]]
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 75%;"| '''Figure 16'''. The displacement of Baijiabao landslide model with the reduction factor 1.06
|}

==6. Conclusions==

By taking MatDEM as an example, a calibration method based on BP neural network is proposed to acquire the micro particle parameters corresponding to the given shear strengths in this paper, and an energy criterion is proposed for determining the FOS of the slope when applying SRM in DEM. And the proposed enhancements are applied in the 3D stability analysis of Baijiabao landslide. The following conclusions can be made:

(1) BP neural network can well represent the complicated relationship between the micro particle parameters and the shear strengths of the particle aggregate in MatDEM. And the proposed BP neural network-based micro particle parameters calibration method is applicable to provide the micro particle parameters when the shear strength is constantly reduced.

(2) The variations of the kinetic energy, the gravitational potential energy and the heat are noticed to be quite different between the stable and instable slope. However, the abrupt in the variation tendency and the magnitude of the kinetic energy is more suitable for determining the FOS of the slope.

(3) When determining the FOS by using the kinetic energy criterion, the displacements of particles in the slope model are permitted, so the FOS actually reveals the slope stability at the global level.

==Acknowledgments==

The work is supported by the Open Research Programme of the Hubei Key Laboratory of Disaster Prevention and Mitigation (China Three Gorges University) (no. 2022KJZ07), the Open Research Fund of Key Laboratory of Geological Hazards on Three Gorges Reservoir Area of Ministry of Education (2020KDZ10), National Natural Science Funds (Project No. 52079070), and the Open Research Fund of State Key Laboratory of Simulation and Regulation of Water Cycle in River Basin (China Institute of Water Resources and Hydropower Research), Grant NO：IWHR-SKL-202020.

==References==

<div class="auto" style="text-align: left;width: auto; margin-left: auto; margin-right: auto;font-size: 85%;">

[1] Sun G., Jiang W., Cheng S.,  Zheng H. Optimization model for determining safety factor and thrust line in landslide assessments. International Journal of Geomechanics, 17(4):04016091, 2017.

[2] Deng D.P., Li L.,  Zhao L.H. LEM for stability analysis of 3D slopes with general-shaped slip surfaces. International Journal of Geomechanics, 17(10):06017017, 2017.

[3] Firincioglu B.S.,  Ercanoglu M. Insights and perspectives into the limit equilibrium method from 2D and 3D analyses. Engineering Geology, 281:105968, 2021.

[4] Lv Q., Liu Y.,  Yang Q. Stability analysis of earthquake-induced rock slope based on back analysis of shear strength parameters of rock mass. Engineering Geology, 228:39-49, 2017.

[5] Zheng H., Sun G.,  Liu D. A practical procedure for searching critical slip surfaces of slopes based on the strength reduction technique. Computers and Geotechnics, 36(1-2):1-5, 2009.

[6] Espada M., Muralha J., Lemos J.V., Jiang Q., Feng X.-T., Fan Q.,  Fan Y. Safety analysis of the left bank excavation slopes of Baihetan arch dam foundation using a discrete element model. Rock Mechanics and Rock Engineering, 51(8):2597-2615, 2018.

[7] Lu, Y., Tan, Y. and Li, X. Stability analyses on slopes of clay-rock mixtures using discrete element method. Engineering Geology, 244: 116-124, 2018.

[8] Nie Z., Zhang Z., Zheng H.,  Lin S. Stability analysis of landslides using BEM and variational inequality based contact model. Computers and Geotechnics, 123:103575, 2020.

[9] Sun G., Lin S., Zheng H., Tan Y., Sui T. The virtual element method strength reduction technique for the stability analysis of stony soil slopes. Computers and Geotechnics, 119:103349, 2020.

[10] Xu D., Wu A., Yang Y., Lu B., Liu F., Zheng H. A new contact potential based three-dimensional discontinuous deformation analysis method. International Journal of Rock Mechanics and Mining Sciences, 127:104206, 2020.

[11] Yang Y., Xia Y., Zheng H.,  Liu Z. Investigation of rock slope stability using a 3D nonlinear strength-reduction numerical manifold method. Engineering Geology, 292:106285, 2021.

[12] Wang B., Vardon P.J.,  Hicks M.A. Rainfall-induced slope collapse with coupled material point method. Engineering Geology, 239:1-12, 2018.

[13] Bui H.H., Fukagawa R., Sako K.,  Wells J.C. Slope stability analysis and discontinuous slope failure simulation by elasto-plastic smoothed particle hydrodynamics (SPH). Géotechnique, 61(7):565-574, 2011.

[14] He L., Tian Q., Zhao Z., Zhao X., Zhang Q., Zhao J. Rock Slope stability and stabilization analysis using the coupled DDA and FEM method: NDDA approach. International Journal of Geomechanics, 18(6):04018044, 2018.

[15] Sun L., Liu Q., Abdelaziz A., Tang X.,  Grasselli G. Simulating the entire progressive failure process of rock slopes using the combined finite-discrete element method. Computers and Geotechnics, 141:104557, 2022.

[16] Bao Y., Sun X., Chen J., Zhang W., Han X., Zhan J. Stability assessment and dynamic analysis of a large iron mine waste dump in Panzhihua, Sichuan, China. Environmental Earth Sciences, 78(2):48, 2019.

[17] Jiang M., Niu M., Zhang F., Wang H., Liao Z. Instability analysis of jointed rock slope subject to rainfall using DEM strength reduction technique. European Journal of Environmental and Civil Engineering, 26(10):4664-4686, 2021.

[18] Shi C., Yang B., Zhang Y., Yang J. Application of discrete-element numerical simulation for calculating the stability of dangerous rock mass: A case study. International Journal of Geomechanics, 20(12):04020231, 2020.

[19] Zienkiewicz O.C., Humpheson C.,  Lewis R.W. Associated and non-associated visco-plasticity and plasticity in soil mechanics. Geotechnique, 25(4):671-689, 1975.

[20] Griffiths D.V.,  Lane P.A. Slope stability analysis by finite elements. Géotechnique, 49(3):387-403, 1999.

[21] Yang B., Jiao Y., Lei S. A study on the effects of microparameters on macroproperties for specimens created by bonded particles. Engineering Computations, 23(6):607-631, 2006.

[22] Zhou Z.H., Wang H.N., Jiang M.J. Macro- and micro-mechanical relationship of the anisotropic behaviour of a bonded ellipsoidal particle assembly in the elastic stage. Acta Geotechnica, 16(12):3899-3921, 2021.

[23] Cheng Y., Wong L.N.Y. A study on mechanical properties and fracturing behavior of Carrara marble with the flat‐jointed model. International Journal for Numerical and Analytical Methods in Geomechanics, 44(6):803-822, 2020.

[24] Liu C., Pollard D.D.,  Shi B. Analytical solutions and numerical tests of elastic and failure behaviors of close-packed lattice for brittle rocks and crystals. Journal of Geophysical Research: Solid Earth, 118(1):71-82, 2013.

[25] Rackl M., Hanley K.J. A methodical calibration procedure for discrete element models. Powder Technology, 307:73-83, 2017.

[26] Hou J., Zhang M., Chen Q., Wang D., Javadi A., Zhang S.  Failure-mode analysis of loose deposit slope in Ya’an-Kangding Expressway under seismic loading using particle flow code. Granular Matter, 21(1):8, 2019.

[27] Su H., Fu Z., Gao A., Wen Z. Numerical simulation of soil levee slope instability using particle-flow code method. Natural Hazards Review, 20(2):04019001, 2019.

[28] Xu G.J., Zhong K.Z., Fan J.W., Zhu Y.J.,  Zhang Y.Q. Stability analysis of cohesive soil embankment slope based on discrete element method. Journal of Central South University, 27(7):1981-1991, 2020.

[29] Wang M.,  Cao P. Calibrating the micromechanical parameters of the PFC2D (3D) models using the improved simulated annealing algorithm. Mathematical Problems in Engineering, 2017:6401835, 2017.

[30] Cheng H., Shuku T., Thoeni K.,  Yamamoto H. Probabilistic calibration of discrete element simulations using the sequential quasi-Monte Carlo filter. Granular Matter, 20(1):11, 2018.

[31] Do H.Q., Aragón A.M., Schott D.L. A calibration framework for discrete element model parameters using genetic algorithms. Advanced Powder Technology, 29(6):1393-1403, 2018.

[32] Benvenuti L., Kloss C.,  Pirker S. Identification of DEM simulation parameters by Artificial Neural Networks and bulk experiments. Powder Technology, 291:456-465, 2016.

[33] He P., Fan Y., Pan B., Zhu Y., Liu J.,  Zhu D. Calibration and verification of dynamic particle flow parameters by the back-propagation neural network based on the genetic algorithm: recycled polyurethane powder. Materials, 12(20):3350, 2019.

[34] Ye F., Wheeler C., Chen B., Hu J., Chen K.,  Chen W. Calibration and verification of DEM parameters for dynamic particle flow conditions using a backpropagation neural network. Advanced Powder Technology, 30(2):292-301, 2019.

[35] Zhou, Y., Wu S., Jiao J.,  Zhang X. Research on mesomechanical parameters of rock and soil mass based on BP neural network. Rock Soil Mech, 32(12):3821-3826, 2011.

[36] Wang H., Zhang B., Mei G., Xu N. A statistics-based discrete element modeling method coupled with the strength reduction method for the stability analysis of jointed rock slopes. Engineering Geology, 264:105247, 2020.

[37] Deluzarche R.,  Cambou B. Discrete numerical modelling of rockfill dams. International Journal for Numerical and Analytical Methods in Geomechanics, 30(11):1075-1096, 2006.

[38] Shen H., Abbas S.M. Rock slope reliability analysis based on distinct element method and random set theory. International Journal of Rock Mechanics and Mining Sciences, 61:15-22, 2013.

[39] Liu C., Xu Q., Shi B., Deng S., Zhu H. Mechanical properties and energy conversion of 3D close-packed lattice model for brittle rocks. Computers & Geosciences, 103:12-20, 2017.

[40] Liu Y., Zhang D., Wang G.-y., Liu C., Zhang Y. Discrete element method-based prediction of areas prone to buried hill-controlled earth fissures. Journal of Zhejiang University-SCIENCE A, 20(10):794-803, 2019.

[41] Luo H., Xing A., Jin K., Xu S., Zhuang Y. Discrete element modeling of the Nayong rock avalanche, Guizhou, China constrained by dynamic parameters from seismic signal inversion. Rock Mechanics and Rock Engineering, 54(4):1629-1645, 2021.

[42] Scaringi G., Fan X., Xu Q., Liu C., Ouyang C., Domènech G., Yang F.,  Dai L. 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). Landslides, 15(7):1359-1375, 2018.

[43] Zhan Q., Wang S., Wang L., Guo F., Zhao D., Yan J. Analysis of failure models and deformation evolution process of geological hazards in Ganzhou city, China. Frontiers in Earth Science, 9:731447, 2021.

[44] Graziani A., Rossini C.,  Rotonda T. Characterization and DEM modeling of shear zones at a large dam foundation. International Journal of Geomechanics, 12(6):648-664, 2012.

[45] Peña J.M., Lozano J.A.,  Larrañaga P. An empirical comparison of four initialization methods for the K-Means algorithm. Pattern Recognition Letters, 20(10):1027-1040, 1999.

[46] Yao W.M., Li C.D., Guo Y.C., Criss R.E., Zuo Q.J.,  Zhan H.B. Short-term deformation characteristics, displacement prediction, and kinematic mechanism of Baijiabao landslide based on updated monitoring data. Bulletin of Engineering Geology and the Environment, 81(9):393, 2022.

[47] Xu S.L., Niu R.Q. Displacement prediction of Baijiabao landslide based on empirical mode decomposition and long short-term memory neural network in Three Gorges area, China. Computers & Geosciences, 111:87-96, 2018.

[48] Yao W.M., Li C.D., Zuo Q.J., Zhan H.B., Criss R.E. Spatiotemporal deformation characteristics and triggering factors of Baijiabao landslide in Three Gorges Reservoir region, China. Geomorphology, 343:34-47, 2019.

[49] Su, H., Yang J., Hu B., Gao X., Ma H. Study of particle size effect of rock model based on particle discrete element method. Rock and Soil Mechanics, 39(12):4642-4650, 2018.