## 2 Methodology
### Governing Equations and Boundary Conditions  
The incident SH-wave motion in the half-space along the z-axis direction can be written in the local coordinate system (x\({}_{1}\)-y\({}_{1}\)-z\({}_{1}\)) as  
\[w_{1}^{i}(r_{1},\theta_{1})=\exp\{i[\omega t-kr_{1}\cos(\theta_{1}-\gamma)]\} \tag{1}\]  
where \(i\) is the imaginary unit and \(\omega\) is the circular frequency of the incident wave.  
The reflected wave will be generated when the incident wave propagates to the ground surface, and the scattered wave will be generated when it encounters the tunnels. Then, the total wavefield in the half-space is the superposition of incident wave, reflected wave, and scattered wave, all of which must satisfy the 2D wave equation.  
\[\frac{\partial^{2}w}{\partial r^{2}}+\frac{1}{r}\frac{\partial w}{\partial r}+ \frac{1}{r^{2}}\frac{\partial^{2}w}{\partial\theta^{2}}=\frac{1}{\beta^{2}} \frac{\partial^{2}w}{\partial t^{2}} \tag{2}\]  
To satisfy the zero-stress boundary condition of the half-space surface, the mirror image method [32] is used to solve the problem. As shown in Figure 2, assuming that there is another incident wave source and two tunnels with an outer diameter, \(a\), in the half-space with the surface as its axis of symmetry, the incident wave is  
\[w_{2}^{i}(r_{3},\theta_{3})=\exp\{i[\omega t-kr_{3}\cos(\theta_{3}-\gamma)]\} \tag{3}\]  
Figure 1: The analytical model of dynamic tunnel-soil-tunnel interaction subjected to out-of-plane SH-wave.  
The scattered wavefield corresponding to different tunnels can be expressed as  
\[w_{j}^{s}(r_{j},\theta_{j})=\sum\limits_{m=0}^{\infty}H_{m}^{(2)}(kr_{j})(A_{j,m }cosm\theta_{j}+B_{j,m}sim\theta_{j})\ j=1,2,3,4 \tag{4}\]  
Here, \(H_{m}^{(2)}(x)\) is the Hankel function of the second kind with argument x and order m. The time factor exp(iut) in Equations (1) and (3) is omitted in the subsequent formulas in this article. In the physical sense, they all represent cylindrical waves propagating outward from their respective centers, oi, satisfying wave Equation (2) and Sommerfeld radiation conditions. From the symmetry of the mirror method, it can be obtained that \(A_{1}\),\(m=A_{3}\),\(m\), \(B_{1}\),\(m=B_{3}\),\(m\), \(A_{2}\),\(m=A_{4}\),\(m\), and \(B_{2}\),\(m=B_{4}\),\(m\).  
The wavefield in the flexible tunnels is  
\[\begin{array}{ll}w_{j}^{f}(r_{j},\theta_{j})&=\sum\limits_{m=0}^{\infty} \left[C_{jm}^{(1)}H_{m}^{(1)}(k_{0}r_{j})+C_{jm}^{(2)}H_{m}^{(2)}(k_{0}r_{j}) \right]\cos m\theta_{j}\\ &+\sum\limits_{m=0}^{\infty}\left[D_{jm}^{(1)}H_{m}^{(1)}(k_{0}r_{j})+D_{jm}^{ (2)}H_{m}^{(2)}(k_{0}r_{j})\right]\sin m\theta_{j};\ j=1,2\end{array} \tag{5}\]  
where \(k_{0}=\omega/\beta_{0}\), and is the Hankel function of the first kind with argument x and order \(m\).  
All the wave motions must satisfy the boundary conditions given by  
\[\sigma_{yz}|_{y_{1}=D}=\mu\frac{\partial w}{\partial y_{1}}|_{y_{1}=D}=0 \tag{6}\]  
\[w|_{r_{1}=a}=w_{1}^{f}|_{r_{1}=a} \tag{7}\]  
\[\mu\frac{\partial w}{\partial r_{1}}|_{r_{1}=a}=\mu_{0}\frac{\partial w_{1}^{ f}}{\partial r_{1}}|  
\[\mu\frac{\partial w}{\partial r_{2}}|_{r_{2}=a}=\mu_{0}\frac{\partial w_{2}^{f}}{ \partial r_{2}}|_{r_{2}=a} \tag{11}\]  
\[\sigma_{\tau z}|_{r_{2}=b}=\mu_{0}\frac{\partial w_{2}^{f}}{\partial r_{2}}|_{r_ {2}=b}=0 \tag{12}\]  
Since the wave functions and boundary conditions are represented in different coordinate systems, the coordinate transformation is required. With the help of [PERSON]'s addition theorem [33; 34] of the oblique coordinate system, coordinate transformation can be carried out between any two coordinates, and the details will not be described again.