## 2 Numerical Simulation with FLAC3D
### Numerical Model
To account for the mechanical influence of the overlying backfill on the top surface of L1, it is simplified as a uniformly distributed load in the numerical model. The uniform load is set to \(q=1\) MPa, which corresponds to the self-weight of a 60-m-high backfill (calculated as \(pgh=1780\) kg/m\({}^{3}\)\(\times\) 10 m/s\({}^{2}\)\(\times\) 60 m = 1.068 MPa). In situ stress measurements have been conducted at the Maoping Lead-Zinc Mine and a linear regression equation describing the variation of ground stress with depth was derived:  
\[ \begin{array}{l}\sigma_{\rm hmax}=-40+0.04172D\\ \sigma_{\rm hmin}=-20+0.01876D\\ \sigma_{\rm v}=-30+0.028D\end{array}  \tag{2}\]  
where \(\sigma_{\rm hmax}\) and \(\sigma_{\rm hmin}\) represent the measured maximum and minimum horizontal stress, respectively; \(\sigma_{\rm v}\) is the vertical stress. The \(\sigma_{\rm hmax}\), \(\sigma_{\rm hmin}\), and \(\sigma_{\rm v}\) are employed to generate the initial values of the stress components \(\sigma_{\rm x}\), \(\sigma_{\rm y}\), and \(\sigma_{\rm z}\) in the numerical model. Specifically, the coefficients 0.04172, 0.01876, and 0.028 in Equation (2) are used to compute the gravity-induced stress gradients of the \(\sigma_{\rm x}\), \(\sigma_{\rm y}\), and \(\sigma_{\rm z}\), and the initial ground stress can form due to the gravity and the applied surcharge on the top boundary of the model. To examine the validity of the formed initial stress based on the in-site data, stress components are extracted along a vertical line connecting the centroids of the model's top and bottom surfaces. These values are then compared with the ground stress obtained from Equation (2). The average interpolation method in FLAC3D is used to extract the stress components and for the subsequent contours, which means the values of zone-based data are averaged to grid points and linearly interpolated to the zone interior [18]. Figure 3 compares the initial ground stress generated by the numerical model with the field-measured values. Their excellent agreement verifies the reliability of the constructed numerical model, the applied boundary conditions, and the formation of the initial stress. After the initial stress is formed in the numerical model, the following steps are used to simulate the process of mining and backfilling:  
(1) Excavate a 0.5-m-thick layer of the rock mass above the stope region and apply a vertical uniformly distributed load \(q=1\) MPa on the top surface of the orebody. Once the numerical model reaches equilibrium, reset the displacements and velocities of all nodes to simulate the displacement release of rock mass after the excavation. \begin{table}
\begin{tabular}{c c c c c c c} \hline \hline
**Material** & **Density** & **Posson’s Ratio** & **Friction Angle** & **Cohesion** & **Young’s Modulus** & **Tensile Strength** \\  & **(kg\(\cdot\)m\({}^{-3}\))** & **Posson’s Ratio** & **(\({}^{\circ}\))** & **(MPa)** & **(GPa)** & **(MPa)** \\ \hline Rock mass & 2800 & 0.3 & 35 & 1.98 & 10 & 1.25 \\ Backfill & 1780 & 0.27 & 33 & / & 0.315 & / \\ \hline \hline \end{tabular}
\end{table} 
Table 1: Parameters of the rock mass and backfill used in numerical simulations. (2) Excavate the primary drifts in one step at the L1 layer and solve the model to an equilibrium state. Subsequently, fill the primary drifts and solve the model again. Repeat the mining and filling process for the secondary drifts at L1. (3) Repeat Step 2 to complete the mining and filling process for the drifts at L2 and L3, respectively. Following the above procedures, numerical studies were conducted to determine the strength requirements of the backfill roof for different drift sizes: (1) 3.5 m \(\times\) 3.0 m (drift width \(\times\) drift height), (2) 4.5 m \(\times\) 4.0 m, and (3) 5.5 m \(\times\) 5.0 m. For each case, a three-dimensional numerical model was established in FLAC3D, and the strength requirements of the backfill roof were determined using the introduced iterative search method.