An Enhanced Discrete Element Modeling Method Considering Spatiotemporal Correlations for Investigating Deformations and Failures of Jointed Rock Slopes

Abstract

The discrete element method (DEM) is commonly employed to analyze the deformations and failures of jointed rock slopes. However, when the iterative calculation process of the DEM modeling should be terminated is still unclear. To solve the above problem, in this paper, a discrete element modeling method based on the energy correlation coefficient is proposed to determine when the iterative calculation process could be terminated, and then applied the proposed method to analyze the deformations and failures of jointed rock slopes. Compared with the existing discrete element modeling method based on the displacement variation coefficient, the proposed method based on the energy correlation coefficient is much more applicable for jointed rock slopes. The main advantage of the proposed method is that there is no need to determine the position of the potential sliding surface, and the displacements of all blocks are no longer counted as statistics, but the spatiotemporal correlations between all blocks are considered. The effectiveness of the proposed method is verified by comparing with the existing method based on the displacement variation coefficient for an abbreviated jointed rock slope. Moreover, the proposed method is successfully applied to analyze a real-world jointed rock slope without an obvious potential sliding surface in which the existing method cannot work.

1. Introduction

The problem associated with the deformations and failures of slopes has long attracted the attention of those in geotechnical engineering fields [1]. Extensive research has been conducted on the deformations and failures of slopes in theory and practice. With the rapid development of computer technology, from the early limit equilibrium method (LEM) to the current numerical simulation methods, numerical calculation methods have been widely adopted to investigate the deformations and failures of slopes. The commonly used numerical calculation methods include the discrete element method (DEM) [2,3,4], finite element method (FEM) [5,6,7,8,9], finite difference method (FDM) [10,11], boundary element method (BEM) [12,13], material point method (MPM) [14,15], discontinuous deformation analysis (DDA) [16,17], and numerical manifold method (NMM) [18,19]. In addition, deep learning is also adopted to analyze typical geological disasters such as landslides, mudslides, and earthquakes [20].

At present, the FEM is the most commonly used method for calculating factors of safety, but it could not be used to reflect the characteristics of structural planes in rock masses, nor is it suitable for modeling large deformations of slopes such as landslides. Specifically, when the FEM is used to model large deformations of rock masses, the computational mesh may be distorted, and this may lead to inaccurate results. In the DEM, the rock blocks are modeled to interact with each other through corners or edges. The interaction of the boundaries between blocks in the DEM could reflect the discontinuity of rock masses and the characteristics of joints. Because the DEM is inherently capable of interpreting the geometric characteristics of jointed rock masses, it could well highlight the characteristics of structural planes in the rock masses; and it is convenient to deal with the failures of rock masses where all nonlinear deformations and failures are concentrated on the jointed surfaces. For the above reasons, DEM is widely used to numerically model mechanical processes such as landslides and rockfall in jointed rock masses, regardless of many blocks or only several blocks being moved.

The most common method used to investigate the deformations and failures of jointed rock slopes in numerical simulations is the DEM. The DEM could better simulate the failure process of discontinuous media, allowing discrete blocks to have limited displacement and rotation, and thus simulating the large displacement, rotation and other deformations of rock blocks. In the calculation process of the DEM, it could automatically identify new contact surfaces, which is suitable for simulating the jointed rock systems or discontinuous block assembly systems under static or dynamic loads. At the same time, the DEM could simulate the calculation after the partial instability failure of the slopes, and then predict how the slopes would be further damaged and the shape of the sliding surface [21].

Cundall [22] first proposed the calculation of the DEM, with later scholars applying UDEC and PFC${}^{2\mathrm{D}}$ to investigate the deformations and failures of two-dimensional slopes. For example, Wang et al. [23] analyzed the stability of a jointed rock slope with a weak interlayer by using UDEC for numerical analysis. The factor of safety (FOS) and failure mode before and after installing steel bolts in the slope are obtained by the strength reduction method. Conrad et al. [24] studied the basic concept of the limit equilibrium method for evaluating footwall slope stability. The numerical model established by UDEC illustrates the potential of the code to analyze these failure mechanisms. Salmi et al. [25] focused on the slope quality grade, used the back analysis method to analyze the slope stability of the karoun trench, and estimated the geo-mechanical parameters. The slope quality grade is used to evaluate the stability of the trench. Then, UDEC is used for numerical simulation to conduct back analysis, so as to predict the mechanical properties of the rock slope. Wolter et al. [26] proposed a comprehensive method to investigate the Madison Canyon landslide induced by earthquakes. Through the creation of engineering geomorphic maps, field investigations, remote ground digital photogrammetry, and preliminary 2D numerical modeling, the purpose is to determine the condition factors, mechanism, motion behavior, and fault evolution.

Numerical simulations have gradually developed from two-dimensional to three-dimensional with the development of computer technology. For example, Bonilla Sierra et al. [27] coupled the discrete fracture network with the DEM to provide a solution to evaluating the mechanical behavior of real three-dimensional structures where fracture persistence could not be assumed. A practical case is studied to show the complete method from photogrammetric data acquisition to numerical simulation of the potential progressive failure process of the rockmass. Wu et al. [28] used the discrete element software 3DEC to capture the basic three-dimensional characteristics of the landslide and form a post-failure structure similar to that observed in the field. Wang et al. [29] used 3DEC to numerically simulate the excavation of a rock slope of the Puli-Xuanwei Expressway in Yunnan Province. By comparing the rock slopes under different stratum angles, directions, rock characteristics, and thicknesses, they obtained the possible failure modes of the slope and the mutative regulation of safety factors under different conditions. Xu et al. [30] used 3DEC to analyze the stratum and surface movement and damage caused by open-pit slope mining.

In addition, Vyazmensky et al. [31] used an FEM/DEM modeling method to study the formation of the basement failure surface in an open-pit slope directly caused by cave mining, and emphasized the importance of rock mass tensile strength and its influence on slope response induced by collapse. Bonilla Sierra et al. [32] proposed the DFN-DEM model combined with photogrammetry technology to evaluate the stability of potentially unstable rock slopes. The model describes the yield and fracture of the intact rock matrix along pre-existing discontinuities. Scholtès et al. [33] proposed a three-dimensional numerical model based on the discrete element method to predict the deformation and failure of rock slope and avoid the problem of over simplification of the limit equilibrium method. Nadi et al. [34] considered the structural anisotropy effect caused by rock joints, blocks and pore water pressure, and studied the static stability of rock slope with an excavation depth of 15 m. Arnold et al. [35] used the particle discrete element program to analyze the dynamic failure of rock slope. Four main failure mechanisms are identified, which lead to three different failure modes, corresponding to the well-documented dynamic rock slope failure modes observed in natural slopes.

Moreover, Wu et al. [36] used 3DEC to study the rock collapse behavior of Xiandu Mountain caused by Typhoon Morakot and simulated the wedge failure of the landslide. Dong et al. [37] established a three-dimensional numerical model by using 3DEC software and combined it with the mechanical properties of the stratum, geological structure, rockmass, and discontinuity to analyze the stress, deformation, and stability of an excavated slope in China. The method used in this paper will be helpful in determining the true three-dimensional stress, deformation, and stability analysis of the rock slope. Regassa et al. [38] performed discrete element numerical simulation using an equivalent discontinuity model to simulate rock movement and failure caused by mining. Liu et al. [39] used the DEM to simulate the failure process of the Xiaogangjian landslide, emphasizing the influence of discrete fracture networks and rainwater infiltration. Weng et al. [40,41] applied the DEM to simulate the deformation behavior of slate slope under the influence of long-term gravity and material metamorphism, used the DEM to simulate the deformation and sliding process of highly active slope elements, and used the model to predict the slope behavior under the condition of wetting geotechnical materials.

However, when to properly terminate the discrete element calculation process remains a difficult challenge. Discrete element calculations are an iterative process. It is necessary to properly terminate according to the specific situation, and determine how many iterative steps to perform in the discrete element calculation process in advance. Currently, specifying the maximum number of steps is used as the most common method. Although this method is simple, it is not suitable in some simple and complex cases, because too few iterative calculation steps would bring about inaccurate numerical calculation results, and too many iterative calculation steps would entail unnecessary calculations.

To date, there has been much research work that attempts to solve the above problem associated with determining when the iterative process could be terminated. For example, Zhang et al. [42] combined the energy-based discrete element method with time-series analysis to determine when the discrete element iterative calculation process could be terminated. Shen et al. [43] analyzed the influence of the discontinuity of rock mass joints on slope stability, and conducted a sensitivity analysis on different input parameters. The global displacement of the rock slope was studied, and the global displacement vector diagram was drawn to judge when the slope was damaged. Although this method is reasonable, it is suitable for only some simple slopes. However, it is difficult to judge the position where the vector map has changed significantly for complex slopes. Marko et al. [44], Mikos et al. [45], and Lin et al. [46] determined when slope failure occurs through the value of the mean and dispersion of displacement.

Recently, Wang et al. [47] proposed a discrete element modeling method based on displacement statistics to investigate the stability of a jointed rock slope with a potential sliding surface. The procedure of the displacement-statistics-based method is as follows.

(1)

The potential sliding surface is roughly determined according to the geological data and conditions;

(2)

The vertices around the determined potential sliding surface are established before the discrete element iterative calculation;

(3)

The displacements of the vertex at the specified iteration interval are obtained; then the mean, standard deviation and the ratio of mean to standard deviation are calculated;

(4)

The convergence of the ratio of mean to standard deviation is analyzed.

However, in the above displacement-based method, it is necessary to determine the position of the potential sliding surface before the discrete element iterative calculation. In fact, determining the position of the potential sliding surface is difficult for some slopes in many cases; the position of the sliding surface would directly affect the final calculation results. In addition, the displacement-based method only considers the displacement of rock blocks around the potential sliding surface, but does not consider the mass or weight of rock blocks around the potential sliding surface.

To solve the above problem, a discrete element method based on energy correlation coefficient is proposed to analyze the deformation and failure of jointed rock slopes. The proposed method is a generalization of the displacement-based discrete element method. This proposed method considers not only the mass or weight of the block, but also the displacement of the block. Compared with the energy-based method [42], the proposed method in this paper considers both the temporal and spatial correlation between all blocks, which is more comprehensive. Meanwhile, the calculation cost and efficiency of the proposed method are basically the same as those of the displacement-based method and the energy-based method. Our contributions can be summarized as follows.

(1)

An energy-based criterion is proposed to terminate the iterative process of discrete element modeling.

(2)

A discrete element modeling method based on the energy correlation coefficient is proposed to analyze the failures and deformations of jointed rock slopes.

(3)

The rationality of the proposed method is verified by a simplified jointed rock slope, and the proposed method is applied to a practical case.

The rest of this paper is organized as follows. Section 2 introduces the proposed method in this paper in detail. Section 3 gives a simplified rock slope model to verify the effectiveness of the proposed method. In Section 4, the proposed method is applied to the deformations and failures analysis of a real rock slope case. Section 5 discusses the novelty and shortcomings of the proposed method. Meanwhile, future work is given in this section. Finally, Section 6 draws several conclusions.

2. Methodology

2.1. Overview

The DEM is essentially suitable for analyzing the failures and deformations of jointed rock slopes, but the discrete element calculation is an iterative process that needs to choose an appropriate time to terminate according to the specified situation; the DEM requires that the number of iterative steps to be performed in the calculation process needs to be determined before the discrete element calculation. Currently, the specified maximum number of steps is the most commonly used method. Although this method is simple, it is not suitable in some simple and complex cases, because too few iterative calculation steps would bring about inaccurate numerical calculation results, and too many iterative calculation steps would entail unnecessary calculations.

In this paper, an enhanced discrete element modeling method is proposed to investigate the deformations and failures of jointed rock slopes through the statistic of the energy correlation coefficient. In this method, the convergence of the energy correlation coefficient is used to determine when the iterative process could be terminated. The proposed method is a generalization of the displacement-based method. The proposed method considers the mass or weight and displacement of the block. Meanwhile, the proposed method in this paper considers the temporal and spatial correlations between all rock blocks.

In short, the basic idea of the proposed discrete element modeling method is to provide a method based on energy analysis to determine when to terminate the DEM iterative process to investigate the deformations and failures of jointed rock slopes. As shown in Figure 1, the proposed discrete element modeling method is roughly composed of the following three steps:

(1)

According to the specific situation of the study area, the geological model and calculation model are established, and the appropriate mechanical parameters are determined;

(2)

The deformation of the jointed rock slope under gravity is calculated by the DEM, and the location and volume of all rock blocks are obtained;

(3)

Correlation analysis is performed to properly terminate the iterative process of discrete element calculation. Specifically, the energy of each rock block under gravity is calculated according to the position and volume of all rock blocks obtained by the second step. Then, the correlation coefficient between the energy at each time and the energy at the initial time is calculated. Finally, the convergence of the energy correlation coefficient sequence is analyzed. If the sequence converges, the discrete element calculation process would be terminated.

2.2. Establishment of the Computational Model

Because MIDAS/GTS software could establish large-scale and complex geological models and could automatically generate high-quality mesh models, the MIDAS/GTS software is used to establish the surface and mesh according to the topographic contours of the study area. In addition, the discrete element modeling software 3DEC is suitable for numerical simulation of jointed rock slopes, but it could not establish the geological model well. Therefore, the combination of MIDAS/GTS and 3DEC is used to establish the final calculation model. The model generated in MIDAS/GTS is imported into 3DEC through a specific interface, by adding multiple groups of key joints in 3DEC, dividing the model into multiple regions, and generating the final calculation model.

After establishing the calculation model, the boundary conditions need to be set. The boundary conditions in the numerical analysis are mainly determined according to the displacement. The boundary condition setting in this paper is that, for the calculation model, there are usually six edge interfaces, including the bottom, top, and surrounding four edges. The displacement of all nodes on the bottom surface is strictly fixed. Instead, all nodes on the top face are free. For the four surrounding faces, the displacement in the normal direction of the face is fixed. For example, the displacement in the X-direction is fixed, and the displacement in the Y-direction and Z-direction is not fixed.

2.3. DEM Modeling

At present, the DEM is widely used to investigate the deformation, failure, and stability of jointed rock slopes, but the discrete element calculation is an iterative process, and when the iterative process could be terminated is uncertain. Currently, specifying the maximum number of iterations is artificially used as the most common method before discrete element calculation.

According to the calculation model established in Section 2.2, different physical and mechanical parameters and joint mechanical parameters are assigned in different regions. The deformations and failures of jointed rock slopes under gravity are analyzed. For a fixed number of iteration steps per interval (for example, every 1000 steps), the position and volume of all rock blocks are obtained in 3DEC software until the specified number of maximum iteration steps or dynamically determined iteration steps is reached.

2.4. Statistical Analysis

In statistics, the most commonly used analysis method is time series analysis. Time series, also known as dynamic series, refers to a numerical sequence that arranges the same statistical indicator values in chronological order. Time series analysis could be approximately divided into three parts: describing the past, analyzing the regular pattern, and predicting the future. The purpose of time series analysis is to predict the changing of this statistical indicator in the future based on the existing real data. Time series analysis is a quantitative forecasting method and is commonly used in forecasting methods in statistics. Currently, time series analysis has been widely used in regional comprehensive development planning, business management, weather forecasting, geological disaster forecasting, crop disease, insect disaster forecasting, environmental pollution control, ecological balance, etc.

In this paper, the method of time series analysis will be adopted to determine when the iterative process could be terminated in the discrete element calculation process. The state of the rock block at different times of the slope is shown in Figure 2. Figure 2 shows that from the initial stable state $T_0$ to the final failure of the slope $T_n$, the slope has experienced $T_0$, $T_1$, $T_2$, …, $T_{n-1}$, $T_n$. Every moment corresponds to a different state of the slope at distinct moments so that the energy of the slope under the action of gravity at each moment could be calculated, and then a complete-time sequence could be formed.

2.4.1. Calculation of the Energy under Gravity

The first step of the statistical analysis is to calculate the energy of all rock blocks of the jointed rock slope under the action of gravity from the initial time $T_0$ to $T_1$, $T_2$, …, $T_{n-1}$, $T_n$.

Based on the positions and volumes of the rock blocks obtained in Section 2.3, the energy of all rock blocks under the action of gravity could be calculated from the initial state to any time thereafter. For a three-dimensional jointed rock slope, when the slope is under the action of gravity for an extended time, some rock blocks on the slope would slide. The change in the position of a rock block on the rock slope from one moment to another is shown in Figure 3. The total displacement of the rock block from $T_0$ to $T_n$ is $\Delta d$. The total displacement $\Delta d$ could be decomposed into $\Delta x$ in the X-direction, $\Delta y$ in the Y-direction, and $\Delta z$ in the Z-direction, so $\Delta d$ could be expressed as Equation (1).

$$\Delta d=\sqrt{{\Delta x}^2+{\Delta y}^2+{\Delta z}^2}$$

(1)

• $\Delta d$—The total displacement (m);
• $\Delta x$—Displacement in the X-direction (m);
• $\Delta y$—Displacement in the Y-direction (m);
• $\Delta z$—Displacement in the Z-direction (m).

The energy could be expressed as the force multiplied by the displacement in the direction of the force. Therefore, the energy under gravity is equal to gravity multiplied by $\Delta z$. $\Delta z$ is equal to the z coordinate of the center position of the rock block at time $T_n$ minus the z coordinate of the center position of the rock block at time $T_0$. The coordinates of the center positions of all rock blocks are obtained in Section 2.3. In addition, the density is a physical and mechanical parameter determined before the iterative calculation of the DEM. Therefore, the energy of all rock blocks under the action of gravity could be calculated according to Equation (2).

$$Energy={\rho }_n{v}_{n}g\Delta z_n$$

(2)

• Energy—Energy under gravity (N·m);
• ${\rho }_n$—Density (kg/m${}^3$);
• $v_n$—Volume (m${}^3$);
• g—Acceleration of gravity (N/kg);
• $\Delta z_n$—Displacement in the Z-direction (m).

The proposed method based on energy correlation coefficient in this paper is an extension of the displacement-based method; this is because the proposed method considers the mass or weight and displacement of each block. The proposed method is suitable for slopes without obvious potential sliding surfaces. The displacement-based method considers the mass or weight of each block to be one, but the weight of each block obtained by the proposed method is the actual value calculated by the volume and density. Meanwhile, the proposed method considers both the temporal correlation and spatial correlation. Therefore, the proposed method based on energy correlation coefficient is more accurate, more comprehensive, and has wider applicability.

2.4.2. Calculation of Correlation Coefficient Based on Energy

The second step of the statistical analysis is to calculate the Spearman correlation coefficient of energy, that is, to calculate the correlation coefficient between the energy under gravity of all rock blocks at any time and the energy under gravity of all rock blocks at $T_1$.

The Spearman correlation coefficient, also known as the Spearman rank correlation coefficient, is a kind of rank correlation coefficient. The Spearman correlation coefficient is a measure of the synergy relationship. It is different from the Kendall correlation coefficient, which describes the synergy relationship between two independent and identically distributed binary distributions. The Spearman correlation coefficient describes one of the three independent and identically distributed binary distributions. In this paper, the Spearman correlation coefficient is used as a statistical indicator due to its wide applicability and monotonic transformation invariance, and it is not restricted by dimensions, and is not affected by outliers. At the same time, it could also describe a certain non-linear relationship. Therefore, the Spearman correlation coefficient is chosen as a statistical indicator. The calculation equation used is as follows.

$$CC=1-\frac{6\sum d_i^2}{n\left(n^2-1\right)}$$

(3)

The calculation process is as follows. By first sorting the data of the two variables $\left(X,Y\right)$, and then writing down the position $\left(X_1,Y_1\right)$ after sorting, the value of $\left(X_1,Y_1\right)$ is called the rank, the difference of the rank is the $d_i$ in the above formula, n is the number of the variable data, and finally the formula could be used to solve the result. The following examples illustrate the meaning of each letter as shown in

Table 1. A simple example of explaining the meaning of $d_i$.

Original Location	Original X	After Sorting	Rank ${\mathit{X}}_1$	Original Y	After Sorting	Rank ${\mathit{Y}}_1$	Square of Rank Difference ${\mathit{d}}_{\mathit{i}}$
1	11	490	5	2	75	6	1
2	490	43	1	75	44	1	0
3	14	30	4	3	42	5	1
4	43	14	2	44	7	2	0
5	30	11	3	7	3	4	1
6	3	3	6	42	2	3	9

.

A relationship curve is drawn between each CC value and the corresponding number of iteration steps. When the number of iterations reaches a certain value, the relationship curve between CC and the number of iterations begins to converge, and the iteration process could be terminated.

The proposed method based on the energy correlation coefficient in this paper considers the mass or weight and displacement of each block. The proposed method considers both the temporal correlation and spatial correlation. Therefore, compared with the displacement-based method, the proposed method could be effectively applied to investigate the deformations and failures of jointed rock slopes in which the potential sliding surface could not be determined easily.

2.4.3. Determination of the Number of Iteration Steps

The third step of statistical analysis is to determine the number of steps to terminate the iterative calculation of the DEM.

According to Section 2.4.1 and Section 2.4.2, the correlation coefficient corresponding to the number of steps in each iteration are obtained. The convergence of the CC curve is analyzed, in which the abscissa is the number of iteration steps and the ordinate is CC, as shown in Figure 4a. Figure 4a shows that, when the iterative process reaches the 8000th step, the curve tends to converge. To make the result more accurate, we calculate the difference corresponding to $T_1-T_0$, $T_2-T_1$, $T_3-T_2$, …, $T_n-T_{n-1}$, and then draw a curve, as shown in Figure 4b. The difference between the 8000th step and the 7000th step is near zero. As the number of iterations increases, the difference is basically equal to zero. At this time, CC remains unchanged. The iteration process could be terminated when the iteration reaches the 7000th step.

The above method for determining the number of iterations is based on the energy correlation coefficient between all rock blocks of the jointed rock slope, and the energy calculation takes the displacement and weight of all blocks into account. Meanwhile, the proposed method considers the temporal and spatial correlation. The proposed method makes the calculation results more accurate. The proposed method is suitable for slopes without obvious potential sliding surfaces. However, the displacement-based method is a special case of the method proposed in this paper. It is only suitable for jointed rock slopes with obvious potential sliding surfaces.

3. Verifications

In this section, a simplified jointed rock slope will be used as an example to verify the deformation and failure of the simple jointed rock slope, and compare it with the discrete element method based on displacement statistics to verify the effectiveness of the proposed method.

The displacement-based discrete element modeling method could be applied to analyze the deformation and failure of jointed rock slopes with a weak interlayer. However, the displacement-based method needs to determine the potential sliding surface before the iterative calculation of the DEM, and the displacement-based method needs to obtain the displacement of rock blocks around the sliding surface, and calculate the average value, standard deviation, and variation coefficient (Mean/SD) of the displacement. Then, according to the convergence of Mean/SD, the number of iterations to terminate the iterative process is obtained.

In this section, a model of a simple jointed rock slope is established, and the displacement-based method and the proposed method are used for calculation, verification and comparison to verify that the method based on the energy correlation coefficient is correct and feasible.

3.1. Computational Model and Parameters

In this section, an abbreviated jointed rock slope is established with a weak interlayer, and the inclination of the weak interlayer is 20${}^{°}$. More details are shown in Figure 5.

MIDAS/GTS is used to build the model, divide the slope model into the tetrahedral mesh, and then the tetrahedral meshes are imported into the discrete element modeling method software 3DEC as a numerical calculation model through a special interface. The quantity of blocks and vertices in the simplified rock slope model are 2491 and 11,455, respectively. This paper will use the Mohr–Coulomb model as the constitutive model, because the Mohr–Coulomb model is the most commonly used model in numerical research to analyze the deformations and stability of slopes [48,49], and it could comprehensively reflect the strength characteristics of the rock. The physical and mechanical parameters of the model are shown in

Table 2. Mechanical parameters of rockmass and weak interlayer.

Type Unit	$\mathit{\rho }$ (kg/m${}^3$)	BULK (GPa)	G (GPa)	C (MPa)	$\mathit{\phi }$ (${}^{°}$)	JKN (GPa/m)	JKS (GPa/m)
Rockmass	2300	3.33	2.0	0.13	32	1.5	0.35
Weak interlayer	1800	0.33	0.2	0.01	19	1.2	0.3

according to the handbook of rock mechanics parameters.

In Table 2, BULK stands for the bulk modulus, G stands for the shear modulus, $\rho$ stands for the density, C stands for the cohesion, $\phi$ stands for the internal friction angle, JKN stands for the normal stiffness of the joint, and JKS stands for the shear stiffness of the joint.

3.2. Computational Results

The simplified rock slope is simulated numerically in Section 3.1 in the discrete element software 3DEC, and the physical and mechanical parameters are assigned in Table 2. The slope is initially unstable, and the initial stress state is mainly caused by geological structure and gravity. In addition, topography, rock mechanical properties, water, and temperature will also affect the initial stress state. As the number of iteration steps gradually increased, the rock blocks on the jointed rock slope began to slide, and some rock blocks began to slide out of the boundary of the model. When the number of iteration steps reaches 14,000, a large number of rock blocks have slipped and accumulated, and the rock blocks slide along the weak interlayer, as shown in Figure 6. Figure 6a shows the overall sliding situation of the simplified jointed rock slope when the number of iteration steps reaches 14,000. With the gradual increase in the number of iterations, the rock blocks on the simplified jointed rock slope will continue to slide, and result in greater displacement. Figure 6b is a diagram of a certain surface when the simplified slope is iteratively calculated to the 14,000th step. Figure 6b shows that the rock blocks slide along the joint surface and the weak interlayer.

The discrete element software 3DEC is used to obtain the positions and volumes of all rock blocks every 1000 steps until 50,000 steps. Then, data cleaning is performed on the acquired data files; the energy under gravity is calculated from the initial state to any state. Then, the Spearman correlation coefficient based on energy is calculated, and the relationship curve between the number of iteration steps and the Spearman correlation coefficient is drawn. To verify the correctness of the proposed method, the result of the proposed method is compared with the displacement-based method. The results obtained using two different methods are shown in Figure 7. We can see that the suggested number of iteration steps in the displacement-based method is 12,000, while the suggested number of iteration steps in the proposed method is 14,000. The inflection points of the two methods are almost the same, and the results are basically consistent. In addition, from Figure 7, it is thought that when the calculation reaches approximately the 14,000th step, the iterative calculation process could be terminated.

4. Application: A Real Case

4.1. Geological and Engineering Background

The Songmugou landslide is located in the middle of the mining area of the Jingu Coal Industry in Guxian County, Shanxi. It belongs to Songmugou Formation, Guyang Village, Guyang Town, Gu County, with coordinates: N: ${112}^{°}{01}^{\prime }{04}^{″}$, E: ${36}^{°}{23}^{\prime }{24}^{″}$. The main terrain is a middle–low mountainous area, with an elevation of 1050 m and a valley cutting depth of 100 m–200 m, and the south side of the point is mostly high mountains and steep slopes with great terrain changes, with an elevation of 1200 m and above, and the slope range of the mountain slope is $50--{60}^{°}$, while the terrain on the north side is relatively flat, with an elevation of 1050 m and below, most of which are the residences of villagers, transportation and administrative centers. The geographic location of the Songmugou landslide is shown in Figure 8, and the geological section of the landslide is shown in Figure 9.

The Songmugou landslide is located on the south side of Songmugou valley. The bottom of the valley has an elevation of 1050 m and a widened “U” shape. The top of the slope on the south side has an elevation of 1200 m, and the original slope is ${55}^{°}$. The slope on the north side is relatively flat, with an elevation distribution of 1050 m–1100 m, of which there are four steep ridges distributed in steps, with a drop of 0.8 m–1.0 m, thus forming multiple loess platforms. The landslide is the unstable deformation of the southern slope, and the sliding body and the front edge of the landslide slide into the valley.

4.2. Computational Model and Mechanical Parameters in the DEM Calculation

This section introduces the establishment of the geological model and calculation model in the discrete element modeling, and determines reasonable physical and mechanical and joint parameters. The calculation model and calculation parameters will be used to analyze the deformation and failure of the Songmugou landslide.

When using 3DEC software to investigate the movement and destruction of a jointed rock slope, one of the most significant steps is to establish a simple but functional calculation model. According to the geological report of the Songmugou landslide and related description of the landslide, reasonable simplification of specific geological conditions is conducted. The rock formations are divided into Lower Shihezi Formation (P1x), Lower Shanxi Formation (P1s), Upper Taiyuan Formation (C3t), bed rock and coal seams. Below the 9# and 10# coal seams excavated is the bedrock. According to the contour lines extracted from the topographic map of the study area, the initial geological model is established in the MIDAS/GTS software, and imported into the 3DEC software, adding simplified joints. Finally, the required calculation model is obtained.

More specifically, the establishment of a calculation model includes two procedures: (1) establishing a geological model without joints; (2) adding joints based on the geological model to establish a final calculation model.

The first procedure is to establish a preliminary geological model, which is completed by using MIDAS/GTS software according to the overview of the study area. The following three main steps are included in the first procedure:

(1)

Extracting the contour data of the ground surface according to the topographic and geological map of the location of the Songmugou landslide, selecting the contour data in a reasonable range where the Songmugou landslide is located, and generating the surface in the MIDAS/GTS software.

(2)

Dividing the surface into triangular meshes in the MIDAS/GTS software.

(3)

Calculating the crossing points of the triangular meshes, and connecting two corresponding triangles in any pair of adjacent curved surface meshes vertically to create an initial geological model that is composed of prisms.

The second procedure is to establish the final calculation model by adding reasonable joints. This procedure was completed using 3DEC software. The following main steps are included in the second procedure:

(1)

A specially developed interface program is used to import the preliminary geological model which is composed of prisms from the MIDAS/GTS software to the 3DEC software.

(2)

The excavated coal seam is cut by multiple planes, and the geological model is divided into multiple different subareas. Then, the required calculation model is established by adding key joints, as shown in Figure 10.

According to the Songmugou landslide geological survey report, the related description of the landslide and handbook of rock mechanics parameters, the physical and mechanical parameters of the rock formation and joint mechanical parameters are shown in

Table 3. The physical and mechanical parameters and the parameters of the joint [47].

Type Unit	$\mathit{\rho }$ (kg/m${}^3$)	BULK (GPa)	G (GPa)	C (MPa)	$\mathit{\phi }$ (${}^{°}$)	JKN (GPa/m)	JKS (GPa/m)
P1x	2560	2.22	1.02	0.73	31	0.20	0.20
P1s	2460	2.86	1.4	2.8	39	0.58	0.25
C3t	2437	4.3	2.8	0.7	30	0.60	0.28
Bed rock	2800	5.57	4.53	11.4	38	0.84	0.36
Coal	1420	0.46	0.19	0.8	20	0.13	0.13

.

Currently, the Mohr–Coulomb model is the most commonly used model for studying the stability of rock slopes in numerical simulations [48,49]. This paper also uses the Mohr–Coulomb model as the constitutive model.

4.3. Computational Results

When using the DEM to simulate the Songmugou landslide, since the distance between the 9# and 10# coal seams is relatively short and the thickness of the coal seams is relatively thin, we will excavate the 9# coal seam, 10# coal seam and the interlayer at one time. The total thickness is 20m. Then, the energy under gravity and the correlation coefficient are calculated corresponding to each iteration step, so as to draw a curve about the relationship between the correlation coefficient and the number of iterations; finally, the number of calculation steps is determined for termination by analyzing the convergence of the curve.

After the simulated excavation of the 9# and 10# coal seams is completed, the Z-displacement diagram of the Songmugou landslide is shown in Figure 11. Figure 11a shows the overall Z-displacement diagram after the excavation is completed, and Figure 11b shows the Z-direction displacement of the cross-section after the excavation is completed. From Figure 11, we find that due to the excavation of the coal seams, the rock mass on the upper part of the coal seam moves downward and slides down along the joint. From Figure 11b, we could see that the split level is approximately 15 m. In addition, the landslide also has a large displacement along the Y direction. As shown in Figure 12, we can see that the sliding displacement is large on the left boundary of the model and the top of the landslide. We compared the numerical simulation results with the on-site survey results and found that they are basically consistent with the actual landslide result, and they all showed basically the same sliding surface, as shown in Figure 13.

In Figure 13, the Songmugou landslide is in the shape of a dustpan and could be roughly divided into three deformation zones.

Zone I: Deformation stability zone. This part of the area does not slide with the main sliding body, but is caused by deformation after the main sliding body slides, and is adjacent to the east pine wood ditch. It is mainly silt and crushed stone that collapses and slides in the direction of the empty surface.

Zone II: Intensive deformation zone. This part of the zone belongs to the most obvious deformation zone affected by landslide action. The upper part is a platform with a large slope. The trees on it were originally the top of the slope and slide together with the landslide. However, due to the sliding force, these trees have been bent and deformed, and due to the difference in the sliding speed of the upper sliding body and the lower sliding body, the lower sliding body is squeezed to produce a bulge and multiple fan-shaped cracks.

Zone III: Unstable deformation development zone. Because the landslide slides slowly, the landslide does not quickly fall to the bottom of the valley, but slowly accumulates along the slope. If coal mining activities proceed further, torrential rain and other factors induce this part. The area may undergo large-scale deformation, and the deposits on the slope will slide into the valley. If the rainfall reaches more than 300 mm/h, it may induce a debris flow.

Due to the mining of the coal seams, the rock mass above the coal seam produces different degrees of slipping. The main slipping is along the Y direction and the Z direction. Compared with the displacement in the Y direction, the displacement in the Z direction is larger and the sliding range is wider. The method proposed in this paper will be used to perform a statistical analysis of the landslide, and the statistical results are shown in Figure 14. Figure 14 shows that when the number of iterative steps reaches the 7000th step, the curves in Figure 14a,b tend to converge. The CC value of the 8000th step is approximately equal to the CC value of the 7000th step. Therefore, when the calculation reaches the 8000th step, the iterative calculation will be terminated.

The above excavation process regards the 9# coal seam, 10# coal seam, and the interlayer. The results of the excavation could verify the feasibility and applicability of the proposed method. However, with increasing excavation depth, a large number of rock blocks accumulate on the surface. In brief, according to the results of the above excavation, it is thought that the proposed method is applicable to more jointed rock slopes that could not determine the position of the potential sliding surface before the iterative calculation of discrete element modeling. In short, the proposed method considers both temporal correlation and spatial correlation between all blocks.

5. Discussion

5.1. The Novelty of the Proposed Method

The novelty of the proposed discrete element modeling method based on the energy correlation coefficient is that it is unnecessary to determine the location of potential sliding surfaces before the iterative calculation of the DEM. The proposed method based on the energy correlation coefficient is a generalization of the method based on the displacement statistics. The proposed method considers the mass or weight and displacement of every block.

In addition, it is no longer to count the displacement of all blocks as a statistic, but to consider the correlation between all blocks. The method considers both time correlation and spatial correlation, and the consideration is more comprehensive. Since the method based on displacement statistics needs to determine the potential sliding surface before the discrete element calculation, it is only suitable for jointed rock slopes that the position of the potential sliding surface is liable to determination. Compared with the displacement-based method, our proposed method is not essential to determine the location of the potential sliding surface in advance. Therefore, the method proposed in this paper has a wider application range and could be applied to general jointed rock slopes that the position of the potential sliding surface is hard to determine.

Specifically, the displacement-variation-coefficient-based method is only suitable for the analysis of jointed rock slopes where the potential sliding surface is easy to determine before the DEM iterative calculation. The method takes the displacements of rock blocks above the potential sliding surface as statistics, not the displacements of all rock blocks. Therefore, the calculation is faster.

The energy-variation-coefficient-based method is a generalization of the above method, and it is no longer necessary to determine the position of the potential sliding surface in advance. The method is suitable for complex jointed rock slopes where the potential sliding surface is difficult to determine. In this method, the energy generated by all rock blocks in the process of movement are taken as statistics. Not only the displacement of rock blocks, but also the mass of rock blocks are considered, so as to eliminate the influence of block size on the calculation results. This method considers the time correlation between rock blocks, but does not consider the spatial correlation.

The comparison of the mentioned methods is listed in

Table 4. The comparison of the mentioned methods.

Method	Comparison	Efficiency
The displacement-variation-coefficient-based method	Needs to determine the location of potential sliding surfaces before the DEM modeling.	Higher
The energy-variation-coefficient-based method	No need to determine the location of potential sliding surfaces before the DEM modeling; Only time correlation is considered.	Lower
The proposed energy-correlation-coefficient-based method	No need to determine the location of potential sliding surfaces before the DEM modeling; Time correlation is considered; Spatial correlation is considered.	Lower

.

The proposed method could be applied in different case studies, such as complex rock slopes, underground mining and so on. In addition, the proposed method is built on the block discrete element method, if similar numerical models such as DDA [16,17] and NMM [18,19] are employed, the proposed method could also be used to analyze the deformations and failures of jointed rock slopes.

5.2. Shortcomings of the Proposed Method

The method proposed in this paper could be adopted to investigate the failures and deformations of simple jointed rock slopes, and could also be used to analyze the failures and deformations of complex jointed rock slopes, regardless of the size or shape of the slopes. However, the method based on the energy correlation coefficient needs to take the mass or weight and displacement of all rock blocks into account, not the mass or weight and displacement of some rock blocks. In this case, it is perhaps computationally a little more expensive. The proposed method in this paper is suitable for jointed rock slopes, but is not suitable for soil slopes. When the deformations of the rock slopes are very small, or when the rock blocks on the rock slopes collapse completely, the applicability of the proposed method in the paper needs to be verified by real-world engineering cases.

5.3. Outlook and Future Work

In the future, the Python scripting language would be used to embed into the discrete element calculation software, so as to realize real-time calculation and automatically terminate the iterative calculation process.

6. Conclusions

In this paper, a new discrete element modeling method based on the energy correlation coefficient is proposed to analyze the failures and deformations of simple and complex jointed rock slopes. The main idea of this method is to determine when the discrete element iteration process could be terminated in investigating the failures and deformations of various jointed rock slopes based on the statistic of the energy correlation coefficient. The proposed method has wider applicability because it is unnecessary to determine the location of the potential sliding surface in advance. The proposed method based on the energy correlation coefficient is a generalization of the method based on displacement statistics; the proposed method considers the mass or weight and displacement of every block. The displacement of all blocks is counted as a statistic; the proposed method considers spatial and temporal correlations, which is more comprehensive. The effectiveness of the proposed method is first verified using an abbreviated jointed rock slope, and then the proposed method is applied to a practical case, which demonstrates that the proposed method is capable of investigating the deformation and failure of jointed rock slopes in which the position of the potential sliding surface is difficult to determine. In the future, we will use the Python scripting language to embed it into the discrete element calculation software to achieve real-time calculation and automatically terminate the iteration.