## 2 Numerical Model of the Tunnel and the Support Schemes  
### Numerical Model  
The dimensions of the model established according to the actual stratigraphic information were 30 m \(\times\) 140 m \(\times\) 30 m (Figure 2). The length of the circular tunnel with a radius of 2.15 m was 120 m. The depth of the tunnel was assumed to be 1095 m. In order to eliminate the influence of the opening on the boundary on stability, the tunnel was excavated inside the model. The model was set up in 3 DEC (Itasca), which is a discrete element method software and addresses two points: (1) allows for finite displacements and rotations of discrete bodies, including complete detachment, and (2) recognizes new contacts automatically as the calculation progresses [25]. The rock mass was assumed to be homogeneous and isotropic, and was arranged so to meet the strain-softening constitutive model suing the Mohr-Coulomb criterion. The estimation of the rock mass properties has been investigated by many researchers. [PERSON] and [PERSON] used a GSI system to estimate the rock mass property values [26]. Based on the laboratory tests, [PERSON] et al. assumed that the mean rock mass modulus was 51% of that of the intact rock, and the rock mass Poisson's ratio was 21% higher than that of the intact rock [27]. The internal friction angle of the rock mass was usually taken to be 0.83 to 0.9 times that of the intact rock ([PERSON] et al. [18]; [PERSON] et al. [28]). [PERSON] et al. used an empirical equation to estimate the rock mass cohesion [29]. According to the above research, the calculation parameters of each lithology were determined and are listed in Table 1.  
The bulk modulus (\(K\)) and shear modulus (\(G\)) given in Table 1 were calculated from deformation modulus (\(E\)) and Poisson's ratio (\(v\)) of the rock mass using Equations (1) and (2).  
\[K=\frac{E}{3(1-2v)} \tag{1}\]  
\[G=\frac{E}{2(1+v)} \tag{2}\]  
In addition, because of the lack of laboratory test data about the strain-softening behavior, some research ([PERSON] et al. [30]; [PERSON] et al. [31]) works were reviewed to estimate the strain-softening property values, as given in Table 2. The top boundary of the model was located at a depth of 1080 m. A stress boundary of vertical stress of 27.54 MPa was applied based on the average weight of 2550 kg/m\({}^{3}\) of the overburden strata. The roller boundary condition (no vertical or displacement) was applied at the bottom. For lateral boundaries, the stress boundary conditions with consideration of the gravity gradient were applied on the corresponding sides in both x- and y-directions. The gravity gradient was also applied to the entire model to simulate real in situ stress. The horizontal stress was calculated by a lateral stress ratio of \(K_{x}=1.35(\sigma_{x}/\sigma_{z})\) and \(K_{y}=0.9(\sigma_{y}/\sigma_{z})\). In this study, the initial stress field was obtained by equilibrating the in situ stress with the applied boundary conditions.  
\begin{table}
\begin{tabular}{c c c c} \hline
**Xtrain-Softening Parameter** & **Fine Sandstone** & **Medium Sandstone** & **Silstone** \\ \hline Residual friction angle, \(\sigma_{r}\), (\({}^{\circ}\)) & 29.38 & 24.91 & 27.79 \\ Residual cohesion, \(C_{r}\), (MPa) & 0.02 & 0.18 & 0.23 \\ Residual tensile strength, \(\sigma_{tr}\), (MPa) & 0.03 & 0.27 & 0.29 \\ Plastic shear strain, \(e^{pts}\) & \(4\times 10^{-3}\) & \(4\times 10^{-3}\) & \(4\times 10^{-3}\) \\ Plastic tensile strain, \(e^{pt}\) & \(6\times 10^{-3}\) & \(6\times 10^{-3}\) & \(6\times 10^{-3}\) \\ \hline \end{tabular}
\end{table}
Table 2: Strain-softening property values of the rock masses.  
Figure 2: 3D numerical model: (**a**) whole model and (**b**) profile view.