Seismic ground response by twin lined tunnels with different cross sections

In this paper, the geometrical effects of shallow twin lined tunnels with different cross sections are investigated to obtain the anti-plane seismic ground motion under vertical/horizontal incident plane SH waves. A model of long two-dimensional lined tunnels is established and embedded in a homogeneous linear elastic half-plane by an applied numerical time-domain boundary element approach. In addition to a brief introduction to the formulation of the method, by considering five tunnel sections including circular, elliptical, horseshoe, square and rectangular, the surface response is sensitized to observe the normalized displacement amplitude/amplification ratio. In this regard, the angle of the incident wave and the frequency of the response are also included in changing the response pattern. To illustrate the results in both time and frequency domains, they are presented as blanket charts, snapshots, and three-/two-dimensional diagrams. The results showed that the seismic response of the surface is extremely affected by the geometric parameters of underground tunnels, which can create different conditions on the ground surface with shifting the direction of the wavefront. Geometrical effect of twin horizontally overlapping lined tunnels. Applying a time-domain half-plane boundary element method. Illustrating the response in time and frequency domains. The effect of depth and distance ratios on the seismic ground motion. Propagating vertical and horizontal incident SH-wave type. Geometrical effect of twin horizontally overlapping lined tunnels. Applying a time-domain half-plane boundary element method. Illustrating the response in time and frequency domains. The effect of depth and distance ratios on the seismic ground motion. Propagating vertical and horizontal incident SH-wave type.


Introduction
Study on the seismic ground motions and damage investigations in the presence of subsurface openings such as transportation tunnels and pipelines has been highlighted among the most important topics of geotechnical earthquake engineering [80]. This issue is more significant for shallow lined tunnels because of their stiffer materials compared to the surrounding medium as well as their proximity effect on the ground surface. These features can change the initial properties of incident waves to create the amplification/de-amplification on different zones of the surface. In this regard, utilizing appropriate approaches for modeling the problems of wave scattering embedded cavities. The dynamic response of twin circular tunnels subjected to incident SH waves was investigated by Balendra et al. [5]. In the following years, Shi et al. [69] studied the propagation of SH waves in the presence of a lined cavity embedded in an anisotropic media. Liang et al. [36] obtained a series solution for surface motion amplification due to underground twin tunnels. Hasheminejad and Kazemirad [22] and Jiang et al. [25] were able to present the seismic response of lined cylindrical cavity in a poroelastic medium. In this period, the effect of underground cavities on earthquake ground surface motion under SH-waves' propagation was investigated by Smerzini et al. [71]. Moreover, Yi et al. [89] analytically solved the scattering of plane SH waves in the presence of a circular lined tunnel with an imperfectly bonded interface. Fu et al. [16] presented an analytical solution for deforming twinparallel tunnels in an elastic half-plane. Gao et al. [18] and Faik Kara [17] analyzed the scattering of SH waves by a horseshoe-and cylindrical-shaped cavities, respectively. Later, Wang et al. [82] obtained analytical solutions of stresses and displacements for deeply buried twin tunnels in a viscoelastic rock. Recently, He et al. [23] analytically modeled the vibration prediction of two parallel tunnels buried in a full space. Although the analytical and semianalytical methods have high accuracy for benchmark purposes, their lack of flexibility in modeling and analysis of arbitrary and sophisticated shaped models with different boundary conditions has forced the researchers to use alternative approaches like numerical methods.
In recent decades, increasing the power of computers besides the development of numerical approaches has made researchers eager to use them for analysis of wave propagation problems as well as predicting the real responses of topographic features more than ever [40]. In this way, one can never claim that the obtained results are completely exact; rather, the main purpose is to move toward accurate responses. Based on the formulation, the numerical methods can be usually divided into two general categories known as the domain and boundary methods. The common domain methods are including the finite element method (FEM) and finite difference method (FDM). These methods require discretization of the whole body including internal parts of the model and its boundaries. A review of technical literature shows that Karakus et al. [32], Yamamoto et al. [88], Mirhabibi and Soroush [41], Zhang et al. [90], Das et al. [14], and Tsinidis [79] were among the researchers who studied the twin tunnel problems by FEM. In the use of FDM can name the studies of Chakeri et al. [9], Do et al. [15], Jamshasb et al. [26], Narayan et al. [44], and Ziaei and Ahangari [91] as well. Some other researchers like [2][3][4]82], Zeng et al. [92], Song et al. [72][73][74] and Song and Rodriguez-Dono [75] presented specific numerical approaches such as the convergence-confinement method (CCM) for tunnels' excavation in rock masses. Although the simplicity of domain methods makes them favorable for seismic analysis of finite media, the models are intricate because of discretizing the whole body and closing the boundaries at a considerable distance from the desired zone.
The boundary methods have overcome the shortcomings of domain methods. In boundary approaches, which are mostly known today as the boundary element method (BEM), due to the concentration of discretization only around the boundary of the desired features, the satisfaction of wave radiation conditions at infinity is achieved automatically. Moreover, the volume of input data/memory seizure and analysis time are remarkably reduced. On the other hand, because of the large contribution of analytical processes in solving various problems, the high accuracy of the obtained results is guaranteed. Therefore, the BEM provides a better manner for modeling infinite/ semi-infinite problems [46]. In this regard, Beskos [8] and Stamos and Beskos [68] were able to present a complete review of the BEM and its applications in modeling of underground structures, respectively. The BEM can be formulated in two general categories including full-plane (FP-BEM) and half-plane (HP-BEM) that each one developed in the transformed (frequency and Laplace) and time domains as well [12]. Among the pioneering researches based on FP-BEM, the studies of Crouch and Starfield [10], Yang and Sterling [85], and Xiao and Carter [84] can be addressed. The effective parameters on the underground tunnel stability were evaluated by Panji et al. [45]. Wu et al. [81] exhibited the stress distribution in subsurface structures. In another study, Panji et al. [48] analyzed the stability of shallow tunnels subjected to eccentric loads. Manolis and Beskos [39], 35], Kattis et al. [31], and Liu and Liu [37] presented their researches in the transformed domain. In the time domain, studies of Takemia and Fujiwara [78] and Kamalian et al. [30] can be addressed as well. In the application of FP-BEM, in addition to truncate the model from a full space, it is required to discretize all the boundaries of the problem including the interfaces, smooth ground surface and enclosing boundaries. This leads to the approximate satisfaction of stress-free conditions on the surface and reduces the accuracy of its results in some cases [1].
In the HP-BEM approach, the discretization of smooth ground surface and definition of fictitious elements for enclosing the boundaries are ignored and the stress-free boundary condition of the surface is satisfied in an exact process. Despite difficult implementation and creating large equations in the HP-BEM compared to the FP-BEM, the mentioned advantages help to make the simple models. The HP-BEM was statically developed by Telles and Brebbia [77] for stress analysis of plane problems and completed by Ye and Sawada [86] who studied the numerical properties of two-dimensional elastostatic problems. In the following years, Panji et al. [49] and Panji and Ansari [50] used static HP-BEM for analyzing the shallow tunnels and pressure pipes embedded in layered soils, respectively. This method was established for analyzing the dynamic problems in the transformed and time domains as well. In the frequency domain, Benites et al. [7] and Yomogida et al. [87] investigated anti-plane/plane strain elastic wave scattering for a system of multiple cavities. Rus and Gallego [67] presented a boundary integral equation (BIE) for sensitivity analysis of inclusion and cavity features in the harmonic elastodynamic medium. Parvanova et al. [47] studied the seismic response of lined tunnels with surface topography. Liu et al. [38] were able to obtain the seismic response of twin lined tunnels subjected to P/SV waves. Utilizing time-domain BEM, Hirai [19] presented the transient response of SH-wave scattering in a half-space. Then, Kamalian et al. [28] obtained the 2D site response of topographic structures. A direct time-domain boundary integral method for a viscoelastic plane with circular holes was suggested by Huang et al. [21]. The transient analysis of wave propagation problems was investigated by Panji et al. [46]. Recently, Panji and Ansari [51] studied the SH-wave scattering by the lined tunnels located in an elastic half-plane. Next, 52, 53 analyzed the seismic behavior of twin and multiple subsurface cavities and obtained the amplification patterns on the surface. The scattering of plane P and SV waves by twin lined tunnels with imperfect interfaces embedded in an elastic half-space was discovered by Huang et al. [24]. In the most recent studies of Panji & Ansari [55], Panji & Habibivand [54], Panji et al. [58,59,60] and Panji and Mojtabazadeh-Hasanlouei [56,57,59,60,61], they investigated the SH-wave dispersion by various topographic features. In another study, Mojtabazadeh-Hasanouei et al. [42] studied the seismic response of multiple subsurface inclusions as well. Recently, Mojtabazadeh-Hasanouei et al. [43] presented a review on SHwave propagation for orthotropic topographic features.
The reviews showed that the twin tunnels are frequently used as Metros, underground transportation tunnels, mining tunnels, and pipelines for water and petroleum conveyance [29]. On the other hand, as far as the authors know, a general seismic investigation was not done in the literature on the geometrical effects of twin lined tunnels embedded in an elastic half-plane. In the present study, a new and fast approach based on a direct half-plane time-domain boundary element method (HP TD-BEM) was applied to analyze the seismic ground response in the presence of twin circular-, elliptical-, horseshoe-, square-, and rectangular-shaped lined tunnels subjected to the propagation of incident SH waves. For this purpose, by using an appropriate sub-structuring process, the models are prepared in two separate parts including a twin-pitted half-plane and a system of twin liners that are assembled on each other. The proposed method is implemented in the general DAS-BEM 1 algorithm [46]. As a numerical scheme by considering some intended parameters, the synthetic seismograms beside two-and three-dimensional amplification patterns of the surface are presented. Moreover, the wave propagation below the surface is illustrated by wave scattering snapshots in different significant times.

Problem statement
As depicted in Fig. 1, a system of twin lined tunnels is embedded in a linear elastic homogeneous isotropic halfplane subjected to the incident out-of-plane SH waves of the Ricker type [65]. The function of Ricker wavelet is defined as Eq. (1): in which f p is the predominant frequency of the wave and the time-shifting is t 0 . By conducting the model in halfplane and satisfying the stress-free boundary conditions on the ground surface (G. S.), free-field displacement ( u ff ) can be obtained as follows [66]: where inc. , ref . , r inc. , and r ref . are obtained from Eqs. (3) and (4): In Fig. 1, Ω is the domain and Γ is the boundary of the body, which are defined separately for the pitted domain and the closed filled liners. b is the radius of tunnels, c is the wave velocity of the half-plane and is the angle of the incident waves. To initiate the formulation procedure, the equation of motion for the anti-plane strain model should be utilized as Eq. (5): In this equation, u j (x, y, t) is out-of-plane displacement and b j (x, y, t) is out-of-plane body force for the j th domain at the point (x, y) for the current time of t. The shear-wave velocity for the j th medium is defined by c j and determined by √ j ∕ j , where j is the shear modulus and j is the mass density for the mentioned domain. Based on the boundary conditions illustrated in Eq. (6), solving Eq. (5) results in a 2D anti-plane semi-infinite medium.
The half-space Green's functions can be achieved by simultaneously taking the singular solution into account for Eqs. (5) and (6)

HP time-domain BEM
To concentrate the meshes exclusively around the boundary of tunnels, the wave source image technique should be used to satisfy the stress-free boundary conditions on the ground surface. This procedure significantly simplifies the model and helps to decrease the volume of input data and time of analysis [46].

BIE
In the first step, the weighted residual integral without considering the boundary condition of Eq. (6) should be applied in Eq. (5). Using boundary methods for eliminating the volumetric integral defined on the domain, ignoring the body forces and contributions of the initial conditions, the direct BIE in the time domain can be achieved as Eq. (7) [6]: In the above equation, the half-space displacement and traction Green's functions of the time domain are shown by u * and q * , respectively. Also, u j and q j are the displacement and traction fields on the boundary; x and ξ are the coordinates of source and receiver, respectively. The Riemann-convolution integrals are shown by u * ⋅ q j , q * ⋅ u j and u ff 1 is the free-field displacement of the surface. The angle of boundary refraction (geometry coefficient) is defined by c( ) [13,20,27].

Numerical implementation
To obtain field variables, the discretization of geometric boundary of the body as well as time axis should be carried out before solving Eq. (7). Therefore, there is no approximation in Eq. (7) before discretizing on the boundaries of the twin-pitted medium and closed filled liners. An analytical process and a numerical procedure are applied in Eq. (7) to obtain its matrix form for temporal and spatial integration, respectively. The details of these processes can be found in Panji and Ansari [51,55].

Applications
The above formulation was implemented in the DAS-BEM algorithm. Some verification examples were solved, which can be found in the study of Panji and Ansari [51,55]. Thus, in the present work this part was ignored. The mentioned program is developed for the seismic analysis of 2D underground arbitrary-shaped structures such as twin lined tunnels embedded in an elastic half-plane.
To study the behavior of twin underground lined tunnels subjected to the seismic SH waves numerically, some key parameters including the shape of tunnels, angle of incident waves ( ), the depth ratio (DR) and the space ratio (SR) of the tunnels are considered. The shapes of circular, elliptical, horseshoe, square and rectangular are considered for the tunnels, respectively. In the following sections, the normalized displacement amplitude (NDA) is the ratio of the Fourier amplitude of the total ground surface motion obtained by BEM to the Fourier amplitude of the incident motion. The flexibility ratio ( F ) is the relative stiffness ratio of tunnels' lining to the surrounding soils and determined as (F = l c l ∕ m c m ) , where l and c l are the mass density and shear-wave velocity of linings; m and c m are the mentioned parameters for surrounding soils, respectively. The dimensionless frequency ( ) is defined as In the last part, some linear graphs were fitted to estimate the maximum responses in different cases for engineering applications (Sect. 6). The numerical procedure is implemented in the MATLAB programming software.  increased. For the circular-, horseshoe-and square-shaped tunnels, the amplitudes in central point of the surface are higher than the elliptical-and rectangular-shaped tunnels. Therefore, the large space occupied by the tunnels intensifies the volume of trapped waves to create a high-risk zone between the tunnels. In the following, Figs. 3 and 4 show the 2D time-domain responses for incident angles ( ) of 0° and 90°, respectively. The results are presented based on three stations marked by 1 to 3 on the ground surface. The stations 1 and 3 are located on the center line of tunnels and station 2 is placed at the central point of the surface between the tunnels. When the incident wave is applied vertically (Fig. 3a), the maximum amplitudes of squareand rectangular-shaped tunnels are recorded equal to 1.3 and 1.7, respectively, which are higher than the other models. In these specific cases, the smooth surfaces of tunnels sections increase the effect of wave trapping between the top boundary of tunnels and ground surface. This effect is weak when the tunnels have rotary cross sections like circular-, elliptical-and horseshoe-shaped models. As can be seen in Fig. 3b, the obtained amplitudes are increased to 2.6 for the second station. Therefore, the middle area between the tunnels is considered as a critical point for surface structures when the incident waves are radiated vertically. In Fig. 4, the horizontal scattering of incident waves increases the maximum amplitudes to 2.8. Therefore, when the incident waves are horizontally propagated, the tunnels directly reflect the waves toward the surface. In this case, the effect of trapped waves is weak, but the severe collision of the reflected waves to the surface creates higher amplitudes.

Snapshots
Presenting the wave scattering snapshots is the best way to see the seismic waves' propagation below the surface in significant times. Therefore, Figs. 5, 6, 7, 8 and 9 are presented to show the diffraction of vertical SH waves in the presence of twin lined tunnels with different shapes. In order for preparation of the snapshots, 14,028 and 10,201 internal points are considered for the domain and each of tunnels, respectively. The surface range is assumed between − 5b and 5b and the depth range of 0 ~ − 5b is considered as well. comparison of the responses in different shapes of tunnels, the snapshots are presented for six specific times including the seconds of 3, 3.5, 4, 5, 5.5 and 6, respectively. The results clearly show the effect of tunnels geometry on waves scattering to form the trapped waves between the boundary of tunnels and ground surface.

Frequency-domain responses
Illustrating the frequency-domain responses is the best way to show the general pattern of ground surface displacements in the presence of twin subsurface lined tunnels with different shapes subjected to seismic waves and Fig. 10 is presented for this purpose. The surface range of − 10b-10b and the dimensionless frequencies ( ) of 0.25-5 are considered, respectively. When the vertical incident waves are applied to the models, the responses are symmetric. The collision of the incident waves to the boundary of lined tunnels separates the direct waves to reflected and refracted phases of waves. A portion of the reflected waves leaves the medium, and another portion deviates from its initial paths to reach the surface. The main volume of the refracted waves is trapped between the tunnels and ground surface where they experience the intermittent reflections. This effect can significantly increase the amplifications between the tunnels. On the other hand, the lowest amplifications are obtained in the location of tunnels where a small volume of waves can enter. This phenomenon is called shadow zone [76], which is because of the stiffer material of linings to the surrounding soils, which creates a safe zone above the tunnels by blocking the passage of the waves. Comparing the responses shows that the maximum amplifications are occurred for the horseshoe-and square-shaped tunnels with the values more than 2.5 in central zone. The horseshoe-shaped tunnels create the large safe zones above the tunnels. On the other hand, when the elliptical-and rectangular-shaped tunnels are considered, the maximum amplifications are distributed in the central and peripheral areas with the lower value of 2.2, which exactly shows the effect of tunnels geometry on the surface displacements. In the following, the 2D frequency-domain responses are presented in Figs. 11 and 12 for the dimensionless frequencies of 0.5, 1 and 3 under vertical and horizontal incident waves. As can be seen in Fig. 11, the fluctuation of the results is increased at higher frequencies and the maximum NDA is recorded in the presence of horseshoe-shaped tunnels with the value of 4.6 at equal to 3 for the central point of ground surface. By applying the horizontal incident waves (Fig. 12), the NDA values are increased for different dimensionless frequencies when compare to Fig. 11. Based on the results (Fig. 12), the maximum displacement on the ground surface is obtained for square-shaped tunnels with the value of 5.4 at equal to 3 in the position of first tunnel. Therefore, when the seismic waves are propagated vertically, the safe zones are located above the tunnels. But, by applying horizontal waves, the upper areas of the tunnels have the highest seismic risk for surface structures.
To complete the responses, Figs. 13 and 14 are presented to show the amplifications versus different dimensionless frequencies at three surface stations. The responses are illustrated for between 0.1 and 3. In Fig. 13a, the maximum amplification of 1.7 in of 1.37 is recorded above the location of tunnels and related to the rectangular-shaped model. For the central point of the surface (Fig. 13b), this value is increased to 2.28 in of 3. Also, in horizontal radiation of incident waves, the maximum amplification is occurred in stations 1 (Fig. 14a) and 2 (Fig. 14b), which is related to the square-and horseshoe-shaped tunnels with the value of 2.5 in the of 3 and 0.9, respectively. Given that a lower volume of waves has reached the third station (Fig. 14c), the maximum value is decreased to 2.3 in of 2.4 for circular tunnels.

Conclusion
The geometrical effect of shallow twin lined tunnels was discovered by a direct half-plane time-domain BEM approach. The tunnels were embedded in a 2D linear elastic half-space subjected to vertical/horizontal incident SH waves. To establish the models, only the interfaces and internal boundaries of the lined tunnels were required to be discretized. A computer program was encoded based on the proposed method [51,55]. Then, by considering some intended parameters such as depth and space ratio of the tunnels and the angle of incident waves, a parametric study was evaluated to present the seismic ground motion in the presence of the twin lined circular-, elliptical-, horseshoe-, square-and rectangular-shaped tunnels, respectively. First, the 3D synthetic seismograms of the surface were illustrated for different models. To complete the results, the 2D responses of time-domain were demonstrated for specific stations on the surface. To show the wave propagation below the surface, some snapshots were presented at significant times. Finally, by converting the time-domain responses to frequency-domain, the displacements of the ground surface and amplification patterns were depicted as well. Some key results obtained from the parametric study are summarized as follows: 1. The HP time-domain BEM was able to prepare the simple models for time-history step-by-step seismic analy-sis of the ground surface in the presence of arbitrarily shaped twin subsurface lined tunnels. 2. Comparing the seismic responses of twin unlined underground cavities [52] with the present results showed that the role of lined tunnels was stronger on the seismic isolation to create the safe area on the ground surface. 3. When the incident wave was applied vertically/horizontally, the maximum amplitudes of 2.6 and 2.8 were recorded on the ground surface, respectively. 4. The maximum amplification occurred for the horseshoe-and square-shaped tunnels with values more than 2.5 in the central zone. 5. By applying the vertical incident waves, the maximum NDA was recorded in the presence of horseshoeshaped tunnels with a value of 4.6 in the central point of the ground surface. 6. In horizontally propagating waves, the maximum displacement was obtained for square-shaped tunnels with the value of 5.4 on the location of the first tunnel.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.