## 2 Numerical Modeling
### Ground Motion Processing  
In the dynamic response analysis of tunnel engineering, selecting seismic motions that align with the seismic environment of the site is critically important. With advancements in technology, artificial ground motions are increasingly being utilized in seismic analyses of tunnel engineering. The generation method for artificial ground motions is described below.  
Using the trigonometric series method, the unsteady seismic acceleration is considered the product of the smooth random process and the outer envelope function considering the unsteady characteristic, as shown in Equation (4):  
\[a(t)=f(t)a_{s}(t) \tag{4}\]  
where \(a(t)\) represents the seismic acceleration time course, \(a_{s}(t)\) is a stationary random process with a (one-sided) power spectral density function and zero mean, and \(f(t)\) is the envelope function.  
The expression for \(f(t)\) is given as follows:  
Figure 4: Mechanical behavior of the contact surface. (**a**) Separation state. (**b**) Contact state. (**c**) Tangential mechanical behavior.  
\[f(t)=\left\{\begin{array}{ll}(t/t_{1})^{2}&t<t_{1}\\ 1&t_{1}\leqslant t<t_{2}\\ e^{-c(t-t_{2})}&t_{2}\leqslant t<t_{3}\\ 0&t_{3}\leqslant t<T\end{array}\right. \tag{5}\]  
where \(t_{1}\), \(t_{2}\), \(t_{3}\), and \(T\) represent the start and end times of the stationary phase, the end time of the attenuation phase, and the total duration of the seismic wave, respectively. \(c\) is a constant controlling the attenuation rate, with a value range of 0.1-1.0.  
The expression for \(a_{s}(t)\) is given as follows:  
\[a_{s}(t)=\sum_{k=1}^{n}C_{k}\cos(\omega_{k}t+\varphi_{k}) \tag{6}\]  
where \(\varphi_{k}\) is a random phase angle uniformly distributed in \((0,2\pi)\) and \(\omega_{k}\) and \(C_{k}\) are the frequency and amplitude of the \(k\)-th spectral component, respectively. \(\omega_{k}\) and \(C_{k}\) are determined based on the generated power spectral density function using the following equation:  
\[\left\{\begin{array}{ll}C_{k}=\sqrt{4S(W_{k})\Delta_{W}}\\ \Delta_{w}=2\pi/T\\ W_{k}=k\Delta_{W}\end{array}\right. \tag{7}\]  
where \(S(\omega_{k})\) represents the given power spectral density function.  
As described above, generating artificial ground motions requires a specified power spectral density function. During the generation process, the power spectral density function is first derived from the design response spectrum specified in China's Code for the Seismic Design of Buildings (GB 50011-2010) [31]. This power spectral density function is then substituted into the trigonometric series cosine model to obtain the stationary random process, which is subsequently multiplied by the envelope function to generate the seismic acceleration time history. In this study, a MATLAB program was developed to implement the above process, generating artificial ground motions suitable for the seismic environment of the study site.  
Figure 5a shows the acceleration time history of the artificial ground motion. The ground motion duration is 40 s, with a time interval of 0.02 s, and the peak acceleration is 0.2 g. Figure 5b illustrates the acceleration response spectrum of the artificial ground motion. As shown in the figure, the acceleration response spectrum of the artificial ground motion closely matches the designed acceleration response spectrum for the site, indicating that the generated artificial ground motion effectively reflects the seismic environment of the tunnel site.  
Figure 5: Input artificial ground motion (PGA = 0.2 g). (**a**) Acceleration time history. (**b**) Acceleration response spectrum.