Seismic Response Analysis of Hydraulic Tunnels Under the Combined Effects of Fault Dislocation and Non-Uniform Seismic Excitation

Abstract

Hydraulic tunnels are prone to pass through faults and high-intensity earthquake areas, which will cause serious damage under fault dislocation and earthquake action. Fault dislocation and seismic excitation are often considered separately in previous studies. For tectonic earthquakes with higher frequency in seismic phenomena, fault dislocation and ground motion are often associated, and fault dislocation is usually the cause of earthquake occurrence, so it is limiting to consider the two separately. Moreover, strong earthquake records show that there will be significant differences in the mainland vibration within 50 m. The uniform ground motion inputs in previous studies are not suitable for long hydraulic tunnels. This paper begins with the simulation of non-uniform stochastic seismic excitations that consider spatial correlation. Based on stochastic vibration theory, multiple multi-point acceleration time-history curves that can reflect traveling wave effects, coherence effects, attenuation effects, and non-stationary characteristics are synthesized. Furthermore, a fault velocity function is introduced to account for the velocity effect of fault dislocation. Finally, numerical analyses of the response patterns of the tunnel lining under four different conditions are conducted based on an actual engineering project. The results indicate the following: (a) the maximum lining response values occur under the combined effects of fault dislocation and non-uniform seismic excitation, indicating its importance in the seismic resistance of the tunnel. (b) Compared to uniform seismic excitation, the peak displacement of the tunnel under non-uniform seismic excitation increases by up to 6.42%, and the peak maximum principal stress increases by up to 28%. Additionally, longer tunnels exhibit a noticeable delay effect in axial deformation during an earthquake. (c) Under non-uniform seismic excitation, the larger the fault dislocation magnitude, the greater the peak displacement and peak maximum principal stress at the monitoring points of the lining. The simulation results show that the extreme response values primarily occur at the crown and haunches of the tunnel, which require special attention. The research can provide valuable references for the seismic design of cross-fault tunnels.

1. Introduction

Underground hydraulic tunnels are an important component of modern social infrastructure [1]. Compared to the surface structure, the underground structure is considered to have good seismic performance because it is constrained by the surrounding rock mass or soil mass. However, seismic damage to underground structures has occurred frequently in recent decades, especially during the 1995 Hanshin earthquake in Japan [2], 1999 Jiji earthquake in Taiwan [3], and 2008 Wenchuan earthquake [4,5]; many underground tunnels have been irrecoverably damaged. Hydraulic tunnels generally have a long line and inevitably pass through fault fracture zones and areas with high seismic intensity. Earthquake damage records indicate that the safety of cross-fault hydraulic tunnels located in seismically active areas is still an important issue.

Over the years, domestic and foreign scholars have mainly carried out three works to solve the above problems. Firstly, studies have been carried out on the response and failure mechanisms of cross-fault tunnels under fault dislocation. Yu et al. [6] combined field investigation and the Finite Element Method (FEM) to evaluate the failure mechanism of the Longxi tunnel crossing the reverse fault. Cui et al. [7] conducted laboratory model tests to compare the effects of several tunnel anti-slip measures under the action of fault viscous slip dislocation. Zhong et al. [8] established a numerical model of a water delivery tunnel crossing multiple strike-slip faults to assess the structural damage. An et al. [9] analyzed the internal forces and deformation characteristics of the tunnel under the action of fault dislocation. Wang et al. [10] explored the deformation and failure mechanism of tunnels when faults were dislocated through large-scale model tests. The second is to study the effect of earthquakes on underground tunnels. Ding et al. [11] conducted dynamic response analysis of an immersed tube tunnel based on the large-scale seismic response simulation method. Anastasopoulos [12] conducted the nonlinear response analysis of a deep tunnel under strong earthquake action. Jiang et al. [13] analyzed the seismic performance of the tunnel model under the condition of uniform seismic input by combining test and numerical simulation. Wang et al. [14] studied the seismic response of the underground tunnel–soil–surface structure interaction system. Sun et al. [15] emphasized the importance of the duration of ground motion for hydraulic tunnels under deep seismic action and found that the longer the duration of ground motion, the higher the risk of hydraulic tunnels collapsing. Thirdly, the seismic response and damage characteristics of cross-fault tunnels are analyzed. Jiao et al. [16] conducted a numerical study on the nonlinear response and failure process of tunnels crossing inactive faults. Zhou et al. [17] investigated the failure modes of shallow buried elliptical tunnels under different geological conditions, finding significant differences between the dynamic responses of tunnels crossing faults and those of tunnels in normal surrounding rock. Sun et al. [18] analyzed the seismic vulnerability of a cross-fault hydraulic tunnel, and then, ref. [19] improved the seismic performance evaluation index of cross-fault hydraulic tunnels based on the partial least squares method.

Currently, the study of tunnel responses under fault dislocation only considers the effect of fault dislocation momentum on the tunnel, neglecting the effect of the fault slip rate on the tunnels. With the increasing demand for water conveyance and the continuous improvement in engineering technology, the axial length of hydraulic tunnels has become increasingly longer. However, the majority of studies on the deformation and force characteristics of tunnels under seismic action have focused on uniform ground motion inputs. Further observations from seismic damage records [20,21] indicate that, unlike shorter tunnels, the damage patterns of longer hydraulic tunnels are unique. The damage locations and types along the tunnel axis are not consistent, exhibiting spatial non-uniformity. These unique and complex damage phenomena are the result of the spatial effects of seismic excitations. Uniform input conditions are insufficient for the seismic analysis of long hydraulic tunnels. Therefore, it is necessary to adopt non-uniform seismic excitation for a more accurate analysis. Some studies have considered oblique incidence [22,23,24,25,26] as an input condition for non-uniform ground motion, but only the travelling wave effect of ground motion is considered in oblique incidence. However, the spatial variability of ground motion that has a significant impact on spatial large-size structures [27,28,29] includes not only the traveling wave effect, but also the attenuation effect, hysteresis coherence effect, and local site effect, and the input conditions of uniform ground motion or oblique incidence are not comprehensive. In recent years, the synthesis method for simulating non-uniform seismic excitations, which can obtain the spectral characteristics at each excitation point and generate multi-point acceleration time histories through a certain rule, has become a focal point. This method not only effectively simulates the spatial variability of seismic excitations but also conveniently models randomness and nonlinearity. Zhang et al. [30] analyze the nonlinear behavior of a concrete gravity dam under near-fault pulse records, and non-pulse records are investigated with consideration of the obliquely incident P waves. Liu et al. [31] proposed a dimensionality reduction method for simulating fully non-stationary seismic random fields by integrating probability density evolution theory, enabling detailed dynamic response analysis of bridge structures. However, these synthesized seismic excitations in previous studies did not account for the weak properties of faults and the effects of fault dislocation. To sum up, the current study only considers the effects of fault dislocation and earthquakes on one side, without considering the effects of both on the tunnel. In reality, earthquakes are often caused by fault activities, and fault dislocation and ground motion usually do not occur separately.

In this paper, we present a three-dimensional numerical model of a hydraulic tunnel grounded in practical applications. Employing the random ground motion synthesis technique, we generate non-uniform ground motions that effectively capture the traveling wave effect, attenuation effect, coherence effect, and non-stationary characteristics of ground motion. We extensively investigate the dynamic response of the hydraulic tunnel when subjected to the combined influence of fault dislocation and non-uniform ground motion. A detailed illustration of the research framework is provided in Figure 1.

2. Spatial Incongruent Ground Motion Synthesis Method and Fault Velocity Function

2.1. Simulation of Non-Uniform Ground Motion Based on Random Vibration Theory

Housner [32] first realized that ground motion is a random ground process and proposed a stationary white noise model of ground motion. However, this model was only discussed from the perspective of pure mathematics and could not reflect the physical nature of ground motion. Japanese scholars [33] proposed a white noise filtered Kanai-Tajimi model that can better describe seismic frequency characteristics, but this model cannot deal with the response problem of low-frequency structures. Many scholars modified the Kanai-Tajimi model to eliminate its low-frequency component. In this paper, the modified Kanai-Tajimi model [34] is adopted as the ground motion self-power spectrum:

(1)$S(\omega )={\displaystyle \frac{{\omega }^6}{{\omega }^6+{\omega }_c^6}}{\displaystyle \frac{1+4{\xi }_g^2\frac{{\omega }^2}{{\omega }_g^2}}{{(1-\frac{{\omega }^2}{{\omega }_g^2})}^2+4{\xi }_g^2\frac{{\omega }^2}{{\omega }_g^2}}}{S}_0$

where ${\omega }_k,{\xi }_g,S_0,{\omega }_c$ are the improved spectral parameters of the Kanai-Tajimi model [35].

The ground motion power spectrum matrix can be constructed based on the random vibration theory. For the ground motion power spectrum matrix S_i(ω_k) with a frequency of ω_k, the equation is as follows:

(2)$S(i{\omega }_k)=\left[\begin{array}{cccc}{S}_1({\omega }_k)& S_{12}(i{\omega }_k)& \cdots & S_{1n}(i{\omega }_k)\\ S_{21}(i{\omega }_k)& S_{22}({\omega }_k)& \dots & S_{2n}(i{\omega }_k)\\ ⋮& ⋮& \ddots & ⋮\\ S_{n1}(i{\omega }_k)& S_{n2}(i{\omega }_k)& \dots & S_n({\omega }_k)\end{array}\right]$

In the matrix, the cross-power spectrum between n and m points on a straight line along the propagation direction of the seismic wave can be expressed by the following equation [35]:

(3)$S_{nm}(i{\omega }_k)=\sqrt{S_n({\omega }_k)S_m({\omega }_k)}{\rho }_{nm}({\omega }_k,d_{nm})e^{-i\omega {\displaystyle \frac{d_{nm}}{v_a({\omega }_k)}}}$

where $S_n({\omega }_k)$ and $S_m({\omega }_k)$ are the auto-power spectra of two points n and m, and d_nm is the distance between two points. v_a(ω_k) is the apparent velocity of the seismic wave, and ${\rho }_{nm}({\omega }_k,d_{nm})$ is the coherence function of two points.

Considering the correlation between the apparent velocity and frequency of seismic waves, $v_a({\omega }_k)$ is solved according to the following formula [36]:

(4)$v_a({\omega }_k)=c_1+c_2\ln ({\displaystyle \frac{{\omega }_k}{2\pi }})$

where c₁ and c₂ are fitting parameters: c₁ = 3344; c₂ = 1095.

The coherence function adopted in this study is the spatial coherence function model proposed by Qu et al. [36], which comprehensively considers multiple coherence models and has a wide range of applicability:

(5)${\rho }_{nm}({\omega }_k,d_{nm})=\left\{\begin{array}{l}{e}^{-a\left({\omega }_k\right)d^{b\left({\omega }_k\right)}}\\ a({\omega }_k)=0.00001678{\omega }_k^2+0.001219\\ b({\omega }_k)=-0.0055{\omega }_k^2+0.7674\end{array}\right.$

The power spectrum matrix $S_i({\omega }_k)$ is a Herrmite matrix and a positive definite matrix. According to matrix theory, $S_i({\omega }_k)$ can be decomposed into the product of a lower triangular matrix and an upper triangular matrix using the Cholesky [31] method:

(6)$S(i{\omega }_k)=\left[\begin{array}{cccc}{l}_{11}({\omega }_k)& 0& \cdots & 0\\ l_{21}(i{\omega }_k)& l_{22}({\omega }_k)& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ l_{n1}(i{\omega }_k)& l_{n2}({\omega }_k)& \dots & l_{nn}({\omega }_k)\end{array}\right]\times \left[\begin{array}{cccc}{l}_{11}({\omega }_k)& {l^{*}}_{21}({\omega }_k)& \cdots & {l^{*}}_{n1}(i{\omega }_k)\\ 0& l_{22}({\omega }_k)& \dots & {l^{*}}_{n2}(i{\omega }_k)\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & l_{nn}(i{\omega }_k)\end{array}\right]$

Using the decomposition term of $S_i\left({\omega }_k\right)$, the amplitude $A_{im}({\omega }_k)$ and phase angle ${\theta }_{im}({\omega }_k)$ under the frequency component ${\omega }_k$ can be obtained, and the ground motion is expressed as the sum of n-term trigonometric series. The synthesis formula of spatially inconsistent ground motion is as follows:

(7)$a_i(t)={\displaystyle \sum _{m=1}^n{\displaystyle \sum _{k=0}^{N-1}{A}_{im}({\omega }_k)\cos [}}{\omega }_{k}t+{\theta }_{im}({\omega }_k)+{\phi }_{mk}]$

where ${\omega }_k=(k-1)\Delta \omega$, $\Delta \omega =2\pi /T_d$, and d is the duration of synthetic ground motion; ω_k is a circular frequency discrete point; ${\phi }_{mk}$ is a random phase angle uniformly distributed over the interval (0, 2π).

Furthermore, the generated non-uniform ground motion is de-stabilized according to the following formula:

(8)${a_i}^{\prime }(t)=f(t)a_i(t)$

where f(t) is the intensity envelope function, which generates the non-stationary amplitude of ground motion. In this paper, a multi-peak intensity envelope model is adopted [37].

(9)$f(t)=\left\{\begin{array}{ccc}{(t/t_1)}^2& & t<t_1\\ 1& & t_1<t\le t_2\\ e^{-c(t-t_2)}& & t_2<t\end{array}\right.$

where c controls the attenuation speed of the attenuation section; t₁ and t₂ control the start and end moments of the stationary section, respectively.

2.2. Fault Velocity Function

The fault dislocation should not only consider the impact of the dislocation momentum and the dislocation form, but also, the impact of the dislocation rate should be fully considered. The sliding velocity function can represent the change in the sliding velocity of the fracture surface with time and space, and the fault dislocation rate is mainly controlled via the sliding velocity function. Commonly used sliding velocity functions mainly include the Brune function, Haskell function, bell function, exponential function, triangle function, etc. [38,39,40]. Considering the balance between the fidelity to actual situations and computational efficiency, this paper chooses the bell function to simulate the viscoelastic slip of the fault.

(10)$v(t)=\left\{\begin{array}{ccc}0& & t<0\\ u(\infty ){\displaystyle \frac{1-\cos \frac{2\pi t}{T}}{T}}& & 0\le t<T\\ 0& & T\le t\end{array}\right.$

where $u(\infty )$ is the final error momentum, and T is the total error time.

3. Numerical Simulation

3.1. Numerical Model

The MZ hydropower station project is located in Yongde County, Lincang City, Yunnan Province. Several developed faults are crossed within the inlet hydraulic tunnel. According to the geological report of the pre-feasibility study, the basic seismic intensity of the site in the MZ tunnel area is estimated to be Ⅷ degree. The geological profile of the right bank hydraulic tunnel of the MZ hydropower station is shown in Figure 2.

Based on the cross-fault underground hydraulic tunnel of the MZ power station project, the three-dimensional numerical model as shown in Figure 3 is established. The horizontal direction parallel to the tunnel axis is selected as the x-axis, the vertical direction perpendicular to the tunnel axis is the y-axis, and the vertical direction is the z-axis. According to the experience of selecting the boundary range of the underground engineering model, the boundary effect affecting the tunnel response can be ignored when the boundary is taken as 3–5 times the diameter of the tunnel [41,42], so the model size is set to 150 m × 60 m × 60 m. The tunnel adopts a circular section with an inner diameter of 7 m, and the lining is made of C30 concrete with a thickness of 0.5 m. In order to meet the requirements of dynamic calculation accuracy [43], the maximum mesh size of the model does not exceed 5 m, and a total of 286,185 grids with 196,731 nodes are divided, and the lining grids are encrypted with 19,861 grids. The Mohr–Coulomb elastoplastic constitutive model was adopted for upper and lower plates and faults, and the linear elastic model was adopted for the lining. The artificial viscous boundaries and free-field boundaries are employed in dynamic analysis [44,45]. The values of the physical and mechanical parameters used in the three-dimensional model are shown in Table 1, which are determined according to the recommended values from the geological report.

3.2. Spatially Incoherent Seismic Acceleration Time-History Curve

Assuming that the input ground motion acceleration time history at the bottom of the model within 50 m is consistent, with a total of three exciting segments (1^#, 2^#, and 3^#), three acceleration time-history curves are generated according to the spatially inconsistent ground motion synthesis method described in Section 2. The reference for the parameters of the self-power spectral model is cited in ref. [35]. Figure 4a shows the acceleration time-history curve generated after correction via displacement reference. It can be seen from the figure that the waveform of each acceleration time-history curve is basically the same. Besides slight differences in amplitude and phase, ground motion also has corresponding time lag, and no displacement deviation can be seen from the displacement time-history curve, as shown in Figure 4b. The development trend of the simulated power spectrum is also basically consistent with that of the target spectrum, as shown in Figure 4c; that is, the generated non-consistent ground motion can well reflect the travelling wave effect, coherence effect, attenuation effect, and non-stationary characteristics of ground motion and meet the requirements of simulation.

3.3. Simulation Process

In order to simulate the dynamic response of the tunnel accurately, the simulation process is divided into material parameter assignment, initial ground stress, static excavation, lining support, and dynamic calculation. The initial ground stress adopts the gravity stress field, and the lateral pressure coefficient is set as k_x = 1.0, k_y = 1.2, and k_z = 1.0. The stress field after support is the initial condition for dynamic calculation. The boundary around the model is a free-field boundary and the bottom boundary is a viscous boundary. Local damping is employed in the model, with a critical damping ratio of 5% and a damping coefficient of 0.1571 [46].

Four typical monitoring surfaces, S1, S2, S3, and S4, are set up in this paper, among which S1 is located in the 1^# excited segment, S2 and S3 are located in the 2^# excited segment, and S4 is located in the 3^# excited segment. The fault intersects with the S3 monitoring section, with an inclination angle of α and a thickness of h. Eight monitoring points are arranged on each monitoring section, with points A, B, C, and D located on the top arch, left haunch, right haunch, and bottom arch of the lining, respectively. Points A’–D’ are located in the surrounding rock. The monitoring scheme is illustrated in Figure 5.

Based on the survey results, the peak horizontal acceleration of seismic motion in the project area is 0.2 g. Considering the seismic intensity and other conditions, the peak acceleration of the seismic waves for the first 20 m, which is the 1^# excitation segment, is set to 2 m/s². The peak accelerations for the 2^# and 3^# excitation segments are adjusted accordingly in proportion. Ground motion for each stimulating section originates from the model’s base, with horizontal vibration input along the x-axis. The 1^# acceleration time-history curve, post-amplitude modulation, serves as input for consistent ground motion. The model’s bottom plate boundary is fixed, and initial downward velocity is applied to the upper plate’s bottom to induce displacement, creating relative motion between the two plates to simulate normal fault dislocation. The proposed final fault dislocation momentum is 0.1 m, with a total dislocation time T of 2 s. To delve deeper into the impact of fault dislocation and non-uniform vibration on tunnel dynamic response, four distinct loading conditions are simulated. Condition 1 involves fault displacement alone, condition 2 solely incorporates inconsistent seismic motion input, condition 3 introduces inconsistent seismic motion following fault displacement, and condition 4 introduces consistent seismic motion post-fault displacement, as depicted in Figure 6.

4. Result Analysis

4.1. Displacement Response Analysis of Lining

Figure 7 depicts the axial distribution curve of z-displacement for each monitoring point along the lining. As illustrated in Figure 7a, the z-displacement responses of each monitoring point exhibit consistency under normal fault dislocation, displaying a distribution akin to the shape of an “S”. Notably, the displacement curve at each monitoring point is most pronounced at the fault, indicating the region with the most significant mutation. Figure 7b shows that after the end of non-uniform ground motion, the displacements of all monitoring points are basically the same along the axis, and the displacements of all monitoring points far from the fault have a small difference, but the displacements of all monitoring points within the fault range have a significant difference, indicating that the fault is the key area for seismic fortification of the lining. Figure 7c shows that under the combined action of fault and incoherent ground motion, the displacement change of the lining located on the upper wall is significantly larger than that of the lining of the lower wall. This is because when the lower wall is fixed and the upper wall is dislocated downward, the influence of fault dislocation on the upper wall is more significant, while the influence of dislocation is weaker the farther away from the upper wall, and the displacement change is mainly affected by incoherent ground motion. Compared with Figure 7b, the axial displacement of the lining changes more significantly, indicating that the fault plays a leading role in the whole loading process. Compared with Figure 7a, the displacement mutation of the overhanging wall lining is more significant; that is, non-uniform vibration aggravates the axial change in the displacement of the lining, indicating that although the fault dislocation plays a leading role, the existence of ground motion still aggravates the deformation of the lining. Compared with Figure 7d, the displacement difference of each monitoring point is more obvious, and the displacement mutation range is also larger, indicating that the displacement response of the lining is inconsistent due to non-uniform vibration, so it is very necessary to consider non-uniform vibration in the seismic design of underground hydraulic tunnels. Under the four working conditions, the displacement changes of the left and right waist arches (monitoring points B and C) basically coincide along the axis.

Figure 8 shows the peak z-displacement and occurrence time at point A on the four monitoring surfaces. Figure 8 illustrates that the displacement response of each monitoring surface is greater under the combined influence of fault displacement and ground motion compared to the effects of fault dislocation or ground motion alone. Therefore, it is more appropriate to consider the combined action of fault and ground motion in the study of a tunnel under complex geological conditions and complex input loads, and the displacement response of lining under the combined action of fault and non-uniform ground motion on each monitoring surface is the largest. It is more realistic to use incoherent ground motion input when studying the effect of fault dislocation and earthquakes on tunnel dynamic responses. The peak displacement under fault dislocation of S1 and S2 monitoring surfaces is larger than that under non-uniform ground motion, because S1 and S2 sections are located on the upper wall and are greatly affected by fault dislocation, while S3 and S4 monitoring surfaces are far away from the upper wall and the displacement response is mainly affected by non-uniform ground motion. The peak displacement times of each monitoring surface in working condition 2 are 14.5 s, 15.3 s, 15.3 s, and 16 s, respectively; the peak displacement times of each monitoring surface in working condition 3 are 4.9 s, 2.9 s, 2.9 s, and 15.9 s, respectively. The inconsistency in the timing of peak displacements at different monitoring sections is due to the traveling wave effect and coherence effect of non-uniform seismic excitation. As the axial distance increases, the time it takes for the lining structure to reach its peak displacement is delayed, which is also the reason for the varying deformations at different parts of the structure. This highlights the necessity of considering the non-uniformity of seismic excitation in the seismic design of large underground geotechnical projects and is also supported by the findings in ref. [47]. While the peak displacement in working condition 4 is about 4.9 s, 4.9 s, 10.2 s, and 10.2 s, and the peak displacement on the monitoring surface of S1 and S2 is earlier than that of S3 and S4. This is because S1 and S2 are located on the upper disk, and the peak displacement time of the S3 and S4 monitoring surface is basically the same due to the non-consistent vibration input. Under the combined influence of fault displacement and incoherent ground motion, the peak displacement of each monitoring surface occurs earlier compared to when subjected solely to incoherent ground motion. This observation suggests that fault dislocation accelerates the deformation of the lining structure.

4.2. Lining Stress Analysis

Figure 9 shows the axial distribution curve of the maximum principal stress at each monitoring point of the lining. As depicted in Figure 9a, the distribution characteristics of the maximum principal stress at each monitoring point along the lining exhibit similarity along the axis under fault dislocation. A phenomenon of stress concentration is observed within the fault range, impacting an area approximately twice the width of the fault. The maximum principal stress shows tension within the fault range, while the area away from the fault is basically compression. The conclusion drawn from Figure 9b is basically the same as that of Figure 9a, but the stress is affected by the fault in a slightly larger range. The comparison reveals that the axial distribution curve of the maximum principal stress in Figure 9a,d is basically a single peak value, and the peak value occurs at the fault, but the maximum principal stress in Figure 9c,d obviously has multiple peaks; apart from faults, the remaining peaks mainly occur in the upper wall. It shows that the stress concentration of the lining structure will occur not only at the fault, but also within the range of the dislocated surrounding rock under the combined action of fault dislocation and ground motion. The further the stress fluctuation of the concrete lining, the more unfavorable it will be to the lining. Therefore, it is necessary to consider the fault dislocation in the seismic response analysis of the cross-fault tunnel.

Figure 10 shows the peak maximum principal stress and occurrence time at point A of the four monitoring surfaces. From Figure 10, it is evident that the peak maximum principal stress response on each monitoring surface is greater under the combined influence of fault displacement and non-uniform ground motion compared to the effects of fault dislocation or ground motion acting alone. The peak maximum principal stress on each monitoring surface follows the following order: operating condition 3 > operating condition 2 > operating condition 4 > operating condition 1. This suggests that ground motion exerts a greater influence on the lining stress response than fault dislocation. Additionally, the stress concentration induced by inconsistent vibration is more pronounced, with larger amplitudes, compared to that caused by consistent vibration. The occurrence time of the peak maximum principal stress of each monitoring surface in working condition 2 is 12.8 s, 16.8 s, 16.8 s, and 17.5 s, respectively. The occurrence time of peak displacement of each monitoring surface in working condition 3 is 0.01 s, 2.84 s, 2.88 s, and 12.06 s, respectively. The time inconsistency of peak maximum principal stress in different monitoring surfaces reflects the traveling wave effect of non-uniform ground motion. The occurrence time of the peak maximum principal stress under the combined action of fault and incoherent ground motion is earlier than that under the single action of incoherent ground motion, which indicates that the fault dislocation will accelerate the stress concentration of the lining structure.

4.3. Influence of Fault Dislocation on Seismic Response of Tunnel

To investigate the impact of fault dislocation on the seismic response of the tunnel, we employ the normal fault dislocation mode with a duration of 2 s. The fault dislocation momentum is varied at 0.1 m, 0.2 m, 0.3 m, and 0.4 m, respectively. The ground motion input comprises non-uniform ground motion. Figure 11 shows the z-direction peak displacement of the monitoring parts of the lining on each monitoring surface under different fault faults. From the figure, it is evident that the peak displacement of the monitoring sections of the lining increases with the fault dislocation momentum. Additionally, the profile located nearer to the upper wall experiences more pronounced changes in peak displacement, while the profile situated farther away from the upper wall shows a slower rate of change in peak displacement. In this simulation, the z-direction peak displacement of each section is top arch > waist arch > bottom arch. Figure 12 shows the peak maximum principal stress at the monitoring parts of the lining on each monitoring surface under different fault momenta. The figure illustrates that as the fault dislocation momentum increases, the peak maximum principal stress at each monitoring section of the lining also increases.

5. Conclusions

Using stochastic process theory, this paper generates ground motion acceleration time-history curves for each excitation section, capturing spatial non-consistency. Additionally, it introduces a sliding velocity function to reflect the velocity effect of fault dislocation. Subsequently, leveraging insights from actual projects, a numerical model of rock surrounding the tunnel lining and fault interaction is established. Four conditions are compared and analyzed: the dynamic response of the tunnel under fault dislocation, non-consistent vibration, fault dislocation combined with non-consistent vibration, and fault dislocation combined with consistent vibration. This study delves deeply into the influence of fault dislocation and non-consistent vibration on hydraulic tunnels. Finally, the preliminary analysis examines the seismic response patterns of tunnels under varying fault dislocation momenta, leading to the following conclusions:

• (1). The z-direction displacement varies significantly within the fault area, and there is a noticeable stress concentration in the maximum principal stress within this region. Therefore, the fault is a critical area that requires special attention in the seismic design of tunnels. Under the combined effects of fault dislocation and non-uniform seismic excitation, the peak z-direction displacement at point D on monitoring section S1 of the tunnel is −9.01 cm, and the maximum principal stress is 6.7 MPa. The percentage increase in displacement due to non-uniform seismic excitation is 5.5%, and the percentage increase in stress is 44.7%. The dynamic response under the combined effects is greater than the response under either fault dislocation or non-uniform seismic excitation alone;

• (2). Compared to uniform seismic excitation, the peak displacement of the tunnel under non-uniform seismic excitation increases by up to 6.42%, and the peak maximum principal stress increases by up to 28%. Under the combined effects of fault dislocation and non-uniform seismic excitation, the maximum time difference for the occurrence of peak displacements at different monitoring sections is 11 s. Under the combined effects of fault dislocation and uniform seismic excitation, the maximum time difference is 5.3 s. The results from non-uniform seismic excitation more closely align with observed damage cases. In longer tunnels, there is a noticeable delay effect in deformation along the axial direction during an earthquake, and fault dislocation shortens the time it takes for the lining structure to reach its peak seismic response;

• (3). Under the non-uniform seismic excitation, the larger the fault dislocation magnitude, the greater the peak displacement and peak maximum principal stress at each monitoring point of the lining. Compared to a fault dislocation magnitude of 0.3 m, when the dislocation magnitude is 0.4 m, the maximum peak displacement of the tunnel increases by 13.04%, and the maximum peak principal stress increases by 40.66%. Moreover, the peak displacement of sections near the upper wall undergoes more significant changes, while those farther away from the upper wall experience slower changes. In simulations, extreme response values mainly manifest at the top arch and lumbar arch, highlighting these critical areas for seismic design considerations.

The proposed approach and results provide a reference for the seismic safety analysis of long underground hydraulic tunnels. However, this study primarily focuses on the impact of the combined effects of fault dislocation and non-uniform seismic excitation on tunnel responses. Future work would analyze the sensitivity of key parameters such as fault dislocation magnitude, fault dislocation velocity, direction, and seismic excitation to expand the research. Additionally, the influence of water within the tunnel will be a critical aspect to consider in future research efforts.

Author Contributions

Conceptualization, H.L. and W.Y.; methodology, W.Y.; validation, Y.C.; investigation, J.F.; data curation, D.L.; writing—original draft preparation, H.L.; writing—review and editing, W.Y.; funding acquisition, H.L. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Data are available from the authors by request.

Conflicts of Interest

Authors Hao Liu and Dexin Li were employed by the company PowerChina Zhongnan Engineering Corporation Limited. The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Table

Figure 1. Research framework.

Figure 2. Geological profile of the right bank hydraulic tunnel of the MZ Project (Unit: m).

Figure 3. Three-dimensional numerical model of hydraulic tunnel.

Figure 4. Non-uniform ground motion acceleration time-history curve, displacement time-history curve, and power spectrum: (a) acceleration history curve, (b) displacement history curve, and (c) power spectrum.

Figure 5. Monitoring scheme of numerical model.

Figure 6. Four working conditions: (a) normal fault dislocation, (b) inconsistent ground motion, (c) incoherent ground motion after normal fault dislocation, and (d) input consistent vibration after normal fault dislocation.

Figure 7. Axial distribution curve of z-direction displacement of monitoring points: (a) working condition 1, (b) working condition 2, (c) working condition 3, and (d) working condition 4.

Figure 8. Peak displacement and time in the z direction at point A of each monitoring section.

Figure 9. Axial distribution curve of maximum principal stress at monitoring points: (a) working condition 1, (b) working condition 2, (c) working condition 3, and (d) working condition 4.

Figure 10. Peak maximum principal stress and time in the z direction at point A of each monitoring section.

Figure 11. Peak displacement in the z direction of monitoring parts of lining on each monitoring surface under different fault dislocation momentum (unit: cm).

Figure 12. Peak maximum principal stress in the z direction of monitoring parts of lining on each monitoring surface under different fault dislocation momentum (unit: 107 Pa).

Table 1

Physical and mechanical properties of materials used in 3D model.

Component	Density/(kg·m−3)	Elastic Modulus/GPa	Poisson Ratio	Internal Friction Angle/°	Cohesion/MPa	Tensile Strength/MPa
Hanging wall	2600	4.0	0.25	33	1.5	0.5
Footwall	2600	4.0	0.25	33	1.5	0.5
Fault	2100	0.68	0.32	33	0.6	0.12
Lining	2500	30	0.20	-	-	-

References

1. Cilingir, U.; Madabhushi, S. A model study on the effects of input motion on the seismic behaviour of tunnels. Soil Dyn. Earthq. Eng.; 2011; 31, pp. 452-462. [DOI: https://dx.doi.org/10.1016/j.soildyn.2010.10.004]

2. Lu, C.; Hwang, J. Nonlinear collapse simulation of Daikai Subway in the 1995 Kobe earthquake: Necessity of dynamic analysis for a shallow tunnel. Tunn. Undergr. Space Technol.; 2019; 87, pp. 78-90. [DOI: https://dx.doi.org/10.1016/j.tust.2019.02.007]

3. Uzarski, J.; Arnold, C. Chi-Chi, Taiwan, earthquake of September 21, 1999, reconnaissance report. Earthq. Spectra; 2001; 17, 190.

4. Chen, Y.; Booth, D. The Wenchuan Earthquake of 2008: Anatomy of a Disaster; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2011.

5. Shen, Y.; Gao, B.; Yang, X.; Tao, S. Seismic damage mechanism and dynamic deformation characteristic analysis of mountain tunnel after Wenchuan earthquake. Eng. Geol.; 2014; 180, pp. 85-98. [DOI: https://dx.doi.org/10.1016/j.enggeo.2014.07.017]

6. Yu, H.; Chen, J.; Bobet, A.; Yuan, Y. Damage observation and assessment of the Longxi tunnel during the Wenchuan earthquake. Tunn. Undergr. Space Technol.; 2016; 54, pp. 102-116. [DOI: https://dx.doi.org/10.1016/j.tust.2016.02.008]

7. Cui, G.; Wang, X.; Wang, Z.; Wang, D. Model tests on the antibreaking countermeasures for tunnel lining across stick-slip faults. Adv. Mater. Sci. Eng.; 2020; 2020, 7937595. [DOI: https://dx.doi.org/10.1155/2020/7937595]

8. Zhong, Z.; Wang, Z.; Zhao, M.; Du, X. Structural damage assessment of mountain tunnels in fault fracture zone subjected to multiple strike-slip fault movement. Tunn. Undergr. Space Technol.; 2020; 104, 103527. [DOI: https://dx.doi.org/10.1016/j.tust.2020.103527]

9. An, S.; Tao, L.; Han, X.; Zhang, Y. Application of two-level design method on subway tunnel crossing active fault: A case study on Urumqi subway tunnel intersected by reverse fault dislocation. Bull. Eng. Geol. Environ.; 2021; 80, pp. 3871-3884. [DOI: https://dx.doi.org/10.1007/s10064-021-02164-y]

10. Wang, T.; Geng, P.; Li, P.; Wang, P.; Wang, L. Deformation and failure of overburden soil subjected to normal fault dislocation and its impact on tunnel. Eng. Fail. Anal.; 2022; 142, 106747. [DOI: https://dx.doi.org/10.1016/j.engfailanal.2022.106747]

11. Ding, J.; Jin, X.; Guo, Y.; Li, G. Numerical simulation for large-scale seismic response analysis of immersed tunnel. Eng. Struct.; 2006; 28, pp. 1367-1377. [DOI: https://dx.doi.org/10.1016/j.engstruct.2006.01.005]

12. Anastasopoulos, I.; Gerolymos, N.; Drosos, V.; Kourkoulis, R.; Georgarakos, T.; Gazetas, G. Nonlinear response of deep immersed tunnel to strong seismic shaking. J. Geotech. Geoenvironmental Eng.; 2007; 133, pp. 1067-1090. [DOI: https://dx.doi.org/10.1061/(ASCE)1090-0241(2007)133:9(1067)]

13. Jiang, L.; Chen, J.; Li, J. Seismic response of underground utility tunnels: Shaking table testing and FEM analysis. Earthq. Eng. Eng. Vib.; 2010; 9, pp. 555-567. [DOI: https://dx.doi.org/10.1007/s11803-010-0037-x]

14. Wang, G.; Yuan, M.; Miao, Y.; Wu, J.; Wang, Y. Experimental study on seismic response of underground tunnel-soil-surface structure interaction system. Tunn. Undergr. Space Technol.; 2018; 76, pp. 145-159. [DOI: https://dx.doi.org/10.1016/j.tust.2018.03.015]

15. Sun, B.; Zhang, S.; Wang, C.; Cui, W. Ground motion duration effect on responses of hydraulic shallow-buried tunnel under SV-waves excitations. Earthq. Eng. Eng. Vib.; 2020; 19, pp. 887-902. [DOI: https://dx.doi.org/10.1007/s11803-020-0602-x]

16. Jiao, H.; Du, X.; Zhao, M.; Huang, J.; Zhao, X. Nonlinear Seismic Response of Rock Tunnels Crossing Inactive Fault under Obliquely Incident Seismic P Waves. J. Earth Sci.; 2021; 32, pp. 1174-1189. [DOI: https://dx.doi.org/10.1007/s12583-021-1483-2]

17. Zhou, H.; He, C.; Wang, X.; Chen, Y.; Li, J. Assessment of the seismic response of shallow buried elliptical tunnels. J. Earthq. Eng.; 2021; 27, pp. 465-487. [DOI: https://dx.doi.org/10.1080/13632469.2021.2009057]

18. Sun, B.; Deng, M.; Zhang, S.; Wang, C.; Cui, W.; Li, Q.; Yan, H. Optimal selection of scalar and vector-valued intensity measures for improved fragility analysis in cross-fault hydraulic tunnels. Tunn. Undergr. Space Technol.; 2023; 132, 104857. [DOI: https://dx.doi.org/10.1016/j.tust.2022.104857]

19. Sun, B.; Liu, W.; Deng, M.; Zhang, S.; Wang, C.; Guo, J.; Wang, J. Compound intensity measures for improved seismic performance assessment in cross-fault hydraulic tunnels using partial least-squares methodology. Tunn. Undergr. Space Technol.; 2023; 132, 104890. [DOI: https://dx.doi.org/10.1016/j.tust.2022.104890]

20. Feng, S.; Dai, S.; Lei, H. Application of an integrated nonuniform seismic random excitation method in tunnel engineering in soft soil areas. Comput. Geotech.; 2023; 153, 105043. [DOI: https://dx.doi.org/10.1016/j.compgeo.2022.105043]

21. Feng, Q.; Hu, Y. Spatial correlation of earthquake motion and its effect on structural response. Proceeding USPRC Bilateral Workshop on Earthquake Engineering; Science Press: Beijing, China, 1982; pp. 5-14.

22. Huang, J.; Shao, W.; Zhao, M.; Han, J.; Zhao, X.; Du, X.; Lv, X. Simplified analytical solution for circular tunnel under obliquely incident SV wave. Soil Dyn. Earthq. Eng.; 2021; 140, 106429. [DOI: https://dx.doi.org/10.1016/j.soildyn.2020.106429]

23. Zhang, G.; Wang, P.; Zhao, M.; Du, X.; Zhao, X. Seismic structure-water-sediment-rock interaction model and its application to immersed tunnel analysis under obliquely incident earthquake. Tunn. Undergr. Space Technol.; 2021; 109, 103758. [DOI: https://dx.doi.org/10.1016/j.tust.2020.103758]

24. Huang, J.; Zhao, M.; Du, X. Non-linear seismic responses of tunnels within normal fault ground under obliquely incident P waves. Tunn. Undergr. Space Technol.; 2017; 61, pp. 26-39. [DOI: https://dx.doi.org/10.1016/j.tust.2016.09.006]

25. Huang, J.; Du, X.; Jin, L.; Zhao, M. Impact of incident angles of P waves on the dynamic responses of long lined tunnels. Earthq. Eng. Struct. Dyn.; 2016; 45, pp. 2435-2454. [DOI: https://dx.doi.org/10.1002/eqe.2772]

26. Huang, J.; Du, X.; Zhao, M.; Zhao, X. Impact of incident angles of earthquake shear (S) waves on 3-D non-linear seismic responses of long lined tunnels. Eng. Geol.; 2017; 222, pp. 168-185. [DOI: https://dx.doi.org/10.1016/j.enggeo.2017.03.017]

27. Li, X.; Li, Z.; Crewe, A. Nonlinear seismic analysis of a high-pier, long-span, continuous RC frame bridge under spatially variable ground motions. Soil Dyn. Earthq. Eng.; 2018; 114, pp. 298-312. [DOI: https://dx.doi.org/10.1016/j.soildyn.2018.07.032]

28. Zerva, A.; Zervas, V. Spatial variation of seismic ground motions: An overview. Appl. Mech. Rev.; 2002; 55, pp. 271-297. [DOI: https://dx.doi.org/10.1115/1.1458013]

29. Jia, H.; Zhang, D.; Zheng, S.; Xie, W.; Pandey, M. Local site effects on a high-pier railway bridge under tridirectional spatial excitations: Nonstationary stochastic analysis. Soil Dyn. Earthq. Eng.; 2013; 52, pp. 55-69. [DOI: https://dx.doi.org/10.1016/j.soildyn.2013.05.001]

30. Zhang, J.; Zhang, M.; Li, M.; Min, Q.; Shi, B.; Song, L. Nonlinear dynamic response of a CC-RCC combined dam structure under oblique incidence of near-fault ground motions. Appl. Sci.; 2020; 10, 885. [DOI: https://dx.doi.org/10.3390/app10030885]

31. Liu, Z.; Liu, Z.; Ai, Q.; Ruan, X.; Tan, C. Global reliability evaluation of a high-pier long-span continuous RC rigid frame bridge subjected to multi-point and multi-component stochastic ground motions. Soil Dyn. Earthq. Eng.; 2023; 164, 107623. [DOI: https://dx.doi.org/10.1016/j.soildyn.2022.107623]

32. Housner, G.; Jennings, P. Generation of artificial earthquakes. J. Eng. Mech. Div.; 1964; 90, pp. 113-150. [DOI: https://dx.doi.org/10.1061/JMCEA3.0000448]

33. Kanai, K. An Empirical Formula for the Spectrum of Strong Earthquake Motions. Bull. Earthq. Res. Inst.; 1961; 39, pp. 85-95.

34. Hu, Y.; Liu, S.; Dong, W. Earthquake Engineering; CRC Press: Boca Raton, FL, USA, 1996.

35. Li, H.; Zhu, B.; Lu, T.; Yang, J.; Xia, X.; Liu, Y. Seismic response analysis of underground rock mass cavern group considering spatial inconsistency of ground motion. Chin. J. Rock Mech. Eng.; 2008; 27, pp. 1757-1766. (In Chinese)

36. Qu, T.; Wang, J.; Wang, Q. A practical model for spatially varying ground motion power spectrum. Acta Seismol. Sin.; 1996; 1, pp. 55-62. (In Chinese)

37. Amin, M.; Ang, A. Nonstationary stochastic models of earthquake motions. J. Eng. Mech. Div.; 1968; 94, pp. 559-584. [DOI: https://dx.doi.org/10.1061/JMCEA3.0000969]

38. Hisada, Y. A Theoretical Omega Square Model Considering the Spatial Variation in Slip and Rupture Velocity. Bull. Seismol. Soc. Am.; 2000; 90, pp. 387-400. [DOI: https://dx.doi.org/10.1785/0119990083]

39. Anderson, J.; Richards, P. Comparison of Strong Ground Motion from Several Dislocation Models. Geophys. J. Int.; 2007; 42, pp. 347-373. [DOI: https://dx.doi.org/10.1111/j.1365-246X.1975.tb05866.x]

40. Brune, J. Tectonic stress and the spectra of seismic shear waves from earthquakes. J. Geophys. Res.; 1970; 75, pp. 4997-5009. [DOI: https://dx.doi.org/10.1029/JB075i026p04997]

41. Su, K.; Zhang, Y.; Chang, Z.; Wu, H.; Wang, T.; Zhou, W. Transverse extent of numerical model for deep buried tunnel excavation. Tunn. Undergr. Space Technol.; 2019; 84, pp. 373-380. [DOI: https://dx.doi.org/10.1016/j.tust.2018.11.034]

42. He, Y.; Sun, X.; Zhang, M. Investigation on the deformation of segment linings in cross-fault tunnel considering the creep behavior of surrounding rock during construction-operation period. Buildings; 2022; 12, 1648. [DOI: https://dx.doi.org/10.3390/buildings12101648]

43. Lysmer, J.; Kuhlemeyer, R. Finite dynamic model for infinite media. J. Eng. Mech. Div.; 1969; 95, pp. 859-877. [DOI: https://dx.doi.org/10.1061/JMCEA3.0001144]

44. Yang, B.Y.; Chen, J.T.; Xiao, M. Seismic response and damage mechanism of lining structures for underground tunnels across fault. Chin. J. Geotech. Eng.; 2020; 42, pp. 2078-2087.

45. Wen, Y.M.; Xin, C.L.; Shen, Y.S.; Huang, Z.M.; Gao, B. The seismic response mechanisms of segmental lining structures applied in fault-crossing mountain tunnel: The numerical investigation and experimental validation. Soil Dyn. Earthq. Eng.; 2021; 151, 107001. [DOI: https://dx.doi.org/10.1016/j.soildyn.2021.107001]

46. Zheng, Y.; Ye, H.; Huang, R. Analysis and discussion on failure mechanism and fracture surface of seismic slope. Chin. J. Rock Mech. Eng.; 2009; 28, pp. 1714-1723. (In Chinese)

47. Yu, H.; Yan, X.; Bobet, A.; Yuan, Y. Multi-point shaking table test of a long tunnel subjected to non-uniform seismic loadings. Bull. Earthq. Eng.; 2018; 16, pp. 1041-1059. [DOI: https://dx.doi.org/10.1007/s10518-017-0223-6]