## 2 Numerical Model
### Mesh Sensitivity Analysis  
The finite element mesh size significantly affects the outcomes of the finite element analysis. It may affect the wave propagation and the tunnel's failure mode. As a result, a mesh convergence analysis was carried out to accurately identify the mesh size in the numerical model.  
[PERSON] [37] found that employing 6-12 elements per wavelength may reduce wave distortion. In this analysis, the longitudinal wave velocity was 4000 m/s, and the major vibration frequency ranged from 10 Hz to 100 Hz. The model predicted a wavelength of 40 m to 400 m. Thus, the study area's grid must have been smaller than 3.3 m to ensure reliable wave propagation. Mesh sizes from 7.5 to 50 cm were used for the convergence test. Analyses were performed with an 8 m/s load with a 100 Hz frequency. Figure 6 shows velocity time history at the distance of 19.4 m for different meshes. Convergence of the solution occurred at a mesh size of 10 cm.  
The failure mode of the concrete tunnel is also dependent on the mesh size. [PERSON] et al. [20] recommended that the mesh size be between 1/12 and 1/8 of the thickness of the slab to achieve the optimal balance between computational efficiency and accuracy. In this study, the lining thickness was 0.3 m. The recommended size was 25-37.5 mm. In this study, we used a 25 to 100 mm mesh size. Figure 7 shows that a large mesh develops less damage while a small mesh size creates more damage. [PERSON] et al. [38] recommended a 25 mm mesh size. Hence, 25 mm was used in this simulation. Further mesh refinement may result in more accurate results, albeit at the cost of increased processing time. The model had meshed into 2,299,420 solid elements, of which 1,552,060 were rock elements, 747,360 were tunnel elements, and 29,134 were beam elements, consisting of 2,328,554 elements in total.  
Figure 6: Effect of mesh size on response of wave propagation.  
where \(V_{p}\) is the wave velocity of the media, \(\xi\) is the small strain damping ratio, \(e\) is the exponential number, and \(f_{p}\) is the predominant frequency. As suggested by [PERSON] et al. (2018), the product of geometrical and material damping may be used to generate the theoretical attenuation equation that accounts for both geometrical and material damping.  
\[V_{R}=V_{r}\Big{(}\frac{r}{R}\Big{)}^{n}e^{-\frac{2\pi f_{p}\xi}{V_{p}}(R-r)} \tag{4}\]  
A separate numerical model for far-field wave propagation was established. It was verified with the theoretical solution shown in Equation (4). The dimensions of the numerical model were \(110\ \times 30\ \times 1\) m, and mesh size was considered as 0.1 m. The blast hole radius was 0.5 m, and the charge length was 1 m. A blast load of 1 m/s with 43 Hz frequency was applied. The simulation was conducted in damped and undamped mediums, and 5% damping was selected for this simulation. The rock characteristics were defined as follows: density = 2500 kg/m\({}^{3}\), Poisson's ratio = 0.25, Young's modulus = 33.3 GPa, for \(V_{p}\) = 4000 m/s. Figure 8a,b show the comparison of numerical and far-field theoretical solutions with damped and undamped media. The far-field solution predicts attenuation at the surface of the rock, and it is reduced by a factor of 2 to predict within the profile attenuation curve. The numerical model predicts PPV well and fits with far-field theoretical solution.