Abstract
The southeast region of Tibet in China is prone to earthquakes at high seismic levels. Since the construction of tunnels in this region faces a seismic hazard, it is necessary to explore the seismic responses of related construction areas there. In this study, a high-precision finite model of a mountain in a construction area was established based on a deep-buried tunnel in Southeast Tibet. Core specimens were obtained through on-site drilling, with the initial tectonic stress in the tunnel area measured by using the hydraulic fracturing method. Dynamic mechanical parameters of the rocks were acquired through in-situ testing. Moreover, indoor tests were conducted to obtain the physical and mechanical properties of the rocks. The tectonic stress, identified through multiple linear regression inversions, was used as the initial condition for dynamic calculation. The dynamic responses of the tunnel-engineering area under seismic action were calculated, followed by an analysis of the dynamic response characteristics. In addition to dramatic fluctuations in mountain topography, complex dynamic response characteristics of the mountain structure in earthquake environments were also revealed. Although the displacement time-history curve of each observation point displayed a consistent fluctuating pattern, the dynamic displacement extremes in different areas of the mountain varied significantly when the time arrived for their appearance. A maximum dynamic displacement of 0.53 m was observed at the observation point with the highest elevation (5,265.8 m). The acceleration amplification effect inside the mountain body was obvious, with the acceleration increasing as the elevation rose. The time for the occurrence of PGA (acceleration extreme) of each observation point in different directions differed substantially. The extreme value of dynamic stress along the tunnel axis conformed to the distribution pattern of the initial tectonic stress. Compared to the initial tectonic stress, the dynamic stress increased significantly at all observation points along the tunnel axis, with the dynamic stress increasing at a larger rate in the areas at a greater buried depth. Horizontal dynamic stress Sx at the maximum buried depth and Sy increased to 13.86 MPa and 10.32 MPa, respectively, and vertical dynamic stress Sz increased to 4.30 MPa. Variations in horizontal stress became more apparent as the elevation rose, in contrast to weakened variations in vertical stress as the elevation dropped. This study aims to explore the dynamic responses of mountains and the variation law of tectonic stress in deep-buried tunnels in earthquake environments in Tibet through qualitative analysis. The results revealed drastic changes in acceleration, dynamic displacement, and ground stress in the tunnel-engineering area. It is hoped that the contents of the study will contribute to the prevention and control of seismic geological hazards in the tunnel-engineering region and render technical support for the design and construction of similar projects in this region.
1. Introduction
Due to active neo-tectonic movement in Southeast Tibet of China and intense seismicity in the Himalayan earthquake fault zone, geological disasters occur frequently and distribute widely in this region. Severe earthquakes trigger a large amount of secondary geological disasters, such as rockfalls and landslides that will drastically threaten the safety of people’s lives and property and the efficiency of engineering constructions in this region [1]. The tectonic stress field and the secondary stress field, which came into being after excavation, are considered the most important factors that affect the static and dynamic stability of tunnels [2]. Because of the construction of many deep-buried tunnels, the tunnelling process is inevitably subjected to the threat of earthquake disasters during the construction cycle that lasts for several years [3]. Therefore, it will be of great significance to analyse the seismic dynamic response characteristics of the mountains in the tunnel-engineering area to predict earthquake disasters.
Tunnels are deeply buried underground. There are two major approaches for the research of geodynamic engineering: the shaking table test and the numerical analysis. With the shaking table test, the motion state and the failure process of the test model under the action of an earthquake are recorded through the arrangement of sensors. Since it reflects the dynamic response characteristics of the rocks under actual earthquake load in an accurate manner, it is widely used in the research study of the dynamic response patterns of rock masses [4–10]. Currently, the vibration test for the structure of rocks is limited to simplifying the rocky slope, so it is unable to measure the limit of tectonic stress. In addition, it is difficult to reproduce this type of test owing to long testing cycles and high costs. On the other hand, numerical simulation has been widely used to investigate the dynamic responses of rocks. Dong et al. [11] calculated and simulated the seismic responses of a tunnel under different conditions, for instance, varying construction environments, buried depths, and cross-sections, to draw the conclusion that the peak acceleration is the key factor affecting the seismic performance of the tunnel. Relying on nonlinear finite element simulations of underground structures, Sun and Dias [12] demonstrated that the tectonic stress environment has a significant effect on the seismic responses of underground structures. With the help of finite difference software, Ding and Liu [13] simulated and calculated the seismic responses of a tunnel under a static condition. The effects of buried depth, diameter, and lateral pressure coefficients of tectonic stress on the stress and deformation of surrounding rocks were also acquired in their experiment. According to the results, the maximum shear stress, the radius of the plastic zone, and the maximum displacement of surrounding rocks increase as the diameter of the tunnel enlarges. Li and Chen [14] concluded that the increase in the buried depth of underground chambers changes the structural resonant frequency and affects the structural dynamic response significantly. In their study, the dynamic stability of the model slope and the influences of ground shaking parameters on earthquake-related dynamic responses were also analysed [15]. Song et al. [16] studied the influences of rock structures on wave propagation characteristics and the dynamic amplification effect by analysing the dynamic acceleration amplification coefficient of rocky slopes. Zhou et al. [17] combined the vibration test and numerical simulation to explore the dynamic responses and damage patterns of tunnels constructed under different geological conditions. The results revealed reduced vibration due to the damping effect of vibration isolators, as well as apparently weaker dynamic responses of the tunnels in fault areas than those of the tunnels with normal surrounding rocks.
Currently, the fundamental research on the dynamic responses of underground engineering mainly focuses on rocky slopes and shallow-buried underground cavern enclosures. Numerical calculation of local tunnel areas is usually performed as a comparative study for the shaking table test [18], which fails to reflect the effects of tectonic stress on dynamic responses. Under the action of high tectonic stress, there is a higher risk of rock burst and large deformation in the surrounding rocks of deep-buried railroad tunnels. Since the tectonic stress effect on shallow-buried areas is weak, the research results in this regard have little reference value for tunnel projects at a greater buried depth. The experimental tunnel is buried deep under a high mountain and is strategically valuable in economic terms. Therefore, it is worthwhile to explore the dynamic seismic response characteristics of this tunnel-engineering area for the prevention and control of seismic hazards. The horizontal principal stress of the boreholes of the tunnel, which is 1,200 m deep and typical of highly representative engineering characteristics, was measured at 42 MPa. In this study, a high-precision finite element model of the mountain was constructed in consideration of the actual conditions of the tunnel project. The tectonic stress boundary was applied to the model based on the data of the inverse performance of the measured structural stress to calculate and analyse the dynamic response law of the tunnel-engineering area under seismic effects. The research study on the dynamic response characteristics of the tunnel under seismic conditions and the law of tectonic stress distribution will be of great significance when it comes to the design of tunnel alignment, control of the stability of surrounding rocks, and the operation and maintenance of tunnels.
2. Project Introduction
2.1. Tunnel
The experimental section of the deep buried tunnel in the southeast region of Tibet is mainly subjected to influences of the warm and humid airflow from the Indian Ocean and the southeast monsoon. The warm and humid airflow from the Indian Ocean goes upstream along the Yarlung Tsangpo River and its tributaries before its entry into the southeast inland region of the Tibetan Plateau, after which the strength of the airflow gradually decreases from downstream to upstream. The elevation of the track surface of the tunnel stands between 2,100 and 2,800 m, with the tunnel cavity falling into the categories of the plateau subtropical climate zone and the plateau temperate climate zone. The multiyear average temperature and precipitation are recorded at 13°C and 1,276.0 mm, respectively. The tunnel, deeply buried in southeast region of Tibet, has a total length of 13.8 km and a maximum buried depth exceeding 1,900 m. Passing through a vast mountainous and canyon area, the tunnel witnesses a fault zone intersecting with the line axis at a large angle. The rock stratum in the experimental area is mainly composed of gneiss, while the fault zone is mainly composed of fault breccia.
Figure 1 presents the time-history curve of seismic horizontal acceleration. Artificial seismic waves were synthesized through the trigonometric series’ superposition method. The characteristic period of the basic ground-motion acceleration response spectrum (0.45 s), the peak acceleration of basic ground (0.2 g), and the acceleration response spectrum and time envelope function of the tunnel-engineering area were provided with the ground shaking parameter zoning map of China as the target reference. The control error of the target response spectrum fitting was below 5%.
2.2. Physical and Mechanical Tests
Cores were drilled in the tunnel cross-hole, followed by the fetch of standard specimens with a diameter of 50 mm and a height of 100 mm from the cores drilled pursuant to the Standards for Test Method of Engineering Rock Masses of China. Rock density ρ was measured via the volumetric method, and the static modulus of elasticity Es and static Poisson’s ratio were detected through the uniaxial compression and deformation test, as shown in Table 1. The tunnel crosses with a fault fracture zone that is more than 80 m in width. Since it is impossible to drill cores in the fault fracture zone, Es was set at 1/4 as that of gneiss and as 0.3 in the fault fracture zone.
Cross-hole wave velocity was used to test the dynamic parameters of the rock mass. The test scheme was carried out as follows: Three holes in the cross-hole were drilled in a straight line at equal spacing, with a hole at the end excited while the other two holes received. Identical elevation of the excitation and reception points was guaranteed, and the test was synchronized from the bottom up.
The measured longitudinal wave velocity of gneiss (VP) ranged between 4,348 and 5,921 m/s with an average value of 5,184 m/s, and the measured transverse wave velocity (VS) ranged between 2,342 and 3,158 m/s with an average value of 2,449 m/s. The dynamic shear modulus (Ed) was calculated based on equation (1), with an average value of 58.5 GPa. The dynamic Poisson’s ratio () was calculated through equation (2), with an average value of 0.26. The Ed of the fault fracture zone was set 1/4 as that of gneiss. Since had no reasonable value for reference, .
2.3. Finite Element Model
For accurate calculation of the seismic dynamic response of the tunnel-engineering area, the construction of a high-precision 3D finite element model of the entire area was required. Specifically, topographic data higher than 5 meters in accuracy were extracted to generate the contour data, followed by the construction of the model using geological modeling software SKUA-GOCAD and the output of S-Gird. As revealed by the geotechnical dynamic calculation, the mesh element size of the finite element model influenced the simulation accuracy significantly, so the element size should not exceed 1/8–1/10 of the wavelength based on the maximum frequency of the input waveform [19]. The basic grid cell size was calculated at 20 m. Figure 2 presents the finite element model, with the green area standing for gneiss and the grey area for the fault.
3. Seismic Calculations and Tectonic Stress Conditions
With the seismic time-history method, the earthquake was treated as the time process and the structure was simplified as a multidegree system. In addition, the seismic response equation of the earthquake wave input was directly calculated through step-by-step integration, so that the seismic response of the structure can be calculated at any time. Compared to the response spectrum method, the calculation results of the time-history method are more accurate, so it stands out for its favorable applicability as a solution to linear and nonlinear problems [20].
In the seismic geotechnical calculation of, the geotechnical hysteresis is realized through damping. In the calculation software, however, the damping is set in the form of Rayleigh damping [21]. The damping ratio ζ stood at 0.025, the first-order vibration frequency ω1 of the mountain model was 0.322, and the second-order formation frequency ω2 was 0.491. Damping coefficient: ; .
3.1. Initial Tectonic Stress Conditions
The tectonic stress in the project area was measured through the hydraulic fracturing method, with the drilling locations presented in Figure 2. The principal stress field in the measured area is dominated by the south, west, north, and east extrusions. The direction of horizontal principal stress, measured through the hydraulic fracturing method, had a horizontal angle with the coordinate system of the calculation model, so it was necessary to convert the data of measured ground stress into tectonic stress. The converted tectonic stress values are detailed in Table 2.
Based on the data of measured tectonic stress, the equation for the tectonic stress field of the tunnel-engineering area was obtained through inverse calculation based on the principle of multiple linear regression [22]. According to the details in (3), the constant term stands at 4.641 MPa, represents the stress per unit dead weight, is the tectonic stress produced by the unit compressive stress boundary under 10 MPa in direction x, and denotes the tectonic stress produced by the unit compressive stress boundary under 10 MPa in direction y. The complex correlation coefficient of the regression equation (R = 0.975) approximates 1.0, suggesting a favorable regression effect as
Figure 3 presents the comparison between the component value of measured tectonic stress and its regression calculation value at each borehole point. As revealed in the figure, the regression calculation values of the stress at various stress points fit well with the values measured. In addition, the general variation trend of the tectonic stress was consistent with the variations in depth, i.e., the value of the stress increased as the buried depth deepened.
(a)
(b)
(c)
(d)
To guarantee closer seismic calculation results to the actual conditions, the tectonic stress was input as the initial stress condition at the beginning, after which the dynamic calculation was performed.
3.2. Boundary Conditions for Dynamic Calculation
Dynamic load was applied at the bottom of the model, with the vertical dynamic load set as 2/3 of the horizontal dynamic load. Seismic load was input in a consistent manner. In the dynamic geotechnical calculation, the radiation damping of the foundation was taken into consideration to realize energy dissipation. The viscous artificial boundary was used for simulation [23, 24]. Since acceleration cannot be applied to the viscous boundary at the bottom of the model as a dynamic load directly, it was necessary to convert it into an equivalent nodal force, which was equal to the product of the stress and the bottom area of the bottom element. The calculation equations for stress time history are listed as follows:where and represent the normal and tangential stress, respectively, is the density of rock, and stand for the average values of the measured wave velocity, and and denote the normal and tangential velocity, respectively. The synthetic velocity time history was obtained by integrating the synthetic acceleration time history, and the linear correction was established to acquire the nondrifting velocity time course curve. Figure 4 presents the tangential velocity time-history curve. The normal velocity was 2/3 of tangential velocity.
4. Analysis of Calculation Results
Since the topography of the tunnel-engineering area fluctuated dramatically, the seismic dynamic response characteristics of the mountain under the influence of earthquakes are complex. The observation points on the longitudinal section of the tunnel axis were selected. Point A was located at the maximum buried depth of the tunnel axis, at an altitude of 2,481.8 m. Point B sat at the mountain apex, corresponding to Point A, with an altitude of 4,435.9 m. Point D, at the shallowest buried depth of the tunnel, was located at the bottom of the valley with an altitude of 2,378.6 m. As the highest point in the calculation model area, Point C had an altitude of 5,265.8 m. The dynamic displacement and acceleration response laws at the above four observation points were analysed in the experiment. Along the vertical line of the maximum buried depth, Points E, F, and G, which had an altitude of 2,052.9 m, 2,867.4 m, and 3,372.7 m, respectively, were selected as the observation points for the research of the dynamic response law. Details of the distribution of the observation points on the experimental section are shown in Figure 5.
4.1. Dynamic Displacement
Figure 6 presents the displacement time-history curves of each observation point in different directions under the action of an earthquake. According to the figure, the variations in displacement time-history curves of each point were consistent with the strengthening and attenuation trends of the input seismic wave, based on which the boundary conditions and damping settings of the calculation model were revealed as reasonable. Tables 3 and 4 and show the maximum positive and negative displacements of the four observation points and the time for their occurrence, respectively. The law of maximum partial displacement was expressed as . Both horizontal partial displacement and vertical partial displacement exhibited an increasing trend as the altitude rose. The maximum partial displacement in the same direction at different observation points differed from that in different directions at the same point. The time-history curves of partial displacement in different directions at varied observation points differed from each other as well. The complexity of the dynamic response of the mountain structure affected by the earthquake was also verified.
(a)
(b)
(c)
(d)
Weak correlations existed in the time-history curves of the horizontal and vertical partial displacements of the four observation points, because no substantial correlations were revealed among the four observation points. Figure 7 presents the time-history curves of the total displacement at each point, revealing the characteristic that the fluctuation of displacement tended to remain consistent. This reflected the unity of the dynamic effects. The total displacement extremum of each point was concentrated at the time node of 15.9 ± 0.05 s. The maximum total displacement of each point was calculated as = −0.44 m, = −0.35 m, = −0.53 m, and = −0.50 m, respectively. Since observation points B and C were on the surface of the mountain, a higher elevation would mean a larger fractional displacement in each direction, revealing the whip-sheath effect caused by the tectonic characteristics of the mountain. Observation points A and D were inside the rock, with a greater buried depth indicating a greater dynamic displacement. This reflected the influence of tectonic stress on dynamic displacement.
4.2. Acceleration
Acceleration is an important physical quantity that reflects the dynamic response characteristics of geotechnical systems. The amplification factor of peak ground acceleration (PGA) was defined as the ratio of PGA at each observation point in the PGA input at the bottom of the model. Figure 8 presents the acceleration time-history curves of the four observation points. As shown in the figure, the fluctuation characteristics of the acceleration time-history curves of the observation points were consistent with the strengthening and attenuation trends of the input vibration load, based on which the boundary conditions and damping settings were revealed as reasonable. The PGA values of observation point A, which had the maximum buried depth, were all smaller than those of the other three observation points, indicating the principal influence of tectonic stress on acceleration response.
(a)
(b)
(c)
(d)
The horizontal PGA input of the model stood at 2 m/s2, while the vertical PGA was 2/3 of the PGA in the horizontal direction. Table 5 shows the PGA values, PGA generation times, and PGA amplification factors of each observation point in different directions. The values of the PGA amplification coefficient of the three observation points A, B, and C exhibited an increasing trend with the rise of the elevation in general. Observation point D had the lowest elevation but a shallow buried depth. Due to the poor constraint of rock vibration and the relatively weak damping effect, the corresponding PGA amplification coefficient was greater than that of observation point A. The input seismic PGA generation time lasted for 7.14 s, but the PGA generation time of the observation points in different directions was all longer than 7.14 s, indicating a lagged dynamic response of the rock. Therefore, it was concluded that the calculation results were consistent with the general rule of seismic dynamic responses.
Observation points A and D had buried depths of 1,954.1 m and 108.0 m, respectively. The time for the occurrence of PGAx in observation point A was 3.4 s; while the time for the occurrence of PGAy and PGAz were 7.95 s and 2.1 s, respectively, all of which were earlier than that of observation point D. Observation point B was located on the top of the mountain at the maximum buried depth (A). The time for the appearance of PGAx at observation point A and the appearance of PGAy was 3.35 s and 8.5 s later, respectively, compared to those at observation point B, while the time for the appearance of PGAz was 1.05 s earlier than that at observation point B. The altitude of observation point C was 800 m higher than that at observation point B. The time for the appearance of PGAx at observation point B and the appearance of PGAy were 4.9 s and 0.05 s later, respectively, compared to that at observation point C, while the time for the appearance of PGAz was 1.41 s earlier than those at observation point C. In general, the variability of the time for the appearance of horizontal PGA was greater than that of vertical PGA, and the time for the appearance of PGA varied widely among the observation points. This indicated complex results from the dynamic response of the experimental mountain under seismic action. Moreover, the influence of the horizontal buried depth of the observation points should be taken into consideration, in addition to the impacts of the elevation and buried depth of the observation points on PGA.
4.3. Dynamic Stress
Tectonic stress, which is considered an important factor associated with the stability of surrounding rocks in the process of tunnelling, was obtained through multiple linear regression and taken as the initial boundary condition for seismic calculation. Figure 9 compares the maximum dynamic stress and the tectonic stress on the tunnel axis in the dynamic time-history calculation and the analysis steps. The horizontal dynamic stresses Sx and Sy and the vertical stress Sz increased prominently compared to the tectonic stress. The dynamic stress curve of the tunnel axis appeared abruptly but varied locally, and the area typical of abrupt variations was correlated to the relief of the terrain, indicating the influence of the mountain structure on dynamic stress. The maximum buried depth of observation point A was situated at tunnel mileage CK6+200. Figure 9 shows the horizontal dynamic stress history curves of Sx, Sy and Sz for the vertical stress of the observation points at different buried depth. The horizontal dynamic stress Sx at the maximum buried depth increased to 13.86 MPa, Sy increased to 10.32 MPa, and the vertical dynamic stress Sz increased to 4.30 MPa. The dynamic stress increased at a larger magnitude at the mileage with the maximum buried depth on the tunnel axis.
(a)
(b)
(c)
The variations in dynamic stress were smaller in the area that was shallow in buried depth, a result consistent with the tectonic stress distribution characteristics. The tunnel intersects with the fault at a comparatively large angle at mileage CK8+600, where the tectonic stress and the vertical dynamic stress experienced abrupt variations while the corresponding horizontal dynamic stress witnessed no obvious abrupt variations, which may be explained by the release of stress.
Figure 10 presents the time curve of dynamic stress variations for the four observation points E, A, F, and G, whose elevations rise sequentially. According to the figure, the fluctuation variation pattern of the curve remained consistent because the four observation points were on a vertical line. As a result, the vertical and horizontal depths of the buried sites varied in a more uniform pattern. Table 6 presents the maximum growth stress values and the time for their occurrence at each observation point. Horizontal stress Sx and Sy increased in the magnitude of maximum growth as the elevation rose, which was attributed to the fact that a higher elevation of the mountain structure causes a smaller lateral stiffness. As a result, the dynamic displacement and the acceleration under the influence of the earthquake exhibited an increasing trend, which would further affect the range of the variations in dynamic stress. The maximum growth of vertical stress Sz decreased as the altitude rose.
(a)
(b)
(c)
5. Results and Discussion
In this study, physical and mechanical parameters of rocks and rock masses were obtained through in-situ tests and indoor tests. A high-precision model of a mountain was constructed based on a deep buried tunnel in the southeast region of Tibet. The initial tectonic stress obtained through multiple linear regression inversion was set as the initial stress condition for dynamic calculation, and the seismic responses of the tunnel-engineering area under the action of a three-way earthquake were calculated and analysed.
The variations in time course curves of dynamic displacement and acceleration at each observation point under seismic action remained consistent with the strengthening and attenuation trends of the input seismic wave. The maximum horizontal displacement of each observation point was larger than the maximum vertical displacement. The fluctuations in the time-course curves of total displacement remained consistent. The total displacement at the observation point with the highest elevation was measured at -0.53 m, and the total displacement increased as the elevation rose in general. However, the river valley areas witnessed larger dynamic displacement due to their weak stiffness.
The PGA of the observation points appears after the input of PGA, reflecting the hysteresis of the dynamic response of the rock. The dynamic response of the mountain under a three-way earthquake is complex. The time for the appearance of PGA in different directions at different points differs substantially. Specifically, the time difference of the acceleration extremum at the highest altitude was less than 1.1 s, while the time difference of the acceleration extremum at the maximum buried depth reached more than 9.5 s. The tectonic stress has a critical influence on acceleration response, and the tectonic action is prominent in the area at a large buried depth. In addition, the corresponding PGA is smaller than that in the hilltop and valley areas typical of small tectonic stresses. The PGA amplification factor increases as the elevation rises. A larger amplification factor of PGA in the valley is attributed to a weaker damping effect.
The dynamic stress extremes in the tunnel axis are consistent with the distribution pattern of the tectonic stress, but are significantly larger than the initial tectonic stress. The increase in stress is larger in the areas at a greater buried depth. The tectonic stress increases by 23.6%, 23.9%, and 7.8%, respectively. The mountain, which is cone-shaped, witnesses decreasing stiffness and initial tectonic stress as the elevation rises. Nevertheless, horizontal dynamic stress extremes SX and SY increase at a larger rate as the elevation rises, while vertical dynamic stress extreme SZ increases at a slower rate as the elevation rises.
Under the action of earthquake, the mountain body generates violent acceleration and dynamic displacement that would cause geological disasters such as landslides, rock falls, and uneven settlement. In addition, the tectonic stress also exhibits violent fluctuations under such circumstances. Apparently, a sharp increase in rock stress will be detrimental to the stability of tunnels during excavation. Therefore, the adverse effects of earthquakes on engineering should be taken into consideration during the construction of tunnels. It is hoped that the results of this study can provide some reference value for the prevention and control of seismic hazards in the same type of tunnels in peripheral areas.
Data Availability
The raw/processed data required to reproduce these findings cannot be shared at this time as the data also form part of an ongoing study.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This study was financially supported by the Foundation of State Key Laboratory of Coal Combustion (FSKLCCA2013) and Natural Science Foundation of Tibet Autonomous Region (XZ202001ZR0025G).