An Energy-Based Discrete Element Modeling Method Coupled with Time-Series Analysis for Investigating Deformations and Failures of Jointed Rock Slopes

: The discrete element method (DEM) can be effectively used in investigations of the deformations and failures of jointed rock slopes. However, when to appropriately terminate the DEM iterative process is not clear. Recently, a displacement-based discrete element modeling method for jointed rock slopes was proposed to determine when the DEM iterative process is terminated, and it considers displacements that come from rock blocks located near the potential sliding surface that needs to be determined before the DEM modeling. In this paper, an energy-based discrete element modeling method combined with time-series analysis is proposed to investigate the deformations and failures of jointed rock slopes. The proposed method deﬁnes an energy-based criterion to determine when to terminate the DEM iterative process in analyzing the deformations and failures of jointed rock slopes. The novelty of the proposed energy-based method is that, it is more applicable than the displacement-based method because it does not need to determine the position of the potential sliding surface before DEM modeling. The proposed energy-based method is a generalized form of the displacement-based discrete element modeling method, and the proposed method considers not only the displacement of each block but also the weight of each block. Moreover, the computational cost of the proposed method is approximately the same as that of the displacement-based discrete element modeling method. To validate that the proposed energy-based method is effective, the proposed method is used to analyze a simple jointed rock slope; the result is compared to that achieved by using the displacement-based method, and the comparative results are basically consistent. The proposed energy-based method can be commonly used to analyze the deformations and failures of general rock slopes where it is difﬁcult to determine the obvious potential sliding surface.


Introduction
The deformations and failures analysis of slopes is a representative problem in geotechnical engineering [1], which have been widely performed in both scientific and practical research. At present, many numerical simulation methods have been used to analyze the deformations and failures of rock slopes, such as the discrete element method (DEM) [2][3][4], finite element method (FEM) [5][6][7], finite difference method (FDM) [8], boundary element method (BEM) [9], material point method (MPM) [10][11][12], general particle dynamics (GPD) [13,14], and numerical manifold method (NMM) [15]. These numerical simulation methods are developed from the primitive limit equilibrium method (LEM). In addition, the smoothed particle hydrodynamics (SPH) method [16] and 2-D/3-D finite element code [17] also have been used to evaluate the stability of slopes. The Mohr-Coulomb failure criteria are the most commonly employed in numerical modeling to investigate the stability of slopes in the most research work [18][19][20][21].
When analyzing the deformations and failures of jointed rock slopes, the most commonly used method is the DEM. Because the DEM is a numerical method based on discontinuities, it can be better applied to discontinuous deformations and failures and large displacements. For example, Wang et al. [22] proposed the displacement-based discrete element modeling method (DSDM) to analyze the stability of jointed rock slopes with a potential sliding surface. Bonilla-Sierra et al. [23] used the DEM to simulate jointed rock masses and carried out numerical research on the progressive failure of rock slopes with non-permanent joints. Shi et al. [24] evaluated the factors affecting the deformations and failures of a dangerous rock mass and the geological conditions of a highway slope in Yunnan Province using the DEM. In addition, a three-dimensional numerical analysis of the strength reduction of safety factors was carried out. Wu et al. [25] used the discrete element simulation 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.
In addition, Xu et al. [26] used the DEM to analyze the strata and surface movement and damage caused by open-pit slope mining. Regassa et al. [27] carried out DEM numerical simulation by using an equivalent discontinuous model (EDM) to simulate rock movement and failure caused by mining. Liu et al. [28] used the DEM to simulate the failure process of the Xiaogangjian landslide and especially considered the influence of the discrete fracture network (DFN) and rainwater infiltration. Weng et al. [29] used the DEM to simulate the deformation behavior of slate slopes under the long-term influence of gravity and material metamorphism. Weng et al. [30] 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 of geotechnical materials. Wu et al. [31] used the discrete element simulation software 3DEC to study the behavior of Xiandu Mountain rock collapse caused by Typhoon Morakot in 2009. The analysis successfully simulated the wedge failure of the landslide. In addition, the 3DEC results also illustrate the local rock sliding phenomenon at the boundary of the source area.
However, determining when the DEM iterative process can be appropriately terminated is still a difficult problem. The DEM calculations involve an iterative process that needs to terminate appropriately according to specific conditions. The number of iterative steps should be determined appropriately before the DEM modeling. Specifying the maximum number of iterative steps is the most common solution. When the specified maximum number of iterative steps is reached, the iterative process would be terminated. The common solution is very simple. However, it is obviously inappropriate in some cases, because too many iterations would be superfluous, while too few iterations result in inaccurate calculation results.
Much research work has been conducted to address the above problems. For example, the overall displacements were investigated, and vector graphics of the overall displacements were plotted. Shen et al. [32], Wang et al. [33], and Salmi et al. [34] suggested that obvious changes in vector graphics could indicate that a slope has been deformed. Although the method is rational, it is only applicable to simple slopes. In other cases, vector graphics would be quite confusing, and the obvious changes in displacements would be difficult to estimate.
Displacement characteristics are commonly used to reflect the deformations and failures of slopes. For example, in collapsed landslides [35,36], the slope is stable at the beginning, compared with the size of the slope, the displacements of rock blocks are quite small. When the slope begins to be unstable, the average displacements around the sliding surface would increase by degrees, and the dispersion degree of the displacements of the rock blocks around the sliding surface would also increase by degrees. In traction landslides [37,38], with the development of landslides, the mean of slope displacements would gradually increase. However, the rock blocks located around the sliding surface have broadly consistent mean of displacements. This means that the dispersion degree of the displacements gradually decreases.
Recently, Wang et al. [22] used the displacement-based discrete element modeling method (DSDM) to investigate the stability of jointed rock slopes. The procedure of the displacement-based method is as follows: (1) According to the geological conditions and the geographical location, the potential sliding surface is determined approximately; (2) In the established calculation model, the vertices around the determined potential sliding surface are determined before the iterative DEM modeling; (3) The iteration interval specified during the DEM calculation (e.g., at the 500th step, the 1000th step, etc.), and the mean, the standard deviation (SD), and mean/SD of the displacements of marked vertices are calculated; (4) The authors analyzed the convergence of the mean/SD values. However, in the aforementioned displacement-based method, there is a problem that it needs to determine an obvious potential sliding surface before DEM modeling. In fact, in many cases, it is difficult to effectively determine an obvious potential sliding surface. The position of the sliding surface would directly affect the final results. Moreover, the displacement-based method only considers the displacements of rock blocks around the potential sliding surface, but the weights of blocks around the potential sliding surface are not considered in the displacement-based method.
To address the above problems, in this paper, an energy-based discrete element modeling method is proposed for analyzing the deformations and failures of jointed rock slopes, which is combined with time-series analysis. The proposed method is a generalized form of the displacement-based discrete element modeling method. The proposed method considers not only the displacement of each block, but also the weight of each block. Moreover, the computational cost of the proposed method is approximately the same as that of the displacement-based discrete element modeling method.
The objective of proposing the energy-based method is to generalize the displacementbased method. As mentioned several times, the most obvious shortcoming of the displacementbased method is that, it needs to previously determine the potential sliding surface before conducting the numerical modeling. In contrast, in the proposed energy-based method, there is no need to determine the potential sliding surface before DEM modeling. In other words, the displacement-based method is only applicable to the cases where the sliding surface can be easily determined. It is not suitable to be applied in the cases where the sliding surface is difficult to be determined. In contrast, the proposed energy-based method is much more applicable than the displacement-based method. This is because it can be applied in the cases where the sliding surface is easy or difficult to be determined. For the case where the sliding surface can be easily determined, the numbers of iterative steps suggested by the two methods are consistent. In short, the proposed energy-based method is a generalized form of the displacement-based method, and it is much more applicable than the displacement-based method.
Contributions of this paper can be summarized as follows: (1) An energy-based criterion is proposed to terminate the DEM iterative process; (2) A new energy-based discrete element modeling method is proposed for analyzing the deformations and failures of jointed rock slopes combined with time-series analysis; (3) The proposed energy-based method is verified through a simple jointed rock slope, and the proposed method is applied in a real case.
The rest of this paper is organized as follows. Section 2 introduces the proposed method in this paper in detail. In Section 3, a simplified jointed rock slope is used to verify that the proposed energy-based method is effective. The proposed method is utilized to investigate the deformation and failure of an authentic slope case in Section 4. In Section 5, the novelty, shortcomings, and the future work of the proposed energy-based method are discussed. Finally, several conclusions are drawn in Section 6.

Overview
The DEM is suitable for the deformations and failures analysis of jointed rock slopes inherently. The DEM calculations involve an iterative process that needs to be terminated appropriately according to specific conditions. However, when to appropriately terminate the DEM iterative process is still a difficult problem. That is, in the applications of discrete element modeling, there are often difficulties determining when to terminate the iterative process. When the number of iterations is too few, the calculation results would be inaccurate, while if the number of iterations is too many, it would be redundant.
In this paper, an energy-based discrete element modeling method combined with time-series analysis is proposed for analyzing the deformations and failures of jointed rock slopes. In the proposed method, the convergence of the SD/mean of energy is used to determine when the iteration terminates. The proposed energy-based method is a generalized form of the displacement-based discrete element modeling method [22], and the proposed energy-based method considers not only the displacement of each block but also the weight of each block.
In short, the most essential idea behind the proposed energy-based discrete element modeling method is that, an energy-based criterion is proposed to terminate the DEM iterative process to analyze the deformations and failures of jointed rock slopes. As illustrated in Figure 1, the proposed energy-based discrete element modeling method is roughly composed of the following three steps.
(1) The geological model and computational model are established according to the details of the study area, and the appropriate parameters for the computational model are determined; (2) The DEM modeling is conducted to calculate the deformation of the jointed rock slope under gravity, and the positions and volumes of all rock blocks are achieved; (3) Time-series analysis is carried out to appropriately terminate the iterative process of DEM modeling. Specifically, on the basis of the achieved positions and volumes of all rock blocks, the energy of each block under gravity is first calculated; then, the sequence of the mean and SD of the energy is calculated. Finally, the convergence of the sequence of the SD/mean values is analyzed. If the sequence of the SD/mean values converges, the DEM calculation would be terminated.

Establishment of the Computational Model
The main characters of the MIDAS/GTS software are that it can be used to establish large-scale and complex geological models visually and generate high-quality meshes automatically. The discrete element modeling software 3DEC is inherently quite suitable for numerically investigating the deformations and failures of jointed rock masses. However, it is difficult to establish large-scale and complex geological models, and also cannot be used to generate high-quality meshes automatically.
Therefore, first, the MIDAS/GTS software is used to establish geological models and generate high-quality mesh models; and then, the 3DEC software is used to numerically investigate the deformations and failures of jointed rock masses on the basis of the geological and mesh models created by the MIDAS/GTS software.
The procedure of establishing the geological and calculation models is as follows. First, according to the geological and topographic map of the study area, the elevation of the ground surface is extracted; the ground surface is generated by MIDAS/GTS software, and the ground surface is divided into triangular surface meshes. Then, the initial geological model is established in the 3DEC software according to the prisms generated by triangular surface meshes. Finally, the initial geological model is divided into regions by cutting multiple planes; and key joints are added to generate the final computational model.
The boundary conditions in the numerical analysis are set mainly in terms of displacement. In this paper, the boundary conditions are as follows. Specifically, for the computational model, there are typically six boundary faces or surfaces, including the bottom, the top, and four sides around. The displacements of all nodes locating on the bottom are strictly fixed, without allowing to move. In contrast, all nodes locating on the top surface are free, i.e., the displacements are not fixed. For the four sides around, the displacements in the normal direction of the surface are fixed. For example, the xdisplacements or y-displacements are fixed, while the displacements in the other two directions are not fixed.

DEM Modeling
The DEM is a powerful numerical modeling technique that it could be used to simulate the deformations and failures and large displacements of discontinuous rock masses. At present, the DEM is widely used to investigate the stability of jointed rock slopes. The DEM modeling is an iterative process. However, when to terminate the iterative process is uncertain. The most common method to terminate the iterative process is to specify the maximum number of iteration steps before conducting the DEM modeling.
Based on the computational model established in Section 2.2, different physical and mechanical parameters and joint parameters are set to different regions; and the deformations and failures of the jointed rock slopes under gravity are analyzed. After a certain number of steps (e.g., for every 1000 steps), the positions and volumes of all rock blocks are obtained in 3DEC software until reaching the specified maximum number of iterative steps or a dynamically determined number of iterative steps.

Time-Series Analysis
Time-series refers to the sequence of values of the same statistical index in the order of their occurrence time. The objective of the time-series analysis is to predict the future values based on the existing historical data. Time-series analysis is a quantitative prediction method, which is widely used as a common prediction method in statistics.
In this paper, the energy-based discrete element modeling method coupled with timeseries analysis is used to determine when to terminate the iterative calculation process in the DEM modeling. The states of rock blocks on the slope at different times are shown in Figure 2. Figure 2 shows that from the initial stable T 0 to the final failure T n , the slope experiences the time-steps T 0 , T 1 , T 2 , . . . , T n−1 , T n , and each time-step corresponds to a different state of the slope. The energy under gravity at each time-step could be calculated, and then forms the time-series. Figure 2. A simple illustration of a rock slope changing from stability to failure.

Calculation of the Energy under Gravity
The first step in the time-series analysis is to calculate the energy of the rock slope under gravity.
More specifically, the energy of all rock blocks under gravity are calculated according to the data achieved in Section 2.3. For a three-dimensional rock slope, when the slope is under gravity, some rock blocks on the slope will slide. Figure 3 illustrates the change in displacement of a rock block from T 0 to T n . As illustrated in Figure 3, the total displacement of the rock block from T 0 to T n is ∆d. The total displacement ∆d could be decomposed into ∆x in the X direction, ∆y in the Y direction, and ∆z in the Z direction, thus ∆d could be expressed as Equation (1) Energy could be expressed as the force times the displacement in the direction of the force. Thus, the energy under gravity is equal to gravity times ∆z. ∆z is equal to the Z coordinate of the block center at T n minus the Z coordinate of the block center at T 0 , and the Z coordinate is obtained in Section 2.3. In addition, density is a mechanical parameter determined before DEM modeling, and the energy of all rock blocks under gravity could be calculated according to Equation (2) Energy-Energy under the gravity (N·m) The energy-based method proposed in this paper is a generalized form of the displacement-based method [22]; this is because the proposed method considers not only the displacement of each block but also the weight of each block. The proposed energy-based method is applicable to slopes without an obvious potential sliding surface. The displacement-based method considers that the weight of each block is one; however, the weight of each block obtained by the proposed energy-based method is the actual value that is calculated by volume and density. Thus, the proposed energy-based method is more applicable.

Calculation of the Coefficient of Variations Based on Energy
The second step in the time-series analysis is to calculate the coefficient of variation (CV) of the energy for all rock blocks of a jointed rock slope.
The CV, which is also known as the dispersion coefficient, is a normalized measure of the degree of dispersion of the probability distribution. CV has no dimension. Its size is affected not only by the degree of dispersion of variable values but also by the average level of variable values. It is defined as the ratio of the SD to the mean. In this paper, the CV based on energy is used as an important index to determine when the iterative process would be terminated In this section, the mean of energy and the SD of energy for all rock blocks are achieved by conducting statistical analysis; and the CV of energy corresponding to each iteration step is obtained. Then, a curve of the relationship between the CV of energy and the number of iterations is plotted. The mean of energy is calculated by Equation (3), the SD of energy is calculated by Equation (4) and the CV of energy is calculated by Equation (5). Every CV of energy and the corresponding number of iterations constitute a time series. When the number of iterations reaches a certain value, the relationship curve between CV and the number of iterations begins to converge, and the iterative process would be terminated.
The CV of energy considers not only the displacement of each block but also the weight of each block. Thus, the proposed energy-based method could effectively be used to analyze the deformations and failures of general cases where it is difficult to determine an obvious potential sliding surface of jointed rock slopes compared to the displacementbased method.

Determination of the Number of Iteration Steps
The third step in the time-series analysis is to determine the number of steps when terminating the iteration.
According to Sections 2.4.1 and 2.4.2, the CV (i.e., the SD/mean) corresponding to each iteration step has been achieved. At this time, the convergence of the SD/mean values is investigated by drawing a curve that the abscissa is the number of iterations, and the ordinate is SD/mean, 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, the difference of SD/mean corresponding to T 1 − T 0 , T 2 − T 1 , T 3 − T 2 , . . . , T n − T n−1 are first calculated , and then a curve is plotted as shown in Figure 4b. The difference between the 9000th step and 8000th step is almost zero. With the increase in the number of iterations, the difference is basically equal to zero. In addition, the value of SD/mean does not change. At this time, it is thought that the iterative process could be terminated when iterating to the 8000th step.
The aforementioned method of determining the number of iterations is based on the statistics, i.e., the SD/mean, for the energy of all rock blocks of a jointed rock slope. The calculation of energy considers the displacements and weights of all blocks. The proposed energy-based method is a generalized form of the displacement-based method, and it is much more applicable than the displacement-based method. In the proposed energy-based method, there is no need to determine the potential sliding surface before conducting the numerical modeling. The proposed energy-based method can be applied in the cases where the sliding surface is easy or difficult to be determined.

Verifications
In this section, analyzing the deformation and failure of a simple jointed rock slope is conducted to verify the effectiveness of the proposed method.
The displacement-based discrete element modeling method [22] is used to investigate the jointed rock slopes with a weak interlayer. However, this method needs to determine a potential sliding surface before DEM modeling. The method needs to obtain the displacements caused by the rock blocks located around the sliding surface sliding, and calculate the mean, SD, and the mean/SD of the displacements, and then obtain the number of iterations to terminate the DEM iterative process according to the convergence of the mean/SD.
In this section, a model of a simple rock slope is established; and the displacementbased method and the proposed energy-based method are used to calculate, compare, and verify the feasibility and correctness of the proposed energy-based method in this paper.

Computational Model and Parameters
In this section, a simple jointed rock slope model with a weak interlayer that the dip angle is 20 • is established, and more details are shown in Figure 5.
The model of the simple slope is divided into tetrahedral meshes; the tetrahedral meshes are imported into 3DEC software by the interface program as the computational model. The numbers of nodes and blocks in the established model are 3377 and 15,802, respectively. The calculation parameters of the model are listed in Table 1.
The Mohr-Coulomb model is the most commonly used in the numerical investigation to analyze the stability of slopes [39,40]; it reflects the strength characteristics of rock comprehensively. In this paper, the Mohr-Coulomb model is employed as constitutive model.
In Table 1, the bulk modulus is expressed as BULK, the shear modulus is expressed as G; the density is expressed as ρ; the cohesion is expressed as c; the internal friction angle is expressed as ϕ; the normal stiffness of the joint is expressed as JKN; and the shear stiffness of the joint is expressed as JKS.

Computational Results
By employing the simplified slope model and mechanical parameters in Section 3.1, the Z-direction displacement is obtained when iterating to the 21,000th step, as shown in Figure 6. The slope initially is unstable; the initial stress state is mainly generated by the geological structure and gravity. In addition, topography, mechanical properties of rock mass, water, and temperature also affect the initial stress state. With increasing of the number of iterations, the slope begins to slide gradually. As a result of the existence of weak interlayer and self gravity, and the relatively small mechanical parameters of weak interlayer, the rock slope slides along the weak interlayer; parts of the rock blocks start to slide out of the model boundary. Figure 6 shows that when the iterative step reaches the 21,000th step, a large number of rock blocks have slipped and accumulated, and the sliding blocks slide along the weak interlayer. The positions and volumes of all rock blocks are obtained every 1000 steps until 60,000 steps. After that, the data of the positions and volumes of every rock block are cleaned; and the energy under gravity is calculated. Then, the mean of energy, the SD of energy, and the SD/mean are calculated; and the curve of the relationship between the number of iterations and the SD/mean of energy is drawn.
To verify the correctness of the proposed method, the result of the displacement-based method is compared for verification. The results obtained by using two different methods are shown in Figure 7. It is thought that the suggested number of iterative steps in the displacement-based method is 18,000, and the suggested number of iterative steps in the energy-based method is 21,000; the results are basically consistent. It is found that the inflection points of the two methods are almost the same. Moreover, from Figure 7, it is roughly thought that when calculating to the 21,000th step, the iterative process could be terminated.

Suggested steps in energy-based method is 21,000
Suggested steps in displacement-based method is 18,000

Geological and Engineering Background
The Yanqianshan iron mine is located in Anshan City, Liaoning Province, China; see Figure 8. The geological section of the slope is shown in Figure 9. It is an old mine, approximately 20 km away from the center of Anshan City. The Yanqianshan iron mine was first constructed in 1960, and the open-pit mining was finished in September 2012. Since 2012, transitional ore body mining has been carried out, and the main ore bodies of the mine are in the west, east, and north sides. By the end of March 2015, the mining of iron ore body in the east and west side will be suspended. The ore body at the north side of the pit is underground-mined for the ore body below the bottom of the pit after the end of the transition period. The mining depth will not exceed -500 m for the Yanqianshan iron mine [26,27,41]. At present, the Yanqianshan iron mine is still active and the daily output is approximately 5500 tons.  The first minelayer The second minelayer The mining area is located in the hilly area of Qianshan Branch, surrounded by mountains. The exposed strata in the mining area include Quaternary metamorphic rocks, the Liaohe Formation, and the Anshan Formation, while the main composition of the rock strata is Anshan metamorphic rocks. The basic structure of this area is a monoclinic structure, with a strike of 270 to 300 degrees and a large dip angle, mostly between 70 and 90 degrees. According to the geological section of the slope, and taking the left and right boundary of iron ore as the boundaries, the study area is divided into three sub-areas (i.e., the South, Middle, and North). Moreover, according to the north-south direction, the part in the south is called the south slope, the part in the north is called the north slope. The most common rocks in the southern sub-area are mixed granite rocks, forming massive and gneissic structures. The rocks in this sub-area are greatly affected by the structure and weathering. Banded magnetite quartzite is common in the middle sub-area. The formation is an Anshan-type iron ore bed (orebody), which is mainly composed of magnetite quartzite and hornblende magnetite quartzite. Carbonaceous phyllite is the most common rock formation in the northern sub-area. The rock mass in this stratum is broken and contains a strong weathering zone.

Computational Model and Mechanical Parameters in the DEM Calculation
In this section, the calculation model in the DEM calculation is established, which is used to analyze the deformation and failure of the open-pit slope.
When using the discrete element modeling software 3DEC to analyze the deformations and failures of jointed rock slopes, the essential step is to establish a simple but practical computational model. According to the project report of the open-pit slope, the actual geological conditions are simplified reasonably. The rocks of the slope are divided into migmatite in the southern subarea, ore in the middle sub-area, and phyllite in the northern sub-area. According to the geological and topographic map of the study area, the altitude of the ground surface is extracted; and the ground surface is generated by MIDAS/GTS software. Then, the initial geological model is established in the 3DEC software according to the prisms generated by triangular surface meshes. Finally, simplified joints are added to obtain the required computational model.
To be more specific, establishing the calculation model consists of two procedures: (1) an original geological model without joints is established; (2) establishing the final calculation model by adding key joints to the original geological model.
The first procedure is to establish the initial geological model of the study area, which is completed by using MIDAS/GTS software. This procedure consisting of three main steps is as follows: (1) According to the geological and topographic map of the Yanqianshan iron mine, the elevation of the ground surface is extracted; a representative section is selected, and the ground surface is generated in MIDAS/GTS software; (2) The ground surface is divided into triangular surface meshes in MIDAS/GTS software; (3) An initial geological model is established according to the prisms generated by triangular surface meshes.
The second procedure is to establish the final calculation model by strictly adding key joints, and the procedure is implemented by the 3DEC software. This procedure consisting of three main steps is as follows: (1) The geological model established according to the prisms generated by triangular surface meshes is transformed from the MIDAS/GTS software to the 3DEC software using a special interface program; (2) The excavated roof and ore body are established by using multiple planes to cut models and to divide regions, see Figure 10a; (3) The required calculation model is established by adding the key joints to the model by calling the command jset; see Figure 10b. The parameters are determined according to the project report of the open-pit to underground mining project of the Yanqianshan iron mine. Table 2 shows the physical and mechanical parameters of the rocks in the three sub-areas; and the mechanical parameters of the joints are shown in Table 3.
Currently, the Mohr-Coulomb model is the most commonly employed for investigating the stability of rock slopes in numerical modeling [39,40]. In this paper, the Mohr-Coulomb model is also employed as constitutive model. In Table 2, E represents the young's modulus; µ represents the possion's radio; ρ represents the density; c represents the cohesion; ϕ represents the internal friction angle.

Computational Results
In the DEM calculation, the open-pit slope was excavated layer by layer to analyze the deformations and failures, and the SD/mean was calculated for each iteration interval for the excavation of every layer. Therefore, the relationship between the number of iterations corresponding to the excavation of every layer and the SD/mean values could be plotted. The Z-direction displacement after the excavation of the roof is shown in Figure 11, and the relationship between the number of iterations and the SD/mean values is shown in Figure 12a. The relationship between the difference of each adjacent two steps and the number of iterations is shown in Figure 12b. The excavation of the roof results in the rocks in the north slope moving along the joints to the middle. When the excavation of the roof is completed, from Figure 11, it is thought that the rock blocks on the north slope begin to slide, but no rock blocks fall on the surface of the first minelayer, and the rock blocks in the south slope are stable. Moreover, the excavation of the roof did not cause horizontal plane cracking or sliding. Figure 12 shows that when the number of iterations approaches the 8000th step, the curve in Figure 12a tends to converge, and the SD/mean of the 9000th step is approximately equal to the SD/mean of the 8000th step. Therefore, it is thought that when the calculation reaches the 8000th step, the iterative calculation would be terminated. To make the calculation convenient and accurate, the 10,000th step is chosen to calculate.
When the excavation of the roof is completed and the iteration is terminated, the excavation of the first minelayer of the ore body is carried out. The Z-direction displacement after excavation of the first minelayer is shown in Figure 13, and the relationship between the number of iterations and the SD/mean values is shown in Figure 14a. The relationship between the difference of each adjacent two steps and the number of iterations is shown in Figure 14b. When the excavation of the first minelayer is completed, Figure 13 shows that the sliding of the north slope is relatively more sliding than that of the south slope, the blocks slide along the joints, and a small number of blocks fall on the iron ore bed. Additionally, a small part of the south slope begins to slide, also along the joint planes. From Figure 14, the curve has appeared inflection point and tends to converge, and it is thought when the number of iteration steps reaches the 7000th step, SD/mean tends to converge, the difference between SD/mean of the 8000th step, and SD/mean of the 7000th step tends to zero. Therefore, it is thought that when the calculation reaches the 7000th step, the iterative calculation would be terminated. To make the calculation convenient and consistent with the excavation of the roof, the 10,000th step is also selected.
When the excavation of the first minelayer is completed and the iteration is terminated, the excavation of the second minelayer of the ore body is carried out. The Z-direction displacement after excavation of the second minelayer is shown in Figure 15, and the relationship between the number of iterations and the SD/mean values is shown in Figure 16a. The relationship between the difference of each adjacent two steps and the number of iterations is shown in Figure 16b. After completing the excavation of the second minelayer, Figure 15 shows that with increasing excavation depth, the sliding blocks on the south slope and the north slope gradually increase, but the sliding blocks on the north slope are still more. All the sliding blocks slide along the joint planes. From Figure 16, it is thought that when the iterative calculation reaches the 10,000th step, the curve tends to converge, and the difference between the SD/mean of the 11,000th step and the SD/mean of the 10,000th step tends to zero. Therefore, it is thought that the iterative calculation would be terminated when another 10,000 steps are calculated.
The above excavation is divided into three times: the first excavation is the roof; the second excavation is the first layer of iron ore; and the third excavation is the second layer of iron ore. The results of the three excavations could verify the feasibility and applicability of the proposed energy-based method. However, with the increase in excavation layers, a large number of rock blocks would accumulate on the surface of the iron mine. The next excavation needs to first remove the sliding blocks, and then carry out the excavation of the underlying layers. In brief, according to the results of the above three excavations, it is thought that the proposed energy-based method in this paper is applicable to more jointed rock slopes that cannot determine the potential sliding surface before DEM modeling.

Discussion
The novelty of the proposed energy-based discrete element modeling method is that, it is more applicable than the displacement-based discrete element modeling method because it does not need to determine the position of the potential sliding surface before DEM modeling.
The proposed energy-based method is a generalized form of the displacement-based method. This is because it can be applied in the cases where the sliding surface is easy or difficult to be determined. The proposed energy-based method considers not only the displacement of each block but also the weight of each block. Moreover, the computational cost of the proposed energy-based method is approximately the same as that of the displacement-based method. In summary, compared with the displacement-based method, the proposed energy-based method does not need to determine the potential sliding surface before DEM modeling. Thus, the proposed energy-based method is much more applicable than the displacement-based method.
The energy-based method is applied to investigate the mining-induced slope deformation in the Yanqianshan iron mine. In the application case, the potential sliding surface of the slope is difficult to be determined; in this case, the displacement-based method is difficult to be applied. The computational results of the Yanqianshan iron mine are sufficient to indicate when to terminate the iterative process. The computational results also indicate that the proposed energy-based method is applicable.
As mentioned above, the energy-based method proposed in this paper could be employed to analyze the deformations and failures of both simple and complicated slopes without considering the shape and size of the slopes. However, the shortcomings of the proposed energy-based method are that it needs to investigate the displacements and weights of all rock blocks, rather than part of the rock blocks.
In the future, the Python script language is planned to be used to integrate into the discrete element calculation software, where it enables the conducting of real-time calculations and automatic termination of the iterative process.

Conclusions
In this paper, an energy-based discrete element modeling method combined with timeseries analysis is proposed to analyze the deformations and failures of jointed rock slopes. The most essential idea behind the proposed method is the definition of an energy-based criterion to determine when to terminate the DEM iterative process in the investigations of the deformations and failures of jointed rock slopes. The proposed energy-based method is more applicable than the displacement-based discrete element modeling method because it does not need to determine the position of the potential sliding surface before DEM modeling. The proposed energy-based method is a generalized form of the displacementbased method; the proposed energy-based method considers not only the displacement of each block, but also the weight of each block, and the computational cost of the proposed energy-based method is approximately the same as that of the displacement-based method. Moreover, the effectiveness and correctness of the proposed energy-based method are verified through analysis of a simplified jointed rock slope by comparing it to the displacement-based method, and then the proposed energy-based method is applied in a real case. In summary, the proposed energy-based method could effectively be used to analyze the deformations and failures of general cases where it is difficult to determine the obvious potential sliding surface of joint rock slopes. In the future, Python script language would be used for integration of the proposed energy-based method with the discrete element calculation software to perform out real-time calculations and automatically terminate the iterative process.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.