Dynamic interaction of twin vertically overlapping lined tunnels in an elastic half space subjected to incident plane waves

The scattering of plane harmonic P and SV waves by a pair of vertically overlapping lined tunnels buried in an elastic half space is solved using a semi-analytic indirect boundary integration equation method. Then the effect of the distance between the two tunnels, the stiffness and density of the lining material, and the incident frequency on the seismic response of the tunnels is investigated. Numerical results demonstrate that the dynamic interaction between the twin tunnels cannot be ignored and the lower tunnel has a significant shielding effect on the upper tunnel for high-frequency incident waves, resulting in great decrease of the dynamic hoop stress in the upper tunnel; for the low-frequency incident waves, in contrast, the lower tunnel can lead to amplification effect on the upper tunnel. It also reveals that the frequency-spectrum characteristics of dynamic stress of the lower tunnel are significantly different from those of the upper tunnel. In addition, for incident P waves in low-frequency region, the soft lining tunnels have significant amplification effect on the surface displacement amplitude, which is slightly larger than that of the corresponding single tunnel.


Introduction
With the rapid developing of underground space and the improvement of underground engineering technology, vertically overlapping tunnels (with vertical alignment) have been widely used for the urban subway and other underground transportation system in many large cities, for instance, the twin metro tunnels in Ankara, Turkey (Karakus et al. 2007), in Tehran, Iran (Chakeri et al. 2011), and in Wuhan, China ). On the other hand, it has been observed that the underground tunnel may suffer from serious damage in great earthquakes, such as Taiwan Chi-Chi (Wang et al. 2001) earthquakes. Hence, in the last decades, the seismic response of underground tunnel has become an attractive research topic in earthquake engineering and has been intensively investigated by numerous researchers analytically or numerically.
The analytical method such as wave function expansion method (WFEM) has been widely used to solve the scattering of seismic waves by a lined tunnel (Lee and Trifunac 1979;Liang and Ji 2006;Liu et al. 2013;Yi et al. 2014). Due to the fact that the analytical methods are usually restricted to simple calculation models, for complex geometrical and material characteristics, it is necessary to develop numerical methods, such as the finite element method (FEM), finite difference method (FDM), the boundary element method (BEM), etc. Kobayashi and Nishimura (1983) used the BEM to solve the dynamic response of an underground tunnel. Stamos and Beskos (1996) solved the three-dimensional seismic response of long lined tunnels in a half space by BEM. Kattis et al. (2003) used the BEM to study the harmonic body waves scattering by lined and unlined tunnels in an infinite poroelastic saturated soil. Rodriguez-Castellanos et al. (2006) studied the scattering of P and SV waves by cracks and underground cavities using the indirect boundary element method (IBEM). Esmaeili et al. (2006) analyzed the dynamic response of plane harmonic waves by a lined circular tunnel using the hybrid boundary and FEM. Yiouta-Mitra et al. (2007) adopted the FDM to study the dynamic response of a circular tunnel in a half space subjected to harmonic SV waves. Yu and Dravinski (2009) investigated the scattering of plane P, SV, and Rayleigh waves by a cavity embedded in an isotropic half space by BEM. Liu et al. (2010) further discussed the diffraction of P, SV waves by a tunnel in an elastic half space using a special BEM. Panji et al. (2013) studied the displacement response of an unlined truncated circular cavity in a homogenous isotropic medium under SH waves by BEM. Recently, Parvanova et al. (2014) investigated the seismic response of a lined tunnel in the half-plane with surfacetraction relief. Pitilakis et al. (2014) presented a series of numerical analysis to investigate the dynamic response of shallow circular tunnels. Alielahi et al. (2015) utilized the BEM to study the seismic ground amplification by unlined tunnels subject to vertically P and SV wave propagation.
Note that the above-mentioned studies are mainly focused on the single-tunnel model. As for twin tunnels or tunnel group, the dynamic interaction between closelyspaced tunnels should be taken into account (Hasheminejad and Avazmohammadi 2007;Smerzini et al. 2009;Chen et al. 2010;Wang et al. 2012;An et al. 2015;Fang et al. 2015;Alielahi and Adampira 2016a). Fotieva (1980) studied the effect of the compressional and the shear waves by twin-parallel tunnels. Balendra et al. (1984) solved the dynamic response of a pair of circular tunnels under SH waves. Okumura et al. (1992) studied the seismic response of the twin circular tunnels by FEM. Moore and Guan (1996) investigated the three-dimensional seismic response of twin lined tunnels in full space using the successive reflection method. Liang et al. (2003Liang et al. ( , 2004) discussed a series of solutions for surface motion amplification of the underground twin tunnels under P, SV waves. Chen et al. (2006) presented the Null-field integral equations for stress field around circular holes under anti-plane shear waves. Zhou et al. (2009) applied a semi-analytical method to discuss the dynamic response of twin-parallel elliptic tunnels embedded in an infinite poroelastic medium. Alielahi and Adampira (2016b) studied the seismic ground response under vertically in-plane waves by the twin-parallel tunnels.
As stated above, there have been several studies on the dynamic models of twin tunnels or tunnels group. However, to the author's best knowledge, for the vertically overlapping tunnels (coincidence of plane projection) shallowly buried in a half space, available theoretical analysis is extremely limited. Note that the investigations of Liang et al. (2003Liang et al. ( , 2004 only considered horizontally twin tunnels by an approximate analytic method. In fact, in Liu and Wang (2012), it has been illustrated that the seismic response of vertically overlapping tunnels is significantly different from that of horizontally flat twin tunnels in a full space. However, it is restricted to deep-buried tunnels. In order to improve the qualitative level of seismic design of vertically overlapping tunnels in a half space, it is necessary to calculate and reveal the dynamic interactions between the upper and the lower tunnels, and the influence of nearby ground surface.
In this paper, we focus on the dynamic interaction between these two vertically overlapping tunnels which are shallowly buried in an elastic half space, based on the indirect boundary integration equation method (IBIEM). It has been demonstrated that this method has several advantages such as reducing dimensions of problems, automatic satisfaction of boundary condition, and high calculation precision (Luco and De Barros 1994). Moreover, the IBIEM does not require element discretization, and it can thus be implemented more efficiently. The rest of this paper is organized as follows. The numerical procedure for IBIEM solution is present in section II. Then, the precision of the method is verified by the satisfaction extent of boundary conditions and the comparison between the degenerated and available solutions. Based on the IBIEM, the effects of key parameters, such as the distance between the two tunnels, the stiffness and density of the lining material, and the incident frequency on dynamic response are investigated in detail through numerical examples. Finally, several conclusions are drawn, which provide some useful insights for the seismic design of underground vertically overlapping tunnels. Due to the semi-analytical feature of the IBIEM, the numerical example in this study can also be regarded as a benchmark scheme for other numerical methods.

Model definition
As shown in Fig. 1, a pair of vertically overlapping lined tunnels are shallowly buried in the elastic half space with the depth d (to the upper tunnel center), the inner and outer radius of the upper tunnel and lower tunnel a 1 and a 2 , a 0 1 and a 0 2 , respectively. The distance between the centers of upper and lower tunnels is denoted as D. The domain of tunnels and the half space are assumed to be elastic, homogenous and isotropic. Let D 0 , D 1 , D 2 denote the domain of the half space, the upper tunnel and the lower tunnel, respectively. S and S 0 denote the outer and inner surface of the upper tunnel; correspondingly, S 0 and S 0 0 denote those of the lower tunnel. The shear modulus, Poisson ratio, and the density in D 0 are l 1 , m 1 , and q 1 , respectively. Then the velocity of the P and SV waves in the half space can be defined by c b1 ¼ ffiffiffiffiffiffiffiffiffiffiffiffi , where, l 2 , m 2 , q 2 , c a2 and c b2 are the material parameters of the two lined tunnels. Suppose that the plane P and SV waves propagate in a half space, then the two-dimensional scattering of plane P and SV waves by the vertically overlapping lined tunnels in a half space needs to be studied. For simplicity, the following calculation is limited to the circular tunnel case, while the IBIEM is well suited for tunnels with arbitrary shapes.

Numerical solutions by IBIEM
IBIEM was first proposed by Wong (1982), since then it has applied extensively, among which Liang and Liu (2009) solved the wave motion problems. Based on the single-layer potential theory, the IBIEM solves wave motion problems in the following way. The whole wave field is divided into free field (without scattering) and scattered field, and the scattered waves are constructed by using a linear combination of fundamental solutions (fictitious wave sources). In order to avoid singularities, the fictitious wave sources are placed at some distance to the physical boundaries of scatters, and the densities of the fictitious wave sources are determined by the boundary conditions. In this paper, the compressional line source and shear line source in a half space are introduced as the fundamental solution.
Assume that the wave potential functions of P and SV waves in the medium are / and w, respectively. The steady-state wave equation under two-dimensional plane strain state can be expressed as where k P , k S are the wavenumber of the P and SV waves, respectively. Define the shear modulus of the medium as l, the relationship between displacement, stress and the wave potential functions can be expressed as (Lamb 1904) where u x , u y denote the horizontal and vertical displacement, respectively, r xx , r yy and r xy are the normal stress and shear stress, respectively.
3.1 Green's functions of compressional and shear wave sources buried in elastic half-space It is known that the wave potential functions of compressional line source and shear line source in the whole space can be expressed by Hankel function of the second kind as 0 ðk S r 2 Þ, with r 2 being the distance between the wave source at ðx S ; y S Þ and the observation point at ðx; yÞ. Note that the time factor expðixtÞ is omitted here and thereafter for simplicity, with x and t being the excitation frequency and time variable, respectively.
According to the boundary condition of the free surface in the half space, and combining with the Fourier transform in the wave-number domain, the potential functions of total wave field can be derived (Lamb 1904).
(1) potential functions of total wave field in a half space under the compressional wave source can be expressed as (2) potential functions of total wave field in a half space under the shear wave source can be expressed as

Wave field construction
Define the inner and outer surfaces of lined upper tunnel are S and S 0 , respectively, then introduce the fictitious wave source surface S 1 around S to construct the scattered field in the half space and the fictitious wave source surfaces S 2 and S 3 to construct the scattered field in the lined tunnel. Similarly, the inner and outer surfaces of the lower tunnel are denoted by S 0 and S 0 0 , and the fictitious wave source surfaces for the lower tunnel are S 0 1 , S 0 2 and S 0 3 . Additionally, for convenience, suppose that the shapes of the fictitious wave source surfaces and the tunnel are identical.
The total wave field in the half space can be obtained as the superposition of the free field in the half space and the scattered field. First, we analyze the free field. Suppose that the plane P and SV waves with excitation frequency x, and with incident angles h a and h b , respectively, then the wave potential function in the orthogonal coordinates can be expressed as u ðiÞ ðx; yÞ ¼ exp½Àik a1 ðx sin h a À y cos h a Þ; ð12Þ Due to the existence of the half space surface, the incident P and SV waves will generate reflected P and SV waves as follows u ðrÞ ðx; yÞ ¼ A 2 exp½Àik a1 ðx sin h a þ y cos h a Þ; ð14Þ where A 2 , B 2 amplitude of the reflect waves, can be referred to Luco and De Barros (1994). Once the lined tunnel exists, scattered field will appear in the half space and the inner of the lined tunnel. Based on the single-layer potential theory, the scattered field in the half space and the tunnel can be constructed by all the compressional line sources and shear line sources. Suppose that the scattered field in the half space is generated by the fictitious wave source surfaces S 1 and S 0 1 , then the displacement and stress can be expressed as are the densities of the compressional line source and the shear line source at x 1 and x 0 1 on fictitious wave source surfaces S 1 and S 0 1 , respectively. G ij;l ðx; x 1 Þ and T ðsÞ ij;l ðx; x 0 1 Þ are the Green's functions for the displacement and the traction in the half space (with the subscripts 1 and 2 corresponding to the compressional line source and the shear line source, respectively), which satisfy the wave equations and surface boundary conditions automatically. Note that subscripts i, j = 1, 2 denote the x, y directions, respectively.
The scattered field in the upper tunnel can be constructed by the superposition of the compressional line sources and shear line sources acted on the fictitious wave source surfaces S 2 and S 3 , which can be expressed as where x 2 D 2 , x 2 2 S 2 , x 3 2 S 3 , dðx 2 Þ, eðx 2 Þ are the densities of the compressional line source and shear line source at x 2 on S 2 , and f ðx 3 Þ, gðx 3 Þ denote those at x 3 on fictitious wave source S 3 ,G ðtÞ i;l , T ðtÞ ij;l are the Green's function for the displacement and the traction in the lined tunnel, respectively.
Similarly, the stress and displacement of the lower lined tunnel can be obtained.

Boundary conditions and the numerical solutions
Due to the adoption of the fundamental solution for the elastic half space, the boundary condition of the free ground surface can be satisfied automatically. Thus, we only need to consider the continuity of displacement and traction on the interface between the lining and surrounding soil, and the zero traction condition on the inner surface of the lining, which are as follows: where the superscripts s, t denote the half space and the tunnel, respectively. For ease of numerical solution, we discrete the inner and outer surface of the tunnels and fictitious wave source surfaces. Suppose that the number of discrete points on the internal and external surface of the tunnels is N, and that of the fictitious wave source surface is N 1 . Then the scattered displacement field and stress field in the half space can be expressed as where b n 1 , c n 1 , and b 0 n 1 , c 0 n 1 are the source densities of P and SV waves for the n 1 -th point on fictitious source surfaces S 1 and S 0 1 , respectively. The scattered displacement field and stress field in the upper tunnel can be expressed as i;2 ðx n ; x 3;n 1 Þ; ð25Þ ij;2 ðx n ; x 3;n 1 Þ; ð26Þ x n 2 D 1 ; x 2;n 1 2 S 2 ; x 3;n 1 2 S 3 ; n ¼ 1; 2; Á Á Á; N; where d n 1 , e n 1 and f n 1 , g n 1 are the source densities of P and SV waves for the n 1 -th point on fictitious source surfaces S 2 and S 3 . Note that here we assume that the discrete numbers of S 2 and S 3 are also N 1 . Similarly, the scattered field in the lower tunnel can be constructed by the discrete wave source on S 0 2 and S 0 3 . From Eqs. 20-22, we can obtain where H 1 , H 2 , H 3 , H 0 11 are the Green's influence matrices relating to the displacements and tractions on the discrete points of the outer surface of the upper tunnel caused by the fictitious wave sources on S 1 , S 2 , S 3 , S 0 1 ; H 0 1 , H 0 2 , H 0 3 , H 11 are the Green's influence matrices for the outer surface of the lower tunnel.T 2 , T 3 are the Green's matrix (stress) relating to the traction on the discrete points of the inner surface of the upper tunnel caused by the fictitious wave sources on S 2 , S 3 and T 0 2 , T 0 3 are the corresponding Green's influence matrices for the internal surface of the lower tunnel. Y 1 , Y 2 , Y 3 , Y 0 1 , Y 0 2 , Y 0 3 are fictitious wave source densities vectors on the fictitious surfaces S 1 ,S 2 , S 3 ,S 0 1 , S 0 2 ,S 0 3 , respectively.F 1 , F 2 are the free field vector related to the displacement and traction on the interfaces.
Equations 27-30 are written as a compact form as HA ¼ B, and this overdetermined equation can be solved by least square (LS) method where H, H T and H T are the conjugate, transpose, and conjugate transpose matrices of H, respectively. After solving the equations about the fictitious wave source densities from (31), we can obtain the scattered field. The total wave field is the superposition of scattered field and Earthq Sci (2016) 29(3): 185-201 189 free field, and then we can calculate the displacement, stress at any location both in the half space and the tunnels.

Verification of accuracy and validation of the numerical solution
Until now, it is still a challenging task to obtain the accurate analytical solution for the scattering of plane P, SV waves by the lined tunnel in an elastic half space due to the difficulty in dealing with the traction-free boundary condition of the half-space surface. Thus, we verify the accuracy by the following steps: (1) test the satisfaction extent of the boundary conditions, (2) examine the numerical stability of the solutions, (3) and degenerate solutions to the single-tunnel case with well-known solutions. Then define the non-dimensional frequency g ¼ xa 1 =pc b1 , with c b1 being the shear wave velocity in a half space medium. To test the boundary conditions, plenty of calculation results show that with the increase of discrete points, the boundary residual value decreases gradually. When the incident frequency g = 2.0, the residual value can reach up to 10 -4 for N = 120 and N 1 = 80. The notations and symbols are summarized in Table 1.
To verify the numerical stability of the solution, Tables 2, 3 and 4 illustrate the convergence of the surface displacement amplitudes and the hoop stress amplitudes on the inner surface of tunnels with the increase of discrete points. Parameters are set a 1 = 1.0, d/a 1 = 4.0, D/a 1 = 3.0, a 1/ a 12 = 0.9, a 0 1 =a 0 2 = 0.9, damping ratio f = 0.001, g = 1.0, v = 0.25, q 2 =q 1 = 5/4, c b2 =c b1 = 5/1(rigid lining), N = 60, 80, and 120, correspondingly N 1 = 40, 60, and 80. It is clearly shown that the surface displacement amplitudes and stress converge quite well with the increase of discrete points. This further validates the excellent numerical stability of this solution.
When the non-dimensional frequency is 0.5, the results indicate that as the depth of the lower tunnel is larger than 30a 1 , the impact to the upper tunnel or the ground surface can be ignored. To degenerate the solution to single-tunnel case, the following parameters are taken: D/a 1 = 200, q 2 =q 1 = 1.0, c b2 =c b1 = 1.0, damping ratio f = 0.001, v = 1/3 and g = 0.5. Figure 2 shows the surface displacement amplitudes of the half space and the dynamic stress concentration factors on the inner surface of the upper tunnel compared with the well-known results of the elastic half space by Luco and De Barros (1994) for vertically incident P waves and the buried depth of d/a 1 = 1.5 and 5.0, respectively. Inner and outer radius of lower tunnel q 1 , q 2 Density of D 0 and tunnels N Discrete points v 1 , v 2 Poisson ratio of D 0 and tunnels g Non-dimensional frequency l 1 , l 2 Shear modulus of D 0 and tunnels f Damping ratio c a1 , c b1 Velocity of the P and SV waves in D 0 h a , h b Incident angle of P and SV waves c a2 , c b2 Velocity of the P and SV waves in tunnels Table 2 Numerical stability verification of surface displacement amplitudes under vertically incident P and SV waves (g ¼ 1:0)   Figure 3 shows the surface displacement amplitudes of the half space compared with the well-known results of the elastic half space by Liu et al. (2013) for vertically incident P waves and the buried depth of d/a 1 = 1.5. The parameters are defined as D/a 1 = 200, a 2 =a 1 ¼ a 0 1 =a 0 1 ¼ 1.1, l 2 =l 1 ¼ 0.8, Poisson ratio v = 1/3 and g = 0.5. It is shown that the results of this study are in good agreement with the references.

Numerical results
A pair of vertically overlapping circular lined tunnels is shallowly buried in the elastic half space, with the buried depth of the upper tunnel d/a 1 = 4.0, and the inner and outer radius ratios a 1 /a 1 = 0.9 and a 0 1 =a 0 2 = 0.9. Considering the variation of the distance between the two tunnels, we choose D/a 1 = 3.0, 4.0, 5.0 (the normal range in practical engineering).
The practical parameters of the lining and the medium in the half-space are defined as follows: the radius and thickness of a real tunnel are a 1 = 3 and 0.33 m,  Luco and Barros (1994) for vertically incident P waves with the buried depth of the upper tunnel (d=a 1 = 1.5 and d=a 1 = 5), distance between the twin tunnels (D=a 1 = 200), frequency (g = 0.5), and Poisson ratio (t = 1/3) respectively. Considering the variation of the stiffness and density of the lining material, take the soft lining, homogenous lining, and rigid lining for parameter analysis. The stiff tunnel (q 2 =q 1 = 5/4, c b2 =c b1 = 5/1) corresponds to the reinforced concrete tunnel in soft soil layer such as in Tianjin, Shanghai of China. The lining of the stiff twin tunnels are made of concrete and the density of tunnels lining is q 2 ¼ 2000 kg/m 3 . The shear wave velocity in the lining is c b2 ¼ 2000 m/s. Poisson ratio of half space and tunnels is v 1 ¼ v 2 ¼ 0.25, and damping ratio is f = 0.001. In contrast, the soft tunnel (q 2 =q 1 = 4/5, c b2 =c b1 = 1/3) denotes the case such that a shotcrete tunnel in granite. Moreover, the unlined tunnel is identical to the case of q 2 =q 1 = 1/1, c b2 =c b1 = 1/1. For the convenience of analysis, the parameters used in follows are relative values, such as the non-dimensional frequency. For simplify, only the vertically incident P and SV waves are considered herein with h a (h b ) = 0°.

Dynamic hoop stress and deformation of vertically overlapping tunnels
First, define the dimensionless dynamic stress amplitudes , where k b1 is the wavenumber of SV waves in an elastic half space. Since the stress of the soft lining is not very large, we only study the hoop stress of the rigid lining in this paper. Figures 4 and 5 show the distribution of the hoop stress along the inner surface of the twin tunnels for vertically incident P and SV waves with distance D/a 1 = 3.0, 4.0, 5.0. For comparison, the case of single tunnel (identical with the upper tunnel) is also presented. The dimensionless frequency takes g = 0.5, 1.0, 2.0, respectively. It shows that the hoop stress under high frequency is quite different from those under low frequency. For low frequency (g = 0.5), the hoop stress of the upper tunnel under different intervals between tunnels (D/a 1 = 3.0, 4.0, 5.0) is similar to those of the single one. However, the hoop stress in the lower tunnel is significantly different among different intervals. For example, the max hoop stress value for D/a 1 = 5.0 is 19.6, while the max hoop stress is 10.7 for D/a 1 = 3.0 under P waves. From Figs. 4 and 5, we can observe that when the frequency is 0.5, the depth of the lower tunnel (the normal distance for practical subways) has little impact on the hoop stress, but the impact to the upper tunnel becomes significant when the frequency is 2.0, and the shielding effect decreases with the increase of the distance between the twin tunnels. The results coincide with the conclusion in Huang and Zhang (2013). At the high frequency (g = 2.0), the effect of the distance between the twin tunnels seems to be not so noticeable, but the hoop stress of the single tunnel is significantly larger than that of the upper tunnel. For example, the max hoop stress value of the single tunnel is 26.88, while that of the twin tunnels is only 9.35. It can be attributed to the strong scattering of highfrequency waves around the lower tunnel which leads to the significant attenuation of incident waves. This phenomenon implies that the lower tunnel may play a protective role for the seismic safety of the upper tunnel for high-frequency waves. Note that Chen and Chen (2008) have investigated the seismic response of vertical doublelayered metro tunnels under near-fault strong ground motion by FEM, and similar conclusions have been obtained. While, the studies of Liang et al. (2003Liang et al. ( , 2004 indicate significant amplification effects due to the dynamic interaction of horizontally twin tunnels, which are substantially different from the model of vertically overlapping tunnels in this study. Figures 6 and 7 illustrate the hoop stress amplitude spectrums on surface of twin tunnels under vertically incident P and SV waves for D/a 1 = 3.0. We select the observation points on the inner surface of the tunnel located at h = 0°, ± 30°, ± 90°for P waves incidence, and h = 0°, ± 30°for SV waves. It is clearly shown that the stress amplitude highly depends on the incident frequency and the space location, and there are many peaks and troughs of the spectrum curve with the peak values usually appearing in low-frequency region g\0.5. As for P waves, the hoop stress amplitude can reach up to 29.5 at h = 0°, but the amplification effects seem not remarkable at the apex of arch (h = 90°). In addition, both the peak value and the spectrum characteristic of the upper tunnel are significantly different from those of the lower tunnel. In general, for the lower tunnel, there are more peak frequencies and the peak value is slightly larger than that of the upper tunnel in this case. Compared with the case of P waves, the dynamic stress concentration seems more significant for SV waves and the peak value appears around h = 30°for low-frequency waves and around h = 0°for high-frequency waves. In addition, the spectrum for incident SV waves oscillates more violently than that for P waves.
Considering the incidence of P waves and the variation of the interval between the tunnels D/a 1 = 3.0, 4.0, 5.0, Figs. 8 and 9 illustrate the hoop stress amplitude spectrums at the points h = 0°and h = 30°on the inner surface of these twin tunnels and the corresponding single tunnel (only the upper tunnel or the lower tunnel exists). It shows that in low-frequency region (g B 0.5), the peak stress amplitude of the upper tunnel is close to the single-tunnel case. However, for the high frequencies, at h = 0°, the hoop stress amplitude of the single tunnel is much larger than the case of twin tunnels. For example, at g = 1.56, the hoop stress amplitudes of the single tunnel and the upper one of twin tunnels with D/a 1 = 3.0 are 13.80 and 7.25, Fig. 4 Hoop stress amplitudes of the inner wall of twin tunnels and the single tunnel (identical with the upper tunnel) subject to vertically incident P waves, with frequency (g = 0.5, 1.0, 2.0), distance between the twin tunnels (D=a 1 = 3, 4, 5), and the rigid lining (q 1 =q 2 = 0.8, c b1 =c b2 = 0.2) Fig. 5 Hoop stress amplitudes of the inner wall of twin tunnels and the single tunnel (identical with the upper tunnel) subject to vertically incident SV waves, with frequency (g = 0.5, 1.0, 2.0), distance between the twin tunnels (D=a 1 = 3, 4, 5), and the rigid lining (q 1 =q 2 = 0.8, c b1 =c b2 = 0.2) Earthq Sci (2016) 29(3):185-201 193 respectively. But at the point h = 30°, the influence of the lower tunnel seems not so noticeable. Figures 10 and 11 illustrate the hoop stress amplitude spectrums at points h = 0°and h = 30°on the inner surface of these twin tunnels and the corresponding single tunnel for incident SV waves. It shows that in low-frequency region (g B 0.5), the peak stress amplitude of the upper tunnel can be larger than the single-tunnel case. For high-frequency waves, in contrast, the shielding effect of the lower tunnel is very pronounced. For example, at g = 1.5, the stress amplitude at h = 0°of the single tunnel is 38.4, while for the upper one of twin tunnels, the amplitudes are 22.7, 20.8, and 25.7 for D/a 1 = 3.0, 4.0 and 5.0 respectively. As for Fig. 6 Hoop stress amplitude spectrums at different points (h = 0°, ±30°, ±90°) on the inner wall of twin tunnels, with the distance between the twin tunnel (D=a 1 = 3), and the rigid lining (q 1 =q 2 = 0.8, c b1 =c b2 = 0.2), subject to vertically incident P waves Fig. 7 Hoop stress amplitude spectrums at different points (h = 0°, ±30°, ±90°) on the inner wall of twin tunnels, with the distance between the twin tunnel (D=a 1 = 3), the rigid lining (q 1 =q 2 = 0.8, c b1 =c b2 = 0.2), subject to vertically incident SV waves Fig. 8 Hoop stress amplitude spectrums at h = 0°on the inner wall of twin tunnels, with the buried depth of the single tunnel (d=a 1 = 4 denotes the single upper tunnel, d=a 1 = 7 denotes the single lower tunnel), the distance between the twin tunnels (D=a 1 = 3, D=a 1 = 4, D=a 1 = 5), and the rigid lining (q 1 =q 2 = 0.8, c b1 =c b2 = 0.2), subject to vertically incident P waves the lower tunnel, the displacement spectrum characteristics become more complicated with the variation of D/ a 1 . As the buried depth increase, the spectrum curves oscillate more rapidly. Comparing the case D/a 1 = 3.0 with the corresponding single-tunnel case, we can observe that, for the high-frequency waves, the stress amplitude in the lower tunnel is decreased under the influence of the upper tunnel. Figures 12 and 13 illustrate the imaginary and the real parts of the deformation of the inner wall of tunnels for vertically incident P and SV waves with D/a 1 = 3.0 and the non-dimensional frequency g = 0.5, 1.0, 2.0. It is Fig. 9 Hoop stress amplitude spectrums at h = 30°on the inner wall of twin tunnels, with the buried depth of the single tunnel (d=a 1 = 4 denotes the single upper tunnel, d=a 1 = 7 denotes the single lower tunnel), the distance between the twin tunnels (D=a 1 = 3, D=a 1 = 4, D=a 1 = 5), and the rigid lining (q 1 =q 2 = 0.8, c b1 =c b2 = 0.2), subject to vertically incident P waves Fig. 10 Hoop stress amplitude spectrums at h = 0°on the inner wall of the twin tunnels and the single tunnel, with the buried depth of the single tunnel (d=a 1 = 4 denotes the single upper tunnel, d=a 1 = 7 denotes the single lower tunnel), distance between the twin tunnels (D=a 1 = 3, D=a 1 = 4, D=a 1 = 5), and the rigid lining (q 1 =q 2 = 0.8, c b1 =c b2 = 0.2), subjected to vertically incident SV waves Fig. 11 Hoop stress amplitude spectrums at h = 30°on the inner wall of the twin tunnels and the single tunnel, with the buried depth of the single tunnel (d=a 1 = 4 denotes the single upper tunnel, d=a 1 = 7 denotes the single lower tunnel), distance between the twin tunnels (D=a 1 = 3, D=a 1 = 4, D=a 1 = 5), and the rigid lining (q 1 =q 2 = 0.8, c b1 =c b2 = 0.2), subject to vertically incident SV waves Earthq Sci (2016) 29(3):185-201 195 Fig. 12 Comparisons of the real part and the imaginary part deformation of the twin tunnels' inner wall by the initial shape, with frequency (g = 0.5, 1.0, 2.0), distance between the twin tunnels (D=a 1 = 3), and the rigid lining (q 1 =q 2 = 0.8, c b1 =c b2 = 0.2), subject to vertically incident P waves Fig. 13 Comparisons of the real part and the imaginary part deformation of the twin tunnels' inner wall by the initial shape, with frequency (g = 0.5, 1.0, 2.0), distance between the twin tunnels (D=a 1 = 3), and the rigid lining (q 1 =q 2 = 0.8, c b1 =c b2 = 0.2), subject to vertically incident SV waves shown that the deformation feature of the upper tunnel may be largely different from that of the bottom one, and for high-frequency waves, the shielding effect of the lower tunnel on the upper one can be also clearly seen for both P and SV waves.
5.2 Ground surface displacement response above the twin tunnels Figures 14 and 15 show the vertical and horizontal surface displacement amplitudes on the ground surface above the tunnels for incident P and SV waves, considering the variation of the stiffness and the distance between these two tunnels. For comparison, the case of single tunnel (identical with the upper tunnel) is also presented. The nondimensional frequency takes g = 0.5, 1.0 and 2.0, respectively. The surface displacement amplitudes are normalized by the displacement amplitudes of incident waves throughout the paper. For simplicity, only the vertically incident seismic waves are considered. It can be seen that the incident frequency has large influence on the surface displacement response. For the soft tunnels, the amplification or deamplification effects strongly depend on the incident frequency. For the rigid tunnels, the shielding effect on nearby ground surface seems relatively more prominent. For example, at g = 2.0 the displacement amplitude at x=a 1 = 0.0 is only 0.27 and 0.42 for P, SV waves incidence respectively, while it is 2.0 for the free field. Furthermore, the displacement amplitudes and spatial distribution above the vertically overlapping tunnels are similar to those of the single one, and the influence of the distance between the twin tunnels on the Fig. 14 The vertical displacement amplitudes of ground surface above the twin tunnels and the case of the single tunnel (identical with the upper tunnel), for soft lining (q 1 =q 2 = 1.25, c b1 =c b2 = 3.0), homogenous lining (q 1 =q 2 = 1.0, c b1 =c b2 = 1.0), rigid lining (q 1 =q 2 = 0.8, c b1 =c b2 = 0.2), with incident frequency (g = 0.5, 1.0, 2.0), distance between the twin tunnels (D=a 1 = 3, D=a 1 = 4, D=a 1 = 5), subject to vertically incident P waves surface displacement amplitudes is not so noticeable, especially for high-frequency cases. Figures 16 and 17 show the surface displacement amplitude spectrums of the points (x=a 1 = 0.0, x=a 1 = 2.0) near the lined tunnels under vertically incident P and SV waves for both the vertically overlapping tunnels and the single tunnel (identical with the upper tunnel). It shows that the lining material of the tunnel has a significant effect on the surface displacement amplitudes and spatial distribution around the tunnels. In the case of soft lining and homogenous lining, the spectrum curves oscillate violently, while the spectra curves oscillate relatively smoothly for the rigid lining tunnel. As for incident P waves at low frequencies, the amplification effect is noticeable. For an example, at g = 0.54, the displacement amplitude at the point x=a 1 = 0.0 for the soft tunnel is up to 3.0 in the case of D=a 1 = 3.0, which is about 50 % larger than that of the free field. However, in the case of rigid tunnels, mainly the shielding effect of the tunnels can be observed, and the surface displacement amplitudes above the single tunnel are larger than those of the twin tunnels. Because under such circumstance, the twin tunnels have more pronounced shielding effect on the seismic waves. Moreover, the effect of the distance between twin tunnels on the surface displacement amplitudes seems not so significant. As for incident SV waves, it seems that the displacement spectrum curves oscillate more quickly, and note that in the case of rigid tunnels, the amplification effect can be observed above the tunnel (x=a 1 = 0.0) with the peak amplitude of horizontal displacement 2.36 for D=a 1 = 3.0 and g = 0.32, which is slightly larger than that of the single-tunnel case. Fig. 15 The horizontal displacement amplitudes of ground surface above the twin tunnels and the case of the single tunnel (identical with the upper tunnel), for soft lining (q 1 =q 2 = 1.25, c b1 =c b2 = 3.0), homogenous lining (q 1 =q 2 = 1.0, c b1 =c b2 = 1.0), rigid lining (q 1 =q 2 = 0.8, c b1 =c b2 = 0.2), with incident frequency (g = 0.5, 1.0, 2.0), distance between the twin tunnels (D=a 1 = 3, D=a 1 = 4, D=a 1 = 5), subject to vertically incident SV waves It should be mentioned that in above examples we only considered the small damping case (the damping ratio is 0.001). We calculated the displacement amplitudes of ground surface with the damping ratio 0.001, 0.01, 0.05.

Conclusions
The scattering of plane harmonic waves by a pair of vertically overlapping lined tunnels shallowly buried in an elastic half space is solved by a high-precision IBIEM. The convergence and numerical stability of the IBIEM for this model are verified. Through detailed numerical analysis, several beneficial conclusions can be drawn.
(1) Numerical results have been shown that the scattering of seismic wave strongly depends on the distance between these twin tunnels, the material properties of the lining and the non-dimensional frequency. The dynamic interaction between these twin tunnels cannot be neglected.
(2) It has been shown that the lower tunnel may play a protective role of isolating the P and SV waves of high frequency, leading to the great decrease of dynamic hoop stress amplitude of the upper tunnel, up to 50 % smaller than the case of the single tunnel at some frequencies. However, for the low-frequency SV waves, on the contrary, additional amplification effect between these twin tunnels can be observed at some locations in the upper tunnel. (3) The soft tunnels may have significant amplification effect on the ground motion above the tunnels, while the rigid tunnels mainly have shielding effect on the nearby surface displacement response. The influence of the vertically overlapping tunnels seems slightly larger than that of the single-tunnel case. (4) The deformation feature and the dynamic stress of the upper tunnel are significantly different from those of the lower tunnel, and there are more peaks and troughs of the spectrum curve of the lower tunnel with large buried depth.
This study was limited to the frequency domain analysis. But from the perspective of engineering, it is more attractive in the time domain. Moreover, the half-space model is suitable for the thick near-surface layer. For more general scenarios, the layered half-space model should be adopted. We leave these limitations for the future study.