Analysis of Flexural Toppling Failure in Rock Slopes Using Discrete Element Method

Flexural toppling failure is a common failure mode of natural and artificial rock slopes, which has caused serious damage to human life and property. In this work, an advanced numerical method called the Universal Distinct Element Code (UDEC) was used to study the mechanism of flexural toppling failure. In total, more than twenty slope models were built and analyzed. Two new parameters (displacement discontinuity and transition coefficient of failure surface) were introduced to present a further understanding of flexural toppling. The results show the failure zone of rock slopes subjected to flexural toppling includes two parts: the first-order instability part (FOIP) and the independent toppling zone (ITZ). The FOIP can be further divided into two subzones: the sliding zone (SZ) and the superimposed toppling zone (STZ). The occurrence of surface deformation discontinuities is the precursor to flexural toppling failure. The first displacement discontinuity occurs on the boundary between the FOIP and the ITZ. The angle, spacing, and angle of the joints, the angle of the slope has a significant influence on the stability of anti-dip bedding rock slopes. However, they do not affect the deformation and failure pattern of the slope.


INTRODUCTION
Flexural toppling failure often occurs in artificial and natural slopes with anticlinal bedding structures (called anti-dip bedding rock slopes, as shown in Figure 1), which has been recognized and studied for several decades. Anti-dip bedding rock slopes are a very common type of slope in nature, accounting for 33% of all landslides (Huang and Li, 2011). Ashby (1971) proposed the concept of "toppling", and De Freitras and Watters (1973) verified the phenomenon through field observations. Then, toppling failure has been considered as an important failure type as sliding. To better analyze it, Goodman and Bray (1976) proposed an exclusive classification for toppling instabilities, and flexural toppling is one of them.
In a two-dimensional model, rock layers in such failure behave like cantilever beams, which bend and deform under the combined action of gravity and external loads. Aydan and Kawamoto (1992) used the solid mechanics of cantilever beams to analyze this type of failure and proved it by base friction tests. Centrifugal tests and physical model tests under excavation at the toe area of the slope were also employed by many researchers (Adhikary et al., 1997;Zheng et al., 2019;Zheng et al., 2018a;Zhu et al., 2020). Based on these experimental results, they made some improvements for the analytical method and gave suggestions for the reinforcement of the slope. Meanwhile, Ning et al. (2019) and Li et al. (2019) used shaking tables to simulate this type of failure under seismic load. The shaking table was also utilized by Fan et al. (2016) to discuss the effect of the weak intercalation on the failure patterns of antidip bedding rock slopes. In addition to seismic load, Ning et al. (2021) and Gu and Huang (2016) suggested that the river downcutting should also be considered as an incentive of flexural toppling failure.
Numerical tools are also often employed by scholars and engineers to discuss and analyze flexural toppling failure, especially when the geological conditions of the slope were very complex. A finite element method was undertaken by Adhikary and Dyskin (2007) to analyze their physical model tests. Based on the distinct lattice spring model, Lian et al. (2018) analyzed flexural toppling failure by incorporating the gravity increase method. Zhang et al. (2020) tried to take the discontinuous deformation analysis (DDA) to simulate the tensile failure of blocks in toppling mode, but this method needs to add contacts at the failure surface in advance. As a commercial software commonly used in geotechnical engineering, the universal distinct element code (UDEC) was also utilized by a lot of researchers. For example, Alzo'ubi et al. (2010) used it to analyze the effect of tensile strength on the stability of slopes. Zheng et al. (2019a) used it to investigate the influence of rock bolts on the failure patterns. Weng et al. (2020) implemented a new failure criterion into this software to predict this failure.
As outlined above, many scholars have studied the mechanism of flexural toppling failure through theoretical analysis, model tests, and numerical simulation since its recognition as an important failure mode. However, the whole process of initiation, development, and penetration of the failure surface of the slope subjected to flexural toppling is still unclear. Compared to model tests, numerical simulations are less expensive and allow easy access to parameters that cannot be measured in model tests Wang et al., 2020;Yang et al., 2020;Meng et al., 2021;Yan et al., 2021). UDEC is ideal in modelling jointed rock masses since it can simulate both the behavior of the intact rock and the joints (Itasca, 2004). Furthermore, UDEC can simulate large deformations, such as tensile cracks between rock layers, which are a key feature of toppling failure. In this work, UDEC is thus used to investigate flexural toppling failure in rock slopes with various geometric and mechanical parameters. Two new parameters (displacement discontinuity and transition coefficient of failure surface) are introduced to present a further understanding of this type of failure. Moreover, based on the results of systematic numerical simulations, the whole evolution process of flexural toppling failure is discussed.

Model Configuration and Boundary Conditions
In order to understand the fundamental mechanism of flexural toppling in rock slopes, the basic calculation model used in this work is an ideal anti-dip bedding rock slope ( Figure 2). The slope   height and angle are 40 m and 63°, respectively, and the joint spacing and angle are 1.6 m and 78°, respectively. The geometric and mechanical parameters of the model are shown in Table 1. Displacement constraints were used in the numerical simulation. Specifically, horizontal and vertical displacements at the base of the model were constrained while only the horizontal displacement was restricted at the lateral boundaries. A series of monitoring points (1 m apart) were laid out on the surface of the slope to monitor the horizontal and vertical displacement. The number of the measuring points is assigned as No.1 to No.46 starting from the toe of the slope to the top of the slope, as illustrated in Figure 2.

Constitutive Model and Model Parameters
The deformation and failure of anti-dip bedding rock slopes involve both intact rock and joints. In this work, two commonly used constitutive models, namely the Mohr-Coulomb model and the Coulomb slip model, were employed to describe the mechanical behavior of the intact rock and the joints, respectively.
Toppling failure mostly occurs in slopes with lithology or lithological combination of thin layered slate, schist, slate, sandstone, and muddy limestone (Tu et al., 2007;Amini et al., 2018;Zheng et al., 2018a;Huang et al., 2018;Zhang et al., 2018;García-Moya et al., 2019;Ning et al., 2021). This work relies on an actual limestone slope facing the Jingzhu Road in Southern China (Zheng et al., 2018a). Therefore, the lithology modelled in this work is limestone with low strength. The parameters used in the calculations are summarized in Table 1.

Influencing Factors Studied in the Simulation
In order to have a comprehensive understanding of the flexural toppling mechanism, four factors were studied in the simulations: angle of the joints, spacing of the joints, strength of the joints, and angle of the slope. Each factor includes 6 to 11 scenarios, and there are 26 models in total. It is important to note that a univariate analysis is used, i.e., when analyzing the influence of one parameter on the flexural toppling, the other parameters keep constant.

RESULTS ANALYSIS
In this section, distribution characteristics of the surface displacement and the plastic zone were studied to access the whole process of formation of the failure surface of the slope, and thus reveal the mechanism of flexural toppling failure. Moreover, the factors that play a controlling role in flexural toppling can be also found. Whether the surface displacement is perpendicular to the joints is a preliminary indication of flexural toppling. If it is perpendicular to the joints, the slope is likely to experience flexural toppling failure, otherwise the slope is experiencing sliding failure.

Effect of the Joint Angle
The prerequisite to the occurrence of toppling failure is that rock layers are misaligned with each other (Goodman and Bray, 1976).
In the shallow part of the slope, the direction of the maximum principal stress is approximately parallel to the surface of the slope, as shown in Figure 3.
Only when the angle between the maximum principal stress and the normal to the joints is greater than the friction angle of the joints (ignoring the joint cohesion), mutual misalignment between adjacent rock layers is possible. It can be written as (Goodman and Bray, 1976): where β c and η are the angle of the slope and the angle of the joints, respectively, and φ j is the friction angle of the joints. The stability of rock slopes in six cases of η 40°, 50°, 60°, 70°, 78°, and 90°were studied in the numerical simulations. Figure 4 shows the distribution of horizontal displacement (xdis), vertical displacement (ydis), and ydis/xdis for each measuring point for different joint angles. The vertical black dashed line indicates the position of the crest of the slope and the vertical red solid line indicates the position where the first displacement discontinuity (FDD) occurs. The horizontal solid line indicates the value of tan β(the angle of the normal to the joints), which is used to determine the direction of the displacement. ydis/xdis tan β means that the direction of the displacement is perpendicular to the joints, which is an important feature of toppling deformation.
From Figure 4, we can find that, as the value of η increases, the slope displacement tends to increase first and then decrease. The measuring point with the largest displacement is not necessarily located at the crest of the slope. When there are displacement discontinuities, the largest displacement tends to occur at the position where the first displacement discontinuity occurs, as shown in Figure 4D, Figure 3F. When η 40°, the value of ydis/ xdis at most measuring points is less than tan β, which means that the direction of the displacement is not perpendicular to the joints. In other words, insignificant toppling deformation of the slope occurs. However, as the value of η exceeds 50°, the displacement direction below the crest of the slope is approximately perpendicular to the joints, except for the FIGURE 3 | Toppling stability analysis in an anti-dip bedding rock slope (Goodman and Bray, 1976 points near the toe of the slope. That is, the slope exhibits typical characteristics of toppling failure. As the value of η increases to 90°, the displacement direction after the crest of the slope is no longer perpendicular to the joints and the toppling deformation becomes unremarkable. Figure 3 indicates that no displacement discontinuities are found at η 40°and 50°. However, the situation changes immediately when the value of η exceeds 50°. Displacement discontinuity means that a tensile crack appears at the corresponding position along the joint, as shown in Figure 4D. The first tensile crack above the toe of the slope forms the trailing boundary of the first-order instability part (FOIP) and controls its range. Thus, the location of the first displacement discontinuity is a key factor in the analysis of  Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 773088 5 toppling failure. To analyze this factor quantitatively, a new parameter named transition coefficient of failure surface (λ t ) is introduced: where m d and m c are the measuring point numbers at the first displacement discontinuity and the crest of the slope, respectively.
The variation in the value of transition coefficient of failure surface (λ t ) as a function of joint angle (η) is plotted in Figure 5A. It is seen that the value of λ t tends to increase overall as the value of η increases, which indicates that the trailing boundary of the FOIP gradually moves upwards. Figures 4B-E show the simulated plastic zone obtained using UDEC for different joint angles (η). Due to the similarity in the distribution of plastic zone at η 60°, 70°, and 78°, only a FIGURE 7 | Surface displacement distribution of slope for different joint spacing (t). N C and N FDD : Measuring points corresponding to the slope crest and first displacement discontinuity, respectively.
Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 773088 6 representative one η 78°is shown here. It can be found that the plastic zone only appears around the toe of the slope at η 40°. The corresponding displacement distribution indicates that no tensile cracks occur along the joints ( Figure 4D). The slope is therefore stable when η 40°, which is consistent with the result obtained from the criteria (Eq. 1) proposed by Goodman and Bray (1976). When the value of η increases to 50°, a through plastic zone is formed from the toe to the top of the slope, which indicates that the slope undergoes a complex shearing-flexural toppling failure ( Figure 5C). Although there is no through plastic zone within the slope at η 78°, conditions for the formation of a landslide are met: the plastic zone forms the sliding surface, and the tensile crack forms the rear boundary of the sliding body ( Figure 5D). Under the action of gravity, rock mass above the plastic zone moves towards the toe of the slope, thus forming the FOIP. After this, due to a large space available, the rock layer adjacent to the free surface undergoes a significant toppling deformation until instability occurs ( Figure 6). In this case, the rock layer works as a cantilever beam and the forces acting on both the left and right sides can are zero. Similar deformation and failure occur in other layers above until the length of the cantilevered section is less than the corresponding self-stabilizing height (Zheng et al., 2018a).

Effect of the Joint Spacing
Joint spacing (t) has been found to be an important factor affecting the stability of anti-dip bedding rock slopes (Adhikary et al., 1997;Adhikary and Dyskin, 2007;Zheng et al., 2020a). Here, six cases of t 0.5, 1.0, 1.6, 2.0, 2.5, and 3.0 m were simulated to further investigate the effect of joint spacing on the deformation and failure patterns of this kind of slope. Figure 7 shows the distribution of horizontal displacement (xdis), vertical displacement (ydis) and ydis/xdis at each measuring point for different joint spacing. Figure 8 shows the failure distribution of the slope for different joint spacing (t).
It can be seen from Figure 7 that the surface displacement tends to increase and then decrease as the joint spacing increases and the maximum value occurs at the first displacement discontinuity. Except for the measuring points after the crest of the slope at t 0.5 m, the direction of the surface displacement is basically the same for different joint spacing. They all show that the displacement at most points is perpendicular to the joints, except for a few points near the toe area of the slope. In other words, the deformation of the slope is dominated by toppling in all cases. These imply that joint spacing has an insignificant effect on the deformation pattern of this kind of slope. Displacement discontinuities (tensile cracks) occur in all cases. From Figure 8A, it can be seen that the first displacement discontinuity moves upwards as the value of t increases. The first displacement discontinuity is located below the crest of the slope when t < 2.5 m, and above the crest after reaching 2.5 m.
Figures 7B,C shows the simulated plastic zone obtained using UDEC for different joint spacing (t). It should be noted that the distribution of the plastic zone within the slope is similar for different joint spacing, so only the results for t 0.5 and 1.0 m are shown here. It can be seen that several rock layers at the toe area undergo shear failure while those above undergo complex tensileshear failure. The range of the FOIP for t 0.5 m is less than that for t 1.0 m. This is because the smaller the joint spacing, the lower the bending resistance of the rock layer, and the smaller the corresponding range of self-stabilizing rock layers at the toe area of the slope (Zheng et al., 2018b;2019b).

Effect of the Joint Strength
According to Eq. (1), to prevent the slope from toppling failure, the friction angle of the joints (φ j ) needs to be greater than a critical value, φ jcr η + β c − 90 o . Taking η 78°and β c 63°into the equation gives φ jcr 51°for the basic model. In this work, the deformation and failure mechanism of slopes were analyzed for eleven cases: φ j 0°, 6°, 12°, 18°, 24°, 30°, 36°, 42°, 48°, 54°, and 60°. Figure 9 shows the distribution of horizontal displacement (xdis), vertical displacement (ydis) and ydis/xdis at each measuring point for different joint friction angles. It should be noted that, for the sake of simplicity, only representative graphics are shown. Figure 10 shows the failure distribution of the slope for different joint friction angles (φ j ) It can be seen that, when φ j 0°, 12°, and 36°, the surface displacements do not differ much, and the maximum values are all located at the first displacement discontinuity. However, after φ j reaches the value of 60°, the surface displacements decrease sharply. Moreover, the first three cases all show discontinuities in surface displacement (i.e., tensile cracks) whereas when φ j 60°, it is continuous. The position of the first displacement discontinuity (trailing boundary of the FOIP) generally moves upwards as the value of φ j increases ( Figure 10A). Figures 9B-E show the simulated plastic zone obtained using UDEC for different joint friction angles (φ j ). It can be seen that, for φ j 0°, 12°, and 36°, the FOIP is formed within each slope, which means that the corresponding slope has failed. However, no plastic zone occurs within the slope for φ j 60°and it is in a stable state. The simulated results demonstrate that it is feasible to estimate the critical friction angle using the equation of φ jcr η + β c − 90 o .

Effect of the Cut Slope Angle
The angle of the slope (β c ) is a key parameter in slope design. The toppling stability of an anti-dip bedding rock slope is also closely related to the value of β c (Eq. 1). Here, six cases of β c 30, 43, 53, 63, 73, and 83°were studied using UDEC. Figure 11 shows the distribution of horizontal displacement (xdis), vertical displacement (ydis) and ydis/xdis at each measuring point for different joint friction angles (only representative graphics are shown). Figure 12 shows the failure distribution of the slope for different slope angles (β c ). From Figure 11, it can be found that the surface displacement increase as the slope angle increases and the position of the maximum displacement point varies with β c . The larger the value of β c , the closer the maximum displacement point is to the toe of the slope. The displacement direction is insensitive to the change of the angle of the slope. Except for a few points near the toe of the slope, the displacement vector of each point is basically perpendicular to the joints (xdis/ydis≈tan β), regardless of the slope angle taken, and all slopes present typical characteristics of toppling deformation. Figure 11 also shows that when β c 30°, no displacement discontinuities occur on the surface of the slope. However, the situation changes immediately as the value of β c exceeds 30°, which implies that the occurrence of slope instability. As the value of β c increases, the overall trend of λ t deceases ( Figure 12A), indicating that the larger the value of β c , the closer the first Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 773088 displacement discontinuity (tensile crack) is to the toe of the slope. The reason for this is that the larger the angle of the slope, the less stable the rock layer is, and thus the smaller the range of the self-stabilizing section of the rock layers. Figures 11B-E shows plots of the simulated plastic zone obtained using UDEC for different slope angles (β c ). When β c 30°, the absence of a through plastic zone within the slope means that the slope is in a stable state. However, the distribution of plastic zones in other cases differs considerably from it. The plastic zone and the tensile crack form the bottom slip surface and the trailing boundary of the sliding body, respectively, which means that the slope is unstable (Figure 12c-12E). The results show that whether a slope undergoes flexural toppling failure or not is closely related to the angle of the slope.

DISCUSSION
Based on the results of the above numerical simulation, we can find that flexural toppling failure starts at the toe of the slope (Figure 13). Under the combined action of gravity and pushing forces from the upper rock layers, the layers at the toe area of the slope are the first to lose stability. The instability mode of each rock layer is controlled by the relative magnitude of its resistance to bending and shearing. If the bending resistance of a rock layer is greater than its shearing resistance, the layer undergoes shearing failure; otherwise, it undergoes flexural toppling failure (Zheng et al., 2018a;2018b).
Due to the small ratio of height to width of the rock layers around the toe of the slope, they have strong bending resistance and are therefore less prone to toppling failure. Under the combined action of gravity and pushing forces from the upper rock layers, the layers at the toe area often undergo shearing failure, forming the sliding zone (SZ). However, as the distance from the toe of the slope increases, the bending resistance of rock layers decreases due to the increase of the layer height. Then, the failure mode of rock layers gradually changes from shearing to flexural toppling, forming the superimposed toppling zone (STZ). Flexural toppling failure extends to the first displacement discontinuity (tensile crack), at which point a through plastic zone appears within the slope ( Figure 5D). The through plastic zone referred to here is a penetration to the tensile crack rather than the surface of the slope, forming the boundaries of the first-order instability part (FOIP) together. After the occurrence of FOIP, the upper rock layers then have enough space for deformation that they undergo flexural toppling failure under the action of gravity, forming the independent toppling zone (ITZ), as shown in Figure 14. The failure surface of an anti-dip bedding rock slope is a complex stepped surface rather than a simple plane ( Figure 13). Thus, finding the failure surface is a hard issue in carrying out stability assessment of this kind of slope. Limit equilibrium methods combined with optimization algorithms such as Genetic Algorithm and Bidirectional Evolutionary Structural Optimization (BESO) method have proved to be effective means of solving this problem (Zheng et al., 2020a(Zheng et al., , 2020bLiu et al., 2021). Taking the first displacement discontinuity as the boundary, the lower rock layers are in close contact with each other despite the occurrence of interlayer slip, and the overall deformation is characterized as 'superimposed cantilever beams'. However, for the upper rock layers, tensile cracks appear between adjacent layers, and the side forces is approximately zero. Each rock layer can be considered as an 'independent cantilever beam' (Figure 14). Thus, separate mechanical models need to be developed for these two parts of the rock layers (Zheng et al., 2018a).

CONCLUSION
In this work, the mechanism of flexural toppling failure in antidip bedding rock slopes is systematically studied using numerical simulation, and the main conclusions obtained are as follows: 1) The failure zone of rock slopes subjected to flexural toppling includes two parts: the first-order instability part (FOIP) and the independent toppling zone (ITZ). The FOIP can be further divided into two subzones: the sliding zone (SZ) and the superimposed toppling zone (STZ). 2) The mechanical behaviors of rock layers in the FOIP and ITZ are different. In the FOIP, rock layers are in close contact with each other and behave as "superimposed cantilever beams".
However, the layers in the ITZ behave as "superimposed cantilever beams" and can deform freely.
3) The occurrence of surface deformation discontinuities can be considered to be a precursor to flexural toppling failure and the first displacement discontinuity is the boundary between the FOIP and the ITZ. 4) The angle of the joints, spacing of the joints, strength of the joints, and angle of the slope have significant influences on the stability of anti-dip bedding rock slopes. However, their variations do not change the deformation and failure pattern of the slope.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

ETHICS STATEMENT
Written informed consent was obtained from the (individual(s) AND/OR minor(s)' legal guardian/next of kin) for the publication of any potentially identifiable images or data included in this article'.