The Time-Dependent FailureMechanism of Rocks and Associated Application in Slope Engineering: An Explanation Based on Numerical Investigation

In this study, a numerical model for long-term deformation and progressive failure of rock slope is presented.,emodel accounts for both rock heterogeneity and the initiation, activation, nucleation, and coalescence of cracks in rock slope through a stochastic local stress field and local rock degradation by using an exponential softening law.,e time-dependent behaviour of rocks is taken as a macroscopic consequence of damage evolution and strength degradation in microstructure. A series of demonstrative slope cases containing preexisting joints are constructed and investigated.,e slope instability occurs at a particular point in time when the rock strength is reduced to a certain value. ,e temporal and spatial evolution of joint linkage structures is numerically obtained, which clearly shows how the local stress field and damage evolution within the joint network contribute to the fracture pattern and the long-term instability.,en, a practical slope case in jointed and layered rock formations in Yunyang city is studied. ,e prevailing failure phenomena of the slope, including gradual surface scaling, sliding collapses, and block falling, are numerically reproduced, with an emphasis placed on the slope failure process and development tendency.,ere is a good agreement on the failure mode and instability time between the numerical simulations and the field observations.


Introduction
In recent years, numerical methods have led to significant developments in the analysis of rock slope stability. A number of numerical models, including PFC [1], FEM [2,3], DDA [4], DEM [5][6][7], coupled FEM/DEM [2], and ELFEN [8], have been developed and applied to model the slope deformation and failure. Based on the previous studies, the following two points of view have been generally accepted.
First, the long-term stability plays an important role in the design of slope engineering of fractured rock mass [9]. Stead et al. [10] analyzed varied failure mechanisms of slope using a combined FEM/DEM method, in which the importance of the consideration of damage, fatigue, and time was emphasized. e physical properties of the slope can be affected by the boundary conditions of the rock mass and subsequent stress redistribution, cyclic loading (e.g., thermal-or hydro-mechanical loading), and near-surface physical, chemical, and biological weathering reactions, which can drive progressive development of fractures through intact rock connecting nonpersistent discontinuities [2,[11][12][13][14]. Such accumulated damage reduces the strength of rock mass and contributes to the failure of slope [15].
Secondly, instabilities of rock slopes are generally related to the preexisting discontinuities. Because in many cases, nonpersistent discontinuities are involved in the rock failure process and the main origin of destabilization of the rock mass is the breakage of rock bridges between those joints [16,17]. e progressive failure mechanisms of rock bridges lead to an entire failure of slope. is progressive failure mechanism explains the time delay between the initial localized deformation and the entire failure [5].
Consequently, the creep (time effect) is one of the unavoidable problems in slope stability analysis. e rock slopes, including the open pit slopes, traffic cutting slopes, and the high slopes in water resources and hydropower engineering including large earthen dams, are susceptible to hazards induced by the time-dependent rock behaviour during and after construction [2,[18][19][20][21][22]. During the entire life of a slope, the strength of rock mass in slope may gradually degrade due to stress redistribution, near-surface physical, chemical, and biological weathering, or cyclic loading (e.g., thermal-or hydromechanical). In addition to the soft rock, even hard rock is susceptible to weathering erosion. For example, the weathering of granitic and gneissic rocks in tropical regions can reach depths of more than 100 m [23,24].Ündül and Tugrul [25] investigated the changes in physical properties of dunite due to weathering. As shown in Figure 1(a), it is found that dunites generally show similar trends to other rock types; i.e., uniaxial compressive strength decreases with weathering process [25,[29][30][31]. e tests on the influences of freeze-thaw cycles on the mechanical properties of biotite granite showed that the elastic modulus reduced exponentially with the number of freeze-thaw cycles, as well as compressive strength and cohesive strength (as shown in Figure 1(b) [28]). Temperature variation has a direct effect on the strength and other physical properties of rocks. As the Arrhenius term exp(− (Q/RT)) has a sound physical background, it is used in most of the creep laws for steady-state creep constitutive equations [32].
In addition, the mechanical behaviour of rock can be strongly affected by water. Water affects rock strength through two important effects, e.g., physicochemical and mechanical effects. e lowering of the effective stress lowers the rock failure strength. e physicochemical effect results in the rock strength being time dependent [20,33,34]. is physicochemical effect is suitable for all kinds of rock, particularly for soft rocks. Zhou et al. [35] found that the mechanical property degradation of silty mudstone with saturated time conforms to an exponential law, as shown in Figure 2.
No matter what causes the damage (damage induced by stress, water, chemical, and biological weathering or cyclic loading), it is believed that there is a close relationship between the macroscopic long-term behaviour and the mesoscopic damage [36][37][38][39]. Based on the characteristics of the long-term strength and damage mechanics, numerous studies on the long-term stability of rock mass have been widely carried out. For example, Korezniowski studied the creep and crack evolution with time and then the macroscopic fracture propagation, which eventually leads to the failure of rock mass from local to global instability because of the deterioration of environmental conditions [40]. Aubertin et al. introduced the damage variable and established a rheological model expressed by the internal state variables for salt rock [41]. e continuum damage mechanics is applied to the rheological analysis of salt rock. Shao et al. believes that the time effect of the materials in macroscopic scale is due to the changing microstructure of materials. A rheological constitutive model which can describe the rheological properties of rock and concrete is established, which is similar to that of the internal damage variable [38]. Kemeny [21] developed a model for the timedependent degradation of rock joint cohesion. Using the model, a probabilistic slope analysis was conducted, and the probability of failure as a function of time was predicted. Eberhardt et al. [2] simulated the influence of damage due to unloading by glacial ice and cohesive strength degradation with time in the stability of the Randa rockslide using finite element modelling. Based on the hypothesis that rock properties would degrade with time because of weathering/ alteration processes, Chemenda et al. [42] employed DEM to reproduce the gravity-driven destabilization of slopes. Utili and Crosta [7] conducted a numerical analysis on slope stability, in which slope instability follows the decrease in rock strength in time and space with the only assumption concerning the weathering distribution within the slope. Based on the theory of mesoscopic damage mechanics, a rheological model is established by Dashnor et al. [43] to describe the long-term characteristics of natural gypsum. In their proposed model, the damage is employed, which can be used to describe the accelerated phase of the rheology. It is proved that there is a close relationship between the macroscopic deformation and the mesoscopic damage [43]. Other studies on the time-dependent deformation and failure of slopes are also be found in the previous literatures [18,19,22,23,44,45].
In this study, as an alternative approach, a rheological model based a strength-degradation law was employed and implemented in RFPA [46]. en, numerical simulations were performed and attempted to gain an insight into the mechanism of time-dependent deformation and damage evolution of jointed rock slopes.

Rheological Model: Rock Failure Process Analysis (RFPA) Code
In microscopic scale, rocks comprise an aggregate of crystals and sometimes amorphous particles of mineral, with sizes in the range of millimeters and centimeters. If the rocks are weathered or stressed, dilatation or local damage will happen. On a larger scale, rocks are mostly continuous and homogeneous, but quite often joints, faults, bedding planes, or different strata appear. e result is that rocks do not behave homogeneously as metals usually do. Under certain conditions, this micro-or macrostructure results in a rather complicated mechanical behaviour. erefore, a rheological model based on micromechanics is worth trying to be built. Particularly, a model to investigate the stability of rock slope should encompass the initiation, activation, nucleation, and coalescence of rock cracks. In this study, the FEM-based code, rock failure process analysis (RFPA), is employed, and the widely accepted damage constitutive law for an element subjected to uniaxial or triaxial stress states is also adopted here, which had been illustrated elaborately by Lemaitre and Desmorat [47] and Li et al. [48]. In RFPA, FEM is used as the basic stress analysis tool, where the four-node isoparametric element is used as the basic element in the finite element mesh. Numerical modelling on the long-term deformation 2 Mathematical Problems in Engineering and failure of jointed rock slopes is conducted based on the following assumptions: (i) e jointed rock slope is assumed to be composed of many mesoscopic FEM elements, as shown in Figure 3. e properties of each element are specified according to Weibull distribution. (ii) e mesoscopic element of rock is assumed to be elastobrittle with a residual strength. e mechanical property of rock is defined by an elastic damage constitutive law, as illustrated in Figure 4, and the residual strain/deformation upon unloading is not taken into account. (iii) e mechanical properties of an element are considered as isotropic when the element size is small enough. e maximum tensile stress criterion is preferred to be used to determine whether an element is damaged. If the element is not damaged in tension, the Mohr-Coulomb criterion is then used. (iv) e time-dependent behaviour of rock is considered as a macroscopic consequence of mesoscopic     [25], Arel and Tugrul [26], and Tugrul [27]); (b) freeze-thaw cyclinginduced degradation of compressive strength of granite (Tan et al. [28]). damage evolution and the stiffness and strength degradation at the elemental scale; i.e., it is assumed that the time-dependent behaviour of rock is considered as a macroscopic consequence of the instantaneous or short-term damage development depending on the local stress level in a certain time interval (Δt) and strength degradation subjected to environmental change in the long-term period.
In this study, based on a general understanding of timedependent behaviour of rocks (e.g., the strength degradation pattern shown in Figures 1 and 2), a rock degradation law, i.e., an exponential relation between time-dependent strength and the time to failure, is introduced when subjected to a constant stress σ i (major stress on this element) smaller than its instantaneous strength σ 0,i , as shown in Figure 5, and expressed as follows: where σ t i denotes the time-dependent strength of the mesoscopic element at the moment t i and σ ∞ denotes the long-term strength at the moment t. At this point, the time is assumed to be infinite, i.e., t i ⟶ ∞; σ 0,i denotes the initial instantaneous strength; and a s denotes the strength degradation coefficient. e damage in the mesoscopic element will happen when the stress σ i satisfies the failure criterion σ 0,i or the time t reaches its failure time t i . en, all the damage variables, stress, strength, and failure times of elements, are updated. Letting σ ∞ /σ 0,i � k, the following equation is obtained: In our numerical model, Young's modulus of the mesoscopic element is also assumed to conform to the degradation law, as shown in equation (2).
In a certain time interval (Δt), the instantaneous damage evolution of each element depends on local stress level and the stress-strain relationship of an element is assumed to be linear elastic. After the given damage threshold is reached, the softening is followed. e damage-induced energy e released from each element can be calculated by e � (σ 2 /2E)V, where σ is the peak strength of the element, E is Young's modulus, and V is the volume of an element. Based on the models above, the acoustic emission (AE) associated with the progressive failure of the rock element can be simulated. A flow chart of the model is presented in Figure 6.

Validation of the Numerical Model
e validity of the numerical model was tested. e sandstone, from Yunyang city in the southwestern China, was chosen in the creep experimental tests by Jiang et al. [49].
e stress levels applied on the rock specimens were 77.8%, 67.8%, and 56.2% of the uniaxial compressive strength of the rock, e.g., 45.13 MPa, 39.34 MPa, and 32.58 MPa, respectively [50], in order to be identical to the study by Jiang et al. [49]. e physical mechanic parameters of tested sandstone specimens are Young's modulus (E 0 ) 6.0 GPa, uniaxial Joint Rock matrix  e ratio of long-term strength to short-term strength was set to be 0.7 according to the experimental data, and the degradation coefficient of the element was chosen with trial and error to be 0.05. e geometry of the model specimen is 100 mm × 50 mm, and a 100 × 50 (5,000 elements) square grid is used, so the element size equals to 1 mm. e shear stress fields and acoustic emission (AE) distributions for the simulation at 45.13 MPa are shown in Figure 7. e grey scale in Figure 7(a) denotes the relative magnitude of the stress, where lighter areas are indicative of higher shear stresses, and vice versa. e black areas in Figure 7(a) represent completely damaged elements. In the modelling, when the ultimate tensile strain or stress of the element is met, the element will be displayed as cracks with black color in the numerical results. One of the main advantages of this model is that there is no need for setting a preexisting crack to simulate the crack initiation and propagation. is approach is similar to a smeared crack model. e advantage is that the mesh topology is kept to be untouched [51][52][53]. Monitoring AE event is a good way to identify the initiation and propagation of cracks of rocks. Figure 7(b) shows the AE evolution with the rock sample failure. In Figure 7(b), each circle is corresponding to an AE event. Shear damage is represented by white circle, and tensile damage is represented by red circle. e diameter of the circles represents the relative magnitude of the AE energy.
Some researchers [54][55][56] suggested that the AE distribution during a rock creep experiment could indicate the initial damage, and the acceleration of the tertiary creep phase could be resulted from microcrack interaction, coalescence, and localization of previous failure. e AE distributions show that the initial state of damage was distributed throughout the sample (Figure 7(b) A). is was due to the wide distribution of weak elements in the numerical model. Some scattered AEs had occurred, resulting in a noninteracting microcrack network. Shear cracks (white circles) dominated the deformation during the initial loading. As damage accumulated, the microcracks became increasingly clustered (Figure 7(b) B)), thus involving (a)  Based on experimental and numerical data, the strain versus time curves are plotted in Figure 8. It is shown that the numerical data are in good agreement with the experimental data. Particularly, three regimes of rock specimen creep, i.e., primary, secondary, and tertiary creep regimes, are clearly visualized in the simulations. In the typical creep curves, it is generally accepted that the secondary creep is due to the stable crack propagation, and the tertiary creep is due to the unstable crack propagation [57].

Numerical Simulation of the Time-dependent Failure of Jointed Rock Slopes
In this section, the model proposed in section 3 is applied to simulate the time-dependent failure of jointed rock slopes.
In the stability analyses of rock slope, a continuous plane or series of interconnected planes is often predefined to determine the failure surface of slope. However, Terzaghi et al.
suggested that the setup of preexisting key discontinuity is more limited in natural rock slopes, and that a complex interaction between natural discontinuities and brittle crack development through intact rock bridges in slope bodies is required to form an evident sliding surface and then induces the macroscopic failure of the slope [17]. In the following sections, several demonstrative examples and a slope case will be numerically simulated to explore the characteristics of time-dependent failure in natural rock slopes.

Numerical Simulation of the Progressive Failure of Jointed
Slopes. As a demonstrative example, a rock slope model, including a pair of underlapping joints with an inclination angle of 45°, as shown in Figure 9, is established. e model geometry is 75 m × 50 m and has been discretised into 60,000 elements. In numerical modelling, the mesh dependency is unavoidable. Many previous studies detailed the aspect related to mesh effect, which showed that if the mesh size is too small, the computational resource will be wasted; in contrast, if the mesh size is too large, an evident calculation error will occur [58,59]. erefore, in this study, we just use an optimal element size to address the local damage evolution and macroscopic failure process of rock slopes. It is expected that this mesh size can meet the requirements of the calculation accuracy. Similar to the equivalent continuous medium model, joint planes can be defined as a set of planes individually. Such sets also can be combined to create a multiple-intersecting, jointed, or blocky system [1,60]. In the modelling with RFPA, the joints are assumed to be preexisting fractures. However, the fractures are not vacant but fully filled with a weak material. Referring to Figure 3, similar to the intact rock matrix, the joint is assumed to be isotropic and elastobrittle with a residual strength. e joint is still represented as the four-node isoparametric element same as those of rock matrix, but it is assumed to have a lower strength and stiffness compared to the rock matrix. Because it is hard to evaluate the accurate strength data for a joint, the mechanical properties of joints are selected referring to the previous studies [61]. e properties of joints are specified according to Weibull distribution. In this respect, only two distinct materials are included in the slope rock mass: joints and matrix.
Gravitational instability of slope is resulted from the interaction of different processes. For an existing slope, besides the dynamic action of earthquake, the most important factors causing the slope failure seem to be the alteration caused by climatic factors, weathering, and fluid circulation within the slope [42,62,63]. As shown in Figure 10, the weathering caused by freezing-thawing cyclic load and seepage flow can evidently soften the slope rock mass. In this study, we impose the following assumption on the demonstration example: the initially stable and inhomogeneous slope is subject to a progressive degradation of effective strength because of rock weathering or other alterations. Ultimately, this slope should undergo gravitydriven damage and macrofracturing increasing with time. Although it is relatively simplistic, the aim of this study is to simulate the slope damage development and investigate the  Figure 2 and other previous studies [35,[64][65][66][67][68], the mechanical properties of the demonstrative slope are chosen and listed in Table 1. Figure 11 shows the strength variation as a function of time in this simulation. In the first 10 months, the degradation of strength is much more rapid. After 10 months, the decrease gradually reaches a constant value in this demonstrative example. Similar to previous studies [2,19,21,42], during simulation, the strength of the slope is incrementally reduced throughout the whole example. Based on the strength degradation law, the spatial and temporal development of the fracture development and linkage are numerically predicted by RFPA, as shown in Figure 12. Figure 12(a) illustrates the fracture development in detail, and the grey scale indicates the relative magnitude of stress. e elements with a lighter shade of grey have relatively higher stresses. e toe of the slope and the rock bridge region show remarkable stress concentrations. Scaled principal stress vectors indicate that the maximum principal     damage first initiates at the upper end of the joint J1. With the strength degradation, joints J1 and J2 are fully fractured and the wing cracks begin to develop at the end of the joints. e orientation of the wing cracks is similar to those theoretically predicted for isolated underlapping joint pairs [2,60].
However, the actual situation is somewhat less symmetric due to the local heterogeneity in rock matrix and the nonuniform stress distribution around the rock bridge. As a result of the local heterogeneity, another two microcracks appear near the rock bridge. All these cracks deterministically select the path of least resistance through the rock bridge, and the random locations of the local inhomogeneity result in an irregular crack trajectory. e variation in the cracking pattern is highly sensitive to the interaction between the isolated cracks and the disorder features of the rock matrix. e crack path is therefore flexural, and the crack surface is rough. In practice, both high-stress and lowstrength failures can occur. e failure often initiates at a high-stress site in homogeneous materials, whereas the failure often initiates at the weak locations, e.g., microfracture, pore, and grain boundary [69]. Local heterogeneity is the main drive of failure at locations with low stress level.
If one carefully observes the rock bridge destruction in local scale, it can be found that the macroscopic shear fractures near the rock bridge are formed through the interaction and coalescence of many microtensile cracks. Due to the local heterogeneity and the stress mode near the rock bridge, microtensile cracks first come into being. With the increase in the number of microtensile cracks, the coalescence of these microtensile cracks controls the final shear fracturing mode. e elastic stress around microtensile cracks promotes a mutual interaction that produces brittle shear planes oriented obliquely to the remote principal stresses. In other words, the "shear" failure surface only generates once tensile fracture has developed to the point where significant cohesion loss has occurred along the path of coalescing fractures and larger displacements emerge through kinematic feasibility and slide mass mobilization.
is failure mode is consistent with the fracturing mode in the previous studies [2,70,71]. Figure 12(b) shows the associated development of the AE events during the fracturing process. AE events are stepwise simulated, and the intensity of the AE events is proportional to the degree of damage (white circles indicate shear failures, and red circles indicate tension failures). e fracture first initiates from the joints, where the element strength is relatively lower and then forward along the joints. During the initial stage, because of the heterogeneity of rock, most of the micro-AE events are discretely located in weak elements of joints. With strength degradation, the AE event activity within the central part of the rock bridge gradually increases.
e AE events indicate that the failure mode of the slope is a combination of tension and shear fracturing, while the tension failure is prevailing. Based on the AE activity in Figure 12(b), one can clearly find the development of microtensional cracks and the formation of the macro-"shear" surface or band in the rock bridge. Simulated horizontal displacement variation at the slope top and associated AE counts during slope failure is shown in Figure 13. It is shown that the isolated damage in early stages cannot induce significant deformation of slope at initial stages, while with the coalescence of isolated damage in a long-term period, a catastrophic failure occurs and the corresponding evident deformation comes into being.
A combination of environmental factors, such as seepage flow, temperature, earthquake, cyclic loading, alternate wetting, and drying, affects the strength of the rock mass and associated strength degradation rate. To investigate the influence of strength degradation rate on the long-term behaviour of slope, another six computational cases are employed in this section. e geometry configuration of these cases is same as that in Figure 12. e coefficient of strength degradation a s for the six cases is selected to be 0.005, 0.01, 0.05, 0.08, 0.1, and 0.2, respectively. e other mechanical parameters are same as those listed in Table 1. Figure 14 is time variation of catastrophic slope failure as a function of the coefficient of strength degradation a s . It is shown that if the slope is in a bad environment, i.e., the coefficient of strength degradation is relatively higher (for the case a s > 0.1 in this particular example), a catastrophic slope failure may occur within a short period of time. In contrast, if the slope surface is protected by effective protective measures and rock mass strength of the slope is less susceptible to erosion attenuation (for the case a s < 0.01), the slope will keep stable in a long-term period.
For the cases closer to the actual situation, another two demonstrative cases, as shown in Figure 15, were considered. e mechanical parameters are the same as those listed in Table 1. To obviate the influence of the randomly distributed joints on the failure pattern, four additional examples are simulated numerically while holding all of the rock parameters and boundary conditions constant (at the values used in case shown in Figure 9). e numerical results are listed in Figures 16 and 17.
Although the time-to-failure of these slopes is different, the general failure surface for all the cases appears to be macroscopic shear mobilization. e macroscopic failure mode is closely related to the local tensile damage, the  (3) integral sliding mobilization along the key discontinuity planes. e whole process could be called as "step-path" failures [5,72]. In most cases, besides a major failure surface that can be  Figure 18: Wangjiang road slope located in Yunyang city (the distribution of unloading fractures is generally characterized by a long strip; the unloading fracture strip is slightly wider near the top of the slopes, while it is relatively narrow near the toe of the slopes; the average width of the unloading fracture strip is about 10 m; the inclination of unloading fractures is generally parallel to the slope surface; these unloading fractures belong to the weak structural planes filled with plastic clay and weathered debris): (a) high slope containing alternating soft and hard rock layers (mudstone and sandstone); (b) the outcrop of cross section near south part of slope.

Surface scaling and collapse
Original slope surface observed, the secondary potential failure surface can be found which has developed from the top of the slope, resulted from internal block displacements. After the first sliding block formed and collapsed along the first failure surface, the second failure surface did not form completely because the degradation process ceased after the first upper block formed. If the slope is further affected by additional stress and strength attenuation, this secondary failure surface may also induce a second macroscopic shear mobilization. e numerical results above clearly indicate the failure process of the slope with discontinuity sets, revealing that the long-term stability of the slope is influenced not only by the quality and quantity of joints but also by the properties of rock bridges. Most of the critical sliding surfaces are formed along the joints by cutting through rock bridges and intact blocks.

Numerical Simulation of the Time-Dependent Failure in a Slope
Case. Yunyang city is located in a subtropical monsoon climate zone in the southwestern China, where the annual average temperature is 18.7°C (extreme minimum temperature is − 4°C, and the highest temperature is 42°C) and the average annual rainfall is 1145.1 mm (the largest annual rainfall is 1752.6 mm, and the minimum annual rainfall is 740.0 mm) [66]. Many high and steep rock slopes situated in Yunyang city were formed with the construction of highway. Sandstone and mudstone are the main lithologies, as shown in the Wangjiang road slope surface in Figure 18(a). Mudstone is likely to be easily weathered due to the water and temperature influence. Another main problem is the existing natural discontinuity sets such as weathered and unloading fractures. Figure 18(b) shows the other end of this slope. Support and protected measures were not implemented in this slope. About two years later, gradual surface scaling, sliding collapse, and block falling occurred in this slope. Figure 19 is a sketch of the prevailing failure phenomena for one section of this slope.
According to the sketch of one typical geological cross section in Figure 20, a two-dimensional planar model was used in the numerical study. Its geometry was 70 m × 45 m and was discretised into 50,400 elements. e rock properties of the numerical model are listed in Table 2.
Varying the coefficient of strength degradation a s , a total of five examples are modelled numerically while holding all other rock parameters and boundary conditions constant. e relationship between catastrophic failure time and coefficient of strength degradation is shown in Figure 21(a).  For the case a s � 0.02, the slope stability can last 26 months, which is consistent with the actual situation in situ observed in Yunyang city. e simulated horizontal displacement variation at the slope top and associated AE activity during the slope failure is shown in Figure 21(b). Figure 22 numerically illustrated the time-dependent failure process of the slope case. Significant stress concentration occurs near region of toe of the slope and the inner rock bridges. Local failures near the surface of slope show the predominant performance at the initial stage. With the degradation of strength, slope deformation develops deeper, and severe slope failure is expected if no support measures are carried out. Particularly, the destruction of mudstone layers is evident, as shown in Figure 22(a) C. Along with the destruction of mudstone layer, the overlying sandstone layer shows an evident tendency of extraversion and falling. e numerically obtained failure modes agree with failure phenomena observed in field, as shown in Figures 18 and 19. During the failure process, different types of rock bridges are observed, which are similar to those in Figures 12, 16, and 17. For example, at Locations A and B, the structure is similar to that for underlapping joints in a contractional geometry. At Locations C and D, the linking structure is similar to that for an overlapping pair of joints in an extensional orientation.
Combining numerical simulation and geological investigation in field, it is found that the following factors are the main factors contributing to the time-dependent instability of this slope case. (1) e presence of unloading fractures: these weak structural planes first reduce the integrity of the slope, and secondly, the strength of the weak structural planes is highly vulnerable to environmental erosion. As a result, the long-term stability of the slope is reduced. (2) e existing of mudstone layer: subjected to the erosion effect in the external environment, the underlying mudstone is easily weathered and its strength is susceptible to be degraded. Particularly under conditions of long-term rainfall, the mudstone is not easy to be drained in a timely manner. As a result, as shown in Figure 18(a), local collapses such as peeling or surface scaling first occur in the mudstone. en, in a longer period of time, a gravity-driven dumping or slider may occur in the overlying rock mass of sandstone.

Conclusions
In this study, a numerical model was developed to study the time-dependent deformation of rock and the associated long-term stability of slopes. In the model, the time-dependent behaviour of rocks is considered as a macroscopic consequence of damage development of microstructure. Comparing the experimental results with numerical simulations, it is found that the proposed model is suitable to investigate the rheological behaviour and time-dependent response of rocks.
Using the proposed model, a series of demonstration slope cases containing preexisting joints are numerically investigated. It is shown that the proposed approach is proven to reproduce the time-dependent failure mechanisms of a jointed rock slope, i.e., local sliding along natural joints, stress concentration around the joint tips, local damage along natural joints, local cutting of rock bridge, and the subsequent damage initiation and propagation in the intact rock. is progressive failure explains the time delay between the initial localized damage and the entire slope collapse. Besides the time-dependent strength degradation of both natural joints and intact rock mass, the distribution of preexisting joint has a significant influence on step-path generation. Excess shearing stresses along a discontinuity cause the development of orthogonal tensile fractures until eventually cross-over fractures form a continuous, but stepped, failure surface and overall failure occurs. e key block is released through the development of cross-over fractures linking the overlapping joints. e phenomenological approach provides supplementary information on the stress distribution and shows the development of the fracture and the interaction of the fractures within the rock bridge in detail that are impossible to be observed in field.
Numerical results for a rock slope case in Yunyang city show that gradual surface scaling, sliding collapses, and block falling are the prevailing failure phenomena of the slope in a long-term period. e degradation of rock properties, hard sandstone interbedded with soft mudstone, and regional structural joints with stress-release fracturing are the predominated factors influencing the long-term stability of the slope case. Slope deformation develops deeper, and severe slope failure is expected if no remedies are carried out. e general fractured pattern and the qualitative description of the failure process, as well as the instability time modelled with the proposed approach, agree with the field observations. Differing from previously adopted macroscopic approaches, the proposed model shows that the complex macroscopic time-dependent behaviour of heterogeneous brittle rocks can result from the small-scale interaction of elements and rock degradation. However, based on the numerical results, many factors associated with the modelling work on the long-term stability of slopes still need be considered in the future. For example, (a) a more realistic degradation law of rock strength need to be further investigated and employed in numerical models because different factors such as stress redistribution, near-surface physical, chemical, and biological weathering, or cyclic loading (e.g., thermal-or hydromechanical) may induce distinct degradation mode and rate of strength, which predominantly determines the long-term response of slopes; (b) to more accurately capture the fracturing pattern, the simulation of failure surface in slope in 3D is needed. Particularly, as an attempt to include more complicated physics, 3D modelling should be encouraged; (c) the characterization of joints in the model should be more clearly presented. Although it is hard for geological engineers to accurately provide these parameters particularly in the case with a large scale in 3D, it is important for the numerical prediction of associated stability of rock slopes.

Data Availability
We claim that the data used to support the findings of this manuscript are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.