## 1 Introduction  
The presence of joints in rocks is a common feature of rock masses. Such joints usually have a significant effect on the rock mass strength and deformation [1, 2]. The joint persistency is non-continuous and is one of the most important geometrical parameters in a rock mass. There are different structural characteristics of rock masses due to differences in joint persistency. The mechanical behavior of a rock mass with non-persistent joints is controlled by complex interactions of joints and intact rock bridges. The rock mass is cut into rock blocks by fully persistent joints, which have discontinuous medium characteristics. In contrast, shear and tensile failure always occur in rock bridges formed by non-persistent joints, which is a typical continuous medium characteristic. The bearing capacity of the non-persistent jointed rock masses mainly depends on the geological characteristics of the joints, such as their strength and the mechanical properties of rock bridges. Thus, estimating the bearing capacity of non-persistent jointed rock masses is one of the most challenging tasks in rock engineering. Many studies have been conducted to analyze the bearing capacity of jointed rock masses using various approaches. Limit equilibrium methods, which have a simple principle and reliable solutions [3, 4, 5, 6], have been widely applied to rock slope stability analysis of persistent joints. However, because of their unreasonable assumptions, these methods are only suitable for analyzing simple failure modes or homogeneous rock masses in which discontinuities play a minor role. With the development of computational mechanics and computational resources, many numerical methods have been proposed to investigate rock masses. Among these methods, the finite element method (FEM) [7, 8, 9] and finite difference method (FDM) [10] have been widely applied to jointed rock-bearing capacity analysis. However, these methods are not suitable for modeling the faults and joints in a rock mass. The distinct element method (DEM) or discontinuous deformation analysis (DDA) can capture the failure modesof jointed rock slopes with different arrangements of joints using discontinued formulations [11; 12; 13]. However, there are still difficulties in modeling a complex network of joints and the complex constitutive relations of rock materials and discontinuity surfaces. Block theory, which is a block identification method based on combinatorial topology concepts [14], was recently used as a new numerical method for considering a jointed rock mass. Block theory can efficiently solve the removability and stability problems of simple non-persistent jointed rock blocks under various engineering conditions [15; 16]. As such, block theory has been developed and widely applied in rock slope engineering [17]. Plastic limit analysis theory is a powerful tool for efficiently analyzing the bearing capacity in soil and rock mechanics with high accuracy. This method focuses on seeking analytical solutions, and as such, it can only treat classical geotechnical problems with simple structures [18; 19], as it is difficult to apply this method to build statically admissible stress fields and kinematically admissible velocity fields for large-scale and complex structures. Thus, a numerical approach to plastic limit analysis was developed. Plastic limit numerical analysis combines the lower bound theorem, the finite element discretization technique, and linear programming [20] and has been widely applied in soil mechanics [21; 22; 23], stability analysis of soil slopes [19; 24; 25; 26], the bearing capacity of jointed rocks [27; 28; 29], and stability analysis of rock slopes [15; 16; 30]. These studies demonstrate that plastic limit numerical analysis is a valid approach to investigating the bearing capacity of geotechnical structures. However, almost all studies have only considered the soil slope, soil foundation, and simple rock slope.