## 3 Numerical Modeling
### Parameters and Numerical Model Setup  
This section employed FLAC3D, a finite-difference-based method, to construct a numerical model for simulating the gypsum mining room excavation. Due to the excessive scope of the gypsum mining area, it is difficult to build a complete large-scale numerical model for the stability analysis, but since the mining area has symmetrical mine rooms, the model was simplified as depicted in Figure 5. It has overall dimensions of 10 m in width and 73 m in height. Given the typical ground conditions of the abandoned mining area, the maximum length of the mining room is about 150 m, which is much longer than the room width, and the model was designed as a planar strain model. The model's bottom and sides were subjected to simple support constraints, while the top was treated as a free boundary. Due to the symmetric nature of one mining room, the simplified numerical model included a mining room sidewall of 3 m \(\times\) 3 m and a mining room of 17 m \(\times\) 3 m. Only the medium-weathered rock layer and pre-existing gypsum rock layer were incorporated into the numerical model, as the gravel and strongly weathered rock layers were confined to the shallow layers of the ground and did not significantly impact the abandoned mining area's stability.  
The mining process inevitably disrupts the initial stress equilibrium within the original rock mass, leading to a redistribution of stress. As stress is rebalanced, the weight of the rock mass above the mining room shifts to the stable adjacent regions, resulting in a stress  
Figure 5: Numerical model setup for current stability.  
concentration along the residual sidewalls. Under the assumption of considering only the self-weight load, the failure of the abandoned mining area is expected to initiate from the sidewalls when the stress exceeds the rock mass's strength limit. To track this, a monitoring point has been established at the top of the sidewall to record vertical stress and displacement, as shown in Figure 5.  
The GSI, mi, and UCS parameters play a crucial role in shaping the peak Hoek-Brown envelope, and these relationships are defined through the equations established by [PERSON] et al. (2002). A significant element of the IMASS model is the residual envelopes and the phenomenon of rock mass softening that occurs after reaching the peak envelope. This phenomenon is extensively elaborated upon in the introduction of \"cohesion weakening frictional Strengthening behavior [20]\". The properties listed in Table 1 represent the essential input properties necessary for the implementation of the IMASS constitutive model. These properties were carefully selected based on geological survey data to ensure the model's accuracy and reliability.  
The critical plastic shear strain for the IMASS model, as indicated in Table 1, employs a multiplier to compute the critical strain, as visually demonstrated in Figure 1. By adjusting this multiplier, the material behavior can be modified to exhibit increased brittleness if the multiplier is set to a value less than 1.0, or heightened ductility if its value is set to be greater than 1.0. This multiplier factor serves to adjust the post-peak mechanical properties of rocks, as rock mass undergoes progressive damage during the loading process. After failure, the rock mass does not completely lose its bearing capacity under certain stress conditions. Therefore, the post-peak mechanical behavior still influences the overall bearing capacity of the model. Specifically, this is manifested by the speed at which the stress-strain curve drops to its residual strength after the peak strength, representing brittleness and plasticity.