Stability Analysis of Three-Dimensional Tunnel Face Considering Linear and Nonlinear Strength in Unsaturated Soil

Abstract

The shear strength of unsaturated soils exhibits significant nonlinearity, while previous studies often simplified it with linear strength models. The objective of this paper is to investigate the distinctions in the stability of three-dimensional (3D) tunnel faces when using linear and nonlinear strength models. A new 3D rotational failure mechanism and an extended form of the Mohr–Coulomb (M-C) failure criterion were integrated into the kinematically limited analysis (KLA) framework to describe the failure characteristics of tunnel faces. Subsequently, the factor of safety (FS) of the 3D tunnel faces was calculated using the strength reduction method (SRM). In the discussion section, the impacts of nonlinear shear strength, matric suction in the unsaturated soils, and the 3D geometric parameters of the tunnel on the stability of the tunnel face were analyzed. The outcomes indicate that, in unsaturated soil conditions, diverse nonlinear strength calculation models and soil types exert disparate influences on the FS of 3D tunnel faces. The main novelty of this study lies in establishing an effective method for assessing the stability of tunnel faces in unsaturated soils.

1. Introduction

A substantial body of the literature includes two-dimensional (2D) analyses concerning the stability of tunnel faces, typically presuming either saturated or completely dry soil conditions [1,2]. Nonetheless, it is essential to recognize that the stability analysis of tunnel faces is intrinsically a three-dimensional (3D) problem. Neglecting spatial effects in the analysis framework of a 2D failure mechanism is evidently impractical in real tunnel engineering scenarios [3]. Furthermore, the soil is often in an unsaturated state due to the influence of matric suction. The strength properties in unsaturated soils differ significantly from those in saturated soils [4,5,6]. Hence, a series of innovative works have been introduced within the academic community to tackle the characteristics of unsaturated soils [7]. To develop a more accurate scheme, it is essential to integrate the influence of unsaturated components into the 3D failure mechanism.

The analysis methods for assessing the stability of geotechnical problems include limit equilibrium (LE) [8,9,10,11], experimental modeling (ME) [12], computational simulation (CS) [13,14], and kinematically limited analysis (KLA) [15,16,17,18]. Among these methods, LE is commonly employed for investigating the stability of tunnel faces by solving static equilibrium equations. For instance, Oreste et al. [19] developed a new triangular prism model based on LE to determine the factor of safety (FS) for tunnel faces reinforced by positioning lines. Anagnostou and Perazzelli [20] utilized a wedge model to assess the stability of tunnel faces under the influence of Clear Span and nonuniform support. Lu et al. [21] focused on calculating the lowest supporting force necessary for tunnel faces experiencing seepage conditions, addressing stability concerns. In contrast to the LE, the KLA has the ability to optimize the failure mechanism and has been widely used in various geotechnical engineering scenarios [22,23,24]. During the course of optimizing the kinematic mechanism, Leca and Domieux [25] were the pioneers involved in proposing a 3D kinematic mechanism involving double blocks, which paved the way for more advanced analyses in tunnel stability assessments. On this basis, Mollon et al. [26] proposed a modified multi-block mechanism to obtain more accurate solutions. The horned failure mechanism has been utilized in the examination of tunnel face safety, and its results are in good agreement with both numerical simulations and experimental observations [27].

In the context of the New Austrian Tunnel Method (NAT), the tunnel face is not subjected to supporting pressure [28,29]. Consequently, evaluating the stability of the tunnel face using the FS is a more appropriate approach than considering supporting pressures. There are two approaches for delineating the FS: the strength reduction method (SRM) and the gravity augmentation method (GAM). In the GAM, the FS is expounded as the quotient of the internal energy dissipation rate to the external energy dissipation rate. Nonetheless, within the SRM paradigm, the FS is articulated as the quotient of the real soil strength to the critical strength in the state of limit equilibrium. This approach yields a more cautious resolution when contrasted with the GAM, consequently rendering SRM a prevalent choice for assessing the stability of Mohr–Coulomb (M-C) materials.

Numerous scholars have delved into the influence of matric suction on shear strength within the domain of geotechnical engineering, particularly concerning the unsaturated constituents of soil strength [30,31,32,33]. To account for unsaturated properties in soils, Fredlund et al. [34] modified the M-C criterion and proposed the concept of matric suction angle. Following this, various estimation formulas were advanced to assess the strengths of unsaturated soils across diverse soil classifications [35,36,37,38,39]. Extensive research has been conducted in the academic realm concerning the analysis of unsaturated soil strength, with a focus on topics such as slope stability and retaining walls. However, there is a noticeable lack of studies addressing the safety of tunnel faces. As an instance, Zhang et al. [40] determined the critical groundwater level depth for retaining walls under varying distributions of matric suction in unsaturated soil, utilizing the temporal equilibrium method. Additionally, Oh and Lu [41] explored the consequences of rainfall-induced failures on engineering slopes under unsaturated conditions. These investigations highlight the significance of matric suction on influencing soil strength and motivate the author to incorporate unsaturated components into tunnel face safety analyses in order to attain more precise solutions.

The primary objective of this study was to investigate the disparities in the stability of three-dimensional (3D) tunnel faces when employing linear and nonlinear strength models. To achieve this, a novel 3D rotational failure mechanism was developed to depict the failure of tunnel faces. In order to account for unsaturation in shear strength, an extended Mohr–Coulomb (M-C) criterion was adopted in this paper. Additionally, the strength reduction method (SRM) was utilized to compute the factor of safety (FS) of tunnel faces. Finally, the study examined the impact of the nonlinear shear strength of the soil, tunnel geometry characteristics, and matric suction in unsaturated soil on the safety of tunnel faces. The key novelty of this study lies in establishing an effective approach to assess the stability of tunnel faces in unsaturated soils.

2. Strength Theory of Unsaturated Soil

Linear and nonlinear forms of the unsaturated soil strength formula can be distinguished [42]. And both forms of equations are employed in this paper to calculate the impact of matric suction on tunnel face safety in unsaturated conditions.

2.1. Linear Shear Strength

In the linear cases, two forms of matric suction distribution are considered: uniform distribution and linear increase with burial depth. Fredlund et al. [34] introduced a formula for calculating soil strength under uniform matric suction based on the M-C failure criterion:

$$\tau =c^{\prime }+({\sigma }_n-u_a)\tan {\phi }^{\prime }+(u_a-u_w)\tan {\phi }^b$$

(1)

where $\tau$ denotes the shear strength of unsaturated soil, while ${\phi }^{\prime }$ and $c^{\prime }$ represent to the soil’s friction angle and cohesion, respectively. $({\sigma }_n-u_a)$ represents the net normal pressure, and $(u_a-u_w)$ represents the matric suction, with $u_w$ representing the pore water pressure and $u_a$ representing the pore gas pressure. ${\phi }^b$ is commonly referred to as the matric suction angle, which reflects the influence of matric suction on shear strength. According to previous studies [43], ${\phi }^b$ is often assumed to have a constant value of 15°.

Equation (1) can also be reformulated in the M-C criterion form as follows:

$$\tau =c+(\sigma -u_a)\tan {\phi }^{\prime }$$

(2)

where the cohesive force, accounting for unsaturated strength components, is redefined as follows:

$$c=c^{\prime }+(u_a-u_w)\tan {\phi }^b$$

(3)

Furthermore, apart from the uniformly distributed matric suction, Zhang et al. [40] introduced an alternative linear distribution form, depicted in Figure 1. In this model, matric suction linearly increases with the buried depth, commencing from 0 at the groundwater level. This relationship can be described as follows:

$${(u_a-u_w)}_z=\rho z$$

(4)

In which $\rho$ represents the constant coefficient for matric suction linearly increasing with depth.

2.2. Nonlinear Shear Strength

The shear strength of soil in unsaturated states is conventionally depicted in a linear fashion; nevertheless, a substantial corpus of empirical evidence attests to the pronounced nonlinearity exhibited by shear strength of soil in unsaturated conditions, particularly when subjected to the effects of matric suction [44,45]. This article posits several nonlinear equations to delineate the inter-relation between the shear strength of unsaturated soil and matric suction. Zhang et al. [46] encapsulated a selection of emblematic nonlinear equations, as shown in

Table 1. Several formulas related to apparent cohesion, as compiled by Zhang et al. (2016) [40].

Author and Year	Apparent Cohesion Induced by Matric Suction
Fredlund et al. (1978) [34]	$c_{app}=(u_a-u_w)\tan {\phi }^b$
Fredlund et al. (1996) [36]	$\begin{array}{l}{c}_{app}=(u_a-u_w){{\Theta }_d}^{\kappa }\tan {\phi }^{\prime }\\ \mathrm{where} \kappa \mathrm{is} \mathrm{a} \mathrm{fitting} \mathrm{parameter}, {\Theta }_d=\frac{{\theta }_w}{{\theta }_s}\end{array}$ ${\theta }_w=\left[1-\frac{\ln \left(1+\frac{u_a-u_w}{{\psi }_r}\right)}{\ln \left(1+\frac{{10}^6}{{\psi }_r}\right)}\right]\frac{{\theta }_s}{{\left\{\ln \left[\exp \left(1\right)+{\left(\frac{u_a-u_w}{a_f}\right)}^{n_f}\right]\right\}}^{m_f}}$
Vanapalli et al. (1996) [35]	$c_{app}=(u_a-u_w)\left[\tan {\phi }^{\prime }\left(\frac{{\theta }_w-{\theta }_r}{{\theta }_s-{\theta }_r}\right)\right]$ $\mathrm{where} {\theta }_w=\left[1-\frac{\ln \left(1+\frac{u_a-u_w}{{\psi }_r}\right)}{\ln \left(1+\frac{{10}^6}{{\psi }_r}\right)}\right]\frac{{\theta }_s}{{\left\{\ln \left[\exp \left(1\right)+{\left(\frac{u_a-u_w}{a_f}\right)}^{n_f}\right]\right\}}^{m_f}}$
Vilar (2006) [39]	$\begin{array}{l}{c}_{app}=\frac{\left(u_a-u_w\right)}{a+b\left(u_a-u_w\right)}\\ \mathrm{where} a=\frac{1}{\tan {\phi }^{\prime }}, b=\frac{1}{c_{ult}-c^{\prime }}\end{array}$
Khalili and Khabbaz (1998) [37]	$\begin{array}{l}{c}_{app}=\left(u_a-u_w\right)\left[{\lambda }^{\prime }\right]\tan {\phi }^{\prime }\\ \mathrm{where} {\lambda }^{\prime }=\left\{\begin{array}{ll}1.0& \mathrm{when} \left(u_a-u_w\right)\le AEV\\ {\left(\frac{u_a-u_w}{AEV}\right)}^{-0.55}& \mathrm{when} \left(u_a-u_w\right)>AEV\end{array}\right.\end{array}$
Bao et al. (1998) [38]	$\begin{array}{l}{c}_{app}=\left(u_a-u_w\right)\left[\zeta \right]\tan {\phi }^{\prime }\\ \mathrm{where} \zeta =\left\{\begin{array}{ll}1.0& \mathrm{if} \left(u_a-u_w\right)\le AEV \\ \frac{\log {\left(u_a-u_w\right)}_r-\log \left(u_a-u_w\right)}{\log {\left(u_a-u_w\right)}_r-\log \left(AEV\right)}& \mathrm{if} \mathrm{AEV}<\left(u_a-u_w\right)<{\psi }_r\\ 0& \mathrm{if} \left(u_a-u_w\right)\ge {\psi }_r\end{array}\right.\end{array}$

. In this context, $S$ signifies saturation, while $S_r$ signifies residual saturation. ${\theta }_w$, ${\theta }_s$, and ${\theta }_r$ correspond to volume moisture content, saturated volume moisture content, and residual volume moisture content, respectively. Furthermore, $c_{ult}$ designates the ultimate undrained soil shear strength in a desiccated state, and ${\psi }_r$ designates residual matric suction. And AEV refers to the Air Entry Value.

Table 2. SWCC parameters by Fredlund and Xing (1994) [47].

SWCC	${\mathit{\theta }}_{\mathit{s}}$	${\mathit{a}}_{\mathit{f}}(\mathbf{kPa})$	${\mathit{n}}_{\mathit{f}}$	${\mathit{m}}_{\mathit{f}}$	${\mathit{\psi }}_{\mathit{r}}$	$\mathbf{AEV}\left(\mathbf{kPa}\right)$	${\mathit{\theta }}_{\mathit{s}}$	$\mathit{\kappa }$	${\mathit{c}}_{\mathit{u}\mathit{l}\mathit{t}}(\mathbf{kPa})$
1	0.4	1	2	1	10	0.5	0.081	1.0	11.4
2	0.4	10	2	1	100	5	0.080	1.8	13.7
3	0.4	100	2	1	1000	50	0.078	2.2	28.4
4	0.4	1000	2	1	10,000	500	0.073	2.5	107.3

enumerates the parameters characterizing the Soil–Water Characteristic Curve (SWCC) for four distinct soil classifications in accordance with the study by Fredlund and Xing [47]. Specifically, NO.1 corresponds to sandy soil, NO.2 corresponds to silt, NO.3 corresponds to clay, and NO.4 corresponds to ultra-fine-grained soil. In Table 2, $a_f$, $n_f$, and $\kappa$ stand for the suction value, slope, and fitting parameter, respectively, which correspond to the inflection point in Figure 2b. It can be observed that various soil types exhibit substantial disparities in their shear strength profiles.

Furthermore, for the aforementioned representative soil specimens, Figure 2a illustrates the nonlinear correlation existing between matric suction and shear strength. Meanwhile, Figure 2b delineates the SWCC of a prototypical soil sample. It is clear from the AEV that shear strength demonstrates a notably linear association with saturated soil conditions. Upon surpassing the critical AEV, the shear strength envelope of the soil shows a nonuniform expansion pattern within the nonlinear domain. As the matric suction approaches the residual zone, sand demonstrates a diminishing tendency with the augmentation of matric suction, whereas plastic soil and silt exhibit a proclivity toward horizontal dispersion.

3. KLA in the 3D Horned Failure Mechanism

Figure 3 serves as a schematic representation of the practical tunnel engineering issue in a water-rich region. A cylindrical tunnel with a diameter denoted as D was excavated above the groundwater level devoid of any supporting pressure. H represents the burial depth of the tunnel, indicating the vertical measurement from the groundwater level to the tunnel’s base. In this paper, we can think of the soil’s unit weight as a fixed amount. And it is hypothesized that homogenous, unsaturated soil clings to an enlarged definition of the M-C yield function, as shown in Equation (3).

3.1. Three-Dimensional Horned Rotational Failure Mechanism

In the investigation of tunnel face safety through the KLA approach, it is imperative to ascertain a rational presumed failure mechanism. This mechanism should adhere to the following assumptions in the KLA: (1) the soil mass complies with pertinent regulations; (2) the soil exhibits ideal plastic behavior. Additionally, there is one need that has to be met, namely the work exerted by outside forces can be no greater than the rate of energy dissipation., as follows:

$${\displaystyle {\int }_V{\sigma }_{ij}{\dot{\epsilon }}_{ij}dV}\ge {\displaystyle {\int }_S{T}_i{v}_{i}dS}+{\displaystyle {\int }_V{X}_i{v}_{i}dV}$$

(5)

where ${\dot{\epsilon }}_{ij}$ and ${\sigma }_{ij}$ indeed represent the strain and stress tensors, respectively. The failure body’s volume and border are precisely denoted by the $V$ and $S$. And $T_i$ appropriately refers to the forces acting on the boundary face, while $X_i$ correctly refers to the body forces acting within the volume $V$. $v_i$ accurately signifies the velocity along the failure face.

Numerous failure mechanisms have been posited within academic sphere as means to emulate tunnel collapse. When using discrete techniques, the horned mechanism differs from these others in that it is more conservative than the multi-block mechanism and easier to design than the rotating mechanism. Michalowski initially introduced the horned failure mechanism, which subsequently evolved into a potent tool for evaluating both slope stability and tunnel face safety. Therefore, considering its efficacy in practical application, this paper opts for the adoption of the horned failure mechanism.

This mechanism encompasses a pair of logarithmic spirals, denoted as AE and BE, sharing the central pivot at O and characterized by the angular velocity of ω. Schematic representations for the horizontal and longitudinal cross-sectional aspects are shown in Figure 3. In the realm of polar coordinates, the expressions for AE and BE are articulated as follows:

$$r=r_{\mathrm{A}}{\mathrm{e}}^{-\left(\theta -{\theta }_{\mathrm{A}}\right)\tan {\phi }^{\prime }}$$

(6)

$$r^{\prime }=r_{\mathrm{B}}{\mathrm{e}}^{\left(\theta -{\theta }_{\mathrm{B}}\right)\tan {\phi }^{\prime }}$$

(7)

In which $(r_{\mathrm{A}}, {\theta }_{\mathrm{A}})$ and $(r_{\mathrm{B}}, {\theta }_{\mathrm{B}})$ denote the polar coordinates of points A and B, respectively. The polar radii, $r_A$ and $r_B$, can be determined using the following equations:

$$r_{\mathrm{A}}=\frac{\sin {\theta }_{\mathrm{B}}}{\sin \left({\theta }_{\mathrm{B}}-{\theta }_{\mathrm{A}}\right)}D$$

(8)

$$r_{\mathrm{B}}=\frac{\sin {\theta }_{\mathrm{A}}}{\sin \left({\theta }_{\mathrm{B}}-{\theta }_{\mathrm{A}}\right)}D$$

(9)

The distance between the rotational center O and the horned mechanism’s center axis is as follows:

$$r_{\mathrm{m}}=\frac{r+r^{\prime }}{2}$$

(10)

Each rotational segment’s radius is as follows:

$$R=\frac{r-r^{\prime }}{2}$$

(11)

The coordinates of the intersection point of the two logarithmic spirals are as follows:

$$\left\{\begin{array}{l}{\theta }_E=\frac{{\theta }_A+{\theta }_B}{2}+\frac{\ln \left(\sin {\theta }_B/\sin {\theta }_A\right)}{2\tan {\phi }^{\prime }}\\ r_E=\frac{\sqrt{\sin {\theta }_A\sin {\theta }_B}}{\sin \left({\theta }_B-{\theta }_A\right)}{e}^{\frac{1}{2}\left({\theta }_A-{\theta }_B\right)\tan {\phi }^{\prime }}D\end{array}\right.$$

(12)

When the tunnel burial depth appears to be shallow, the horned mechanism might emerge above the ground face. With such instances, the ground face CD, the tunnel face AB, and the two logarithmic curves AC and BD serve to define the bounds for this mechanism. The parameters ${\theta }_C$ and ${\theta }_D$ might be used separately for determining the coordinates of points C and D, where this mechanism intersects the surface of the ground, in the following manner:

$$C-r_B\cos {\theta }_B+r_A\cos {\theta }_C{e}^{-\left({\theta }_C-{\theta }_A\right)\tan {\phi }^{\prime }}=0$$

(13)

$$C-r_B\cos {\theta }_B+r_B\cos {\theta }_D{e}^{\left({\theta }_D-{\theta }_B\right)\tan {\phi }^{\prime }}=0$$

(14)

As depicted in Figure 3, when the horn faces, the failure mechanism primarily involves three constituents, whereas if the horn remains submerged, it solely comprises two elements. Two autonomous variables, ${\theta }_A$ and ${\theta }_B$, determine the configuration of the horned failure mechanism, and this configuration is then refined and resolved based on these two parameters.

3.2. Balance Equation for Work and Energy

In accordance with the KLA method, the balance equation for work and energy arises by equating the rate of energy dissipation to the rate of external work. And the self-weight of soil is a primary contributor to the rate of outside work. Meanwhile, the internal energy dissipation consists of two constituents: effective internal cohesion $c^{\prime }$ and matric suction-induced apparent internal cohesion $c_{app}$. The formulation of the balance equation for work and energy is defined as follows:

$${\dot{D}}_{\mathrm{c}}+{\dot{D}}_{\mathrm{s}}={\dot{W}}_{\gamma }$$

(15)

where ${\dot{D}}_{\mathrm{c}}$ and ${\dot{D}}_{\mathrm{s}}$ correspondingly represent the rates of energy dissipation attributed to effective internal cohesion $c^{\prime }$ and matric suction-induced apparent internal cohesion $c_{app}$. ${\dot{W}}_{\gamma }$ signifies the external work generated by one’s own weight. And the rate of outside work equals the summation of the gravitational work rate for each component, as illustrated below:

$$\begin{array}{cc}\hfill {\dot{W}}_{\gamma }& =2\omega \gamma \left[{\displaystyle {\int }_{{\theta }_{\mathrm{A}}}^{{\theta }_{\mathrm{B}}}{\displaystyle {\int }_{d_1}^R{\displaystyle {\int }_0^{\sqrt{R^2-y^2}}{\left(r_{\mathrm{m}}+y\right)}^2\sin \theta \mathrm{d}x\mathrm{d}y\mathrm{d}\theta }}}+{\displaystyle {\int }_{{\theta }_{\mathrm{B}}}^{{\theta }_{\mathrm{C}}}{\displaystyle {\int }_{-R}^R{\displaystyle {\int }_0^{\sqrt{R^2-y^2}}{\left(r_{\mathrm{m}}+y\right)}^2\sin \theta \mathrm{d}x\mathrm{d}y\mathrm{d}\theta }}}\right.\hfill \\ & \left.+2\omega \gamma {\displaystyle {\int }_{{\theta }_{\mathrm{C}}}^{{\theta }_{\mathrm{D}}}{\displaystyle {\int }_{-R}^{d_2}{\displaystyle {\int }_0^{\sqrt{R^2-y^2}}{\left(r_{\mathrm{m}}+y\right)}^2\sin \theta \mathrm{d}x\mathrm{d}y\mathrm{d}\theta }}}\right]\hfill \end{array}$$

(16)

It is imperative to establish the basic premise that these collapsing blocks exhibit rigidity in order to determine the rate of internal energy dissipation. This premise enables the integration for the rate of internal energy dissipation as follows:

$$\dot{D}={\displaystyle {\int }_{S_{\mathrm{t}}}cv\cos {\phi }^{\prime }}\mathrm{d}{S}_{\mathrm{t}}$$

(17)

Here, $\mathrm{c}$ signifies the cohesive strength of the soil, and $S_{\mathrm{t}}$ signifies the failure face. Using Equation (17), it is possible to deduce the energy dissipation resulting from the effective cohesion of the soil $c^{\prime }$.

$$\begin{array}{cc}\hfill {\dot{D}}_{\mathrm{c}}& =2\omega c^{\prime }\left[{\displaystyle {\int }_{{\theta }_{\mathrm{A}}}^{{\theta }_{\mathrm{B}}}{\displaystyle {\int }_{d_1}^R\frac{R{\left(r_{\mathrm{m}}+y\right)}^2}{\sqrt{R^2-y^2}}\mathrm{d}y\mathrm{d}\theta }}\right.+{\displaystyle {\int }_{{\theta }_{\mathrm{B}}}^{{\theta }_{\mathrm{C}}}{\displaystyle {\int }_{-R}^R\frac{R{\left(r_{\mathrm{m}}+y\right)}^2}{\sqrt{R^2-y^2}}\mathrm{d}y\mathrm{d}\theta }}\hfill \\ & +\left.{\displaystyle {\int }_{{\theta }_{\mathrm{C}}}^{{\theta }_{\mathrm{D}}}{\displaystyle {\int }_{-R}^{d_2}\frac{R{\left(r_{\mathrm{m}}+y\right)}^2}{\sqrt{R^2-y^2}}\mathrm{d}y\mathrm{d}\theta }}\right]\hfill \end{array}$$

(18)

where $d_1$ and $d_2$ can be ascertained through geometric relationships:

$$d_1=r_{\mathrm{A}}\frac{\sin {\theta }_{\mathrm{A}}}{\sin \theta }-r_{\mathrm{m}}$$

(19)

$$d_2=r_{\mathrm{C}}\frac{\cos {\theta }_{\mathrm{C}}}{\cos \theta }-r_{\mathrm{m}}$$

(20)

Similarly, the energy dissipation brought about by apparent cohesion $c_{app}$ is as follows:

$$\begin{array}{cc}\hfill {\dot{D}}_{\mathrm{s}}& =2\omega \left[{\displaystyle {\int }_{{\theta }_{\mathrm{A}}}^{{\theta }_{\mathrm{B}}}{\displaystyle {\int }_{d_1}^R{c}_{app}\frac{R{\left(r_{\mathrm{m}}+y\right)}^2}{\sqrt{R^2-y^2}}\mathrm{d}y\mathrm{d}\theta }}\right.+{\displaystyle {\int }_{{\theta }_{\mathrm{B}}}^{{\theta }_{\mathrm{C}}}{\displaystyle {\int }_{-R}^R{c}_{app}\frac{R{\left(r_{\mathrm{m}}+y\right)}^2}{\sqrt{R^2-y^2}}\mathrm{d}y\mathrm{d}\theta }}\hfill \\ & +\left.{\displaystyle {\int }_{{\theta }_{\mathrm{C}}}^{{\theta }_{\mathrm{D}}}{\displaystyle {\int }_{-R}^{d_2}{c}_{app}\frac{R{\left(r_{\mathrm{m}}+y\right)}^2}{\sqrt{R^2-y^2}}\mathrm{d}y\mathrm{d}\theta }}\right]\hfill \end{array}$$

(21)

In which the apparent cohesion $c_{app}$ induced via matric suction can be formulated as follows:

$$c_{app}=(u_a-u_w)\tan {\phi }_b$$

(22)

$$c_{app}=(u_a-u_w)\tan {\phi }_b=\rho z\tan {\phi }_b$$

(23)

where $z$ represents the vertical distance from the groundwater level to a point on the tunnel failure face, and it can be expressed as follows:

$$z=r_{\mathrm{A}}\cos {\theta }_{\mathrm{A}}-\left(r_{\mathrm{m}}+y\right)\cos \theta +z_0$$

(24)

3.3. Implicit Solution for the FS

The computation of the safety factor serves as an efficacious method for assessing the stability of tunnels excavated without applied support pressure. In this paper, the SRM is employed to determine the FS, a method that has found extensive application in M-C materials [48,49]. Within the framework of the SRM, the parameters, namely cohesion and the angle of internal friction, are progressively reduced until failure occurs. The FS is defined as the ratio of the current strength of the soil to the soil strength at a reduced state that leads to ultimate failure. Consequently, the FS can be written the following way:

$$FS=\frac{\tan {\phi }^{\prime }}{\tan {\phi }_{\mathrm{d}}^{\prime }}=\frac{c^{\prime }}{c_{\mathrm{d}}^{\prime }}$$

(25)

As the soil strength, characterized by cohesion $c_{\mathrm{d}}^{\prime }$ and the friction angle ${\phi }_{\mathrm{d}}^{\prime }$, is progressively diminished, the tunnel face approaches the ultimate failure state. Please be aware that the computation of the apparent cohesion $c_{app}$ induced via matric suction is interconnected with the internal friction angle and should not be further reduced subsequent to the reduction in the internal friction angle.

The SRM is utilized to formulate an implicit equation [50]. In this paper, the FS is derived through a binary search approach. The flowchart depicted in Figure 4 elucidates the principles of this method. A parameter without dimensions has been included to facilitate computation, as outlined below:

$$\Delta =\max \left\{\frac{1}{\omega \gamma r_{\mathrm{A}}^4}\left({\dot{W}}_{\gamma }-{\dot{D}}_{\mathrm{c}}-{\dot{D}}_{\mathrm{t}}\right)\right\}$$

(26)

For a specific FS, the variables ${\theta }_A$ and ${\theta }_B$ are utilized in an optimization procedure to arrive at the highest value $\Delta$ for Equation (26). Initially, $F{S}_1$ and $F{S}_2$, two initial amounts for the FS, are given, where $F{S}_1$ corresponds to the upper limit of the FS for a stable state, while $F{S}_2$ corresponds to the lower limit of the FS for an unstable situation. And the FS will be computed as the arithmetic average for $F{S}_1$ and $F{S}_2$. When $\Delta <0$, signifying the stability of the tunnel face, the mean FS is employed instead of the lower limit $F{S}_1$. Conversely, when $\Delta >0$, indicating the instability of the tunnel face, the upper limit $F{S}_2$ is substituted with the FS. This iterative process continues until either the absolute magnitude of $\Delta$ falls below 10⁻⁶ or the disparity between $F{S}_1$ and $F{S}_2$ narrows to a level below 10⁻³. MATLAB (R2019b) Software is employed to implement the entire optimization procedure. The design idea of the procedure is formulated on the basis of the enumeration method, which means that all possible results are calculated one by one and the maximum result is extracted from them. And the enumeration precision and range of the corresponding variables are as follows: precision 0.5°, range of 2°~86° for ${\theta }_{\mathrm{A}}$; precision 0.5°, range of 4°~88° for ${\theta }_{\mathrm{B}}$. Although this method seems to be more computationally intensive, it ensures the accuracy and repeatability of the data.

4. Results and Discussion

4.1. Comparisons

A numerical model based on FLAC 3D is developed to validate the correctness and effectiveness of the method proposed in this study. Due to the symmetry of the tunnel section, only half of the model is considered, as shown in Figure 5. The model boundary conditions are set at the bottom and both sides of the tunnel. The formula for apparent cohesion is selected from the results of Fredlund et al. (1996), as presented in Table 1. The unsaturated soil properties are chosen from SWCC NO. 3 in Table 2, and deformation characteristics are set as follows: Poisson’s ratio $\mathrm{v}=0.35$ and elastic modulus $E=1.0 GPa$. The remaining parameters are consistent with those of the KLA method employed in this study.

Two types of suction distribution are considered: uniform distribution and linearly increasing distribution. The comparison results are displayed in Figure 6. It can be observed that when suction is uniformly distributed ($u_a-u_w=50 \mathrm{kPa}$), the results from this method closely match those of the numerical simulation. In the case of linearly increasing suction ($u_a-u_w=\rho z, \rho =5$), the solutions from this method are slightly larger than the numerical simulation results, but the maximum difference does not exceed 9%. This discrepancy is attributed to the need for grid cell division in numerical simulation, leading to discontinuities in soil property parameters. Furthermore, the displacement contour in Figure 5 demonstrates that the 3D horned failure mechanism employed in this study aligns well with the simulated shape, providing additional evidence of the method’s correctness.

4.2. Parameters Analysis

Through parameter analysis, this study investigates the impact of the nonlinear shear strength, geometric attributes of the tunnel, and matric suction on the safety of the tunnel cross-sections. Here are some of the parameter configurations: $\gamma =20 \mathrm{kN}/{\mathrm{m}}^3$, ${\gamma }_{\mathrm{w}}=10 \mathrm{kN}/{\mathrm{m}}^3$, $c^{\prime }/\gamma D=0.1$ and $z_0=0$.

4.2.1. Case 1—Uniformly Distributed Suction

Figure 7 depicts the FS of the tunnel face, considering uniform matric suction, where ${\phi }_b$ represents a constant. Figure 7a delineates the correlation between the FS and the burial depth ratio C/D for scenarios with D = 10 m, $u_a-u_w=$ 0, 25, and 50 kPa. Figure 7b and Figure 7c, respectively, present the FS for different values of C/D and D. Figure 7d demonstrates the influence of the internal friction angle of the soil on the FS at varying matric suction levels.

As depicted in Figure 7, it is evident that matric suction has a favorable impact on tunnel face safety. Figure 7a,b indicate that as C/D increases, the FS initially improves before reaching a plateau. Once the burial depth ratio C/D surpasses a specific threshold, the horned failure mechanism will become completely submerged underground, and the stability of the tunnel face will no longer be influenced by variations in C/D. In contrast to the scenario where suction is not taken into account, the maximum difference, following the consideration of uniform suction of 50 kPa, can reach up to 33.4%. Figure 7c,d indicate that, under the condition of C/D = 1, as the suction and pressure increase, the FS exhibits nearly linear growth. Additionally, smaller values of D prove advantageous for the stability of the tunnel face, and changes in D have a relatively limited effect on the influence of the FS.

4.2.2. Case 2—Linearly Increasing Suction

This section presents the safety factor of the tunnel face under the condition of linearly increasing matric suction. Figure 8a,b illustrate the trends in how the FS is affected by C/D and D across various linear increase coefficients. Figure 8c,d delineate the impacts of variations in $\rho$ and $\phi$ on the FS.

It is evident that an increase in $\rho$ can significantly enhance the stability of the tunnel cross-sections. As illustrated in Figure 8a, the FS increases as C/D rises up until C/D surpasses a certain threshold. Under conditions of smaller $\rho$, the critical C/D value at which the FS transitions to a stable state is relatively low. From Figure 8b, it is evident that variations in D do not affect the FS. Referring to Figure 8c, it can be observed that when C/D is small, the FS linearly increases with $\rho$, whereas when C/D is large, the pattern of the FS variation with respect to $\rho$ begins to exhibit nonlinearity. This suggests that when C/D is large, the impact of p on the FS becomes more pronounced. In Figure 8d, an increase in $\phi$ can bolster the FS until $\rho$ and $\phi$ attain specific larger values.

4.2.3. Case 3—Nonlinear Shear Strength

The impact of soil strength nonlinearity on tunnel cross-section stability, based on the nonlinear strength equations from Table 1 and various soil characteristics from Table 2, is depicted in Figure 9, in which Figure 9a–d represent the nonlinear shear strengths of sandy soil, silt, clay, and ultra-fine grained soil, respectively. In this section, specific parameters were preselected: $\phi =34°$, D = 5 m, and C/D = 1.

As depicted in Figure 9, varying assessing equations for nonlinear strength and soil kinds result in various impacts for matric suction in soil on the FS of tunnel face, as shown in Figure 9. In the case of sandy soil, the curves for Vanapalli et al. (1996) [35] and Bao et al. (1998) [38] exhibit a flat and overlapping trend, indicating that these two models can be considered for saturated conditions. Regarding the other three models, the pattern of the FS varies nonlinearly within the range of $u_a-u_w=$ 40 kPa, and subsequently, it linearly increases with the augmentation of matric suction. In comparison to the case where the suction angle is constant (${\phi }^b=15°$), under the nonlinear equations, suction has a relatively minor effect on the stability of sandy soil tunnel faces.

For fine-grained soils, the curves of the Bao et al. (1998) [38] and Vanapalli et al. (1996) [35] models exhibit a similar pattern, both increasing nonlinearly at low suction, before gradually decreasing and ultimately degenerating into a horizontal and overlapping configuration. The curves for Fredlund (1996) [36] and Vilar (2006) [39] indeed exhibit a close resemblance, both eventually converging to stability. Additionally, the curve for Khalili and Khabbaz (1998) [37] is notably higher than the other nonlinear curves. When compared to sandy soil, the nonlinearity of fine-grained soil has a more pronounced effect on the stability of the tunnel face.

In the case of clay, when matric suction remains within 40 kPa, the trends of variation for all nonlinear curves, except for the Vilar model, are largely consistent with ${\phi }^b={\phi }^{\prime }=34°$. Subsequently, the growth trends of the curves become more gradual. Only the curve for Vilar et al. (2006) [39] falls below ${\phi }^b=15°$ and ultimately levels off. For ultra-fine grained soil, apart from the Vilar model, the safety factor curves are essentially consistent with the safety factor curve at ${\phi }^b={\phi }^{\prime }=34°$, and the nonlinearity can be considered negligible. This aligns with the findings of Zhang et al. (2014) [46] ${\phi }^b={\phi }^{\prime }=34°$ regarding larger AEV values. These curves are all significantly higher than the curve at ${\phi }^b=15°$. From the comparison in Figure 9, it can be concluded that matric suction has the most significant impact on ultra-fine grained soil, followed by clay and silt, while its effect on sandy soil is the least pronounced.

5. Conclusions

The novelty of this study lies in proposing an effective method for assessing the stability of 3D tunnel faces in unsaturated soils. The study investigated the differences in the stability of 3D tunnel faces when using linear and nonlinear strength models through the application of the kinematically limited analysis (KLA) and strength reduction method (SRM). This investigation encompassed diverse tunnel geometric factors, distinct soil kinds, and varying nonlinear strength models and groundwater suctions. The impact of matric suction on soil strength was modeled through the extended M-C failure criterion. By incorporating the horned failure mechanism, the external work rate and internal dissipation work rate were computed, and an implicit solution for the FS was derived using the SRM. Subsequently, a comparison with the numerical results and a parameter analysis were performed. Based on the aforementioned findings and discussions, the following conclusions can be drawn:

(1)

In the case of a linearly increasing matric suction distribution, the suction stress distribution initially increases linearly with the distance from the groundwater level and then gradually becomes less steep, transitioning to a nonlinear pattern. The suction stress curve derived from the nonlinear estimation equation for extremely fine-grained soil will degrade into a linear form. The suction stress distribution in fine-grained soil is more complex than in other soils, exhibiting a decreasing trend in suction stress when the distance from the groundwater level is sufficiently large.

(2)

Considering the influence of suction stress is advantageous for enhancing the stability of tunnel faces. This impact varies across different soils, with the least effect on sandy soil, followed by fine-grained soil (silt), and the greatest impact on extremely fine-grained soil. When matric suction reaches a sufficiently high level, the increase in suction stress may slow down or even decrease, as observed in fine-grained soils, leading to the stabilization of or a reduction in the stability of tunnel faces. For extremely fine-grained soil, the contribution of suction stress to the stability of tunnel faces is nearly a linear function of ${\phi }^b={\phi }^{\prime }$.

(3)

Compared to the uniform matric suction distribution condition, the stability of tunnel faces under the condition of linearly increasing matric suction distribution is more influenced by the burial depth ratio C/D. For sandy and fine-grained soils, the stability of tunnel faces is less affected by C/D, and the difference is not significant under saturated conditions. For extremely fine-grained soil and clay, the factor of safety increases linearly with C/D until it stabilizes. The stability of the tunnel faces is significantly improved compared with that under saturated conditions. These two types of soil can significantly reduce the economic costs of practical tunnel construction.