A simplified approach to reliability evaluation of deep rock tunnel deformation using First-Order Reliability Method and Monte Carlo simulations

: The Monte Carlo simulation (MCS) and First-Order Reliability Method (FORM) provide a reliability analysis in axisymmetric deep tunnels driven in elastoplastic rocks. The Convergence-Confinement method (CV-CF) and Mohr-Coulomb (M-C) criterion are used to model the mechanical interaction between the shotcrete lining and ground through deterministic parameters and random variables. Numerical models synchronize tunnel analytical models and reliability methods, whereas the limit state functions control the failure probability in both ground plastic zone and shotcrete lining. The results showed that a low dispersion of random variables affects the plastic zone's reliability analysis in unsupported tunnels. Moreover, the support pressure generates a significant reduction in the plastic zone's failure, whereas the increase of shotcrete thickness results in great reduction of the lining collapse probability.


INTRODUCTION
Several geotechnical and structural parameters are involved in deep tunnels and underground excavation analysis, which induce high or undefined structural risks. The reliability analysis evaluates structural failure probabilities in tunnels by applying random and deterministic parameters in analytical solutions or numerical methods.
The first serious discussions and analyses of tunnel reliability emerged during the 1990s with [1]- [6]. Low and Tang [5] focused on identifying and evaluating the reliability and risk parameters through numerical procedures, synchronizing reliability codes developed in software like Microsoft Excel, and analytical methods for axisymmetric tunnels. This work used the First Order Reliability Method (FORM) and Monte Carlo Simulation (MCS) in the reliability study. Low and Tang [7], [8] developed and improved the FORM numerical procedures. Li and Low [9] conducted a series of reliability analyses of an axisymmetric tunnel, considering normal and non-normal random variables of the rock mass and the Low and Tang [5] FORM algorithm. The internal pressure produces a reduction of risk; moreover, the authors checked the reliability results using MCS. Lü and Low [10] and Song et al. [11] studied the reliability of axisymmetric tunnels through the FORM and Second-Order Reliability Method (SORM), synchronizing these methods with analytical solutions based on Mohr-Coulomb (M-C) and Hoek Brown criteria.
Bjureland et al. [12] applied the MCS associated with the Convergence-Confinement method (CV-CF method) to analyze the failure probability in a shotcrete lining of a circular tunnel. The M-C criterion describes the plasticity behavior for the excavated rock mass, and the Observational method checks the failure results, verifying the necessity of putting into action security measures.
The present paper analyses the reliability methods impact on analytical tunneling methods, examining two different tunnels with statistical and deterministic distinct parameters. This study aims to investigate the reliability results of the plastic zone excavation and the shotcrete lining for both tunnels, establishing interesting comparisons between them. For this, the reliability procedures associate the use of the FORM Low and Tang method, developed in Low and Tang [5], and the MCS. The reliability algorithms, developed in VBA and MATLAB codes, employ two analytical methodologies for axisymmetric tunnels: the M-C plasticity criterion for excavated rock masses and the CV-CF method to analyze the interaction between the shotcrete lining and rock mass.

Reliability Index: The Low and Tang FORM Interpretation
Hasofer and Lind [13] defined the second-moment reliability index β as a better approach in civil engineering design in comparison with usual safety factors [7]. Equation 1, proposed by Low and Tang [7], represents the matrix formulation of the Hasofer-Lind reliability index for correlated normal variables: where i x is the set of random variables of the vector x, i µ is the set of mean values of the vector μ , and i σ is the set of standard deviations of the vector σ . The domain F represents the failure regions, whereas [ ] is the correlation matrix, expressed by [14]: where i, j ρ is the correlation coefficient between two random variables i x and j For the non-normal distributions, the basic idea is to calculate the equivalent normal mean N μ and equivalent standard deviation N σ , which replaces the values of μ and σ in Equation 1. The Rackwitz-Fiessler transformation defines the following analytical solution for N μ and N σ [15]: In a classical explanation, β is the smallest distance from the mean value point to the limit state surface in units of directional standard deviations [7]. Low and Tang [5], [7] proposed an alternative interpretation of the reliability index β , which is a simple method of defining the index, in the original space of the random variables, through the perspective of an expanding ellipsoid. Figure 1 shows this interpretation. As shown in Figure 1, the design point is a tangency point of the expanding dispersion ellipsoid with the limit state surface, which separates safe values from unsafe values (failure domain). The reliability index β is the axis ratio / R r of the ellipse that touches the limit state surface [8]. The dispersion ellipsoid expands from the mean-value point according to the following probability density function of a multivariate normal distribution [8]: where to minimize β , defined by Equation 1, is to maximize the value of the probability density function defined by Equation 5. Hence, it is possible to define the smallest dispersion ellipsoid tangent to the limit state surface (see Figure 1) and the most probable failure point (design point). The index β is particularly useful in evaluating the failure probability through the application of the CDF of β , ( ) Φ β , and considering standard normal variables [9]: Low and Tang [8] performed a new efficient algorithm using FORM, where the first and third terms under the square root in Equation 1 are equivalent standard vectors n     . Thus, Equation 1 is defined as: The optimization routine Solver of Microsoft Excel solves Equation 7. Autonomously, the routine changes the values of components i n (computing the parameters N μ and N σ , when necessary) and subjecting the limit state function to a zero-value restriction. The components in the original space, i x , are program-calculated from: ( )

Monte Carlo Simulation
The Monte Carlo Simulation (MCS) consists of developing an analytical and numerical model and using many simulations or cycles to check the real behavior of a system. Each cycle provides output variables (results) based on random input variables, using statistical methods and probability distributions [16].
Considering a random variable X, with a CDF ( ) i X F x , to generate sample i x values included in X, the MCS code needs to accomplish some steps. First, MATLAB numerical procedures generate a sample value of i u between 0 and 1 for a random variable with uniform distribution. Next, the inverse of X F computes a sample value of i x ( ( ) For two correlated normal random variables, the numerical procedure generates two random numbers 1 u and 2 u , used to compute the random variables x and y from the following equations [17]: 1 2 2 ln cos(2 ) where the mean and standard deviation of the variable y, which are dependent on x value, are evaluated through the following equations: Each MCS cycle consists of introducing random input variables in an appropriated mathematical model: the limit state or performance functions ( )  Figure 3 shows the axisymmetric tunnel characteristics, excavated in a homogeneous elastic-plastic rock mass and with a radius i R . The cylindrical coordinate system r e , e θ , and z e represents the stress and displacement components.  In the fully elastic regime, the radial displacement ( ) ie u r at any a point of the rock mass located at distance r from the center of the excavation is given by ( [18], [19]):

Elastoplastic fields in the surrounding rock mass
where ν and E are the Poisson's rate and Young's modulus, respectively. In the elastoplastic regime, the mechanical analysis of the tunnel deformation indicates that two regions should be distinguished, as illustrated in Figure   The analytical methods use the M-C criterion concepts for the elastoplastic behavior, defined according to Salençon [18] and Corbetta [19]. The plastic radius p R is: where: tan c a ϕ = (18) in which c is the cohesion and ϕ is the friction angle. The limit pressure lim P defines a rock mass pressure when p r R = as follows: Referring to Figure 4, the radial displacement ip u in the elastoplastic region i where Equation 21 defines the dilatancy coefficient b K , ψ representing the dilatancy angle, whereas Equations 22, 23 and 24 characterize the parameters A, B, and C.

Convergence-confinement method
The CV-CF method allows modeling in a simplified manner the evolution of normalized radial displacement i U of the rock mass and support, or the convergence at the tunnel wall, associated with the internal pressure i P in a specific cross-section of the tunnel wall. The convergence i U is given by: (elastoplastic behavior). The principle of the CV-CF method is schematized in Figure 5.  Figure 5 shows the Ground response curve, or Convergence curve (CV), which describes only the mechanical behavior of the unsupported excavation wall, while the Support curve, or Confinement curve (CF), describes the mechanical behavior of the elastic shotcrete lining.
When the excavation front is close to a specific cross-section, changes occur in the ground stress state. At some distance, before the excavation reaches the cross-section, there are small changes of the stress state, and elastic convergence ie U occurs as i P decreases. The internal pressure i P continues to decrease where, eventually, the decrement of stresses in the rock mass will reach a limit lim P , and the stresses cause plastic behavior of the ground in a zone with a radius p R surrounding the tunnel periphery [12]. Thus, considering an interval i lim P P P ∞ ≥ ≥ , the CV curve shows an elastic behavior (blue curve in Figure 5), and for the interval lim i 0 P P ≥ ≥ , the CV curve shows an elastoplastic behavior. Equations 16 and 19 defines p R and lim .

P
The CF curve avoids the ground confinement loss, which causes normalized radial displacements of the rock mass. According to the red curve in Figure 5, a ground confinement loss generates a growing internal pressure i P and convergence i U on the cross-section lining. Moreover, numerical or analytical approaches estimate the convergence 0 U at the beginning of the support installation. Panet and Guenot [20] proposed an analytical approximation for 0 U : where i, U ∞ is the convergence for elastoplastic behavior, at a larger distance between the excavation face and the crosssection (Equation 25 with i 0 P = ). The pressure at the shotcrete lining is evaluated as [21]: where s t is the shotcrete lining thickness, while s E and s ν are the Young modulus and the Poisson Coefficient of the support, respectively. The maximum support pressure max P is defined by the following Equation: By using Equation 16 one can calculate p R with i eq P P = . Then Equation 30 is solved numerically by the bisection method inside the MCS MATLAB code.

PERFORMANCE FUNCTIONS AND ADOPTED PARAMETERS
The limit state function, or performance function, ( ) The ratios L and max The geotechnical random variables were obtained in the research of Hoek [6] (tunnel A) and Brantmark [24] (tunnel B). Tables 1 and 2 show the geotechnical statistical parameters for tunnels A and B.  for tunnel B. It should be emphasized that the correlation between some of the design parameters significantly affects the reliability predictions. In that respect, the assumption of negatively correlated shear strength parameters was found to be conservative with respect to uncorrelated variables [25]. The analysis developed, for instance, in Laso et al. [4] has included the correlation between some of the relevant problem parameters, such as correlation between the ground Young modulus and ground density, shotcrete average resistance and shotcrete layer thickness, as well as between tunnel radius and shotcrete layer thickness. Although the subject of correlation is not fully addressed in the present reliability study, this important issue will be foreseen in the continuation line of this research.
The rock mass failure criterion features a non-associated flow rule in both tunnels, with the following deterministic parameters: for tunnels A and B, respectively. For both tunnels, a radius of i 4.5m R = is adopted. Table 3 shows the statistical parameters for the shotcrete lining. Both tunnels use the same parameters, which are statistically independent and normally distributed. The Poisson coefficient is c 0.25 ν = . Li and Low [9] suggested applying the lognormal or beta distribution instead of the normal distribution when the coefficient of variation a random variable (ratio between standard deviation and mean) is 0.25 or higher. This is necessary to avoid negative geotechnical parameters, which is irrational. All coefficients of variation are smaller than 0.25 in tunnels A and B, except for the cohesion c in tunnel A, where, alternatively, the cohesion is constrained to be greater than zero.

RELIABILITY ANALYSIS REGARDING THE EXTENT OF PLASTIC ZONE
This section analyzes the reliability of tunnels A and B, contemplating the performance function of Equation 31.
The first reliability analysis performs the characterization of the fixed resistances in ( ) 1 g x regarding unsupported tunnels, whereas the second analysis of the reliability methods will contemplate the shotcrete lining. The 2007 FORM algorithm of Low and Tang and MCS compute the reliability indexes and failure probabilities, employing the analytical theories explained in section 2.
In summary, the Low and Tang algorithm is a code developed in Visual Basic (VBA) language and Microsoft Excel. A constrained optimization routine, denominated Solver, assists in the FORM procedure development, where the equations and parameters described in sections 2.1 and 4 have been employed directly in the Excel worksheets. Low and Tang [7], [8], and Li and Low [9] give more details about the FORM algorithm of Low and Tang.
First, it is necessary to verify the results obtained by FORM and MCS to check the algorithms' proper operation. Li and Low [9] analyzed the plastic zone reliability, employing the FORM Low and Tang code, the performance function ( ) 1 g x , and tunnel A parameters. Thus, Table 4 shows a comparison between results of this work and Li and Low [9], with the following deterministic values: the parameters of tunnel A, 2.50MPa P ∞ =  The failure probabilities evaluated in Tab. 4 have shown a good agreement with the Li and Low [9] results, considering both FORM and MCS results of this work. Thus, considering random variables with normal distribution, the reliability algorithms used in this work have shown efficiency when applied in unsupported and supported tunnels analysis.
The first analysis intends to verify the influence of resistance L in the reliability index β , on the plastic zone reliability in tunnels A and B. Figure 6   Comparing the curves in Figure 6, a great sensibility on the β variation with the increase of L occurs for tunnel B (red curve). For example, the index in tunnel B increases from  Table 5 illustrates the values of β and f p according to the curves of Figure 6.  In contrast, the reliability index β increases from 0.14 to 0.65 for tunnel A and from 0.88 to 5.64 for tunnel B, with L variation. The greater variation of reliability results in tunnel B is a consequence of the correlation between the random variables and the plastic zone results p i / R R , performed by MCS. . Due to this, the reliability results of tunnel A, in Figure 6 and Table 5, shows a low variation when adopting 2. In an overall analysis, these results say that the designers must be careful when adopting elevated values for L, in reliability analysis that involves highly correlated parameters. For this purpose, this study adopted the average resistance values in the performance function ( )    The MCS code provides the plastic zone histograms in Figures 8a and 8b, where the MATLAB application Distribution Fitter adjusts probability distribution functions with these histograms. To summarize, the i P increasing causes the reliability index increasing and failure probability reduction during the reliability analyses.

RELIABILITY ANALYSIS REGARDING THE TUNNEL CONVERGENCE
In this section, the reliability analysis evaluates the probability of eq P and eq U exceeding the allowed parameters ). The shotcrete lining parameters are those described in section 4 and Table 3.
As evaluated in the previous section, Equation 14 checks the MCS accuracy considering some number of simulations and a shotcrete lining thickness c 10 cm t = . This work adopts 5E+06 simulations, which generate the following results for tunnels A and B: f 4.06% These results say that tunnel B needs more shotcrete than tunnel A to get acceptable reliability levels due to the far-field pressure in tunnel B is much higher than in tunnel A.   The results in Table 8 Figures 10a and 10b show the support failure region, between 6.00 and 7.00 mm approximately, where ( ) (tunnel B) have close values. This occurs due to eq u standard deviations since the eq u histogram in Figure 10a is sparser than the corresponding histogram in Figure 10b. As well as in Figure 10, it is easy to evaluate the failure probability through ( )    Figures 12a and 12b show an analysis through the CV-CF method, using the mean values of tunnels A and B to verify the accordance with the results defined in Table 8. In these cases, it has been considered c 13 cm t = and c 20 cm t = for tunnels A and B, respectively. The results in Figure 12 must be interpreted with caution due to eq U referring to the shotcrete lining displacement, which starts from the displacement 0 U . There are similarities between the design parameters defined in Figures 12a and 12b and the design parameters of Table 8, where it can be seen that the interaction between the Convergence and Confinement curves has been influencing the reliability analysis. The Confinement curve of tunnel A tends to perform better than tunnel B due to the internal pressure i P being higher in tunnel B.
These results suggest that the shotcrete lining of tunnel A, considering c 13 cm t = , produces acceptable stiffness results s K , which increases the inclination of the Convergence curve in Figure 12. Tunnel B, likewise, produces acceptable s K results using c 20 cm t = .
The designer must consider adopting support parameters higher than the before adopted, aiming to increase s K , thus keeping the failure probability low.

SUMMARY AND CONCLUSIONS
This study performed a reliability analysis in deep axisymmetric tunnels, integrating FORM and MCS methods with tunnel analytical methodologies. VBA and MATLAB codes verified the reliability in rock masses and shotcrete lining, examining two different tunnels.
The first aim was to check the plastic zone reliability through VBA FORM and MATLAB MCS algorithms. The reliability results verification has shown an excellent agreement between the reliability indexes and failure probabilities, obtained by both algorithms, and with the Li and Low [9] results.
The reliability study of the variation in the resistance L, considering unsupported tunnels, suggests that, in general, it is necessary to pay attention to the adopted L value in rock masses with low dispersion. This work recommends adopting a resistance value L around the mean value to avoid higher f p values. The increase of internal pressure i P produces great changes in the plastic zone reliability. Both tunnels develop significant f p reductions when i 0.3MPa for tunnel B. Moreover, the internal design pressure eq P for these thicknesses is greater than 0.3MPa , which keeps plastic zone failures close to zero.
This paper contributes to check safety or failure probabilities in tunnel design stages. The reliability codes provided statistical data for different rock mass parameters. So, designers and engineers can use these methodologies to ensure security in underground design.
A natural progression of this work is to link two-dimensional and three-dimensional numerical models with the reliability codes. Moreover, this study suggests using reliability analysis in different cross-section geometries and supports (steel lining, precast concrete lining, bolts).