Numerical study on the characteristics of the seismic response of subway shield tunnels beneath rock mountains: A case study in China

In this study, we investigated the characteristics of the seismic response of a subway shield tunnel beneath a rock mountain through the finite-element method based on a practical subway tunnel project in China. We discuss the influence of the spectral nonstationarity of seismic motion on tunnel responses. First, a set of quasistationary (QS) and fully nonstationary (FN) bedrock motions were artificially simulated. Second, the background information of the referenced project was briefly discussed, and correspondingly, several refined numerical models were established. Finally, the mutual relationship between a subway tunnel and the overlying rock mountains under QS bedrock motions was studied, and the results under the QS and FN bedrock motions were compared to clarify the effect of the spectral nonstationarity of ground motions. We found that i) an overlying rock mountain can suppress the acceleration response of an underground tunnel, and combined with the seismic motion, high confining pressure owing to the existence of an overlying mountain aggravates the radial shear-stress response and damage to the tunnel; ii) the spectral nonstationarity of input motions results in higher acceleration responses at the mountaintop and on the tunnel and causes greater radial shear stress and damage to the tunnel, even though the input PGAs and waveforms for QS and FN motions are the same. Therefore, artificial seismic motions exhibiting spectral nonstationarity may provide conservative results.


Introduction
Tunnels have become one of the most important lifeline engineering in modern cities. It has long been considered that underground structures, such as subway tunnels, can sustain earthquakes. However, recent events, such as the Kobe Earthquake in Japan (1995), Duzce Earthquake in Turkey (1999), Chi-chi Earthquake in Taiwan (1999), Bam Earthquake in Iran (2003), and Wenchuan Earthquake in China (2008), have shown that tunnels are susceptible to severe damage due to seismic loading (Cilingir and Madabhushi, 2011). Thus, underlying earthquake hazard constitutes a severe threat to the safe operation of tunnels built in seismically active areas.
Many researchers have conducted seismic analyses of tunnels. Hashash et al. (2001) provided an important indication of the seismic design and analysis of underground structures. Cilingir and Madabhushi (2011) investigated the effects of input motion characteristics on the seismic behavior of circular and square tunnels using experimental and numerical methods. Pakbaz and Yareevand (2005) studied the interaction between the ground and tunnel lining during earthquake excitation by finite difference analysis. Chen et al. (2018) conducted several dynamic centrifuge tests to explore the effect of isolation layer and cross-section dimensions on shield tunnels. Fabozzi et al. (2018) employed 3D full dynamic analysis to investigate the effect of ground motion asynchronism on the distribution of internal forces in tunnel linings. Miao et al. (2018aMiao et al. ( , 2018b) studied the influence of ground motion spatial variation on the seismic response of subway tunnels using a large-scale finite-element model. Focusing on the tunnel embedment ratio, soil-tunnel interface conditions, lining thickness, tunnel shape, and input ground motion, Patil et al. (2018) conducted parametric studies and discovered that the distortion of lining depends on the depth of embedment and flexibility ratio of the tunnel. Yan  Mountain tunnels are assumed to be seismic resistant as they are situated deep in rock layers (Towhata et al., 2008). Recently, however, surveys on earthquake damage have given disproved this widely accepted view . Consequently, an increasing number of studies on the seismic performance of mountain tunnels have been reported. Shen et al. (2020), Guan et al. (2019), and Xu et al. (2016) conducted shaking table tests to investigate the seismic behavior and damage mechanism of mountain tunnels. Andreotti and Lai (2019) proposed a numerical method for developing fragility curves for mountain tunnels subjected to transversal seismic loading. Significant efforts on the collection, organization, and classification of seismic damage to mountain tunnels due to earthquakes have been made by many researchers, including Dowding and Rozan (1978), Asakura and Sato (1996) affecting the response of mountain tunnels. Although the seismic response of mountain tunnels has been extensively studied, the mutual relationship between subway shield tunnel and overlying rock mountain under seismic loading has received little attention.
Earthquake ground motions are nonstationary in both time and frequency domains. Temporal nonstationarity refers to the variation in the intensity of ground motion in time, whereas spectral nonstationarity is the variation in the frequency content of the motion in time (the predominant frequency of earthquake ground motion decreases gradually within the duration) (Rezaeian and Der Kiureghian, 2010). Fully nonstationary (FN) ground motion exhibits both temporal and spectral nonstationarity, whereas quasistationary (QS) motion exhibits temporal nonstationarity only. The variation of the frequency of ground motion may have an underlying influence on the characteristics of the seismic response of underground tunnels and the interaction between tunnels and adjacent structures or mountains. However, the impact of the spectral nonstationarity of seismic motion on the seismic response and damage characteristics of tunnels and underground-tunneloverlying-rock-mountain systems has not been explored.
Herein, to investigate the mutual relationship between subway tunnels and overlying mountains and the impact of spectral nonstationarity of ground motion on the seismic response of underground tunnels, we conducted a numerical study based on a practical subway tunnel project in China. FN and QS ground motions were simulated. Several refined numerical models were established, considering the dynamic behavior of soil and the degradation of concrete due to damage. Parametric analyses were conducted to clarify the characteristics of the dynamic interaction between underground tunnels and overlying mountains and the difference between the results under FN and QS motions.

Engineering background
The 361.314-m interval tunnel project between the Simenkou Station and Pengliuyang Road Station of railway transit line 5 in Wuhan, China, was investigated as a case study. The external and inner diameters of the lining are 6 and 5.4 m, respectively. C60 concrete was used to construct the lining. The location of the interval tunnel is shown in Fig. 1. The interval tunnel passes under several buildings with multiple floors, as well as the Snake Mountain where the Yellow Crane Tower is located. Fig. 2 shows the geological conditions of the region. Double tunnels are below the horizon, surrounded by hard rock, and directly beneath the Snake Mountain. The excavation was performed using millisecond-delay blasting and static crush methods.

Simulation of bedrock motion
Because there is hardly a ground motion record near the engineering site, artificial bedrock motion was simulated. Besides, to demonstrate the influence of the spectral nonstationarity of seismic motion on the dynamic response of tunnels and rock mountains, QS and FN bedrock motions were generated, respectively. Recently, artificial synthesis technology of ground motion has developed to a large extent (Deodatis, 1996a(Deodatis, , 1996b 1,2,... where   denotes the bandwidth, and N is the frequency interval. ( , )  S t is the evolutionary power spectral density function (PSDF) at the base rock, u  is the upper cutoff frequency, and  y is a random phase angle uniformly distributed in [0, 2π]. The evolutionary Clough-Penzien acceleration spectrum (Clough and Penzien, 1993) was employed to model the PSDF of FN seismic motion, which has the following form: Here, t1, t2, and c are set to 6, 10, and 0.5, respectively; max a denotes the PGA value;  is the peak factor and set to 2.8 (Seya et al., 1993). The sampling frequency and upper cutoff frequency are set to 100 and 25 Hz, respectively. Based on the above-mentioned method and models, FN seismic motion could be simulated. Notably, the simulation of QS bedrock motion is a special case of FN bedrock motion, and QS bedrock motion can be derived by combining PSDF that is invariant with time with the classical SRM. Besides, the same intensity envelope function () at , PGA value max a , and random phase angle  y were employed to simulate FN and QS seismic motions. The duration of the simulated motions is set to 20 s, and the time interval to 0.01 s.

Setup of study cases and corresponding models
To observe the dynamic interaction between a tunnel and an overlying rock mountain, the following three models were established: i) Free-field (FF) model: rock mountain is included in this model, but tunnels are not. This model was set as a control to study the influence of the existence of tunnels on the seismic response of a rock mountain. It is depicted in Fig. 3(a).
ii) Tunnel-rock-mountain interaction system (TRM model): both rock mountain and underground tunnels are considered. This model is used to reveal the mutual relationship between underground tunnels and rock mountains. It is depicted in Fig. 3 iii) Tunnel-rock system (TR model): this consists of tunnels and surrounding rock only (removing rock mountain from TRM model). This model was set as a control to study the influence of overlying rock mountains on the seismic response of underground tunnels. It is depicted in Fig. 3   The above-mentioned models were established using the Abaqus/Explicit platform. The external and inner diameters of the lining are 6 and 5.4 m, respectively. C60 concrete was used to construct the lining. The site model is 700 m long and 68 m thick. According to the geological survey report, four soil or rock types were considered (Fig. 3). Four-node quadrilateral elements (CPE4R) were employed to discrete the model (including the tunnel and site). Each tunnel model consists of two layers of elements and 64 elements. According to Karabalis and Beskos (1997), Hatzigeorgiou and Beskos (2010), and Kuhlemeyer and Lysmer (1973), eight to ten finite elements per wavelength are adequate to accurately describe the elastic wave propagation in solids. To minimize the computational cost, eight elements per wavelength were assumed and employed (i.e., lmax ≤ Vs/8fmax, where lmax is the maximum element length, Vs denotes the shear wave velocity, and fmax is the maximum frequency that can be accurately transmitted, and it was set to 25 Hz in this study). Besides, a maximum analysis time step of 1 × 10 −5 s was set. Dynamic contact (Zhuang et al., 2015) was incorporated into the model to account for the slip, gapping, and contact on the interface between rock and lining. In the normal direction of the interface, the normal contact compressive stress mutually transfers through the contact constraint. The element nodes on the surfaces satisfy Hooke's law and the harmonized condition of displacement. If the interaction surfaces of rock and lining are separated, the contact constraint vanishes, and the contact boundary condition on the interface reduces to a common free boundary. The tangential contact shear stress (TCSS) can also be transferred on the interface. If TCSS is greater than the critical shear stress crit  , slip occurs. Coulomb's friction law is employed to describe the tangential mechanics, expressed as follows: where c is the friction coefficient of the interface between the soil and the concrete (0.4 in this paper) and P is the normal contact stress at the interface.

Constitutive model
The Davidenkov viscoelastic dynamic constitutive model, developed by Martin and Seed (1982) and improved by Zhao et al. (2017), is employed herein to describe the dynamic characteristics of soil, expressed as follows: where τ and γ are the shear stress and strain, respectively; Gmax is the initial shear modulus; A, B, and γ0 are dimensionless fitting parameters; na is the coefficient controlling the scale of the hysteresis loop; τc and γc denote the shear stress and strain at the last stress reversal point, respectively; τex and γex are the shear stress and strain at the last extreme value point, respectively. The fitting parameters were obtained through resonant column tests on soil samples in the engineering site. The fundamental material properties and fitting parameters for soils or rocks are listed in Table 2. The concrete damage-plasticity constitutive model proposed by Jeeho and Fenves (1998) was adopted to describe the crack damage and degradation of concrete. Material property parameters for the lining are listed in Tables 3 and 4. Rayleigh damping was employed, and the damping ratio  was set to 0.05 for the whole model. Then, the natural frequencies of the tunnel model were calculated, and the first-and third-order natural frequencies were used to calculate the mass and stiffness proportional damping factors and , respectively, based on the following equation: where 1 s  denotes the fundamental frequency of the site.

Boundary conditions and loading
The viscous-elastic artificial boundary ( where K denotes the spring stiffness, and C denotes damping. Subscripts N and T denote the normal and tangential directions of the boundary plane, respectively, and subscript l denotes a certain element node on the boundary; G, ρm, cs, and cp are the mechanical parameters of the medium, which denote the shear modulus, mass density, shear wave velocity, and compression wave velocity, respectively; R is the distance between the scatting wave source and the boundary nodes; Al is the total truncated boundary length of all elements containing boundary node l; X and Y are modified coefficients (herein, 0.8 and 1.1, respectively). One set of FN and QS bedrock motions with PGA = 0.3 g was simulated, as shown in Fig. 4, and the corresponding displacement time-histories, cumulative aria intensity (CAI), and normalized timefrequency spectra (the normalization was conducted by dividing the Fourier amplitudes at one time instant by the maximum value at the time instant) derived using short-time Fourier transformation are plotted in the same figure. The waveforms for both FN and QS motions are similar because the same initial values of PSDF, intensity envelop, and random phase angles were used to simulate the two motions. However, the lower frequency content of the FN bedrock motion tends to be dominant as time goes by, whereas that of the QS one is invariant with time. Therefore, the discrepancy between the displacement time-histories and CAI for the two motions gradually increases with time. These bedrock motions are inputted as SV waves. The equivalent node force derived from these motions was finally inputted into the model to calculate the seismic response (Zhang et al., 2016).

Results and discussion
The above-mentioned numerical model was used to simulate the seismic response of the tunnelrock-mountain dynamic system, and the results are discussed herein. The influence of underground tunnels on the seismic response at a mountaintop is presented in Section 5.1, and the location of the observation point is shown in Fig. 5(a). In Section 5.2, the impact of overlying mountains on the seismic response of underground tunnels is discussed. The right-line tunnel is selected as the study target, and the location of the observation points is shown in Fig. 5(b). Besides, the effect of the spectral nonstationarity of bedrock motion is discussed.

Effect of underground tunnels on the overlying mountain
The acceleration responses at observation point F in the case with the FF model are compared with those of the TRM model (Fig. 6). For all cases, the acceleration responses for both models are approximately the same, indicating that the presence of an underground tunnel does not affect the seismic response at the mountaintop. This might be because the wavelength of seismic motion is much greater than the dimension of the tunnel, and thus, the input wave is not sensitive to the tunnel. More importantly, for all cases regarding the two models, the acceleration response under FN motion is higher than that under QS motion, even though the input PGAs and waveforms for both motions are nearly the same, which has not been considered before. This is attributed to the difference in the prescribed PSDF between the two motions. The PSDF of FN motion varies with both time and frequency, whereas that of the QS motion varies with only frequency. Due to the decrease in the predominant frequency over time, FN motion exhibits a more significant low-frequency content and higher energy (Fig. 4(d)).

Effect of overlying mountains on underground tunnels
The  Figs. 10-14 compare the radial shear-stress response of the tunnel lining for TR and TRM models at observation points A-E, respectively. The initial values of the radial shear stress in the two models are different. More importantly, for all observation points and both input motions, the amplitude of variation of radial shear-stress responses in the TRM model is higher than that in the TR model under seismic waves. Moreover, the peak shear stress increases faster in the TEM model than in the TR model as input PGA increases, especially at observation points B (shoulder) and D (hance) (Fig. 15). The difference in the peak and initial values for each study case is plotted to eliminate the influence of different initial values among cases and facilitate comparison. The tunnel in the TRM model is subjected to a higher confining pressure than that in the TR model, and higher external confining pressure may cause higher shear force on the cross-section of the lining  [10][11][12][13][14]. Thus, under identical input motions, the shear stress for the TRM model is higher. Moreover, under earthquakes, the ovaling deformation and vibration of tunnels cause damage and cracks, which lower the cross-section area bearing the shear force, and the shear stress on the tunnel induced by confining pressure increases. When a seismic wave causes cracks, high confining pressure facilitates crack propagation, which further lowers the cross-sectional area bearing the shear force. Higher input PGA may cause more and deeper cracks and a higher shear-stress response of the tunnel. Meanwhile, the effect of confining pressure on crack propagation is more notable. In other words, the confining pressure gives rise to a second-order effect on crack propagation, increasing shear stress: the more severe the cracks caused by seismic waves (in other words, the higher the input PGA), the more significant the impact of confining pressure on the crack extension. In this process, higher confining pressure may cause a more remarkable second-order effect. Therefore, owing to higher confining pressure, the peak shear stress in the TRM model increases 'faster' than that in the TR model as PGA increases. In addition, the stress response of a tunnel under FN motion is higher than that of the tunnel under QS motion. This is explained in Section 5.1.
The results provide a reference for practical underground tunnel projects beneath rock mountains or similar projects. First, after the occurrence of cracks or damage, high confining pressure may significantly affect the tunnel; thus, adequate anti-crack measures should be taken for underground tunnels beneath rock mountains. Second, the shear strength of an underground tunnel (especially the shoulder and hance) beneath a rock mountain should be adequately reinforced to resist severe earthquakes. Moreover, the code for seismic design of buildings (GB 50011-2010) specifies that both earthquake records and artificial seismic motions should be considered in the time-history analysis of structures. Therefore, considering the results herein, FN seismic motions should be employed in the design of tunnel projects when artificial motions are used to facilitate the design of the tunnel.
The tensile damage factors (Dt) at observation points A-E for TR and TRM models under a PGA of 0.3 g are plotted in Fig. 16. In general, damage to the tunnel lining is not severe due to the surrounding hard rock. Dt is low at points A, C, and E but relatively high at points B and D, indicating that the shoulder and hance of a tunnel are more susceptible to earthquakes. More importantly, Dt in the TRM model is higher than that in the TR model, attributed to the combined effect of confining pressure and seismic loading. Moreover, spectral nonstationarity causes greater damage, and this tendency becomes more remarkable with time. As shown in Fig. 4, the discrepancies in the displacement amplitude and CAI between QS and FN motions increase gradually with time. Compared with QS motion, FN motion gradually increases the deformation and strain in the tunnel. Therefore, damage under FN motion is higher than that under QS motion.

Discussion
To justify the reliability of the obtained results, several samples of FN and QS motions (PGA = 0.2g) with the same phase angle were further fabricated and compared, and the relative displacement (herein, the displacement under FN motion minus that under QS motion) between FN and QS motions were observed. For conciseness, herein, nine sets of FN and QS motions are given, and the relative displacement for each set is plotted in Fig. 17. For the nine sets of motions, displacement under FN motion is higher than that under QS motion. The trend shown in Fig. 17 is consistent with that in Fig.  4(c). Besides, the seismic response of underground structures highly correlates with the displacement and deformation of the surrounding soil. Thus, FN motion would give rise to a higher displacement response and, thus, greater deformation and shear stress in a tunnel. The difference in the displacement amplitudes between FN and QS motions increases with time. Thus, damage in a tunnel under FN motion would be greater than that under QS motion, which is consistent with the tendency shown in Fig. 16. Therefore, the results of this study can effectively reflect the general regularity.

Conclusions
In this study, taking a practical subway tunnel project as a case study, we investigated the seismic response characteristics of underground-tunnel-overlying-rock-mountain interaction systems. QS and FN bedrock motions were considered. Several numerical models were established, and several study cases were set to investigate the mutual relationship between an overlying rock mountain and an underground tunnel under seismic loading and evaluate the influence of the spectral nonstationarity of seismic motion on the seismic response of a mountaintop and a tunnel. The main conclusions are as follows: i) An overlying rock mountain can suppress the acceleration response of an underground tunnel. Combined with the seismic motion, high confining pressure owing to the existence of an overlying mountain increases the radial shear-stress response and damage of the tunnel. An overlying rock mountain can also change the frequency content of the acceleration response of the tunnel lining.
ii) The spectral nonstationarity of input motions can result in a higher acceleration response to the mountaintop and tunnel and cause greater radial shear stress and damage to the tunnel lining, though the input PGAs and waveforms for QS and FN motions are the same. Therefore, artificial seismic motions exhibiting spectral nonstationarity may provide conservative results.