Dynamic Stability Analysis of Subsea Tunnel Crossing Active Fault Zone: A Case Study

Abstract

The rock strength in an active fault zone is low and the surrounding rock is fractured and has poor stability, making any subsea tunnel crossing the active fault zone extremely susceptible to disasters such as tunnel collapse, sudden water ingress, and mud inrush. This poses a potential threat to the construction project, making the dynamic stability analysis of a subsea tunnel crossing an active fault zone of great significance. This study takes the second subsea tunnel crossing the Cangkou Fault in Jiaozhou Bay as the engineering background and conducts numerical simulations by employing different lining stiffnesses for tunnel excavation, as well as applying dynamic loads. The dynamic stability of the subsea tunnel crossing the active fault zone is evaluated by comparing and analyzing the lining’s displacement, peak acceleration, and stress characteristics. This study explores the disaster-causing mechanisms of active fractures, determining that the hazard of orthogonal misalignment in an active fault zone is the least severe, while the hazard of opposite misalignment is the most severe. This research provides a basis for disaster prevention and mitigation in active fracture zones.

1. Introduction

Subsea tunnels are often located at the edges of continental frames and islands where seismic activity is intense, crossing various geological strata and sometimes intersecting with active fault zones. When the construction of tunnels inevitably crosses active fault zones, the activity patterns and deformation characteristics of these zones during the effective period of the project will affect the safety and stability of the engineering work [1,2]. Potential earthquakes and misalignments of fracture zones can lead to the opening and closing of rock mass structural surfaces, slippage, and the initiation–expansion–coalescence of cracks at these active fault zones crossed by subsea tunnels. This reduces the stability of the rock mass and, in severe cases, can even cause structure damage and deformation to the tunnel, destroy the tunnel’s waterproofing facilities, and eventually trigger a chain of secondary geological disasters, such as large deformations, collapses, sudden water ingress, and mud inrush, posing a serious threat to the long-term safety and maintenance of subsea tunnels [3,4].

Due to the principle of avoiding active faults in tunnel engineering design, there are few engineering cases of tunnels crossing active fault zones. The characteristics of internal force distribution, deformation patterns, and failure modes within the rock mass of fault zones under fault displacement and seismic actions are closely related to the dynamic response of submarine tunnels to fault zones [5,6,7]. The seismic damage to fault-crossing tunnels is considered to have two main causes: one is the vibrational damage to the tunnel caused by earthquakes; the other is the shear damage to the tunnel caused by the fault. When the engineering geological conditions where the tunnel is located are good, the tunnel itself has better seismic performance. Tunnel seismic damage generally occurs more frequently in areas with complex geological conditions, such as poor rock mass quality and significant variations in stratum conditions. He et al. [8] found that tunnels crossing fault zones are prone to damage in earthquakes, and the fault zone sections are the controlling areas for the tunnel’s resistance to seismic damage. The main forms of seismic damage are displacement and various types of cracking. Currently, the main theoretical analysis method for tunnels’ seismic resistance is the analytical method, with the pseudo-static method being the most widely used within the analytical methods. However, due to constraints that differ from above-ground structures, the results of this method tend to be conservative. Research on the dynamic response and stability of tunnels can also be conducted through model experiments and numerical simulations. Fang et al. [9], based on a shaking table, studied the dynamic response patterns and characteristics of fault-crossing tunnels under different vibration directions. They found that faults could weaken the propagation of seismic waves in the rock mass to a certain extent. Shahidi et al. [10] studied the longitudinal seismic dynamic response characteristics of the Koohrang-III water conveyance tunnel in Iran through numerical simulation.

The impact of active fault zones on submarine tunnels mainly has three aspects: one is earthquakes induced by sudden stick–slip deformation, the second is the effect of long-term creep deformation on the structure, and the third is the impact of distant earthquakes on active fault zones [6,7,11,12]. Among these, the creep deformation of the fault zone under tectonic stress is long-term, and it is necessary to design corresponding tunnel structures to prevent its harm to the tunnel structure in advance. Sudden stick–slip deformation, characterized by large displacement, a high displacement gradient, and the potential to induce high-magnitude earthquakes, is mainly controlled by limiting the range of damage in engineering projects. In response to the impact of active fault displacement on tunnel structures, Kiani et al. [13] proposed a new method for assessing tunnel vulnerability curves based on centrifuge experiments and studied the deformation and force patterns of fault-crossing tunnels under the action of normal faults. Ma et al. [5] revealed the damage mechanisms of tunnels under the displacement effect of strike–slip and thrust faults through numerical simulation, from aspects such as lining deformation, stress distribution, and the distribution of plastic zones, finding that lining damage is significantly influenced by the displacement of active faults. Liu et al. [14], based on model experiments, studied the effect of fault dip angle on the mechanical response characteristics of tunnel structures, discovering that as the fault dip angle decreases, the tunnel structure’s mode of failure shifts from bending to shear damage. Due to the complex force characteristics of submarine tunnels in complex seabed environments under seismic actions, their damage mechanisms are diversified under the effect of earthquakes. Influenced by the interaction between marine engineering geology and tunnel structures, as well as the stability of marine engineering geological bodies, the issue of seismic disasters becomes even more complicated [15,16,17].

In this study, a series of numerical simulations are conducted to analyze the influence of the dynamic seismic action and the motion of active fault zones on the dynamic stability of a subsea tunnel crossing an active fault zone. Section 2 briefly introduces the engineering background of the second subsea tunnel in Jiaozhou Bay, Qingdao. Section 3 shows the detailed numerical simulation and study of the dynamic stability of the subsea tunnel under different support conditions. Section 4 shows the results of the numerical simulation of the dynamic stability of the subsea tunnel under seismic effects and fault motion. Section 5 is the discussion and limitations section of this paper. The conclusions are given in Section 6.

2. Engineering Background

The second subsea tunnel in Jiaozhou Bay will connect Qingdao and Huangdao and be a transportation corridor and critical infrastructure asset linking various blue economic zones. With a total length of 14.36 km, it will be Asia’s longest subsea road tunnel and the world’s most giant subsea concealed dredge road tunnel [18].

The coast of Jiaozhou Bay has typical mountainous bedrock, and the nearby topography and geomorphology are influenced by rock lithology and geological structure, and external forces influence the geomorphology. On the whole, Jiaozhou Bay has a large belly and small mouth, deep water in the southeast, sloping topography from north to south, and a flat seabed in the north of the bay. According to the preliminary survey results, the stratum of the Jiaozhou Bay area is mainly granite, with local development of gabbro veins, and the surface layer is covered with a small amount of residual soil and a weathering layer. The geotechnical type is relatively singular and uniform, except for the possible development of a local tectonic fracture zone. There is no other particular geotechnical distribution. The area where the project is located crosses two first-class geotectonic units of North China, and the Yangzi fault block, and the project line is basically located in the area of the North China fault block where the seismic structure is developed. Among them, the north–northeast-oriented fractures are large in scale and have a long history of activity, while the northwest-oriented fractures show a Neozoic nature and control the Late Cenozoic basin. The main source of earthquakes in the project area is the strong seismic activity of the Tanlu Fault Zone and Bohai–Penglai Fault Zone. Combined with the regional geological background, the near-field area of the project has the conditions for earthquakes, and earthquakes of magnitude 6 to 7 may occur within one hundred years.

The second subsea tunnel in Jiaozhou Bay crosses the Cangkou Fault (F52 in Figure 1), which has more fractured rocks and develops several cross-sections. During the construction process, many disaster problems may emerge, such as issues causing local relative displacement of the tunnel and changing the force and surrounding rock stress of the tunnel structure, which will cause incalculable damage to the subsea tunnel.

3. Numerical Simulation

3.1. Model Geometry

This study uses the second subsea tunnel in Jiaozhou Bay, crossing the Cangkou Fault, as an example. The first consideration is to build a model of a suitable size. The mesh division needs to strike a balance between computational cost and accuracy. Excessive density incurs higher computational times, while insufficient density compromises accuracy. Balancing comprehensive site conditions and model simplification, this paper determines the subsea tunnel model, for the left and right directions, should extend 40 m on either side. The upward direction, or buried depth, is set at 45 m, downward as 25 m, and along the tunnel axis as 120 m. The tunnel radius is 6.5 m. After dividing the grid, there are 109,739 cells and 116,684 nodes. The rock mass in the numerical simulation is modeled using the Mohr–Coulomb model, and the excavated section is set as a null model.

The model (shown in Figure 2) is divided into eight groups. The active fault zone is 10 m thick, with an inclination of 62.9°, and is located between 20–65 m in the Y-direction. The overburdened miscellaneous fill is 15 m thick. The blue part represents the hanging wall of the active fault zone, and the brown part represents the footwall of the active fault zone. The blue and brown areas are slightly weathered granite, divided into two by the active fault zone. The tunnel inside the middle of the model represents the tunnel. The initial support, the secondary lining, and the grouting reinforcement are set around the tunnel one by one.

To control the accuracy and make the dynamic calculation more clear and definite, the cell size must meet specific requirements. For seismic waves to be accurately described in the model, the dimensions of the dividing cell must satisfy:

$$\Delta l\le \frac{\lambda }{10}\sim \frac{\lambda }{8}$$

(1)

where λ is the wavelength corresponding to the highest frequency.

The allowable cell size is 5.3 m, while the maximum cell size in this model is 3.2 m, which meets the accuracy requirement.

Along the tunnel, a total of nine typical cross-sections are axially selected at 10 m, 20 m, 30 m, 40 m (within the fault), 50 m, 60 m, 70 m, 90 m, and 110 m for the placement of monitoring points. For the layout of these monitoring points, four are placed in each tunnel section at the top of the arch, arch waist, arch foot, and bottom of the invert, respectively. The monitoring points’ layout and specific locations are shown in Figure 3.

3.2. Material Parameter

For the parameters of geotechnical units, based on the existing engineering cases and the geological data of Jiaozhou Bay, Qingdao, the following necessary stratigraphic mechanical parameters are determined, as shown in

Table 1. Mechanical parameters of the stratum.

Material	Volumetric Mass/kg/m3	Bulk Modulus/MPa	Shear Modulus/MPa	Cohesion/kPa	Internal Friction Angle/°	Strength of Tension/kPa
Seafloor cover	1900	17.9	1.86	5	18.0	0
Active fault zone	2200	18	4530	300	45	100
Slightly weathered granite	2720	7550	5660	6050	58	4257

.

Three contact surfaces are set before, after, and at the top of the active fault zone, and two more contact surfaces are set to separate the overlying soil and rock layers. The parameters of the contact face are that the tangential stiffness is 4.66 × 10¹⁰ N/m, the cohesion is 400 kPa, and the internal friction angle is 24°.

3.3. Boundary Conditions

The displacement boundary conditions of the model in this paper are constrained to be 0 in the X and Y positive and negative semi-axis directions and at the bottom, with the bottom being a static boundary, while the top is free in the Z-direction without constraint. The displacement boundary conditions are achieved by specifying the velocity component of each boundary as 0. The dynamic boundary condition is set after the ground stress equilibrium, and each stress value returns the initial stress state to 0.

The form of seismic wave used in this paper is a synthetic wave, acting at the bottom, applied along the tunnel’s axial Y direction and perpendicular to the axial X direction, but discounted by 50% for the vertical tunnel axial X direction, according to the ratio of 0.5:1 with the seismic wave along the axial tunnel input. Considering the calculation time cost and accuracy control, seismic waves up to 5 s are selected in this paper. The dynamic analysis is solved by concentrating the real mass of the nodes in the grid density, and the critical calculation time steps during the calculation are as follows:

$$\Delta t_{critical}=\min \left\{\frac{V}{C_p{A}_{\max }^f}\right\}$$

(2)

where V is the volume of the tetrahedral subunit; Cp is the p-wave velocity; and A^f_max is the maximum surface area associated with the tetrahedral subcell.

The primary seismic intensity of the proposed project area is 7 degrees, the design primary seismic acceleration value is 0.10 g, the third group is selected in the seismic design grouping, and the site category is I₁ according to each geotechnical layer and overburden in the site area, and its characteristic period value is 0.35 s. This paper uses a form of local damping, and its value is taken as 0.2.

Underground rock formations are subject to fault motion after geological tectonic movements or extrusion by large ground stresses, producing significant displacements along the active fault zone. Due to the complexity of the stratigraphic dislocation, the dislocation is imposed after the tunnel excavation and support, and in this paper, three cases are considered for the direction factor.

(1)

Applying velocities for Y of 3 × 10⁻⁵ m/s and Z of 2 × 10⁻⁶ m/s to the Y = 0 surface.

(2)

Applying a velocity for Y towards 3 × 10⁻⁵ m/s to the Y = 0 surface.

(3)

A set of velocities for Y toward 3 × 10⁻⁵ m/s and 3 × 10⁻⁵ m/s are respectively applied to the two faces of the model, Y = 0 and Y = 120 m, in a squeezed manner.

3.4. Numerical Simulation Scheme

In this study, three major groups of eight working conditions for the simulation schemes are designed for the seismic dynamic and fault motion effects of the subsea tunnel crossing the active fault zone, and the tunnel lining stiffness is shown in

Table 2. Simulation schemes considering different supporting structures, dynamic effects, and fault motions.

No.	Excavation Method	Factor I	Factor II	Factor III
1	Up and down step method	Lining II	/	/
2		Lining I	Seismic dynamic	/
3		Lining II
4		Lining III
5		Lining IV
6		Lining II	Seismic dynamic	Orthogonal stagger (Y and Z directions)
7				Unidirectional stagger (Y-direction)
8				Extrusion stagger (Y-direction)

.

4. Results

Through the numerical simulation of multiple control variables, the stability analysis of the face of the subsea shield tunnel under different geological and construction conditions is carried out, and the influence of different factors on the face instability is studied.

4.1. Basic Features of Undersea Tunnel Excavation by up and down Step Method

The up-and-down step method excavates the upper and lower sections successively and advances simultaneously, which is conducive to the stability of the tunnel face. It has advantages over the full-section method in terms of simplicity of operation, low disturbance, and construction process safety. The displacement clouds after the upper and lower step method of excavation are shown in Figure 4, the acceleration clouds in Figure 5, and the stress clouds in Figure 6.

Following excavation by the up-and-down step method, it can be obtained that the displacements generated by the tunnel are mainly settlements. The highest acceleration peaks are found at the top and waist of the tunnel. Based on the trend of the data from the measured points arranged in the tunnel section, it can be determined that the fault caused the greatest impact at Y = 40 m.

4.2. Dynamic Stability of Tunnel Structures with Different Support Parameters under Seismic Effects

Further, the dynamic stability is analyzed under different liner stiffnesses in comparison, and four groups of liner stiffnesses are set in increasing order, with the bulk moduli of 3.89 GPa, 7.78 GPa, 15.56 GPa, and 38.9 GPa; and shear moduli of 2.915 GPa, 5.83 GPa, 11.66 GPa, and 29.15 GPa, respectively.

In the total displacement cloud of the section at Y = 40 m (Figure 7), the total deformation trend is the same under different liner stiffnesses, and the total displacement decreases with the increase of liner stiffness. The results of each displacement component at Y = 40 m profile with different liner stiffnesses are shown in

Table 3. Maximum displacement of surrounding rock under different lining stiffness in the section at Y = 40 m.

Liner Type	X-Displacement/ cm	Y-Displacement/ cm	Z-Displacement/ cm	Displacement Magnitude/cm
Liner I	1.2035	2.0786	−40.701	40.767
Liner II	1.1865	2.0785	−40.624	40.690
Liner III	1.1608	2.0785	−40.481	40.548
Liner IV	1.1620	2.0785	−40.242	40.309

. The X-direction and Z-direction displacements decrease with increasing liner stiffness. Although the difference is not significant, the trend is consistent with the displacement magnitude. The peak stresses of different liner stiffnesses are demonstrated in Figure 8.

In the case of specific initial support, the X- and Z-directional (reflected as settlement) displacements become smaller as the liner stiffness increases. The absolute values of maximum and minimum principal stress and maximum shear stress all increase with the increase of liner stiffness under particular initial support, indicating that its stress pattern changes and takes on more stress.

4.3. Dynamic Stability of Tunnel Structures under the Effects of Fault Motion

The forms of staggering in different ways are compared. From Figure 9, we can see that the peak displacement of the model lining structure is extrusion stagger > unidirectional stagger > orthogonal stagger, and the results show that the most unfavorable case for the lining structure should be the third case. The stress indicators (Figure 10) show that the maximum and minimum principal stresses and maximum shear stresses under orthogonal stagger and unidirectional stagger are tiny. In contrast, the stresses under phase stagger are much larger than the first two, reflecting the most unfavorable stresses on the lining structure. In the orthogonal and unidirectional staggering, compressive stresses appear in the parts of the maximum principal stress cloud where no velocity is applied, and tensile stresses appear at the ends of the tunnel in the minimum principal stress cloud. In the phase staggering, compressive stresses appear at the ends of the tunnel in the maximum principal stress cloud, and tensile stresses appear in the middle in the minimum principal stress cloud. This indicates that the fault motion greatly impacts the lining structure, and it also means that the lining in the subsea tunnel interval near the active fault zone has higher requirements.

Considering the factors of fault motion and seismic activity, since the direct cause of fault motion is the displacement and motion of the two plates, the direct-acting surface is the most critical part. The results also confirm that the maximum and minimum principal stresses and maximum shear stresses are on the direct-acting surface, and both sides are maximally under the effect of fault motion and bear most of the force, while the lining at the farther side bears less force. It should be that the active fault zone is the defense’s focus.

5. Discussion and Limitations

When the active fault zone starts to misshape, the rock layers on both sides are subjected to shear or tensile compression, and relative displacement occurs between the upper and lower plates [9,14]. Furthermore, it contacts the seafloor overburden, transmitting shear stress to the far side. Lateral bending may occur, leading to damage on the seafloor surface. The activity of the active fault zone causes displacement and mechanical effects only in a specific range on both sides of the active fault zone, while it has less influence on the more distant rock masses. For different ways of movement and displacement of the fracture zone, its tension and compression states do not change. Mainly, the changes and differences are of specific values also affect the range and peak of its force-deformation area. When the pressure on both sides of the active fault zone is very high, the phenomenon of tension damage will also occur, further causing the separation and gap between the upper and lower plates of the active fault zone. The width and angle of the active fault zone do not change the tunnel structure’s deformation, displacement, tension, and shear. In practice, this is often accompanied by earthquakes, cracking, depressions, or bulges in the ground that release energy from the deeper rocks.

Only when fault motion occurs are the tunnel and the surrounding rock subjected to tremendous stress values, and the peak acceleration and relative displacement at this time are also enormous. The cross-sea tunnel through the active fault zone is mainly affected by shear damage, which causes relative displacement and misalignment. This can cause displacement, deformation, and obvious dislocation of the tunnel structure, or even collapse and destruction of the whole, which poses a considerable risk of harm during construction and regular use. Therefore, special attention should be paid to exploring the specific location of the detailed active fault zone and its detailed parameters in the project’s construction, and to strengthening the advance forecasts, real-time monitoring, and risk assessment to better cope with various risk hazards [19,20,21].

Liu et al. [2] showed the danger of the subsea tunnel face when crossing the Cangkou fault, and that the deformation of the tunnel face can be controlled by increasing the tunnel burial depth (but not exceeding the limit), reducing the tunnel diameter, and increasing the surrounding rock strength.

The active fault zone set in this study is tilted in a single direction after some simplification, while in practice, it is mainly tilted in two directions. In this study, the active fault zone is considered to suffer from the effects of seismic dynamics and fault motion under the separate working conditions. However, the more complex dual actions of simultaneous effects are not considered.

Based on this study, adding methods such as fluid–solid coupling and flow coupling to simulate the pore water pressure and rock force will bring the model closer to the actual situation. Designs including a flexible lining and damping layer can also be added to offset the large deformation caused by fault motion or dynamic seismic effects.

6. Conclusions

In this study, we use the second subsea tunnel in Jiaozhou Bay as the background and establish a finite difference model to analyze the force and deformation of the tunnel lining structure under the action of seismic dynamics and fault motion and obtain the following conclusions.

(1)

In the case of specific initial support, both X- and Z-directional (reflected as settlement) displacements become smaller with the increase of lining stiffness, but the Y-directional displacement remains constant. The absolute values of the maximum and minimum principal stresses and maximum shear stress all increase with the increase of liner stiffness.

(2)

The active fault zone undergoes motion, and it can be known that larger displacements occur when specific external influences are applied. It is concluded that the most unfavorable situation will be caused when phase extrusion is applied. The multi-directional disturbance and external forces applied will make the rock formation more prone to motion, but the displacement values are not significant due to the presence of the lining structure.