IBEM Simulation of Seismic Wave Scattering by a 3D Tunnel Mountain

-e seismic wave scattering by a 3D tunnel mountain is investigated by the indirect boundary element method (IBEM). Without loss of generality, the 3D physical model of hemispherical tunnel mountain in an elastic half-space is established, and the influence of the incidence frequency and angle of P or SV wave on the mountain surface displacements is mainly examined. It is shown that there exists quite a difference between the spatial distribution of displacement amplitude under the incident P wave and the one under SV wave and that the incidence frequency and angle of wave, especially the existence of tunnel excavated in the mountain, have a great effect on the surface displacements of mountain; the presence of the tunnel in the mountain may cause the greater amplification of surface displacement, which is unfavorable to the mountain projects. In addition, it should be noted that the tunnel may suffer the more severe damage under the incident SV wave.


Introduction
With the rapid development of economy, the highway, railway, and other infrastructures gradually extend to the mountainous area, and the mountain tunnel accounts for an increasing proportion in the line, for example, 70.2% in China Sichuan-Tibet railway and 86% in Japan central Shinkansen. A large number of theoretical research and earthquake damage investigations [1][2][3][4][5] have illustrated that the mountain/hill topography has a significant effect on the seismic surface motions, and, due to the scattering and interference of seismic waves inside the mountain, there exists a strong interaction between the tunnel and the mountain, which is of great importance for the seismic design of mountain tunnels. e existence of mountain or tunnel will significantly change the spatial distribution characteristics of ground motion and affect the damage degree and the distribution of local ground motions [6][7][8][9][10][11][12][13][14][15]. At present, the seismic response analysis methods commonly used can be divided into the analytical methods and the numerical methods. In terms of the analytical methods, Mow and Pao [16] adopted the wave function expansion method to study the dynamic stress concentration of a single cavity in infinite space under the incidence of elastic wave. Lee and Trifunac [17,18] further derived the analytical solution of the underground circular cavity to the incident SH wave. Yuan and Liao [19] gave a closed-form solution of two-dimensional scattering of plane SH waves by a cylindrical hill of circular-arc cross section in a half-space by the wave function expansion method and studied the effect of the height-to-width ratio of hill on surface ground motion. Smerzini et al. [11] conducted the effect of underground cavities on surface earthquake ground motions by considering the factors such as the cavity size, its embedment depth, the excitation frequency, and the incidence angle. Liang et al. [20][21][22][23] derived the analytical solution of the incident plane P and SV and SH waves by a semicylindrical hill, discussed the effect of incidence frequency and incidence angle on the surface motion of hill, and concluded that a hill greatly amplifies the incident plane wave. It is worth noting that the analytical methods mainly focus on the two-dimensional scattering of hill/tunnel, and there are relatively few studies on the seismic response of tunnel mountain.
Compared with the analytical methods, the numerical methods are more suitable for complex models, including the finite difference method (FDM) [24][25][26], the finite element method (FEM) [27,28], and the boundary element method (BEM) [29][30][31][32]. By the explicit finite element method, Li and Lu [33] established the model of underground caverns and studied the dynamic response of caverns, which is mainly determined by factors such as the characteristics of incident seismic wave, mountain topography, joint surface, and the caverns location. Zhu et al. [34] used the explicit finite difference method to analyze the influence of mountain topography and large geological structure on the seismic response of underground cavern and surrounding rock. Hao and Zhang [35] adopted the explicit finite element analysis method combined with the artificial transmitting boundary theory to evaluate the adjacent terrain effects on ground motion and found that the amplification of seismic waves by multimountain shape in a specific frequency band was more obvious than that by a single mountain. In terms of boundary element methods, Zhou and Chen [36,37] studied the scattering effect of P and SV and SH waves on irregular terrain and analyzed the amplification effect of terrain on ground motion qualitatively and quantitatively. Taking the large underground cavern as an example, Liu et al. [38] studied the dynamic interaction among multiple adjacent mountains by the dynamic indirect boundary element method and evaluated the amplification effect of ground motion for the incident P and SV waves.
It can be seen that the above-mentioned studies mainly focus on the 2D analysis or simple 3D topography, and the research on the amplification effect of 3D tunnel mountain is still poor. It is known that, compared with the finite element method and the finite difference method, the IBEM can accurately satisfy the radiation conditions at infinity without introducing the artificial boundary, and the boundary discretization is only needed. In addition, the identical precision can be reached for both stress and displacement, and no numerical dispersion at high frequency occurs. erefore, in this paper, the indirect boundary element method (IBEM) is used to solve the 3D scattering of P and SV waves by tunnel mountain; on the basis of the precision validation, the influence of the incidence frequency and incidence angle of wave on 3D tunnel mountain surface displacement amplitude is discussed in detail for a typical calculation example, and the conclusions have important theoretical and practical significance for revealing the physical mechanism of wave propagation in the tunnel mountain and determining the seismic precautionary requirements of engineering structures on the tunnel mountain more reasonably.

Physical Model
e model of a 3D tunnel mountain in an elastic half-space is shown in Figure 1, in which L1, L2, and L3 refer to the horizontal surface boundary, the mountain boundary, and the tunnel boundary, respectively; a and a 1 are the radii of hemispherical mountain and semicircular tunnel, respectively. e seismic wave is incident in x-z plane. Due to the spatial attenuation effect of the scattered wave, 5 times of the wavelength is taken to meet the calculation accuracy at the discrete range of surface.

Representation of Integral Equation.
In the absence of physical force, the equation of steady-state motion in the isotropic elastic solid medium is where λ and μ are the Lame constants, u is the displacement vector, ∇ is the Hamilton operator, ∇ 2 is the Laplace operator, ρ is mass density, and ω is the circular frequency. e time factor exp(iωt) will be omitted hereafter. By the single layer potential theory, in three dimensions, the displacement u i (x) and stress t i (x) of the scattered wave field at position x in the i direction can be expressed as the integration on a continuous surface S [39], respectively, and are shown as follows: where G ij (x, y) and T ij (x, y) represent Green's functions of displacement and traction, respectively, and ϕ j (y) can be regarded as the load density at y in the j direction on the boundary surface. e potential function is defined as follows: where q is a given constant (h is the compressional wave number, and k is the shear wave number) and r � |x − y|, in which x and y represent the positions of the field point and the source point, respectively. Green's functions of displacement and stress can be expressed as follows [40]: where i, j � 1, 2, 3 which correspond to the x, y, z direction and n i are the direction cosines of the unit normal vector. When x ≠ y, equations (2) and (3) can be calculated directly by the Gauss integral method:

2
Shock and Vibration where S y is the element area. When x � y, the integrals in equations (2) and (3) are singular. A discrete element is replaced by a disk with identical area, and a virtual uniform load is applied on the disk. Green's function expansion is adopted to solve the problem. , where R e is the radius of equivalent disk of the discrete elements.

Wave Field
Analysis. e total wave field can be composed of the half-space free field and the scattering field. e free field is the wave field of elastic wave incident on the halfspace bedrock. In the half-space, the P wave and the SV wave with circular frequency ω are incident at angles θ α and θ β , respectively. e wave potential functions can be expressed as follows: e incident P and SV waves will reflect at the half-space surface, and the reflection wave functions are as follows: where the reflection coefficients A 2 and B 2 were given in the literature [39]. h and k are the wavenumbers of the compressional wave and the shear wave, respectively.
According to the Huygens principle, the scattering field is generated by applying three orthogonal virtual uniform loads on the surface free boundaries near both the mountain and the tunnel. e total wave field is obtained by the superposition of the free field and the scattering field, and it can be expressed as follows: e three-dimensional horizontal surface L1, the hemispherical mountain surface L2, and the tunnel boundary L3 are all the surface free boundaries, at which the tractions are zero. According to equation (11), the boundary condition can be expressed as follows: Substituting equation (3) into equation (12) yields the following equation: In order to solve equation (13) numerically, the quadrilateral element is chosen to discretize the mountain surface and the nearby half-space surface. It is assumed that the number of discrete elements of the entire free boundaries is N, and the virtual uniform loads are applied on each element. For each boundary element with unknown constant load, a system of algebraic equations can be obtained as follows:

Validation of the Method
In order to validate the accuracy and efficiency of the IBE algorithm for the site response analysis, the ratio of tunnel radius to the mountain radius is 0.01. e numerical results for vertically incident plane P and SV waves by the IBEM are compared with those of [40,41]. e geometry of the 3D Gaussian shaped mountain with circular plane projection is expressed as follows: where a z � 0.5a is the height of mountain and a is the maximum horizontal dimension of the mountain. e P and SV waves are incident along z-axis, respectively. e parameters are given as follows: the shear-wave velocity of half-space c s � 400 m/s, density ρ � 2000 kg/m 3 , Poisson's ratio ] � 1/3, the viscous damping ratio of material is 0.01, and the convergence accuracy is 0.001. e dimensionless frequency is defined as η � 2a/λ � ωa/πc s � 1.0, in which λ is the wavelength of shear wave. e displacement amplitude of the mountain surface and the ground surface along x-axis (y � 0) are shown in Figure 2. It can be seen that the results by the IBEM have quite good agreement with those of [40,41]. erefore, the IBEM used for solving the seismic response of 3D tunnel mountain is suitable in this paper.

Results of Numerical Analysis
e 3D model of tunnel mountain is shown in Figure 1. e following parameters will be used for numerical analysis: the mountain radius a � 50 m, the tunnel radius a 1 � 5 m, Poisson's ratio ] � 0.25, the shear velocity of half-space c s � 1000 m/s, the density ρ � 2000 kg/m 3 , and the viscous damping coefficient � 0.01. Figures 3-6 show the displacement amplitude of the mountain surface for the vertically incident P wave with individual frequency, respectively. It can be seen that the peak displacements mostly appear on the mountainsides and near the tunnel entrance, and the higher the incidence frequency of wave, the more significant the displacement amplification effect in y direction. When η is 0.5, 1, 2, and 5, the corresponding peak displacement in y direction is 0.6, 2.1, 2.2, and 2.5 (Figures 3(b), 4(b), 5(b), and 6(b)), respectively. e peak displacement in the principal direction z is different from the ones in the secondary direction, and the maximum displacement appears on the mountain top. When η is 0.5, 1, 2, and 5, the peak displacement in the principal direction is 4.3, 5.2, 4.5, and 4.2 (see Figures 3(c), 4(c), 5(c), and 6(c)), respectively, and presents the ring distribution. With the increase of the incidence frequency of wave, the focal area of displacement amplification is reduced gradually.

Surface Displacement of the 3D Tunnel Mountain under the Incident P Wave.
For the case of incident wave with high frequency, the scattered waves interfere strongly with each other, which leads to the multiple peaks of displacement in the principal direction; for example, for η � 5, this phenomenon can be seen obviously in Figure 6(c). However, no analogous phenomenon occurs for the incident wave with lower frequency; as for the case in z direction, the only single peak displacement appears in the middle, and the values decrease gradually from the middle to the surrounding. As can be seen in Figure 3(c), no obvious interference occurs for η � 0.5. e surface displacement amplitudes of the 3D tunnel mountain for the incident P wave at 30°are shown in Figures 7(a)-7(c)-10(a)-10(c). It can be found that the distribution of the peak displacements in y direction is analogous to the one for the vertical incidence of P wave, and the location of peak displacement in x direction is different from that in the case of vertical incidence; that is, the peak displacements mostly appear on one side of the mountainside and near the tunnel entrance, and the amplification is more significant near the tunnel entrance; for example, for η � 0.5, the peak displacement at the tunnel entrance is 3.1 (see Figure 7(a)), while it is 2.7 for the case of vertical incidence (see Figure 3(a)). Meanwhile, the peak displacement in z direction is not on the mountain top but moves to the side. As shown in Figure 7(c), for the lower frequency incidence (η � 0.5), the location of single peak in z direction gradually moves from the mountain top (see Figure 3(c)) to the side. As shown in Figure 10(c), analogous to the case of vertical incidence, multiple peaks also appear for the high frequency incidence, which is caused by the interference of scattering waves. Figures 11-14 show the displacement amplitude of the mountain surface for the vertically incident SV wave with individual frequency, respectively. It can be seen that, for the lower frequency incidence, due to the existence of tunnel, the peak displacement of the tunnel 4 Shock and Vibration         mountain in x direction appears near the tunnel entrance; for example, for η � 1, the peak displacement near the tunnel entrance reaches 12.7 (see Figure 12(a)).

Surface Displacement of the 3D Tunnel Mountain under the Incident SV Wave.
For the case of lower frequency incidence of SV wave, e peak displacements in the secondary directions y and z mostly appear near the interface between the horizontal surface and the mountain, as illustrated in Figure 11(b); for η � 0.5, the maximum displacement in y direction appears near 45°in each quadrant of the interface between the mountain and the ground, and its value is about 1.6. Considering the high frequency incidence (see Figure 14), the scattering waves interfere strongly with each other, the phenomenon of multiple displacement peaks occurs, and, in particular, the location of peak displacements in x direction appears at the tunnel entrance, which is unfavorable to the earthquake resistance of tunnels.
When the SV wave is incident at α � 30°, the surface displacement amplitudes of the 3D tunnel mountain at individual frequency are shown in Figures 15-18. When the incidence frequency is lower, the amplification of the seismic wave in both the principal direction x and the secondary directions y and z is greater than that for the case of the vertical incidence and is particularly obvious in x direction; for example, for η � 0.5, the value of peak displacement is 3.5 for the vertical incidence (see Figure 11(a)); however, the value of peak displacement reaches 7 (twice as much as the vertical incidence) for the incidence at α � 30°(see Figure 15(a)), and the location of peak displacement is at one side of the tunnel entrance.
rough the above analysis, it is well known that the incidence direction of wave has a significant effect on the distribution of peak displacement. erefore, the incidence direction of the seismic wave should be considered in the seismic fortification and analysis. In addition, the comprehensive evaluation should be conducted for the urban plan and the earthquake resistance of structure in mountainous areas.

Conclusions
In this paper, based on the elastic wave theory combined with the idea of "zonal fit," the wave scattering by the 3D tunnel mountain under the incident P and SV waves is examined by the indirect boundary element method (IBEM). e following conclusions can be drawn through the numerical analysis: (1) e amplification of surface displacement amplitude of the 3D mountain is strongly influenced by the incidence frequency of wave, the incidence direction of wave, the wave type, and especially the existence of the tunnel excavated in the mountain. (2) For the vertically incident P wave, the maximum displacement in z direction appears on the mountain top, and the amplification of surface displacement near the tunnel entrance is also significant. For the oblique incidence, the peak location moves from the mountain top to one side of the mountain. As the incidence frequency of wave increases, the amplification area is gradually reduced, and the displacement extremisms appear alternately along the tunnel extension. (3) For the incident SV wave, the amplification of surface displacement in x direction near the tunnel entrance is also significant. For the vertical incidence, the value of peak displacement can reach 12.4, and the spatial distribution is symmetrical about the tunnel axis. In addition, the amplification at the oblique incidence is greater than that at the vertical incidence and is particularly obvious along the tunnel axis. (4) e presence of the tunnel in the mountain may cause the greater amplification, which is unfavorable to not only the mountain projects but also the tunnel. e mountain projects may suffer more severe damage under the incident P wave, while the tunnel suffers more severe damage under the incident SV wave. erefore, the amplification should be paid enough attention to in the earthquake resistance of tunnel mountain engineering.

Data Availability
e test data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.