## 3 Methodology
### Numerical Evaluation
#### 3.2.3 Numerical Evaluation of the Slope Stability  
In the present study, FOS was used to investigate the influence of different factors on slope stability. The method used to model the weak interlayers significantly influenced the slope stability analysis, and the most widely used were Weak Zone (WZ) [43; 44; 45], Ubiquitous-Joint (UJ) [46; 47; 48], and Interface (IF) [14; 49]. In a recent study [50], [PERSON] et al. validated each method and compared the sensitivity of each method to modeling parameters. Their study revealed that the method used to model weak interlayers (WZ, UJ, and IF) significantly influenced the slope stability analysis. It is still challenging to determine the most suitable modeling method, so the study conducted numerical evaluations using all these methods:  
1. The WZ method considers the weak interlayer as a 'weak zone' with isotropic materials and employs the Mohr-Coulomb Model [41] for modeling weak interlayers. This method treats potential sliding surfaces as isotropic materials with lower mechanical parameters in all directions, making them more conservative than other methods.  
Figure 6: (**a**) Model dimensions and area of layered excavation. (**b**) Model 1 for Block 8 (with anchorage structure and weak interlayer). (**c**) Model 2 for Block 9 (with anchorage structure and weak interlayer).  
2. The UJ method regards the weak interlayer as a 'weak zone' with anisotropic materials and uses the built-in Ubiquitous-joint constitutive model for modeling weak interlayers. This model accounts for the presence of an orientation of weakness (weak plane) in a Mohr-Coulomb model. The criterion for failure on the plane, whose orientation is given, consists of a composite Mohr-Coulomb envelope with tension cutoff. The position of a stress point on the latter envelope is controlled again by a non-associated flow rule for shear failure and an associated rule for tension failure [41]. This method maintains anisotropy, with lower mechanical properties along the bedding plane and consistent properties with other rock masses in other directions.
3. The IF method considers the potential sliding surface as a two-dimensional plane unit and takes into account properties such as friction angle and cohesion. This approach does not involve the concept of interface thickness, making it closer to real-world conditions. However, the values for shear stiffness and normal stiffness of the interface need to be determined artificially. These two parameters were estimated using the recommended method in the FLAC3D manual [41]: \[k_{n}=\frac{EE_{r}}{s(E_{r}-E)}\] (1) \[k_{S}=\frac{GG_{r}}{s(G_{r}-G)}\] (2) \[G_{r}=\frac{G_{r}}{2(1+\mu_{r})}\] (3) \[G=\frac{E}{2(1+\mu)}\] (4)  
where \(k_{n}\) is the normal stiffness, \(k_{s}\) is the shear stiffness, \(s\) is the joint spacing (\(s=0.5\) m), \(E_{r}\) is the [PERSON]'s modulus of the intact rock, \(E\) is the [PERSON]'s modulus of the rock mass, \(G_{r}\) is the shear modulus of the intact rock, \(G\) is the shear modulus of the rock mass, \(\mu_{r}\) is the [PERSON]'s ratio of the intact rock, \(\mu\) is the [PERSON]'s ratio of the rock mass (\(\mu_{r}=\mu=0.3\)). This estimation method requires [PERSON]'s modulus and shear modulus of the rock mass. However, [PERSON]'s modulus of the rock mass is difficult to determine in practical applications. It is assumed that [PERSON]'s modulus of the rock mass is consistent with the bedding plane. Although this consideration lacks experimental verification, it is conservative. The normal stiffness \(K_{n}\) and shear stiffness \(K_{\mathrm{s}}\) of the interface are \(1.09\times 10^{10}\) Pa/m (K\({}_{\mathrm{n}}\)) and \(4.23\times 10^{9}\) Pa/m (\(K_{\mathrm{s}}\)) in the present study.