## 1 Introduction  
Slope failures are natural phenomena induced either by (1) natural factors, such as rainfall, changes in the groundwater level [1], loss of vegetation, and earthquakes, or by (2) human activities, such as additional loading and construction work at the top of slopes, excavation at the toe, and underground mining [2]. The methods of underground mining include the room and pillar mining method, caving method, and joint mining method. As one of the caving methods, sublevel caving using the filling of post-mining space to control geo-stress is a highly efficient, commonly used, top-down, underground bulk mining method. During underground mining, such as sublevel caving mining, progressive slope failure may be triggered. There are a variety of landslide classifications. Some types of slopes remain dormant after failure, while others show progressive or delayed failure in response to a disturbance. The progressive slope failure induced by sublevel caving mining will occur many times with the progress of the excavation. The slope experiencesa cyclic process that cycles between a stable-unstable-stable-unstable condition; that is, intermittent instability. To reduce the damage of slope failures on infrastructure and human life, numerical modeling methods are crucial for predicting and analyzing slope failure. These numerical methods typically fall into two categories: (1) continuum-based methods, such as the finite element method (FEM) and finite difference method (FDM), and (2) discontinuum-based methods, such as the material point method (MPM) and discrete element method (DEM). The DEM is a powerful numerical modeling technique used to simulate discontinuous rock mass movement and failure [3,4]. [PERSON] and [PERSON] [5] used the DEM for block systems to study the failure and characteristics of jointed rock slopes. The MPM has also been widely used to analyze and simulate the deformation and failure of slopes [6, 7, 8]. However, the DEM is time-consuming and the performance of DEM is usually limited by computing capacity. When analyzing the deformation and failure of a large-scale slope, the DEM is quite computationally expensive. The MPM discretizes the continuum into a group of Lagrange particles with information, such as mass, velocity, position, and stress, and the processes of solving momentum equations and computing spatial derivatives are conducted on the Euler background mesh. When the particles pass through the background mesh, the gradient of the interpolation function will not be close to zero, resulting in a strong disturbance of the particle stress and cell-crossing noise. Therefore, for the modeling of large-scale rock slopes, the MPM is numerically unstable due to a large number of particles passing through the background mesh. The FDM and FEM are also commonly used to analyze slope failure. Compared with the DEM, they are more computationally efficient. However, both the FEM and FDM are mesh-dependent. When meshes undergo a large deformation, mesh-dependent numerical simulation methods experience problems, such as a distorted mesh, making the calculation impossible [9]. To address this problem, [PERSON] et al. [10] and [PERSON] et al. [11] used a hybrid FEM-DEM approach to analyze and evaluate the stability of a slope in open-pit and underground mining. [PERSON] and [PERSON] [2] used the h-adaptive mesh refinement technique to simulate the deformation and damage of a slope by frequently updating the slope geometry. [PERSON] et al. [12] combined the Hoek-Brown strength reduction technique and the FEM to evaluate and analyze the slope stability. On the other hand, the limit equilibrium method and strength reduction method (SRM) [13] are also widely used to evaluate the slope stability and determine the slip surface. However, for the progressive slope failure induced by sublevel caving mining, this strategy needs to reconstruct the model repeatedly according to the terrain conditions after the landslide, which is inefficient.