## 4 Computational Analysis and Discussion
### Factor of Safety
The results from the diagrams show that at the same burial depth of 80 m, the factor of safety calculated by [PERSON]'s method is significantly lower than that calculated in this paper, particularly evident in silo-type caverns, which corroborates analysis. Interestingly, the factor of safety for a 60 m tunnel depth calculated by [PERSON]'s method is very close to that for an 80 m depth calculated in this paper, suggesting that there is about a 20 m discrepancy between the two methods when calculating tunnel depths. Overall, the factor of safety decreases with increasing internal pressure, and this rate of decrease also slows down with increasing internal pressure. From the figures, it is evident that with an increase of 20 m in burial depth, the factor of safety for silo-type caverns significantly improves. This suggests that burial depth has a more significant impact on the stability of silo-type caverns, which is attributed to the substantial increase in the volume of the three-dimensional overlying conical rock mass as the burial depth increases. Assuming that the standard safety factor is set to 2.0 and the burial depth is 60 m, the ultimate pressures of tunnel type and silo type corresponding to the rock mass of Class I are roughly 30 MPa and 60 MPa, respectively; the ultimate pressures corresponding to the rock mass of Class II are 26 MPa and 49 MPa; and the ultimate pressures corresponding to the rock mass of Class III are 13 MPa and 30 MPa. It can be seen that the decrease in the ultimate pressure of the rock body of Class III is larger, which is due to the greatly  
Figure 13: Variationin potential fracture surfaces with mesh and computational step for silo-type caverns overlying rock masses under Class III rock mass (\(K=0.5\)): (**a**) the number of mesh elements; (**b**) the computational step length  
finding potential fracture surfaces into an initial value problem of ordinary differential equations, based on the limit stress field and the Mohr-Coulomb failure criterion. Following this, a mechanical analysis is conducted on the uplifted rock mass using the identified fracture surfaces, converting the equilibrium equations into a form of stress integration. By introducing interpolation equations containing unknown coefficients, the number of equations is made equal to the number of unknowns, allowing for a solution. This method is applicable to both two-dimensional and three-dimensional problems, corresponding to tunnel-type plane strain issues and silo-type three-dimensional issues, respectively. Using this method, potential fracture surfaces were solved for both tunnel-type and silo-type LRCs under three rock mass quality grades. The analysis indicates that the orientation of the fracture surfaces is significantly influenced by the horizontal geostress coefficient; the greater the horizontal geostress coefficient, the more the fracture surfaces tend to be horizontal. This implies that sites with a geostress coefficient greater than 1 should be prioritized when selecting locations for LRCs. The fracture surfaces are also affected by the internal friction angle, and this effect is related to the geostress coefficient. When the horizontal geostress coefficient is less than 1, a smaller internal friction angle causes the fracture surfaces to be more horizontal, while the opposite is true when the horizontal geostress coefficient is greater than 1. The impact of depth and cavern type on the fracture surfaces is relatively minor. When the horizontal geostress coefficient equals 1, the fracture surfaces of silo-type LRCs are nearly 45 degrees from the horizontal direction, while for tunnel-type LRCs, they are close to 52 degrees. Based on the analysis of the fracture surface and assuming that the geostress coefficient is 1, the factor of safety of two types of LRCs with the same volume is calculated and analyzed for different depths of burial and different rock quality grades.