Parametric Finite Element Investigation of Slope–Tunnel–Foundation Interaction in Layered Soils Under Hydro-Mechanical Coupling

Traditional geotechnical practice ensures stability by confirming that the resisting capacity of the ground sufficiently outweighs the driving forces through a predefined safety margin. In contrast, when utilizing numerical frameworks like the FEM or FDM, safety is not an inherent output, as the soil's strength is modeled as a fundamental material property integrated with its stress-strain response.

In contrast, the strength reduction method (SRM) has become a standard technique in numerical modeling for determining the factor of safety (F_s). This procedure involves the systematic degradation of the soil’s shear strength components—specifically the cohesion (c) and the internal friction angle (φ)—until the system reaches a state of instability. The degree of this reduction is governed by a strength reduction factor (SRF), which is mathematically expressed as follows (Eq. (1)):

$\sum M s f=F_s=\frac{\tan g \phi_{\text {input }}}{\tan g \phi_{\text {reduced }}}=\frac{c_{\text {input }}}{c_{\text {reduced }}}$      (1)

In this formulation, the subscript input denotes the original shear strength values, while reduced signifies the factored parameters applied during the simulation. The procedure begins with an F_s of unity, which is then progressively scaled up. This iterative process continues until the numerical model fails to reach equilibrium—characterized by non-convergence or disproportionate strain—which signifies the onset of collapse. The specific magnitude of F_s recorded at this limit state is identified as the overall safety factor of the slope.

4.1 Before the construction of the tunnel

The numerical domain was defined with dimensions of 40 m in width and 12 m in height for the 2D simulations, while the 3D model incorporated a longitudinal length of 8 m. The stratigraphy is composed of three distinct horizontal strata: a 3 m thick basal clay layer, a 5 m intermediate peat layer, and a 4 m superficial sand layer. The slope geometry is characterized by a height (H) of 4 m and an 8 m horizontal reach. To represent the soil’s constitutive response, a Mohr–Coulomb elastic-perfectly plastic framework was utilized. Standard boundary conditions were enforced, including fixed constraints and symmetry considerations where applicable. For this specific scenario, the influence of the phreatic surface was not considered. The geomechanical parameters for each layer are detailed in Table 1.

Table 1. The characteristics of the soils used

Before tunnel excavation, the F_s calculated using the PLAXIS 3D computational suite was 1.44, whereas the PLAXIS 2D simulation yielded a value of 1.41. This reflects a marginal increase of approximately 1.62% in the three-dimensional analysis compared to the 2D results. The close agreement between these two numerical approaches indicates that the dimensionality of the model has a negligible impact on the initial stability assessment, a finding that aligns with previous research in the field [14].

Figures 2 and 3 show the results obtained by both PLAXIS 2D and PLAXIS 3DT finite element analyses, respectively, in terms of increments of total displacement. It can be observed that the same breaking circle is obtained by both software, but the displacements are a little different. The total displacements obtained using PLAXIS 2D equal 2.94 mm (Figure 2), while the total displacements obtained using PLAXIS 3DT (PLAXIS 3D Tunnel module) equal 10.5 mm (Figure 3). The comparison shows that 3D modeling provides a more realistic representation of deformation mechanisms due to three-dimensional stress redistribution, without implying superior accuracy.

Figure 2. The total displacement increments before the construction of the tunnel using PLAXIS 2D finite element

Figure 3. The total displacement increments before the construction of the tunnel using PLAXIS 3DT (PLAXIS 3D Tunnel module), finite element

To ensure the reliability of the numerical results, the obtained findings were compared with previously published studies. Similar trends in slope stability reduction due to tunnel excavation were reported by Zhou et al. [10] and Srivastav et al. [15], confirming the influence of underground excavation on stress redistribution and deformation patterns.

The small difference in safety factor between 2D and 3D analyses is consistent with results reported in the literature, where global stability is not significantly affected by dimensionality, while deformation fields are strongly influenced by 3D effects such as arching and confinement.

Therefore, PLAXIS 3DT software is used in all the following work.

4.2 After the excavation of the tunnel

A linear elastic constitutive law was adopted to simulate the structural response. Figure 4 illustrates the discretized geometry employed for the finite element mesh. To reflect typical TBM excavation depths, the tunnel axis is situated 8 m below the ground surface. The support system consists of a 0.35 m thick (e) reinforced concrete lining, defined by a unit weight of 2500 kg/m³, a stiffness modulus of 30 GPa, and a Poisson’s coefficient of 0.15.

Figure 4. The total displacement increments after the excavation of the tunnel

The safety coefficients F_s obtained before the construction of the tunnel are F_s = 1.44, and after the construction of the tunnel is F_s = 1.45. The safety factor remains unchanged after the construction of the tunnel because the slope is stable before and after the construction of the tunnel, and because the tunnel is far from the rupture circle (reference model). Figures 3 and 4 show the increments of total settlements obtained before and after the construction of the tunnel of 6.60 mm and 10.5 mm, respectively. It can be observed that the settlements are more significant in the areas closest to the tunnel, and that the significant settlements are at the ground surface and propagate towards the tunnel excavation. This indicates that the displacements increase after tunnel excavation due to stress redistribution and ground relaxation.

4.3 Location of the tunnel

Figure 5 demonstrates that F_sremain stable (approximately 1.43–1.45) when tunnels are placed at distances 0, 2D, or 4D, similar to the pre-construction baseline of 1.44. However, positioning the tunnel at 5D, between the crest and mid-slope, significantly impacts stability, leading to failure as shown in Figure 5(c). The proximity of the tunnel to the critical failure mechanism dictates the slope's stability, with greater distances ensuring higher F_s.

Figures 6(a) and (b) illustrate that while the safety factor experiences a marginal reduction as tunnel position increases up to X_t = 4D, a significant drop in stability occurs beyond this point. This critical threshold results in a sharp increase in total settlement, reaching 200 mm, and signifies complete slope failure.

X_t = 0 represents the least detrimental position for the studied configuration to reduce the coupling effect between the tunnel and the slope. Therefore, this tunnel position was adopted for the subsequent analyses.

(a) X_t = 2D

(b) X_t = 4D

(c) X_t = 5D

Figure 5. Effect of tunnel position on the failure mechanism of the slope

(a) Safety factor 

(b) Total displacements

Figure 6. Variation of safety factor (F_s) and total displacements for different tunnel locations

4.4 Effect of slope angle (β)

The sensitivity of the system's stability to the gradient of the slope β during the excavation phase was assessed by testing a range of inclinations: (11.31°, 18.43°, 26.56°, 33.39°, and 38.66°) as presented in Figure 7(a), there is an inverse relationship between the slope angle and the safety margin (F_s); a steeper profile leads to a progressive reduction in stability, culminating in a threshold of failure. For example, a configuration with a slope of 39° yields an F_s of 1.39. Furthermore, Figure 7(b) demonstrates the correlation between stability levels and kinematic response, showing that total displacements are amplified as the safety factor diminishes in response to more aggressive slope geometries.

(a) Safety factor 

(b) Total displacements

Figure 7. Variation safety factor (F_s) and total displacements for different slope inclination (β)

4.5 Effect of loading increments (q)

The influence of foundation-induced surcharges on the stability of the integrated slope-tunnel model was examined by applying a series of pressure increments, ranging from 10 to 60 kPa. The numerical simulation incorporates several key hydraulic assumptions: a steady-state flow regime is maintained within the porous medium, and the soil strata exhibit isotropic permeability (k_x = k_y). Furthermore, the tunnel boundary is modeled as fully pervious; it is kept at internal atmospheric pressure to facilitate the instantaneous removal of seepage. This drainage condition ensures that the structural lining remains free from internal hydrostatic thrust.

The foundation is made of reinforced concrete with a thickness of 0.5 m, and modeled using PLAXIS 3D by defining the "plate" element with a flexural stiffness EI = 2.92 × 10⁵ kN·m²/m and an axial stiffness EA = 1.4 × 10⁷ kN/m.

The impact of external loading on the stability of the slope is investigated through two main configurations: the pre-construction phase and the post-tunneling phase. Within the second framework, two distinct scenarios are compared: one ignoring hydrological influences and the other incorporating the groundwater table. The performance of the system is evaluated based on total and effective stress distributions, cumulative displacements, and global stability margins.

4.5.1 Analysis of total stresses and effective stresses

To assess how foundation surcharges modify the stress state within the soil, a specific pressure of q = 60 kPa was applied. This magnitude corresponds to the ultimate limit of the soil's bearing capacity. Data extraction points were strategically positioned at the tunnel-soil interface and beneath the foundation structure to monitor these variations.

In the initial state, the phreatic surface is assumed to coincide with the base of the clay stratum. Following the tunnel's construction, the analysis differentiates between two distinct hydrological states: (i) a condition where groundwater effects are omitted and (ii) a scenario involving active water seepage.

(a)

(b)

(c)

Figure 8. Groundwater (pore-water) pressure distribution under three conditions: (a) initial state (pre-excavation); (b) post-excavation without groundwater inflow; (c) post-excavation with groundwater inflow

Within a homogeneous and isotropic domain (characterized by a uniform permeability k), the water table in the unconfined aquifer is initially horizontal at a depth of 8 m, as depicted in Figures 8(a) and (b). The distribution of pore water pressure remains identical between the pre-excavation phase in Figure 8(a) and the post-excavation phase without inflow in Figure 8(b). This stability is due to the lack of hydraulic head variation, as no drainage occurs through the tunnel lining.

In contrast, when groundwater inflow is permitted after excavation, a localized depression in the water table is observed. This results in a non-linear phreatic surface that curves downward toward the tunnel crown, as illustrated in Figure 8(c).

While tunnel excavation typically triggers a redistribution of stresses within the surrounding ground [16], the effective stress levels in this study remain nearly constant between configurations in Figures 9(a) and (b), hovering around −85 kPa. This consistency is attributed to the tunnel's location, which lies beyond the influence zone of the slope's critical failure surface.

The data presented in Figure 9 confirm that, under a surface surcharge of q = 60 kPa, the peak effective stress reaches approximately −85 kPa for both the pre-excavation state and the post-excavation phase without hydraulic effects. In contrast, the inclusion of groundwater seepage (Figure 9(c)) leads to a significant reduction in effective stress, dropping to roughly −52 kPa.

This decline is primarily driven by the buildup of pore water pressure resulting from infiltration toward the tunnel. Consequently, the discrepancy in effective stress between Figures 9(b) and (c) is governed by hydraulic processes rather than mechanical excavation itself, as the tunnel geometry remains constant in both scenarios.

(a)

(b)

(c)

Figure 9. Variation of effective stresses (σ’_yy, in kPa) of the slope–tunnel system due to overload before and after tunnel excavation

(a)

(b)

Figure 10. Variation of the mean total stress (σ_yy, in kPa) for the different cases of tunnel excavation

Figure 11. Variation of mean total stress (σ_yy, in kPa) according to the overload, for the different cases of tunnel construction

The effect of foundation-induced loads on the total stress distribution is illustrated in Figure 10. Before tunnel excavation, the applied loads are progressively distributed with depth, following the classical stress diffusion pattern in soils. The maximum total stress is observed for a surcharge of q = 60 kPa, reaching a value of σ_yy = −80 kPa (Figure 10(a)).

After tunnel excavation, and when groundwater effects are taken into account, the stress distribution is significantly modified. The maximum total stress increases to approximately σ_yy = −130 kPa (Figure 10(b)), indicating a concentration of stresses around the tunnel and beneath the loaded area.

In addition, localized stress relief (decompression zones) is observed at the tunnel invert and along the tunnel contour. These decompression zones may weaken the surrounding soil structure and can contribute to instability mechanisms. In particular, they may promote tunnel collapse in saturated sandy soils, as reported in previous studies [15].

Figure 11 illustrates the variation of total stresses under different loading and excavation conditions.

The results show that before tunnel excavation and in the absence of groundwater effects, the stress distribution follows a typical gradual increase with depth, corresponding to classical stress diffusion in soils.

However, after tunnel excavation combined with groundwater infiltration, a significant increase in total stress is observed, particularly around the tunnel and beneath the loaded area. This behavior can be attributed to:

• Stress redistribution induced by tunnel excavation,
• Reduction in effective stress due to increased pore water pressure,
• Concentration of stresses around the tunnel structure.

Additionally, localized decompression zones appear near the tunnel boundary, indicating a loss of confinement that may contribute to instability mechanisms.

The results indicate that groundwater infiltration significantly amplifies stress redistribution, leading to increased total stress concentrations and the development of localized decompression zones around the tunnel, which may promote instability.

4.5.2 Field of displacement and settlement

To assess how surface surcharges influence excavation-induced ground movements, the soil's peak bearing pressure of q = 60 kPa was applied as a reference load. As established in the preceding stress analysis, the tunnelling process by itself does not drastically alter the overall displacement or settlement patterns. Nevertheless, the study examines two primary configurations before and after excavation, with a specific emphasis on the role of hydraulic seepage and its impact on deformation.

(a)

(b)

(c)

Figure 12. The values of total displacements (mm) of the tunnel slope assembly under the effect of load q = 60 kPa

Figure 12 illustrates the cumulative displacement patterns across the various simulated scenarios. When hydraulic influences are omitted, the peak settlement declines to roughly 280 mm, as depicted in Figures 12(a) and (b). This discrepancy emphasizes the critical role of water seepage in exacerbating ground deformations and intensifying settlement patterns within the integrated slope–tunnel framework.

The observed increase in settlement under groundwater conditions is also in agreement with findings reported by Rodríguez et al. [17], where pore pressure variations significantly impact soil deformation.

The results show that peak settlement occurs during tunnel excavation when hydraulic effects are included, reaching a magnitude of nearly 320 mm (Figure 12(c)). These deformations are primarily localized under the foundation area and follow the slope's critical slip surface. Displacement values diminish progressively with depth, reaching roughly −20 mm at the tunnel invert.

This reversal in displacement sign highlights the development of localized heave or decompression zones in the vicinity of the tunnel, which may be indicative of excavation-induced instability mechanisms [18].

Figure 13. Distribution of maximum settlements according to the loads applied to the surface of a slope due to tunnel excavation

Figure 14. Variation of factor of safety (F_s) as a function of the loads applied to the surface of a slope due to tunnel excavation

Figure 13 presents the variation of maximum settlement as a function of applied surface loading. The results clearly show that the settlement increases with increasing load (expected behavior), the presence of groundwater significantly amplifies settlement values of ~280 mm without groundwater, and of ~320 mm with groundwater.

This increase is mainly due to the reduction in effective stress caused by pore water pressure, increased compressibility of saturated soil, and stress release and relaxation induced by tunnel excavation.

The settlement profile also exhibits a nonlinear trend, reflecting the complex interaction between loading, soil properties, and hydraulic conditions.

The settlement evolution exhibits a nonlinear behavior, strongly influenced by the combined effects of surface loading, tunnel excavation, and groundwater infiltration, with the latter significantly increasing deformation.

The results indicate the safety factor decreases progressively with increasing load (Figure 14). The reduction is more pronounced when groundwater infiltration is considered. This behavior is explained by:

• Increased driving forces due to surcharge loading,
• Reduction in shear strength caused by decreasing effective stress under saturated conditions.

The safety factor decreases progressively with increasing surcharge, with a more significant reduction observed under groundwater conditions due to the decrease in effective stress and soil shear strength.

The slope remains stable under a 30 kPa foundation load and a 26.56° angle, yielding a 1.347 safety factor that exceeds the required 1.3 threshold [19-21]. This indicates a satisfactory margin of stability under these specific conditions.

The combined analysis of stress distribution, settlement, and safety factor highlights the dominant role of hydro-mechanical coupling in slope–tunnel systems. While mechanical loading alone induces moderate changes, the presence of groundwater significantly amplifies both deformation and instability risks. These findings demonstrate that neglecting hydraulic effects may lead to a substantial underestimation of failure potential in practical engineering applications.