Probabilistic Stability Analysis of Deep Rock Tunnel Excavated by Mechanized Tunneling Considering Anisotropic Initial Stresses

Abstract

Due to the complex geological conditions, the construction of underground tunnels in the rocks exhibiting creep behaviors is a great challenge and the structures’ long-term stability is of main importance. The excavation by tunnel boring machines (TBM) of a tunnel that is supported by a double flexible/concrete liner is considered a potential technical solution. However, the efficiency of such a support system on the stability at long-term of a tunnel must be verified in the real condition by considering both the effects of uncertainties of the time-dependent rock behavior and anisotropy of initial stress state. In this paper, the subset simulation and the Sobol global sensitivity analysis are chosen to estimate the failure probability of the tunnel supports and quantify the importance of different random parameters, respectively. The results reveal that considering the anisotropy of initial stresses increases the failure probability, especially in the concrete support elements. In addition, the parameters of initial stresses are the ones of main importance for both liners according to the Sobol indices. Therefore, the anisotropy of initial stresses and the related uncertainties should be considered for reliable tunnel designs.

1. Introduction

Advanced techniques have been largely developed to ensure the stability at the short- and long-terms of tunnels. For example, the innovative support system using the double flexible/concrete liner of tunnel drilled by TBMs in the complex geological conditions has been proposed in tunneling projects [1]. However, the complexity of the geological formations’ behavior combined with their high uncertainty remains a great challenge on the stability prediction of tunnel that is usually excavated in the anisotropic condition of initial stress state [2,3,4,5,6].

In the last two decades, the probabilistic approach has been largely accepted as the most performant method to account for explicitly the uncertainty effect instead of using different factors of safety as in the traditional method. This stochastic approach allows quantifying the uncertainty propagation from input parameters (e.g., geological formations’ behaviors) to a structure’s responses (e.g., maximum stress), or estimating the structure’s failure probability (Pf) with respect to one or several failure modes. Some well-known methods, such as the Monte Carlo simulation (MCS), the first/second order reliability method (FORM/SORM), or the response surface method (SRM) have become standard and have been widely used in different tunnel projects [3,4,7,8].

So far, the reliability analysis of tunnels was mainly performed in the geological formations whose behaviors are assumed to be time-independent and characterized by an elastoplastic law. For example, ref. [4] proposed an approach to assess the system reliability of a tunnel verifying three failure modes (i.e., support capacity, tunnel convergence, and rock bolt length) in which the uncertainty of the Hoek–Brown parameters was accounted for. In many other contributions, scholars considered the uncertainty of the elastic perfectly plastic soils/rocks using the Mohr Coulomb model to analyse the stability of underground structures [3,8] and of the tunnel face [7]. The investigations conducted in these studies elucidated the important effect of the elastic moduli and strength of geological formations on the failure probability of constructed undergrounds. However, the adopted assumption of time-independent behavior can underestimate the obtained results, especially in the tunnel projects excavated in the host rock exhibiting a significant creep behavior. It is worth to note that in some cases, the time-dependent behavior of rock can contribute for more than 70% of the total convergence whilst the evolution in time of convergence (i.e., the well-known time delay effect) can significantly affect the stability of tunnel supports [9,10,11].

Thus, the stability of tunnels must also be evaluated at the long-term (i.e., at the end of their service life of about 100 years), especially for deeply buried tunnels constructed in the rock masses owing a high time-dependent behavior. This issue was extensively considered in the literature, but the uncertainty of creep behavior of the geological formations is usually neglected. For example, adopting some simple rheological models to describe the creep phenomenon of viscoelastic rocks, different scholars derived closed-form solutions to estimate displacements in rock masses and stress states in tunnel supports [12,13,14,15,16]. These solutions provided an interesting tool for quick analyses and parametric studies of tunnel responses at the long-term, which are usually expensive and time-consuming in numerical simulations. From their studies, these authors demonstrated the significant effects of creep rocks when an evolution in time of convergence on the tunnel surface induces a quick increase of compressive stresses in liners. Thus, failure of tunnel supports can occur during their service life as a consequence of the time delay effect when the stress state exceeds the constituent materials strength of each liner, even if these liners were well designed to verify the stability condition at the short-term (e.g., at the first months/years after the excavation).

The consideration of uncertainties of time-dependent behaviors for rock masses on the stability of deep tunnels was firstly conducted in the recent contribution of [17]. In their work, the time-dependent behavior of rock mass is characterized by the viscoelastic Burger model represented by four parameters: two springs and two dashpots of the connected Maxwell and Kelvin elements. The reliability analysis of a circular tunnel excavated in this kind of creep rock was undertaken using the direct MCS. More precisely, under the hydrostatic initial stress, the probabilistic response of the tunnel was assessed by MCS thanks to a developed closed-form solution. Their vast numerical investigations showed the important effects of Burger’s parameters on the evolution in time of stress states and hence the stability of liners at long term of 100 years.

In this work, the probabilistic analysis of a deep tunnel excavated in viscoelastic Burger rocks is reconsidered by accounting for the uncertainty of the initial stresses at far field and particularly of their anisotropy. For the sake of clarity, it is supposed that the position at depth H of deep tunnel is largely higher than the diameter of tunnel D (i.e., H/D > 10) that the variation of initial stress state around tunnel can be neglected. In addition, the studied deep tunnel is supported by a system of double elastic liners: an outer flexible liner and an inner concrete liner; both are installed at the same time. Owing a low elastic modulus, the outer flexible liner can absorb the high evolution in time of the rock convergence and reduce the transmitted loading to the concrete liner. Furthermore, a firmly bonded assumption of the outer liner with the rock mass and the inner concrete layer is adopted.

In comparison with the case of isotropic initial stresses as considered in [17], the behavior of the tunnel supports is more critical when the moment induced by anisotropic initial stress increases the concentration of compressive stress in liners as shown in the deterministic analyses performed by different authors [2,5,6,18,19]. From a statistical point of view, it is more rational to consider the anisotropy of the initial stresses and the related uncertainties compared to the assumption of isotropic initial stresses using constant parameters.

The present study is among the first attempts to consider the uncertainty and anisotropy of the initial stress for the probabilistic assessment of a deep tunnel drilled in Burger rocks and supported by double liners. Additionally, a probabilistic analysis procedure for the above-mentioned problem is proposed in the paper. It combines an effective tunnel stability analysis tool with the subset simulation and the Sobol method. The procedure allows evaluating the reliability of a deep tunnel from different aspects at a first design stage, as different results can be provided by the procedure within a short computational time.

2. Problem Statement and Chosen Methodologies

2.1. Problem Statement

The considered problem consists of analysing the stability at the long-term (100 years) of a circular tunnel accounting for the uncertainty of both the creep rock behaviors and the anisotropic initial stresses. The time-dependent behavior of the viscoelastic Burger rock is characterized by two springs ($G_M, G_K$) and two dashpots (${\eta }_M, {\eta }_K$) of the connected Maxwell and Kelvin elements. The tunnel is supported by a system of double elastic liners which consists of an outer flexible liner and an inner concrete liner. The outer one is assumed to be firmly bonded with the rock mass and with the inner liner; both liners are installed at the same time $t_0$. Figure 1 presents the considered tunnel constructed in the viscoelastic Burger rock with double supports.

Whilst the deterministic analysis of deep tunnels constructed in viscoelastic materials has been largely performed [5,6,18,19], the time-dependent behavior’s uncertainty of rock masses has not yet received a lot of attention. Recently, Do et al. [17] investigated the stability of the tunnel supports using reliability analyses. Their stochastic analyses by MCS demonstrated the significant effects of Burger parameters’ uncertainty on the stress state in liners and hence their stability at the long-term. More precisely, the $Pf$ of these elements was estimated at 100 years using the following limit state function (LSF):

$${\mathrm{LSF}}_1={\sigma }_{e\_max}- {\sigma }_e\left(\mathit{X}\right)$$

(1)

where ${\sigma }_{e\_max}$ is the allowable value of the equivalent stress ${\sigma }_e\left(\mathit{X}\right)$ in the support element; this latter depends on the random input parameters gathered in the vector $\mathit{X}$. An analytical solution of this equivalent stress in each support element with respect to $\mathit{X}$ was presented in the work of [17]. However, this last contribution was limited to the case of isotropic stress state at far field and only the uncertainty of the rheological rock was considered (i.e., the vector of random variables consists of four parameters $\mathit{X}=\left[G_M, G_k, {\eta }_M, {\eta }_K\right]$ in [17]).

In comparison with [17], the uncertainty of the anisotropic stress state at far field is involved in the present study. Noting $\lambda$ and $P_0$ as the anisotropic ratio and the vertical principal value of stress state at far field, the vector of random variables is now characterized by six parameters $\mathit{X}=\left[G_M, G_k, {\eta }_M, {\eta }_K, \lambda , P_0\right]$. Furthermore, a volumetric strain based LSF for the outer flexible liner is proposed as follows:

$${\mathrm{LSF}}_2={ϵ}_{v\_max}- {ϵ}_v\left(\mathit{X}\right)$$

(2)

According to Equation (2), the elastic behavior of the outer flexible liner is limited to an allowable value ${ϵ}_{v\_max}$ of the volumetric strain before attaining the non-linear behavior stage due to the pore collapse mechanism.

The calculation of the equivalent stress ${\sigma }_e\left(\mathit{X}\right)$ in the concrete liner and the volumetric strain ${ϵ}_v\left(\mathit{X}\right)$ in the flexible liner is based on a semi-analytical solution presented in the work of [14]. To derive this solution of a circular tunnel under an anisotropic stress state of a viscoelastic Burger rock, Song et al. [14] combined the Laplace transform technique and the complex potential theory. More precisely, the authors expressed the complex potentials in the form of the truncated Laurent series whose coefficients are numerically determined from the inverse Laplace transform. For details of this solution, the interested readers can refer to [14].

2.2. Probabilistic Methods

This section briefly presents the methods selected in this work to probabilistically assess the considered tunnel, including the subset simulation (SS) for the reliability analysis and the Sobol method for the global sensitivity analysis (GSA).

2.2.1. Subset Simulation

The SS, a variant of MCS, improves the efficiency of MCS in estimating low $Pf$ values by solving a series of simpler reliability problems [20]. In the context of SS, the original event $E_{SS}$ is firstly decomposed into a sequence of intermediate ones $\left[E_{SS\_1}, E_{SS\_2},\dots , E_{SS\_m}\right]$ and the target $Pf$ can be expressed as [21]:

$$Pf=P\left(E_{SS}\right)=P\left(E_{SS\_1}\right){\displaystyle \prod }_{i=2}^{m}P(E_{SS\_i}|E_{SS\_i-1})$$

(3)

where $P(E_{SS\_i}|E_{SS\_i-1})$ is the probability of the event $i$ conditioned on the event $i-1$ and $P\left(E_{SS\_1}\right)$ is the probability of the first event. The value of $m$, representing how many intermediate events are needed, is determined by gradually generating conditional samples until reaching the failure domain.

In this paper, the conditional samples of each event are generated using the modified Metropolis–Hastings Markov Chain Monte Carlo sampling [22]. The threshold and the sample number of each event is set, respectively, as 0.1 and 5000 following the recommendations given in [22,23]. The SS has been widely used in the literature to provide quick and reliable $Pf$ estimates of geotechnical structures [24,25,26].

2.2.2. Sobol Global Sensitivity Analysis

The global sensitivity analysis (GSA) aims to quantify the sensitivity of a quantity of interest (QoI) with respect to each random variable (RV) over its entire varying range. Among many methods, the Sobol indices have received much attention since they can give accurate results for most models [27], and are adopted in the present work.

The Sobol GSA is based on the variance decomposition of a random model output ($Y$). The first-order Sobol index of a variable $x_i$ can be given as:

$$S_i=\frac{Var\left[\mathbb{E}\left(Y\parallel x_i\right)\right]}{V_t}$$

(4)

where $V_t$ is the total variance of the model output. For the $Var\left[\mathbb{E}\left(\mathrm{Y}\parallel x_i\right)\right]$, the inner expectation operator $\mathbb{E}\left(\cdot \right)$ is the mean of $Y$ considering all possible $x_{~i}$ values while keeping $x_i$ constant; the outer variance $Var\left(\cdot \right)$ is taken over all possible values of $x_i$. The $S_i$ is a measure about the effect of the variable $x_i$ alone.

There is also the total Sobol index for a variable $x_i$, which represents the sum of all the Sobol indices involving this variable. The definition can be written as:

$$S_{Ti}=S_i+{\displaystyle \sum }_{i\ne j}{S}_{ij}+\cdots S_{1,\dots ,M}$$

(5)

where $M$ is the number of all the random variables.

Estimation of these Sobol indices (first-order and total order) can be achieved by constructing a triplet of matrices [28]; each matrix has $M$ columns and $N_s$ rows. In this work, the $N_s$ is set as 10,000 following the suggestion of [29]. Application of the Sobol GSA to geotechnical engineering can be found in [30,31,32].

3. Numerical Results

This section defines firstly a reference case by giving the geometry and mechanical parameters of the tunnel to be assessed. Then, the numerical results of the tunnel liners’ safety evaluation at 100 years using deterministic, reliability, and sensitivity analyses are presented.

3.1. Definition of the Reference Case

The reference case is described in

Table 1. Input parameters for the analytical model.

	Parameter	Symbol	Value
Initial stress	Vertical stress	$P_0$	6.8 MPa
	Anisotropic coefficient	$\lambda$	0.66
Burger model	Viscosity of Kelvin	${\eta }_K$	2. 07 × 107 GPa.s
	Viscosity of Maxwell	${\eta }_M$	4. 14 × 109 GPa.s
	Shear modulus of Kelvin	$G_K$	0.34 GPa
	Shear modulus of Maxwell	$G_M$	3.45 GPa
Geometry	Tunnel radius	$R$	4.5
	Thickness of 1st liner	$l_1$	0.2
	Thickness of 2nd liner	$l_2$	0.45
Two liners	Elastic modulus of 1st liner	$E_1$	0.1 GPa
	Elastic modulus of 2nd liner	$E_2$	30 GPa
	Poisson coefficient of 1st liner	${\nu }_1$	0.001
	Poisson coefficient of 2nd liner	${\nu }_2$	0.2
Others	Installation time of the two liners	$t_0$	0.001 day
	Excavation rate	$v_l$	0.75 m/day
	1st parameter of deconfinement rate function	$m_1$	0.7
	2nd parameter of deconfinement rate function	$m_2$	1

which includes all the necessary input parameters. They are quite similar as the ones in [17] except for the anisotropic stress state which was not considered. For a better understanding of the Burger model, excavation and convergence parameters, and their value selection, readers are referred to [12,13,17]. A quasi-instantaneous installation of both liners is assumed by profiting the advantage of the presence of the outer flexible liner. In fact, each liner stress state in the tunnel support system depends strongly on their rigidity as well as their installation time. The discussions in [17] highlighted that an early installation of the liners induces the higher compressive stress, particularly in the inner concrete liner which is also amplified by the rigidity (i.e., higher elastic moduli) of the outer elastic liner. However, in the case of flexible liner, the low elastic moduli of this support allow it to absorb the increase in time of the rock convergence and hence reduce the transmitted load to the concrete liner. This advantage allows to install at the same time two liners that facilitates the excavation and installation of a tunnel support system by TBM.

For the reliability analysis, two failure modes will be considered for analyzing the tunnel support system: the exceeded volumetric strain in the outer flexible liner (Equation (2)) and exceeded equivalent stress in the inner concrete liner (Equation (1)).

3.2. Deterministic Analyses

The two QoIs (i.e., the volumetric strain in the outer liner and the equivalent stress in the concrete layer) are firstly studied in a deterministic framework with the parameters of Table 1.

Although in most cases, the initial stress state in the horizontal direction is smaller than the vertical stress (i.e., $\lambda <1$), it is also stated the higher horizontal stress in some cases [33]. Thus, in order to better investigate the anisotropic effect of initial far field stress, two other values of the anisotropic ratio ($\lambda$= 1 and $\lambda$= 1.2) are considered and their results are compared with the reference case. It is noted that $\lambda$= 1 refers to isotropic initial stresses. As illustrated in Figure 2, the QoIs, at the intrados of the corresponding liner, are calculated as a function of the inclined angle ranging from 0° (springline) to 90° (crown). As expected, for the case of $\lambda$ = 1, one observes uniform stresses and deformations from the springline to the crown in each liner. An increase of the anisotropic degree of initial stresses (case of $\lambda$ = 1.2) induces a higher compressive orthoradial stress at the tunnel roof and hence the equivalent stress in the concrete liner is maximum at this position. Inversely, the volumetric strain in the flexible liner is highest at the springline. However, at the same position, this deformation decreases quickly and becomes lowest for $\lambda$ = 0.66. The same trend can be stated for the equivalent stress at the concrete liner roof. However, for both cases of anisotropic stress states, the maximum equivalent stress in the principal support element of tunnel (i.e., concrete liner) is higher than the one of the isotropic stresses. This expected observation confirms the discussions in the literature [5,6,18,19]. Following that, the anisotropic initial stress generates a non-uniform deformation distribution in the outer liner and of the transmitted load/pressure at the extrados of the concrete liner. This non-uniform pressure induces bending moments in the concrete liner and results in, as a consequence, higher compressive hoop stresses in this last support element.

An investigation with different values of anisotropic ratio in the range $\lambda \in \left[0.66, 1.2\right]$ as depicted in Figure 3 is then performed. Only the results at the springline and crown are captured in this figure, which confirms the above-mentioned observations. Regarding the maximum equivalent stress in the concrete liner (Figure 3b), these results demonstrate the important role of the anisotropy of initial stress state on the stability of this support element. The maximum equivalent stress in this tunnel support element is the lowest in the isotropic stress case.

3.3. Probabilistic Analysis of the Reference Case

This section evaluates probabilistically the tunnel liners by following the procedure proposed in this paper. A flow chart of the procedure is given in Figure 4. It starts by preparing a computational model and defining the input uncertainties. The analytical model introduced in [17] is suggested for a preliminary study of deep tunnels due to its high computational efficiency which is an attractive feature for probabilistic analyses. It is followed by estimating the $Pf$ of the double liners by using the SS which allows having a first idea of the tunnel’s reliability. Then, the samples generated in the first subset of the performed SS are further explored to quantify the variation such as a confidence interval of the QoIs (i.e., volumetric strain and equivalent stress). Lastly, the importance of each stochastic parameter is quantified via a Sobol sensitivity analysis so that all the parameters can be ranked according to their importance. The advantage of the procedure lies in the fact that it can provide many valuable results for understanding and evaluating the tunnel reliability which is favorable for a rational decision.

The parameters of the viscoelastic rock behaviors and the initial stresses are modelled by Lognormal-distributed random variables as shown in

Table 2. Uncertainty modelling of the Burger and initial-stress parameters.

Parameter	Symbol	Unit	Distribution	Mean	CoV
Vertical stress	$P_0$	MPa	Lognormal	6.8	15%
Anisotropic coefficient	$\lambda$	/	Lognormal	0.66	15%
Viscosity of Kelvin	${\eta }_K$	GPa.s	Lognormal	2.07 × 107	15%
Viscosity of Maxwell	${\eta }_M$	GPa.s	Lognormal	4.14 × 109	15%
Shear modulus of Kelvin	$G_K$	GPa	Lognormal	0.35	15%
Shear modulus of Maxwell	$G_M$	GPa	Lognormal	3.45	15%

. Their mean values are taken from Table 1 and a coefficient of variation (CoV = 15%) is assumed due to the lack of relevant information in the literature. Other parameters in Table 1 are kept constant. Although the anisotropic ratio $\lambda >1$ can be found in some cases as mentioned above, in the following reliability analysis, it is supposed that this parameter $\lambda$ is bounded at the upper value 1 to be convenient with the usually observed cases.

3.3.1. $Pf$ Estimation

The $Pf$ estimation by SS in each support element is conducted using the corresponding LSF defined in Equations (1) and (2). The values of ${ϵ}_{v\_max}=45\%$ and ${\sigma }_{e\_max}=45\mathrm{MPa}$ are chosen correspondingly as the threshold of the volumetric strain in the flexible liner and the maximum equivalent stress in the concrete liner.

The estimated $Pf$ values of the tunnel liners at different locations (every 15° from 0° to 90° with respect to the springline) are presented in Figure 5. Due to symmetry, only a quarter of the tunnel is presented. Consistent with the previous deterministic analyses in the case $\lambda <1$, the failure probability in the outer flexible liner increases from 0° to 90° with a maximum value of 3.2% found at the roof. Inversely, the $Pf$ in the concrete liner decreases from the springline to the crown with the highest value of 17.5%. Under the configuration of the selected thresholds, the adopted parameters in Table 1, as well as the uncertainty modelling, these results show that the more critical failure mode is in the concrete liner.

3.3.2. Variation of the Two QoIs

In the SS-based reliability analyses, a relatively large sample size (i.e., 5000) is considered for each subset, which allows building a distribution of the concerned QoI using the 5000 deterministic results of the first subset as its samples are generated over the entire input space.

As an illustration, Figure 6 presents the probability density function (PDF) of both QoIs at three positions in the support elements. For the volumetric deformation in the outer flexible liner, the distribution moves to higher values and tends to be more spread out from the springline (0°) to the crown (90°), represented by the higher mean value as well as the higher variation. Regarding the equivalent stress in the concrete liner, an inverse trend can be stated when the distribution spread out from the crown to the springline. The results can help to explain statistically the $Pf$ trends observed in Figure 5. Additionally, further exploiting the 5000 samples allows to obtain more statistics’ estimates of the QoIs such as the mean, standard deviation, and confidence bounds [34] as shown in

Table 3. Statistics of the volumetric strain and equivalent stress.

Angle	Volumetric Strain (%)				Equivalent Stress (MPa)
	Mean	Std (1)	Low (2)	Up (3)	Mean	Std	Low	Up
0	22.69	5.19	14.00	34.31	37.96	7.89	24.14	55.22
45	27.19	5.57	17.23	39.28	27.39	5.61	17.71	39.61
90	31.88	6.58	20.35	46.50	16.70	5.62	7.71	29.47

Notes: ⁽¹⁾ standard deviation; ⁽²⁾ lower or ⁽³⁾ upper bound of the 95% confidence interval.

.

3.3.3. Sensitivity Analysis

The Sobol GSA, as presented in Section 2.2.2, is carried out in this part to quantify the sensitivity of the QoI in each liner with respect to each parameter of the vector $\mathit{X}=\left[G_M, G_k, {\eta }_M, {\eta }_K, \lambda , P_0\right]$. For the sake of clarity, only the critical positions of the liners are considered: flexible liner at the crown and concrete liner at the springline. As observed in Figure 7, among the four parameters of the Burger model that characterize the creep behavior of the viscoelastic rock, the Kelvin’ spring element $G_K$ affects the most the stability of both liners. The influence of the Maxwell spring element $G_M$ is higher than the one of the two dashpots. These results are consistent with the observations provided from the parametric studies in [17]. However, to complete the previous contribution, the resulted obtained in this work highlighted that the initial stress $P_0$ with a highest value of Sobol index is the most affecting factor on the failure probability in each support element. With respect to the absolute value of initial stresses, the sensitivity analysis exhibits a moderate effect of anisotropic ratio $\lambda$ on the stability of tunnel supports.

3.4. λ Effects with Different Thicknesses of the Liners

To provide more insights about the effect of the anisotropic initial stresses (i.e., anisotropic ratio λ), the failure probabilities of the tunnel supports calculated with different cases are compared: isotropic stress with constant value λ = 1 (case 1), anisotropic stress with constant value $\lambda$= 0.66 (case 2), and anisotropic stress with random variable λ as shown in Table 2 (case 3). In Figure 8, these failure probabilities at the springline of the concrete liner and at the crown of the outer flexible liner are calculated versus different liners’ thickness. The results are consistent when the failure probability in each liner decreases with respect to its higher thickness. In comparison with the isotropic stress state, the anisotropic cases (i.e., cases 2 and 3) reveal the higher failure probability in the concrete liner for the whole range of considered thickness. Regarding the previous deterministic analysis, this result can be expected when the anisotropy of initial stresses induces in this concrete support element higher equivalent stress with respect to the case of hydrostatic stress state. Concerning the $Pf$ of the outer liner, as a function of the thickness, the trend is not clear. More precisely, the lower failure probability in the anisotropic cases at small thickness can become higher than the one of the isotropic cases when the thickness of this flexible support element overpass a value of about 25 cm. Between the two anisotropic cases, no significant difference of failure probability is observed for both liners. To reinforce this statement, in the last numerical investigations, $Pf$ was calculated as a function of the mean value and CoV of the anisotropic parameter $\lambda$ (Figure 9), keeping the same liners thickness and other random parameters as in Table 1 and Table 2. Whilst the lower mean value of $\lambda$ induces a significant increase of the $Pf$ concrete liner, the failure probability of the outer liner is slightly reduced (Figure 9a). The results also reveal the negligible effect of the uncertainty of $\lambda$ on the failure probability of tunnel supports (Figure 9b), especially in the range CoV ∈ [5–25%].

4. Conclusions

The long-term stability of a deep tunnel constructed in viscoelastic rocks is investigated by considering the uncertainty of rock time-dependent behaviors as well as the anisotropic initial stresses. Two support elements of the circular tunnel are supposed to be installed immediately by profiting the capacity of the outer flexible liner that reduces the transmitted loading to the inner concrete liner. Based on the closed-form solution presented in the literature, the subset simulation is conducted to estimate the failure probability at one hundred years of the volumetric deformation in the flexible liner and of the equivalent stress in the concrete liner. Comparing with the hydrostatic stress state, the anisotropy of initial stresses induces higher equivalent stresses and as consequence higher failure probability in the concrete liner. The Sobol’s global sensitivity analysis shows that the initial stress state uncertainty is the most influential factor while the uncertainty of the anisotropic stress state ratio can be considered as negligible. This study confirms the necessity to consider the initial stresses anisotropy and the related initial stress uncertainties for the long-term stability analysis of underground structures constructed in creep rocks. The proposed probabilistic analysis procedure could be useful for tunnel reliability evaluation in practice especially for a preliminary stage at it is able to provide a variety of informative results within a short time.