Stability of a Connecting Tunnel in a Shaft–Tunnel System under High Hydraulic Gradient and Staged Excavation: Implication from Numerical Modelling

1  Introduction

Hydropower development has become a key approach for meeting public and industrial energy demands [1–3]. It is also consistent with low-carbon and environmentally friendly energy strategies. To achieve water diversion and conveyance, hydraulic tunnels are widely used in hydropower projects [4–6]. Shaft and tunnel combined structures are commonly adopted in these systems. Vertical shafts regulate hydraulic head differences. Tunnels enable long-distance water conveyance. This structural combination offers high hydraulic efficiency and flexible engineering layouts. As a result, it has been extensively applied in underground engineering [7,8]. Despite these advantages, shaft and tunnel combined structures are relatively vulnerable components of underground systems. The combined effects of excavation unloading and construction disturbance are prone to inducing stress concentration and deformation incompatibility in the surrounding rock mass, resulting in a high sensitivity to instability during the construction stage [9–11]. Consequently, understanding the mechanical response and stability evolution of shaft and tunnel connection zones is a critical engineering issue in hydraulic tunnel construction.

Current research on the mechanical response and stability evolution of intersecting tunnel structures mainly involves field investigations [12–14], physical model tests [15,16], and numerical simulations [17–19]. Field monitoring provides important evidence for revealing deformation and damage characteristics of intersecting tunnel structures under real engineering conditions. For example, Ding et al., [20] established a microseismic (MS) monitoring system in the intersecting-tunnel zone of the tailrace tunnel at the SJK Hydropower Station. Based on the spatiotemporal activity of MS events induced by excavation unloading, they showed that structural weaknesses and multi-cavern effects intensify rock-mass damage under excavation–unloading disturbance, thereby enabling an assessment of surrounding-rock stability performance during construction. However, due to constraints in monitoring layout, construction conditions, and the irreversibility of the engineering process, field monitoring results are generally limited to evaluating the stability performance of a predefined design and providing references for construction adjustments and feedback-based optimization of design parameters, rather than offering quantitative support for excavation scheme development at the design stage.

Physical model tests play an important role in revealing failure mechanisms and deformation characteristics of complex tunnel structures and enable direct observation of crack propagation and instability processes under controlled conditions [21–24]. Hydraulic tunnel engineering is a typical form of underground construction, where the design stage is often constrained by incomplete geological information and uncertainties in surrounding-rock parameters [25]. Due to the high cost and experimental complexity, physical modelling approaches are limited in conducting systematic parametric comparison analyses. With advances in testing equipment and computational performance, numerical simulation methods are capable of realistically reproducing the mechanical behaviour of surrounding rock and construction-induced disturbances while maintaining engineering scale and complex boundary conditions [26–28]. Wang et al. [29] developed a refined three-dimensional numerical model to systematically analyzed the construction response and interaction mechanism of double-track tunnels under different layout forms. However, in the construction of the connection between shafts and tunnels, step-by-step excavation methods such as the step method are often adopted, resulting in obvious irregularity and three-dimensional face effect in the excavation section [30]. Xiao et al. [31] proposed a multiple Levels of Detail (multi-LoD) Building Information Modelling (BIM)-integrated framework that links parametric tunnel modelling with three-dimensional numerical analysis to assess structural performance and sustainability under internal water pressure. Sun et al. [25] evaluated the stability of a hydraulic tunnel subjected to partial lining removal and staged variable cross-section excavation through three-dimensional numerical modelling. The study quantified stress redistribution, plastic zone development, and safety factors during sequential construction and validated the results with in situ monitoring data. Zhang et al. [32] established a thermal–hydraulic–mechanical coupled numerical model to simulate the failure process of water-resistant rock mass in water-rich fault tunnels under high-temperature and high-pressure conditions, and validated the model through self-developed tunnel water-inrush model tests. Their results showed that coupled hydraulic, thermal and mechanical effects can significantly influence the instability process, failure kinematics and stress path of water-resistant rock mass.

In underground excavation under high water levels or water-rich strata, groundwater seepage is one of the key factors affecting surrounding-rock stability [33–35]. The presence of groundwater not only reduces the effective stress and shear strength of the rock mass, but may also impose significant hydraulic pressure on the surrounding rock under high-head conditions [36–38]. During the reconstruction of existing hydraulic tunnels or the construction of new connecting tunnels in hydropower expansion projects, a special construction condition may occur during the flood season [39]. When the reservoir water level rises and overtops the construction cofferdam, the intake foundation pit and the excavated shaft may become hydraulically connected to the reservoir. Under this condition, the unexcavated rock mass between the shaft and the advancing connecting tunnel acts as a temporary blocking section. On the shaft side, the blocking section is subjected to high external water pressure, whereas on the tunnel side, staged bench excavation and drainage create a relatively depressurized boundary. The resulting asymmetric hydraulic boundary, combined with excavation unloading and three-dimensional face effects, may promote plastic-zone expansion, local strength degradation and progressive instability of the remaining rock mass. If the retained length of the blocking section is insufficient, plastic zones may become connected and water inrush or local collapse may occur. Although existing studies have improved the understanding of tunnel stability under excavation disturbance and seepage conditions, the stability of the unexcavated blocking section in a shaft–connecting tunnel system under reservoir–shaft hydraulic connection has received limited attention. In particular, the hydraulic boundary in this situation is different from that of ordinary water-rich tunnel faces or conventional pressurized hydraulic tunnels. The blocking section is controlled by the combined effects of asymmetric pore-pressure evolution, staged bench excavation, excavation–support interaction and three-dimensional face effects. In addition, the rational retained length of the blocking section is still often determined by engineering experience or a single safety index, while the relationship among plastic-zone connectivity, local strength reserve and global strength-reduction convergence remains insufficiently clarified.

In this study, the stability of the blocking section in a shaft–tunnel connection under high-head hydraulic conditions is investigated using a three-dimensional hydro-mechanical coupled numerical model established in FLAC3D. The model explicitly incorporates staged bench excavation, excavation–support interaction, pore water pressure evolution, and asymmetric hydraulic boundary conditions induced by reservoir–shaft hydraulic connection. An integrated stability assessment method is proposed by jointly using plastic-zone connectivity, point safety factor, and the global strength reduction factor of safety to determine a rational retained length that satisfies engineering safety requirements. Section 2 defines the engineering background and outlines the shaft-tunnel hydraulic-connection condition. Section 3 describes the 3D modelling strategy and the staged excavation-support implementation. Section 4 analyses the surrounding-rock response characteristics induced by the pre-excavation stage of the connecting tunnel. Section 5 presents the integrated criteria to identify the stability threshold and recommends the rational unexcavated blocking section length for flood-season construction control. Section 6 summarizes the main findings and highlights the implications most relevant to engineering practice.

2  Engineering Background

This project is part of a large-scale hydropower capacity expansion scheme, in which a new intake is constructed based on the existing T5 tunnel and hydraulically connected to the original tunnel system through a shaft and a connecting tunnel (Fig. 1a). For detailed engineering information, please refer to [25]. The cofferdam constructed in the intake area is designed for dry-season conditions (Fig. 1b), with a crest elevation of EL. 458 m.

Figure 1: (a) Engineering layout; (b) construction cofferdam; (c) geological profile; (d) longitudinal profile of the connecting tunnel.

During the flood season, the reservoir operates at an elevation of approximately EL. 472.5 m; therefore, from July to November each year, the excavated intake foundation pit and shaft structure become directly connected to the reservoir. As of mid-June 2025, the shaft excavation had reached EL. 381.5 m. The secondary concrete lining above EL. 393.63 m has been completed, with the shaft crown elevation at EL. 414 m and a final lining diameter of 13.716 m. Due to the approaching flood season, excavation activities at the shaft have been temporarily suspended (Fig. 1d). The connecting tunnel is excavated from the existing T5 tunnel after removal of the original lining. The excavation and support works of the horizontal curved section have been completed, and the current heading has advanced to chainage 0+21 m. The tunnel is constructed using a three-bench excavation method. The upper bench has been excavated to chainage 21, the middle bench lags behind by 6–10 m, and the lower bench will be excavated after breakthrough.

According to the exposed geological conditions (Fig. 1c), the surrounding rock of the connecting tunnel mainly consists of argillaceous limestone, with local dolomitic limestone. The rock mass is relatively poor in integrity, with an average Geological Strength Index (GSI) value of approximately 30. When the reservoir water level remains below EL. 442 m, the measured seepage rate at the tunnel face is approximately 10–15 L/min, indicating a certain degree of permeability in the surrounding rock. Under high-water-level operation, once the shaft becomes hydraulically connected to the reservoir, the unexcavated rock mass ahead of the connecting tunnel will be subjected to a high-head seepage environment, and its stability characteristics are therefore of significant engineering concern.

During the construction stage of this project, the shaft may become hydraulically connected to the reservoir during the flood season, and the rock mass ahead of the connecting tunnel will be subjected to significant high-head seepage. The connecting tunnel is excavated using a staged bench method, and the surrounding rock undergoes continuous unloading and stress redistribution. Under the combined effects of hydro-mechanical coupling and staged excavation disturbance, the stability of the unexcavated rock mass ahead of the tunnel face becomes a critical factor governing construction safety. This unexcavated rock mass serves as a transitional barrier between the reservoir and the working face during construction. Its mechanical state is controlled not only by rock mass strength parameters, but also by hydraulic head conditions and excavation advance distance. Insufficient stability reserve may lead to plastic zone expansion or even penetration, thereby reducing overall seepage resistance and bearing capacity, and increasing the risk of water inrush and local instability. Therefore, under shaft-tunnel hydraulic connection conditions, it is necessary to perform plastic zone penetration analysis, point safety factor evaluation, and strength reduction calculations for the remaining unexcavated plug section of the connecting tunnel.

3  Numerical Modelling Framework

3.1 Numerical Model Establishment

Fig. 2 shows the overall three-dimensional Computer-Aided Design (CAD) model of the headworks of the hydropower expansion project. A unified spatial coordinate system was established for numerical analysis. The X-axis is perpendicular to the T5 tunnel axis, the Y-axis is aligned with the T5 tunnel axis, and the Z-axis is oriented vertically upward. The coordinate origin is located at the intersection between the T5 tunnel and the new connecting tunnel, and the Z-axis elevation is assigned according to the actual engineering datum. The computational domain extends from −70 to 190 m in the X-direction, from −110 to 110 m in the Y-direction, and from 320 to 490 m in the Z-direction. The model includes the existing T5 tunnel, the horizontal curved section, the new connecting tunnel, the shaft, and the vertical curved section. The geometric configuration incorporates the tunnel excavation face, staged bench excavation sequence, and the interfaces between initial shotcrete support and final concrete lining.

Figure 2: The overall 3D model of the project.

The overall three-dimensional numerical model consists of more than 1.5 million elements, predominantly high-precision hexahedral elements, ensuring reliable engineering simulation and sufficient computational accuracy. All major structural components, including the T5 tunnel, horizontal curved section, connecting tunnel, shaft, and vertical curved section, were discretized into independent zones and blocks (Fig. 3).

Figure 3: (a) FLAC3D mesh model; (b) groupings of shaft-tunnel system; (c) excavation groupings of bench method; (d) support structure and element schematic diagram.

The connecting tunnel was strictly subdivided according to the staged bench excavation method into upper, middle, and lower benches, comprising a total of 20 blocks, including eight backfilled concrete lining blocks. This discretization enables step-by-step simulation consistent with the actual construction sequence. In the numerical model, the concrete lining was simulated using solid elements. Rock bolts and cable anchors were modelled using cable elements. Shotcrete layers and steel arches were represented using liner elements with equivalent stiffness verification. The computational parameters of this project mainly refer to the original design report. Meanwhile, based on the on-site feedback information, it is known that the surrounding rock of the tunnel is mainly argillaceous limestone, and partly dolomitic limestone, with an average GSI of around 30. Therefore, the parameters adopted in the numerical calculation are shown in Table 1.

3.2 Initial Stress and Hydraulic Conditions

The maximum burial depth of the key research sections, including the T5 tunnel, horizontal curved section, shaft, and connecting tunnel, is approximately 100 m. Deterministic structural planes are not well developed in the surrounding rock. Therefore, the initial stress field was mainly determined by the self-weight stress of the rock mass, together with the corresponding hydraulic load.

Under the self-weight stress condition, the vertical stress was estimated from the unit weight of the rock mass and the overburden depth. Since this expression is well known, it is not listed separately as an equation. The horizontal in-situ stress was estimated using the lateral pressure coefficient:

σh=K0σv(1)

where σh is the horizontal in-situ stress, σv is the vertical in-situ stress calculated from the overburden weight, and K0 is the lateral pressure coefficient. For an isotropic elastic rock mass under one-dimensional strain conditions, K0 can be approximated as:

K0=v1−v(2)

where v is Poisson’s ratio of the rock mass. In this study, v = 0.25, and thus K0 = 0.33. Therefore, the horizontal stress was initialized as approximately one third of the vertical overburden stress. The self-weight stress field was further adjusted according to the actual topography and then superimposed with the hydraulic load. As shown in Fig. 4a, the burial depth of the T5 tunnel and the connecting tunnel is generally 40–100 m, and the initial maximum principal stress is approximately 1.5–3.5 MPa, increasing with burial depth.

Figure 4: (a) Initial stress distribution; (b) initial pore water pressure.

In this project study, the initial water level elevation of the reservoir was set at 442 m above sea level. Fig. 4b illustrates the distribution characteristics of pore water pressure in the rock mass. It can be seen that the initial pore pressure of the surrounding rock of the tunnel in the connection section is approximately 0.8 MPa. In addition, in the numerical model, a surface force equivalent to the hydraulic gradient was simultaneously applied to the surface of the flooded area to ensure the balance of internal and external forces.

3.3 Excavation and Support Sequence

Based on the zoned excavation simulation of the T5 tunnel, horizontal curved section, shaft, and vertical curved section, a refined excavation process simulation was further implemented for the connecting tunnel, which is the focus of this study. As shown in Fig. 5, the connecting tunnel was constructed using a staged bench method, divided into upper, middle, and lower benches. The upper bench advanced first, the middle bench lagged by 6–10 m, and the lower bench was excavated after tunnel breakthrough. Support measures, including initial shotcrete, rock bolts, and concrete lining, were incorporated according to the actual construction sequence.

Figure 5: The overall 3D model of the project.

To realistically simulate excavation-induced disturbance, the zone relax technique in FLAC3D was adopted. This approach progressively reduces the stress, stiffness, and density of elements within the excavation zone, enabling automatic stepwise relaxation during block-by-block excavation. The method gradually diminishes reaction forces along the excavation boundary, thereby reproducing the face support effect and the gradual stress release process during excavation cycles. It also mitigates the amplification of surrounding-rock disturbance caused by quasi-inertial effects inherent in explicit numerical solutions. When the relaxation coefficient decreases to zero, the elements are assigned a NULL constitutive model and subsequently removed from the computational domain, completing the excavation simulation.

The timing of support installation has a significant influence on surrounding-rock deformation, plastic zone development, and internal forces within the support system [22]. Therefore, accurate simulation of the support sequence is essential for the reliability of numerical results and the assessment of surrounding-rock stability. Based on the design report, previous numerical modelling experience, and back-analysis of field monitoring data, the support installation schedule was determined in this study. Initial shotcrete, rock bolts, and steel arches were installed when approximately 30% of stress release had occurred. Concrete lining was applied when the stress release reached approximately 70%. The stress release ratio was controlled by adjusting the zone relax reduction factor in FLAC3D, allowing progressive unloading of excavation elements to predefined levels before support installation.

3.4 Numerical Implementation and Convergence Criteria

The numerical analysis was performed using FLAC3D 9.0. The model contains more than 1.5 million zones, predominantly high-precision hexahedral zones. Local mesh refinement was applied around the shaft–tunnel connection, connecting tunnel, T5 bifurcation zone, and retained blocking section, where stress redistribution, pore-pressure gradient, deformation concentration, and plastic-zone development were expected. The connecting tunnel was subdivided into upper, middle, and lower benches, including 20 excavation/support blocks and eight backfilled concrete lining blocks, to reproduce the staged construction sequence. Concrete lining was simulated using solid zones, rock bolts and cable anchors were represented by cable structural elements, and shotcrete together with steel arches was represented by liner structural elements. The model was first solved to initial mechanical and hydraulic equilibrium under gravity stress and prescribed hydraulic boundaries. Staged excavation was simulated using the zone relax technique, in which the stress, stiffness, and density of excavation zones were progressively reduced until the zones were assigned a null model and removed. Initial support was activated at approximately 30% stress release, and concrete lining was activated at approximately 70% stress release. For each excavation stage, convergence was judged by the decrease of unbalanced force ratio to a stable low level and by the stabilization of displacement increments at representative monitoring points. In the strength reduction analysis, global instability was identified from the combined evidence of loss of convergence, continuously increasing displacement, and connected plastic-zone development.

4  Hydro-Mechanical Response during the Initial Excavation Stage

4.1 Deformation Response in the Shaft-Tunnel System

The displacement contours in Fig. 6a indicate that the maximum excavation-induced deformation of the T5 tunnel reaches approximately 100–120 mm. As the original T5 tunnel has been in service for several decades, its historical deformation is not relevant to the present analysis. The displacement field associated with the original excavation was therefore reset to zero in FLAC3D to eliminate pre-existing deformation effects. After partial removal of the concrete lining at the T5C bifurcation portal and excavation of the first two 10 m segments of the horizontal curved section, most of the surrounding rock deformation is concentrated within approximately 30–50 mm. However, localized deformation near the newly excavated tunnel boundary and face is significantly larger, with peak values reaching approximately 70–95 mm, as indicated by the red zones in Fig. 6b. Therefore, the 30–50 mm range represents the dominant deformation level rather than the local maximum value.

Figure 6: Distribution characteristics of surrounding rock deformation. (a) Excavation of T5 tunnel; (b) excavation of horizontal section tunnel; (c) excavation of the upper and middle sections of the vertical shaft and the initial section of the connecting tunnel.

Due to the relatively shallow burial depth of the shaft and the comparatively low in-situ stress level, excavation-induced deformation in this region remains generally below 50 mm. In addition, the concrete lining above EL. 393.63 m has already been completed, further restraining deformation development. However, a locally overhanging rock mass is observed near EL. 385 m on the side adjacent to the connecting tunnel. In this area, localized displacement reaches 120–200 mm, indicating structural asymmetry and stress concentration effects.

In contrast, the initial section of the connecting tunnel is located at a burial depth of approximately 100 m and is subjected to the highest pore water pressure. Fig. 6c reveals significantly larger excavation-induced deformation, typically ranging from 50 to 150 mm. This result reflects the combined influence of deeper overburden stress and elevated hydraulic loading.

4.2 Stress Distribution Characteristics in the Shaft-Tunnel System

After excavation of the T5 tunnel, no significant stress concentration is observed in the surrounding rock, and the local maximum principal stress does not exceed 5 MPa (Fig. 7a). Fig. 7d indicates that stress relaxation occurs around the tunnel periphery. Localized tensile stress develops in the shallow zone, with a maximum tensile stress of approximately 0.25 MPa. These results suggest that, under the current burial depth and stress level, the T5 tunnel excavation does not induce severe stress amplification but leads to moderate tensile stress development near the excavation boundary.

Figure 7: (a) The maximum principal stress distribution of the surrounding rock in the T5 tunnel; (b) the maximum principal stress distribution of the surrounding rock in the horizontal curved section of the tunnel; (c) the distribution of the maximum principal stress in the upper part of the shaft and the initial section of the connecting tunnel’s surrounding rock; (d) the minimum principal stress distribution of the surrounding rock in the T5 tunnel; (e) the minimum principal stress distribution of the surrounding rock in the horizontal curved section of the tunnel; (f) the distribution of the minimum principal stress in the upper part of the shaft and the initial section of the connecting tunnel’s surrounding rock.

After removing half of the concrete lining at the T5C bifurcation portal and excavating the first two 10 m sections of the horizontal curved segment, the stress field changes significantly, as illustrated in Fig. 7b,e. The maximum principal stress distribution (Fig. 7b) shows localized stress redistribution around the newly excavated boundary. The minimum principal stress contours (Fig. 7e) reveal enlarged tensile stress zones near the shallow surrounding rock, where the maximum tensile stress increases to approximately 0.3 MPa. Following installation of shotcrete and lining support, the tensile stress region is substantially reduced, indicating effective control of stress relaxation.

The stress response of the shaft and the initial section of the connecting tunnel is presented in Fig. 7c,f. Due to the shallow burial depth of the shaft, no evident stress concentration is observed around the shaft wall. However, tensile stress develops locally near the overhanging rock mass in the middle portion of the shaft. In contrast, the initial section of the connecting tunnel exhibits more pronounced stress concentration, with local maximum principal stress reaching 4–5 MPa. This is attributed to its greater burial depth and higher pore water pressure. Fig. 7f further indicates tensile stress development at the upper and middle bench faces of the connecting tunnel, with peak tensile stress ranging from 0.2 to 0.3 MPa. After support installation, both stress concentration and tensile stress zones are effectively suppressed.

4.3 Plastic Development and Pore Pressure Evolution

Fig. 8a indicates that, excavation of the T5 tunnel induces a plastic zone with a depth of approximately 10~16 m, which exceeds the effective support range of the systematic rock bolts; therefore, the surrounding-rock load must be jointly carried by shotcrete and the concrete lining. To eliminate the influence of the initial equilibrium state and the plasticity inherited from the historical excavation of the original T5 tunnel, the yield state was reset in FLAC3D. On this basis, after partial removal of the concrete at the T5C bifurcation portal and excavation of the first two 10 m sections of the bifurcation segment, the plastic zone depth reaches about 8–15 m (Fig. 8b), which also exceeds the bolt-supported depth; thus, the structural load transfer relies on the participation of the concrete lining. The shaft is characterised by a shallow burial depth and a relatively low initial in-situ stress level; consequently, plasticity around the shaft is generally limited, with a plastic zone depth mostly less than 5 m and a local maximum of about 8 m near the invert.

Figure 8: (a) The plastic zone of the surrounding rock in the T5 tunnel; (b) the plastic zone of the surrounding rock in the horizontal curved section of the tunnel; (c) the plastic zone in the upper part of the shaft and the initial section of the connecting tunnel’s surrounding rock; (d) the pore water pressure distribution of the surrounding rock in the T5 tunnel; (e) the pore water pressure distribution of the surrounding rock in the horizontal curved section of the tunnel; (f) the distribution of pore water pressure in the upper part of the shaft and the initial section of the connecting tunnel’s surrounding rock.

As illustrated in Fig. 8c, the initial section of the connecting tunnel is located at a greater burial depth and experiences a higher pore pressure level; after excavation, the plastic zone depth can reach 7~14 m, clearly exceeding the support range of the systematic bolts and requiring the concrete lining to share the load.

Fig. 8d presents the pore pressure distribution after excavation of the T5 tunnel and completion of the lining. During excavation, pore pressure around the tunnel decreases markedly due to drainage; after lining construction, a partial rebound in local pore pressure is observed because the lining behaves as a relatively impermeable boundary. Fig. 8e shows the pore pressure distribution after partial removal of the concrete at the T5C portal and excavation of the first two 10 m sections of the horizontal curved segment. A pronounced decrease in pore pressure is observed in the surrounding rock near the newly excavated boundary. In engineering practice, effective drainage measures play an important role in reducing pore pressure and enhancing surrounding-rock stability. Fig. 8f depicts the pore pressure distribution after excavation of the upper–middle part of the shaft and the initial section of the connecting tunnel. Pore pressure around both the shaft and the tunnel decreases significantly. Given the relatively high initial hydraulic head in this area, pore pressure evolution is critical to surrounding-rock stability. Therefore, integrated grouting and drainage measures should be adopted in practice to ensure long-term stability of the surrounding rock.

5  Stability Evolution during the Subsequent Advancing Excavation

5.1 Evolution of Pore Water Pressure during Connecting Tunnel Excavation

According to the long-term reservoir water level records, the reservoir level rises to approximately EL. 472.5 m during July to November each year. The cofferdam in the intake area is a dry-season cofferdam with a crest elevation of EL. 458 m. During the flood season, the reservoir water overtops the cofferdam and establishes direct hydraulic connection between the intake foundation pit, the excavated shaft, and the reservoir.

Fig. 9 illustrates the pore pressure distributions in the shaft and connecting tunnel under different hydraulic conditions. As shown in Fig. 9a, when the reservoir level rises and becomes hydraulically connected to the shaft, pore pressure around the shaft increases significantly. The elevated hydraulic head propagates downward along the shaft, leading to higher pore pressure at the shaft bottom and adjacent rock mass. Fig. 9b presents the pore pressure field before further excavation of the connecting tunnel. Although pore pressure around the shaft remains elevated, the unexcavated rock mass ahead of the tunnel still provides hydraulic isolation. The pore pressure distribution around the connecting tunnel remains relatively uniform at this stage.

Figure 9: Characteristics of pore water pressure variation. (a) Reservoir water overflows the cofferdam; (b) reservoir water level rises; (c) subsequent connecting tunnel excavation.

Fig. 9c shows the pore pressure evolution after the connecting tunnel advances forward. As excavation progresses, drainage pathways are formed, resulting in a significant reduction in pore pressure around the tunnel boundary. However, under shaft–reservoir hydraulic connection, the pore pressure on the shaft side remains high. This condition generates a pronounced hydraulic gradient across the remaining unexcavated rock mass. The increased hydraulic gradient adversely affects the seepage resistance and stability of the rock mass. Therefore, quantitative evaluation of hydraulic gradient distribution and stability characteristics under different reserved lengths is necessary in subsequent analyses. The simulated pore-pressure increase of approximately 0.1–0.3 MPa is consistent with the hydraulic-head variation during the flood season. The initial reservoir level was set at EL. 442 m, whereas the flood-season reservoir level rises to approximately EL. 472.5 m. The corresponding head difference is about 30.5 m, which gives a theoretical hydrostatic pressure increment of approximately 0.30 MPa according to Δp = ρgΔh. Therefore, the upper bound of the simulated pore-pressure increment is physically reasonable. The lower values, approximately 0.1 MPa, occur because the pore-pressure response is affected by the drainage boundary around the excavated tunnel, the permeability of the surrounding rock, and the spatial distance from the shaft-side high-head boundary.

5.2 Surrounding Rock Response during Connection Tunnel Excavation

Fig. 10 presents the displacement evolution characteristics of the surrounding rock in the connecting tunnel and shaft region during different excavation stages. As the connecting tunnel advances from the initial section, the overall deformation level of the surrounding rock shows a decreasing trend because the burial depth, in-situ stress level, and pore-pressure level gradually decrease. The dominant deformation of most surrounding rock is generally within 40–80 mm. It should be emphasized that this range represents the main deformation level rather than the absolute local maximum. In Fig. 10a, limited red zones are observed near the excavation boundary and tunnel face, indicating localized displacement concentration where the peak displacement may exceed 80 mm. This local amplification is mainly caused by excavation unloading and three-dimensional face effects. When excavation proceeds further, the displacement distribution remains generally stable, and the deformation is mainly concentrated around the tunnel boundary. However, when the excavation face approaches the shaft, the pore-pressure distribution shown in Fig. 9c indicates that a pronounced hydraulic gradient develops across the remaining blocking section. Under this condition, local deformation increases significantly, and non-convergent deformation may occur in some areas, as shown in Fig. 10d.

Figure 10: Deformation distribution of surrounding rock in connected tunnels. (a) Chainage 0+40 m; (b) Chainage 0+58 m; (c) Chainage 0+76 m; (d) Chainage 0+94 m.

As shown in Figs. 11 and 12, during the subsequent advance of the connecting tunnel excavation, the burial depth gradually decreases, accompanied by a reduction in the in-situ stress level and pore water pressure. Consequently, stress concentration in the surrounding rock becomes less pronounced, and the overall stress field tends to be more uniform. Meanwhile, evident stress relaxation develops around the connecting tunnel, with localized tensile stresses reaching approximately 0.2–0.3 MPa. After the installation of shotcrete and concrete lining support, the stress relaxation is effectively restrained, and the stress distribution becomes more stable.

Figure 11: The distribution of the maximum principal stress in the surrounding rock of the connecting tunnel. (a) Chainage 0+40 m; (b) Chainage 0+58 m; (c) Chainage 0+76 m; (d) Chainage 0+94.

Figure 12: The minimum principal stress distribution of the surrounding rock of the connecting tunnel. (a) Chainage 0+40 m; (b) Chainage 0+58 m; (c) Chainage 0+76 m; (d) Chainage 0+94.

When the excavation face approaches the shaft location, significant stress relaxation occurs within the reserved blocking section. Under the combined effects of external water pressure and the absence of strong lining support, local collapse or instability may develop in the blocking rock mass.

As illustrated in Fig. 13, during the continued excavation of the connecting tunnel, the depth of the plastic zone slightly decreases with the reduction of burial depth, in-situ stress, and pore water pressure, generally ranging from 6 to 12 m. However, when the excavation advances close to the shaft, the remaining unexcavated blocking section enters an overall yielding state, and the plastic zone tends to become continuous. Under the action of external water pressure and without sufficient lining support, this blocking section is prone to instability and potential failure.

Figure 13: The distribution of plastic zones in the surrounding rock of connecting tunnels. (a) Chainage 0+40 m; (b) Chainage 0+58 m; (c) Chainage 0+76 m; (d) Chainage 0+94 m.

5.3 Stability-Based Length Analysis of the Reserved Rock Plug

The excavation stages were selected at two different resolutions according to the purpose of the analysis. For the overall evolution analysis shown in Figs. 11–13, representative chainages of Chainage 0+40, 0+58, 0+76, and 0+94 m were adopted at an interval of approximately 18 m. This interval corresponds to the characteristic length of the retained blocking section and provides a global view of stress redistribution, deformation evolution, and plastic-zone development during progressive excavation.

In contrast, Figs. 14–16 focus on the stability threshold of the remaining unexcavated rock-plug. Because the critical transition was expected to occur near Chainage 0+76–0+82 m, a refined interval of 6 m was adopted for Chainage 0+70, 0+76, 0+82, and 0+88 m. This finer interval is consistent with the staged bench excavation control scale and enables more accurate identification of the transition from a discontinuous plastic zone to a nearly connected plastic zone and from convergent to non-convergent behavior under strength reduction.

Figure 14: The distribution of plastic zones in the unexcavated blocking section. (a) Chainage 0+70 m; (b) Chainage 0+76 m; (c) Chainage 0+82 m; (d) Chainage 0+88 m.

Figure 15: The distribution of safety factors in the unexcavated blocking section. (a) Chainage 0+70 m; (b) Chainage 0+76 m; (c) Chainage 0+82 m; (d) Chainage 0+88 m.

Figure 16: The distribution of overall safety factor of the unexcavated blocking section. (a) Chainage 0+70 m; (b) Chainage 0+76 m; (c) Chainage 0+82 m; (d) Chainage 0+88 m.

It should be noted that the stability assessment in this section focuses on the upper-bench advancing stage because this stage controls the minimum retained length of the unexcavated blocking section before breakthrough. According to the actual construction sequence, the upper bench advances first, the middle bench lags behind by approximately 6–10 m, and the lower bench is excavated after breakthrough. Therefore, before breakthrough under flood-season high-water-level conditions, the lower bench is not excavated and does not directly control the retained rock-plug length.

Considering that the connecting tunnel is excavated using the benching method, the actual excavation profile is irregular and cannot be accurately characterized by an equivalent circular diameter. In addition, significant three-dimensional face effects exist. Therefore, it is difficult to calculate the safety factor of the blocking section using design codes or simplified analytical solutions. In this study, a three-dimensional numerical simulation approach was adopted. The analysis comprehensively accounts for stress redistribution induced by staged excavation, pore water pressure evolution, and excavation–support interaction effects. The stability of the remaining blocking section between the shaft and the connecting tunnel was evaluated based on the plastic zone penetration pattern, point safety factor reserve, and non-convergence behavior observed in strength reduction analysis. On this basis, a rational recommendation for the blocking section length was proposed.

Fig. 14a–d illustrates the evolution of the plastic zone in the surrounding rock of the shaft and connecting tunnel as excavation advances from Chainage 0+70 to 0+88 m. It can be observed that when the upper bench reaches Chainage 0+82 m, under the combined effects of external water pressure acting on the shaft and the absence of strong lining support in the newly excavated section, the plastic zones within the unexcavated blocking section between the shaft and the connecting tunnel nearly become continuous. This condition indicates a pronounced reduction in stability and a high risk of instability.

Fig. 15a–d presents the distribution of point safety factors and stability zoning characteristics during the same excavation stages. The results show that after lining support is installed, the point safety factors of the shaft and tunnel surrounding rock increase significantly, and overall stability is improved. However, when the upper bench advances to chainage 0+82 m, under the combined influence of shaft external water pressure and insufficient lining support in the newly excavated section, zones with point safety factors lower than 1.2 within the blocking section between the shaft and the connecting tunnel become nearly continuous. This indicates insufficient strength reserve and a markedly reduced safety margin.

Fig. 16a–d presents the distributions of element unbalanced forces and the displacement convergence characteristics of monitoring points in the surrounding rock of the shaft and connecting tunnel during excavation from Chainage 0+70 to 0+88 m, under the strength reduction condition of K = 3.0.

It can be observed that when the upper bench of the connecting tunnel advances to Chainage 0+82 m, under the combined adverse conditions of external water pressure acting on the shaft and the absence of strong lining support in the newly excavated section, the unexcavated blocking section between the shaft and the connecting tunnel exhibits global non-convergence in both stress and deformation. This is manifested by continuously increasing element unbalanced forces and displacement monitoring curves that fail to stabilize, ultimately resulting in overall instability of the blocking section.

Therefore, to ensure that the global safety factor of the blocking section is not less than 3.0, the maximum advance of the upper bench of the connecting tunnel shall not exceed Chainage 0+76 m. Considering comprehensively the plastic zone penetration characteristics, the reserve of point safety factors, and the results of strength reduction analysis for the remaining unexcavated blocking section, it is recommended that during the flood season the advance of the upper bench of the connecting tunnel be limited to within Chainage 0+76 m.

6  Discussions

Under high-water-level operating conditions, the stability of the remaining blocking rock mass in a shaft–connecting tunnel system is not directly governed by the excavation advance distance or changes in overburden depth. Instead, it is primarily controlled by the pore water pressure–effective-stress evolution path driven by the hydraulic gradient. After the reservoir level rises and becomes hydraulically connected to the shaft, the pore water pressure in the surrounding rock on the shaft side increases markedly. Meanwhile, staged excavation and drainage on the connecting-tunnel side create a relatively depressurized zone in the vicinity of the tunnel face, establishing a pronounced hydraulic head difference within the blocking section and thereby amplifying the hydraulic gradient. As the excavation face progressively approaches the shaft, the superposition of the external high-head recharge and the internal excavation-induced unloading promotes gradual plastic-zone expansion within the rock wall, with an increasing tendency toward connectivity (Fig. 13). Notably, during this stage the maximum principal stress does not increase synchronously (Fig. 11); it may even decrease locally as the overburden reduces, while the plastic zone continues to develop. This indicates that instability of the blocking section is better interpreted as a progressive hydro-mechanically coupled failure process rather than an instantaneous failure dominated by conventional stress concentration.

The progressive hydro-mechanically coupled failure mechanism identified in this study is also consistent with experimental findings on the water-induced degradation of limestone. Zhu et al. [40] investigated the mechanical properties and nonlinear energy evolution of holed limestone under hydro-mechanical coupling through triaxial compression tests. Their results indicated that hydro-mechanical coupling significantly affects the deformation, strength degradation and energy dissipation characteristics of limestone. Fan et al. [41] further studied the microstructure damage and mechanical-property evolution of limestone subjected to high-pressure water, and demonstrated that high-pressure water can induce microstructural deterioration and mechanical degradation of limestone. These material-scale observations provide a useful explanation for the engineering-scale response observed in the present shaft–tunnel system. Under reservoir–shaft hydraulic connection, the retained blocking section is subjected to an asymmetric hydraulic boundary: high pore pressure is maintained on the shaft side, while staged excavation and drainage reduce pore pressure on the tunnel side. This hydraulic gradient reduces the effective stress and promotes progressive plastic-zone expansion in the argillaceous limestone surrounding rock. Therefore, the nearly connected plastic zone and strength-reserve degradation observed when the retained length decreases to approximately 12 m can be interpreted as the engineering-scale manifestation of water-induced weakening and hydro-mechanical damage accumulation in limestone.

The three criteria: plastic-zone connectivity tendency, point safety factor, and the convergence behavior of the strength reduction method, are highly consistent at the critical stage. This multi-criteria framework captures the progressive degradation and critical transition behavior more reliably than a single factor-of-safety metric and improves the robustness and interpretability of threshold identification for underwater rock-plug problems. The numerical results are also sensitive to pore water pressure conditions; therefore, the recommended consolidation grouting range of approximately 10 m around the tunnel perimeter was determined based on the simulated plastic-zone depth during staged excavation. As shown in Fig. 13, the plastic zone around the connecting tunnel generally develops to a depth of approximately 6–12 m before the blocking section becomes nearly connected. Therefore, a grouting reinforcement depth of about 10 m was adopted as an engineering control value, covering the dominant plastic-damage zone while remaining consistent with practical reinforcement feasibility. In sections where local plastic-zone depth approaches or exceeds 12 m, the grouting depth should be locally increased according to geological exposure and monitoring feedback.

7  Conclusion

This study investigated the hydro-mechanical response and stability evolution of the unexcavated blocking section under high reservoir water level conditions. A three-dimensional numerical model was established to simulate staged bench excavation, pore pressure redistribution, and excavation–support interaction, and to evaluate the stability of the blocking section between the shaft and the connecting tunnel. It should be noted that the analysis was based on available engineering geological data and simplified representative hydraulic boundary conditions; therefore, the effects of spatial variability in rock mass properties and transient reservoir-level fluctuations were not fully considered. The following conclusions were drawn.

(1)   Hydraulic connection increases the local pore pressure in the surrounding rock of the connecting tunnel by approximately 0.1–0.3 MPa, which is consistent with the hydraulic-head rise from EL. 442 to EL. 472.5 m and is further affected by excavation-induced drainage. Failure of the blocking section follows a progressive hydro-mechanically coupled mode governed by asymmetric hydraulic boundaries and pronounced three-dimensional face effects.

(2)   The blocking section exhibits a typical progressive deterioration behavior. When the upper bench advances to Chainage 0+82 m (12 m remaining), the plastic zone and the low-strength-reserve band tend to become continuous, and the analysis becomes globally non-convergent at K = 3.0, indicating a pronounced increase in hazard level and a substantially elevated risk of progressive instability.

(3)   An integrated criterion is proposed based on plastic-zone connectivity, local strength reserve, and SRM convergence. Using this criterion, a remaining rock-plug length ≥18 m (Chainage 0+76 m) is recommended as the construction control limit. When the remaining length decreases to approximately 12 m (Chainage 0+82 m), the system enters a markedly higher-risk regime; further advance should be avoided and mitigation measures should be implemented promptly.

Acknowledgement: Not applicable.

Funding Statement: This research was funded by the National Natural Science Foundation of China of Funder, grant number U23A2060, 52474150 and 42277461.

Author Contributions: Conceptualization, Feng Gao, Yuepeng Sun and Nuwen Xu; methodology, Guotao Meng and Yuepeng Sun; software, Guotao Meng; validation, Feng Gao, Yuepeng Sun and Nuwen Xu; formal analysis, Guotao Meng; investigation, Heyi Yang and Xianglin Huang; resources, Feng Gao; data curation, Heyi Yang and Xianglin Huang; writing—original draft preparation, Yuepeng Sun; writing—review and editing, Nuwen Xu; visualization, Guotao Meng; supervision, Nuwen Xu; project administration, Feng Gao and Heyi Yang; funding acquisition, Nuwen Xu. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: The data that support the findings of this study are available from the corresponding author upon reasonable request.

Ethics Approval: Not applicable.

Conflicts of Interest: The authors declare no conflicts of interest.

References

1. Dhaubanjar S, Lutz AF, Pradhananga S, Smolenaars W, Khanal S, Biemans H, et al. From theoretical to sustainable potential for run-of-river hydropower development in the upper Indus basin. Appl Energy. 2024;357(3):122372. doi:10.1016/j.apenergy.2023.122372. [Google Scholar] [CrossRef]

2. Guan YH, Post S, Zhao D, Zhang SL, Becker S. Overview of the application status and development trends of hydropower and geothermal power in New Zealand. Energy Built Environ. 2025;6(3):564–84. doi:10.1016/j.enbenv.2024.06.003. [Google Scholar] [CrossRef]

3. Jia JS. A technical review of hydro-project development in China. Engineering. 2016;2(3):302–12. doi:10.1016/J.ENG.2016.03.008. [Google Scholar] [CrossRef]

4. Sun YP, Su HJ, Xiao PW, Li P, Li B, Zhou X, et al. Visualization and early warning analysis of damage degree of surrounding rock mass in underground powerhouse. Int J Min Sci Technol. 2023;33(6):717–31. doi:10.1016/j.ijmst.2022.12.011. [Google Scholar] [CrossRef]

5. Zhang YX, Zhang JX, Wang CS, Ren XH. An integrated framework for improving the efficiency and safety of hydraulic tunnel construction. Tunn Undergr Space Technol. 2023;131(2):104836. doi:10.1016/j.tust.2022.104836. [Google Scholar] [CrossRef]

6. Li B, Xu NW, Xiao PW, Xia Y, Zhou X, Gu GK, et al. Microseismic monitoring and forecasting of dynamic disasters in underground hydropower projects in southwest China: a review. J Rock Mech Geotech Eng. 2023;15(8):2158–77. doi:10.1016/j.jrmge.2022.10.017. [Google Scholar] [CrossRef]

7. Duan KM, Zhang GF, Sun H. Construction practice of water conveyance tunnel among complex geotechnical conditions: a case study. Sci Rep. 2023;13:15037. doi:10.1038/s41598-023-42192. [Google Scholar] [CrossRef]

8. Zhang JH, Sun MQ, Fan ZX, Yuan Y, Sun Q. Shaking table test on shaft-tunnel junction in liquefiable sand site. Tunn Undergr Space Technol. 2024;154(6):106097. doi:10.1016/j.tust.2024.106097. [Google Scholar] [CrossRef]

9. Wang LB, Niu XK, Zhao YH, Li W, Song W, Zhang CP. Study on the mechanism and rapid treatment method of leakage disease at the junction between the shaft and shield tunnel. Eng Fail Anal. 2023;154(4):107639. doi:10.1016/j.engfailanal.2023.107639. [Google Scholar] [CrossRef]

10. Shrestha BP, Panthi KK, Basnet CB. Evaluation of squeezing mechanisms in a headrace tunnel passing through phyllitic slate rock mass: a case study from the Himalaya. Results Eng. 2026;29(7):108900. doi:10.1016/j.rineng.2025.108900. [Google Scholar] [CrossRef]

11. Zhang XM, Deng ZY, Liu BG, Ma X. An integrated analytical approach for predicting structural performance and cracking behavior in composite linings of deep hydraulic tunnels. Tunn Undergr Space Technol. 2026;171:107429. doi:10.1016/j.tust.2025.107429. [Google Scholar] [CrossRef]

12. Yang R, Luo HJ, Sun WT, Xin J, Lu HP, Yang T. Stability study of deep-buried tunnels crossing fractured zones based on the mechanical behavior of surrounding rock. Appl Sci. 2026;16(7):3473. doi:10.3390/app16073473. [Google Scholar] [CrossRef]

13. Tian XX, Song ZP, Zhang YW. Monitoring and reinforcement of landslide induced by tunnel excavation: a case study from xiamaixi tunnel. Tunn Undergr Space Technol. 2021;110(7):103796. doi:10.1016/j.tust.2020.103796. [Google Scholar] [CrossRef]

14. Li X, Xu NW, Mao HY, Xiao PW, Zhou X, Li B, et al. Deformation characteristics and damage evolution analysis of weak interlayer zone in fractured underground cavern. Tunn Undergr Space Technol. 2024;147(1):105686. doi:10.1016/j.tust.2024.105686. [Google Scholar] [CrossRef]

15. Liu SY, Mao YD, Du ZF, Gao L, Cao SD, He MC, et al. Investigation of deformation mechanisms in intersecting sections of interchange tunnels using the negative Poisson’s ratio anchor excavation compensation method. Struct Civ Eng. 2025;19(5):742–59. doi:10.1007/s11709-025-1179-2. [Google Scholar] [CrossRef]

16. Cui XP, Li PF, Ma ZH, Xu Q, Jia F, Ge ZG, et al. Model test study on instability mechanism during shield under-crossing existing pipeline. Eng Fail Anal. 2025;167(1):109041. doi:10.1016/j.engfailanal.2024.109041. [Google Scholar] [CrossRef]

17. Lei HY, Liu YN, Hu Y, Jia R, Zhang Y. Active stability of the shield tunneling face crossing an adjacent existing tunnel: transparent clay model test and DEM simulation. Can Geotech J. 2023;60(6):864–84. doi:10.1139/cgj-2022-0294. [Google Scholar] [CrossRef]

18. Lu Y, He BG, Li Q, Fu H. Amplification and attenuation mechanisms of blast-induced vibration adjacent to cross passage in deep tunnels. Transp Geotech. 2026;56(55–65):101747. doi:10.1016/j.trgeo.2025.101747. [Google Scholar] [CrossRef]

19. Sun YP, Xu NW, Xiao PW, Sun ZQ, Li HL, Liu J, et al. Characterizing large deformation of soft rock tunnel using microseismic monitoring and numerical simulation. J Rock Mech Geotech Eng. 2025;17(1):309–22. doi:10.1016/j.jrmge.2024.05.029. [Google Scholar] [CrossRef]

20. Ding QF, Li B, Su HJ, Xu NW, Li XH, Deng XY. Damage mechanism and stability analysis of rock mass in the high geo-stress tunnel subjected to excavation. Geomat Nat Hazards Risk. 2022;13(1):75–93. doi:10.1080/19475705.2021.2015458. [Google Scholar] [CrossRef]

21. Huang F, Wang ZW, Zhang M, Wang YT, Luo YL, Pan QJ. Model test of the collapse mode of surrounding rock induced by tunnel excavation in a karst area based on PIV. J Cent South Univ. 2026;33(3):1332–47. doi:10.1007/s11771-026-6222-5. [Google Scholar] [CrossRef]

22. Du JM, He C, Xu GW, Wang B, Chen X, Ma GY. Mechanism of the inclined and advanced rockbolt support system and dynamic evaluation of its reinforcement range in mechanized tunneling for high-speed railways. Transp Geotech. 2026;56(08):101739. doi:10.1016/j.trgeo.2025.101739. [Google Scholar] [CrossRef]

23. Wu J, Yao AJ, Wang C, Qin SW. Hydro-mechanical performance and stability of tunnel faces excavated entirely within confined aquifers: physical model and numerical validation. Symmetry. 2026;18(3):507. doi:10.3390/sym18030507. [Google Scholar] [CrossRef]

24. Shao YJ, Ma G, Zhang WQ. Large-scale 3D printed model test on seepage distribution in water diversion tunnel and surrounding fractured rock. Tunn Undergr Space Technol. 2025;158:106389. doi:10.1016/j.tust.2025.106389. [Google Scholar] [CrossRef]

25. Sun YP, Meng GT, Huang XL, Yang HY, Gao F, Hua C, et al. Stability evaluation of a hydraulic tunnel subjected to partial lining removal and variable cross-section excavation: implication from numerical modelling. Int J Numer Anal Methods Geomech. 2026;50(5):2133–50. doi:10.1002/nag.70220. [Google Scholar] [CrossRef]

26. Xiao PW, Yang XG, Li B, Zhou X, Sun YP, Ding XC, et al. Roof arch collapse of underground cavern in fractured rock mass: in situ monitoring and numerical modeling. J Rock Mech Geotech Eng. 2025;17(5):2778–92. doi:10.1016/j.jrmge.2024.05.018. [Google Scholar] [CrossRef]

27. Lin XT, Wu QJ, Cao L, Zhao ZY, Su D, Chen X. Experimental study on failure mechanism of tunnel face during shield tunneling beneath existing underground rectangular structure. Tunn Undergr Space Technol. 2026;167:107072. doi:10.1016/j.tust.2025.107072. [Google Scholar] [CrossRef]

28. Zhu XY, Zhang ZQ, Yu H, Ge H. Water inrush disasters behavior in weakly cemented red-bed tunnels using 3DEC-flow coupled model. Tunn Undergr Space Technol. 2026;170(3):107380. doi:10.1016/j.tust.2025.107380. [Google Scholar] [CrossRef]

29. Wang CW, Liu XL, Hu N, Wang EZ, Yao WL, He ZH. Construction responses and interaction of mechanized twin tunnels considering layout forms: insights from refined three-dimensional numerical modeling. Comput Geotech. 2025;183(3):107203. doi:10.1016/j.compgeo.2025.107203. [Google Scholar] [CrossRef]

30. Tu QZ, Wu W, Wang Y, Dong BS, Zhu HH. Automated numerical simulation-driven dynamic design method for complex geological tunnels. Transp Geotech. 2026;56:101784. doi:10.1016/j.trgeo.2025.101784. [Google Scholar] [CrossRef]

31. Xiao F, Chen XL, Zhu YM, Xie P, Salimzadeh S, Zhang QB. Multi-LoD BIM integrated design framework for pressurised tunnel: hydro-mechanical coupling simulation and sustainability assessment. Tunn Undergr Space Technol. 2025;158(21–26):106404. doi:10.1016/j.tust.2025.106404. [Google Scholar] [CrossRef]

32. Zhang CP, Tu SQ, Liu TX, Li W. Stability of water-resistant rock mass of fault tunnels under high-temperature and high-pressure conditions. Transp Geotech. 2026;57:101893. doi:10.1016/j.trgeo.2026.101893. [Google Scholar] [CrossRef]

33. Luo Y. Influence of water on mechanical behavior of surrounding rock in hard-rock tunnels: an experimental simulation. Eng Geol. 2020;277(B7):105816. doi:10.1016/j.enggeo.2020.105816. [Google Scholar] [CrossRef]

34. Chen ZQ, Zhan HX, Xu GW, Zou XH, Chen KP, He C. Investigation on tunnel face stability and failure mechanism in soft rockmass under coupled high geostress and seepage flow conditions. Tunn Undergr Space Technol. 2026;169(2):107279. doi:10.1016/j.tust.2025.107279. [Google Scholar] [CrossRef]

35. Pujades E, Jurado A. Groundwater-related aspects during the development of deep excavations below the water table: a short review. Undergr Space. 2021;6(1):35–45. doi:10.1016/j.undsp.2019.10.002. [Google Scholar] [CrossRef]

36. Fan HB, Zhu ZG, Song YX, Zhang SY, Zhu YQ, Gao XQ, et al. Water pressure evolution and structural failure characteristics of tunnel lining under hydrodynamic pressure. Eng Fail Anal. 2021;130(2):105747. doi:10.1016/j.engfailanal.2021.105747. [Google Scholar] [CrossRef]

37. Liu X, Fu BM, Wang DC, Jiang AN, Fang Q, Li JY. Analytical study of non-linear seepage field in underwater tunnels with different computational models and boundary conditions. Transp Geotech. 2025;55(3):101674. doi:10.1016/j.trgeo.2025.101674. [Google Scholar] [CrossRef]

38. Ding YC, Zhu HY, Lai HP, Dai ZR. Study on the failure characteristics of tunnel lining structures under water infiltration in loess-mudstone composite stratum. Eng Fail Anal. 2026;186(1):110497. doi:10.1016/j.engfailanal.2025.110497. [Google Scholar] [CrossRef]

39. Niu XQ, Zhang CJ. Some key technical issues on construction of ultra-long deep-buried water conveyance tunnel under complex geological conditions. Tunn Constr. 2019;39(4):523–36. (In Chinese). [Google Scholar]

40. Zhu KX, Zhong ZL, Wang NY, Li ZZ. Mechanical properties and nonlinear energy evolution of holed limestone under hydro-mechanical coupling. Rock Mech Rock Eng. 2025;58(1):807–26. doi:10.1007/s00603-024-04190-0. [Google Scholar] [CrossRef]

41. Fan SY, Song ZP, Xu T, Zhang YW. Investigation of the microstructure damage and mechanical properties evolution of limestone subjected to high-pressure water. Constr Build Mater. 2022;316(5):125871. doi:10.1016/j.conbuildmat.2021.125871. [Google Scholar] [CrossRef]