## 4 Numerical Simulations
### Numerical Modeling  
An FLAC3D numerical model for the test-scale tunnel was developed. The dimensions of the numerical model were the same as the test-scale model, as shown in Figure 12. Some critical aspects of the establishment of a numerical model were as follows:  
(a) Meshing determines the computational efficiency and reliability of the simulation. Based on the literature [55], the size of the element (\(\Delta l\)) is:  
\[\Delta l\leq\frac{\lambda}{10},\lambda=\frac{V_{S}}{f_{max}} \tag{1}\]  
where \(\lambda\) is the wavelength of seismic wave propagation in the model and \(f_{max}\) is the maximum frequency of the input motions with a value of 40 Hz according to Figure 7. \(V_{S}\) is the shear wave velocity of the model mass surrounding the model test. \(V_{S}\) has a value of 30 m/s with a 2.4 m box height and 0.08 s time difference in peak acceleration response between the model bottom and top in the working case considering an El-Centro wave with z-direction input. Thus, the maximum size of an element in the numerical model in this paper was less than 0.075 cm. A total of 139,840 elements were applied to the tunnel numerical model. Figure 12a shows the established numerical model. (b) The Mohr-Coulomb Tension Crack model with an inner friction angle (\(\varphi\)) and cohesion (\(c\)) was selected to simulate the tunnel lining. The Mohr-Coulomb model was adopted for the surrounding rock mass and fault zone. The elastic model was adopted for the isolation layer and the mechanical parameters of the isolation layer of the composite material were determined by the trial-and-error method. An interface element was set  
Figure 12: Setup of the numerical model used in the current study: (**a**) dimensional division of the model, (**b**) sketch of the simulation loading process, (**c**) sketch map of the boundary in dynamic analysis, (**d**) schematic diagram of fault contact. to simulate the interface between the surrounding mass and the fault (Figure 12b). The normal stiffness (\(K_{n}\)) and shear stiffness (\(K_{s}\)) of the interface element were calculated by:  
\[K_{s}=K_{n}=10 max\left[\frac{\left(K+\frac{4}{5}G\right)}{\Delta Z_{min}}\right] \tag{2}\]  
where \(K\) is the bulk modulus, \(G\) is the shear modulus, and \(\Delta Z_{min}\) is the smallest width of an adjoining zone in the normal direction. The parameters of these above element types are tabulated in Table 6. The isolation layer has an elastic modulus with a range of 5 MPa to 20 MPa. (c) The bottom boundary of the numerical model was assumed as the rigid boundary [56] (Figure 12c). The lateral boundaries of the model are applied by the built-in absorbing boundary in FLAC3D, i.e., the free field (Figure 12c), expressed as  
\[F_{i}=-\rho C_{P}\left(v_{i}^{m}-v_{i}^{ff}\right)A+F_{i}^{ff}i=x,y,z \tag{3}\]  
where \(\rho\) is the density of the rock mass, \(C_{P}\) is the P-wave velocity at the model boundary, \(C_{S}\) is the S-wave velocity at the model boundary, \(A\) is the area of influence of a free-field grid point, \(v_{i}^{m}\) is the i-drection velocity of a grid point in the main grid at the model boundary, \(v_{i}^{ff}\) is the \(i\)-direction velocity of a grid point in the side free field, and \(F_{i}^{ff}\) is the free-field grid point force in the \(i\)-direction. (d) After the application of gravity, the vertical sliding support at the vertical sides, the fixed support at the bottom of the model and the free boundary at the top of the model, and the initial stress equilibrium of the numerical model were obtained. Then, creep-slip displacement in the z-direction (Figure 11d) was achieved by moving the hanging wall with a small normal velocity, which kept the whole creep-slip process quasi-static to avoid dynamic effects. Finally, the lateral and bottom surfaces were applied to the free-field and rigid boundaries, respectively. The subsequent seismic excitation was simulated by inputting the acceleration-time histories (Figure 7) at the bottom of the model.