3D Numerical Investigation of Face Stability in Tunnels With Unsupported Face

The paper studies the stability of unsupported tunnel faces by analyzing the results of a large number of 3D numerical analyses, in various ground conditions and overburden depths. The analyses calculate the average face extrusion (Uh) by averaging the axial displacement over the tunnel face. Limiting face stability occurs when the average face extrusion becomes very large and numerical convergence becomes problematic. Using the results of the numerical analyses, a dimensionless “face stability parameter” is defined, which depends on a suitable combination of ground strength, overburden depth and tunnel size. The face stability parameter correlates well with critical tunnel face parameters, like the safety factor against face instability, the average face extrusion, the radial convergence of the tunnel wall at the excavation face, the volume loss and the deconfinement ratio at the tunnel face. Thus, semi-empirical formulae are proposed for the calculation of these parameters in terms of the face stability parameter which is obtained from basic tunnel and ground parameters. Thus, useful conclusions can be drawn for the safety factor against face instability, the volume loss and the deconfinement ratio at the excavation face.

Control of face stability is very important in tunnelling, as incidents of face instability are frequent, severely affect the cost and construction schedule of tunnels and can damage surface structures and utilities in shallow urban tunnels. In mechanized tunnelling with active face pressure (e.g. EPB and Slurry TBMs), the risk of face instability is controlled by the applied face pressure (e.g. Litsas et al. 2017), which is often adjusted empirically from past performance in ''similar'' conditions (e.g. behaviour in previously excavated tunnel sections). In typical tunnelling projects using conventional techniques (SCL / NATM), face stability is often assessed empirically or by simplified limit equilibrium analyses (e.g. Leca and Dormieux 1990;Anagnostou and Kovari 1996;Kim and Tonon 2010), while numerical analyses calibrated with systematic measurements of face movement are rather sparse. When the risk of face instability is considered unacceptable in SCL / NATM tunnels, the size of the excavation face is reduced or active face support measures are applied, such as fiber-glass (FG) nailing, forepoling, or even leaving a ground wedge to provide some lateral pressure on the excavation face. The main reason of the extensive empiricism in assessing face stability, is the lack of systematic measurements of face deformations and that quantitative assessment of face instability requires the definition of a suitable ''safety factor'', and its calculation using complex three-dimensional (3D) numerical analyses with realistic constitutive models and suitably measured/ estimated ground parameters. Although seemingly trivial, even the definition of a ''safety factor'' for face stability analyses is not always straight forward, let alone its numerical calculation. In mechanized tunnelling, the safety factor against face instability is usually defined as the ratio of the applied face pressure to the minimum face pressure required for stability. Calculation of the safety factor requires determination of that minimum face pressure for stability, which is usually achieved by a variance of the increased external load method (Zienkiewicz et al. 1975) via a series of numerical analyses with gradually decreasing face pressure until the tunnel face becomes unstable (i.e., until the numerical model ceases to converge, or face displacements start to increase rapidly). Similar definition of the safety factor and calculation techniques can be used in conventional tunnelling with supported excavation faces, e.g. in cases where the tunnel face is supported by a grid of fiber-glass nails providing an ''equivalent'' face pressure.
In the very common case of conventional tunnelling with unsupported excavation face, the ''safety factor'' against face instability can be defined (and calculated) as the ratio of some ''strength'' over a corresponding ''applied shear stress''. In continuum mechanics analyses (i.e., excluding structurally-controlled instabilities), this definition is straight forward when ground strength is modelled via perfect plasticity with the Mohr-Coulomb (MC) failure criterion, i.e., in: 1. Analytical methods (e.g. Horn 1961, Atkinson andMair 1981;Panet 1995;Anagnostou and Kovari 1996), which calculate the safety factor of the tunnel face on a suitably selected potential failure surface, by some form of limit equilibrium of a critical ground wedge at the excavation face. These methods are widely used in practice, although they include simplifying assumptions about the selected wedge for a complex 3D problem as the tunnel excavation face. 2. Numerical methods, where the ''safety factor'' is usually defined and calculated by the Strength Reduction Method (Zienkiewicz et al. 1975), i.e., by performing a series of analyses with gradually reducing ground strength, until the tunnel face becomes unstable (i.e., until the numerical algorithm ceases to converge, or face displacements start to increase rapidly). In such analyses, the safety factor is the inverse of the strength reduction factor causing face instability. Useful design charts are often produced for the safety factor versus ground strength, tunnel depth and size (e.g. Kavvadas et al. 2009, Prountzopoulos 2012. The above methods exploit the property of the MC failure criterion that ground strength is a linear combination of cohesion (c) and friction angle (tanu). Thus, face instability can be achieved numerically by applying the same ''strength reduction factor'' to both components of strength; the safety factor is the inverse of that factor.
When ground behaviour is modelled more realistically than Mohr-Coulomb perfect plasticity, investigation of face stability requires the use of numerical analyses. Published literature on numerical analyses of face stability in tunnels with unsupported face using such constitutive laws (e.g. based on the Hoek-Brown failure criterion and/or hardening/softening plasticity) is very sparse, because such analyses are usually problem-specific, i.e., they check if a specific tunnel face is stable, by testing the convergence of the numerical model for given ground and geometrical parameters, but are difficult to generalise in other cases. Furthermore, in continuum numerical analyses, ''stable'' faces correspond to relatively small face deformations, while face instability is often related to a non-converging analysis (i.e., large deformations). In ''stable'' faces, it is not easy to define the available safety factor or calculate the margin from face instability. The reason of this difficulty is that, in constitutive laws other than Mohr-Coulomb perfect plasticity, ground strength is controlled by non-linear combinations of model parameters, rendering the strength reduction method inapplicable. Furthermore, other analogous techniques (like increasing suitable external loads until failure) cannot be applied in tunnel excavation with an unsupported face, because a stable tunnel face does not have any external load (face pressure is zero). This common problem becomes evident in designs attempting to apply the ''partial factor method'' in Ultimate Limit State (ULS) analyses of stability problems with ground failure controlled by criteria other than Mohr-Coulomb perfect plasticity and/or cases where ground failure is not caused by external loads as in bearing capacity of footings (see e.g. Frank et al. 2004;Franzen et al. 2019).
In conclusion, although numerical analyses can be performed to check if a tunnel face is stable for specific ground and geometrical parameters, there is lack of guidance in assessing the available safety factor of unsupported tunnel faces and difficulty in using measurements to predict upcoming face instability, because systematic deformation and ground stiffness measurements at the excavation face are very sparse due to technical difficulties in a continuously advancing face. Even in cases where complex 3D numerical analyses are performed to study specific tunnel conditions, it is useful to have guidance on the effects of varying ground conditions and/or tunnel depth on face stability, without having to perform additional analyses for each case. Furthermore, it is useful to have guidance in optimally selecting the required analyses of face stability, among the usually wide range of ground conditions and tunnel depths in practical tunnelling problems. In such cases, it is useful to have guidance from results of full 3D numerical models, which are more accurate than axisymmetric tunnel models commonly used in face stability analyses (e.g. Bernaud and Rousset 1996;Graziani et al. 2005), especially in shallow tunnels where the effect of gravity is more pronounced and conditions of face instability are more frequent and catastrophic.
The present paper attempts to fill that gap, by providing a semi-empirical expression of an equivalent safety factor against face instability in tunnels with unsupported face, in terms of dimensionless quantities of ground strength, tunnel depth and diameter. The safety factor of face stability is obtained from a dimensionless ''face stability parameter'' (K f ) which is found (numerically) to control the average ''face extrusion'' (U h = average axial displacement of the excavation face) for a wide range of ground strengths, failure modes, tunnel depths and sizes. The paper also proposes semi-empirical expressions to calculate the average face extrusion (U h ) and the degree of deconfinement (k) at the tunnel face, in terms of the controlling face stability parameter (K f ). These expressions are derived from the results of a large set of three-dimensional (3D) numerical analyses of the excavation of shallow (H/D = 2.5 -5) and deep tunnels (H/D = 10 -20) with unsupported face in various ground conditions, using the Mohr-Coulomb failure criterion in shallow tunnels (where stiff soils are predominant) and the Hoek-Brown failure criterion in deep tunnels (where rockmass is usually encountered). The analyses focus on the behaviour of the excavation face, by calculating the average face extrusion (U h ), suitably normalized to give a dimensionless ''face extrusion parameter'' (X f ). The results of the analyses show that the face extrusion parameter (X f ) is correlated well with the ''face stability parameter'' (K f ) which depends on ground strength, tunnel depth and size. It is shown that, as K f decreases and approaches unity, the face extrusion parameter (X f ) starts to increase rapidly, indicating incipient face instability. This allows to define and calculate the safety factor of an unsupported tunnel face (SF f ) by the face stability parameter (i.e., SF f = K f ) and establish a relationship among ground strength, tunnel depth and size at limiting face instability (when K f = 1). The proposed semi-empirical relationship between K f and X f can be used to calculate the average face extrusion (U h ), the radial wall convergence (U R ) and deconfinement ratio (k) for various combinations of ground strength, tunnel depth and size. The proposed relationship can be used in preliminary calculations of the safety factor and the degree of deconfinement of the tunnel excavation face, in conventionally excavated tunnels with unsupported face.
The numerical analyses used to obtain the proposed correlations have not been verified by actual data from real tunnels, because systematic tunnel face deformation and ground stiffness measurements (required to calibrate the 3D numerical models) are very sparse. Furthermore, although the analyses include both shallow and deep tunnels in a wide range of ground conditions, the produced correlations do not include very shallow (H/D \ 2.5) or very deep (H/D [ 20) tunnels, where face instability mechanisms may be different. Finally, the results do not include cases of face instability controlled by structural discontinuities, because the numerical analyses treat ground as a continuum, and thus are applicable in cases of stiff soils and very weathered or heavily fractured rockmasses.

Numerical Analyses
A large set of three-dimensional (3D) numerical analyses were performed, using the commercial Finite Element Code Simulia Abaqus, for the excavation of shallow (H/D = 2.5 to 5) and deep tunnels (H/D = 10 to 20) with unsupported face and a wide range of ground properties and tunnel depths. Typical ovalshaped tunnel sections were studied, with width D = 10 m and 6 m (Fig. 1).

Shallow Tunnels
In shallow tunnels, the overburden depth (H), measured from the tunnel axis up to the ground surface, varied in the range H = 15 to 30 m, with examined cases: H/D = 2.5, 3.5 and 5. More shallow tunnels were not examined, as the mechanisms controlling face stability are different (chimney effects) and ground conditions in the very shallow overburden are rarely uniform. Eight-node hexahedral finite elements with full integration were used in the analysis (Fig. 2). Following a sensitivity analysis, the extent of the finite element mesh was sufficiently large to minimize boundary effects in all directions. The finite element mesh included the left half of the tunnel, because the tunnel section is symmetrical with respect to the vertical axis. The tunnel was excavated in a single phase (full face excavation) with excavation steps of 1 m (equal to the size of the elements in the axial direction). In each excavation step, a relatively stiff, 30 cm thick, shotcrete liner was installed on the tunnel wall (full ring) at distance 1 m behind the excavation face. The shotcrete liner was modelled by 4-noded shell elements, as linearly elastic with a relatively low concrete E modulus equal to 10 GPa to account for concrete setting time. To eliminate end effects, the conditions at the excavation face were studied when tunnel excavation had advanced to a distance L = 4-6 D from one end boundary, leaving a clear distance 5-7 D from the other end boundary The study included relatively stiff ground conditions with unit weight c = 20 kN/m 3 , horizontal geostatic stress coefficient K o = 0.5 and 1.0 and linearly elastic-perfectly plastic behaviour, with elastic modulus (E) and yielding according to the Mohr-Coulomb criterion (c = cohesion, u = friction angle). Table 1 shows the sets of ground parameters used in the parametric analyses.
In all cases, the elastic Poisson ratio was v = 0.33. According to the Mohr-Coulomb failure criterion, the ''ground strength'' (r cm ), equivalent to the Uniaxial Compressive Strength, was calculated by the formula: The total number of numerical analyses for the shallow tunnels was 72 (two tunnel sizes, three tunnel depths, two K o values, and six sets of material parameters).

Deep Tunnels
In deep tunnels, the overburden depth was H = 100, 150 and 200 m. The finite element mesh was similar to   Fig. 2, with higher overburden depth. Tunnel excavation and liner construction followed the same procedure as for the shallow tunnels.
The study included weak fractured rock with unit weight c = 25 kN/m 3 , horizontal geostatic stress coefficient K o = 0.5 and 1.0, intact rock properties r ci = 10 MPa and E i = 2000 MPa, Poisson ratio v = 0.33 and Geological Strength Index (GSI) in the range 15 to 45. The rockmass was assumed linearly elastic-perfectly plastic, yielding according to the Generalised Hoek-Brown failure criterion (Hoek et al. 2002) with various parameters (m b , s, a). Table 2 shows the sets of rockmass parameters used in the parametric analyses.
The ''rockmass strength'' (r cm ) and ''rockmass modulus'' (E m ) for the various GSI values were calculated by the following empirical formulae (Hoek and Diederichs 2006): The total number of numerical analyses for the deep tunnels was 48 (two tunnel sizes, three tunnel depths, two K o values, and four sets of material parameters).

Face Extrusion
Each of the numerical analyses calculates the axial displacement (face extrusion) at all integration points on the tunnel face when tunnel excavation has advanced far from the side boundaries. These values are averaged over the tunnel face to give an ''average face extrusion'' (U h ) which is then normalized by the tunnel width (D) and a modulus-to-depth factor (E / p o ) to give the dimensionless ''face extrusion parameter'' (X f ): where E is the elastic Young modulus of the ground (soil or rockmass) and p o = 0.5 (1 ? K o ) c H is the average overburden pressure at the tunnel axis (average of vertical and horizontal geostatic stresses). From the 72 (shallow) ? 48 (deep) = 120 numerical analyses, the present database includes the results of 83 analyses (51 shallow and 32 deep tunnels), as the remaining 37 (21 ? 16) analyses failed to converge, as the combination of ground strength and tunnel depth (ground stress) produced uncontrollable face extrusions (too low strength for the tunnel depth). Each calculated value of the face extrusion parameter (X f ) was then correlated with the corresponding values of various forms of strength-to-stress ratios, with the objective to select the optimal form (giving the best correlation). Figure 3a plots the calculated face extrusion parameter (X f ) versus the corresponding value of the classical ratio of rockmass strength (r cm ) to average overburden pressure (p o ), often used to describe tunnel behaviour (e.g. Hoek 2000). The correlation of the two parameters is poor, especially at low strength-to-stress values (r cm /p o \ 0.5), where face stability problems are expected, indicating that r cm // p o is not optimal for face stability analysis. Figure 3b optimizes the correlations of Fig. 3a, by plotting the calculated face extrusion parameter (X f ) with the semi-empirical dimensionless ''face stability parameter'' (K f ), an optimal expression of the strength-to-stress ratio combining ground strength (r cm ) and overburden stress (cH) with the depth-tosize ratio (H/D) and the K o parameter, by the formula: Despite the possibly different face stability mechanisms of shallow and deep tunnels, the calculated face extrusions fit in a narrow band. This is probably due to the fact that the numerical analyses have not studied very shallow tunnels (H/D \ 2.5) where chimney-type failure mechanisms become more pronounced. The best fit curve of the data points shown in Fig. 3b is expressed by the formula: This formula can estimate the face extrusion parameter (X f ) and, via Eq. 2, calculate the average face extrusion (U h ) for given K f , i.e., a shallow or deep tunnel with of size (D), overburden depth (H) in ground with strength (r cm ). Control analyses have shown that the above formula can also be used in tunnel shapes different than those shown in the present study ( Fig. 1), including excavation of the top heading of a tunnel, via an equivalent tunnel size: D = 1.15 sqrt(A), where A is the section area of the tunnel or phase.
In Fig. 3b, large (K f ) values correspond to good face stability conditions, i.e., large ground strength, and/or relatively shallow and small size tunnels, with degrading face stability conditions (i.e., increasing face extrusion) as (K f ) decreases. The scaling factor ''3.8'' in the definition of (K f ) (Eq. 3) was selected such that, when K f = 1, the rate of the face extrusion parameter (X f ) increases rapidly, indicating that the tunnel face approaches limiting face stability, although the numerical model converged in about 50% of the models with K f \ 1 giving large face extrusions. Based on this remark, the condition K f = 1 provides limiting face stability, while tunnel faces with K f \ 1 are considered unstable. At limiting face stability, Eqs. (2), (3) and (4) give the limiting face extrusion and limiting ground strength: where: (r cm ) lim is the lowest ground strength to ensure limiting face stability for a given tunnel size (D) and overburden depth (H). For a specific ground type at the excavation face, the available ground strength can be calculated from Eq. (1a) for soils and Eq. (1b) for rockmasses and compared to the limiting value (Eq. 5b) to assess whether the tunnel face is stable or not. Figure 4 plots the above limiting ground strength (r cm ) lim versus (H/D) for two values of the horizontal geostatic stress coefficient K o = 0.5 and 1.0. For example, in a tunnel with overburden depth H = 4 D, the limiting ground strength for face stability is: The above definition of the face stability parameter (K f ) can assist in the calculation of the safety factor (SF f ) of tunnel faces against instability, by defining the safety factor as the ratio of the available ground strength to the corresponding limiting ground strength, and using Eqs. (5b) and (3): Thus, the safety factor of the tunnel against face instability (SF f ) is equal to the face stability parameter (K f ), i.e.: Combining Eqs.
(2) and (4), the average face extrusion (U h ) for given safety factor (SF f ) is given by the formula: Figure 6 plots the predicted average face extrusion (U h ) versus the safety factor of the tunnel face, for typical values of the modulus-to-strength ratio E/p o-= 75, 100 and 150. For SF f \ 1, the average face extrusion increases rapidly. At limiting face stability (SF = 1), the average face extrusion is equal to 1-2% of the tunnel size (D). Figure 7 shows the profile, along the tunnel axis, of the radial convergence (U R ) of the tunnel wall and the distribution of the face extrusion (U h ) on the tunnel face. Radial wall convergence occurs and in the tunnel core, ahead of the tunnel face. The volume loss (VL) is defined as the reduction (DV) of the volume (V) of the core ahead of the tunnel face per unit volume of the core. Volume loss is caused by the radial convergence of the tunnel wall in the core which ''squeezes'' the core giving face extrusion. Assuming that the profile of the radial wall convergence in the core (length L, tunnel section area A) is approximately linear, with maximum value of the radial wall convergence at the excavation face equal to (U R ), then:

Radial Wall Convergence, Volume Loss and Deconfinement
and the volume loss (VL) is: The radial convergence (U R ) of the tunnel wall at the excavation face can be obtained from the average face extrusion (U h ) (calculated via Eqs. 7 and 6) by assuming that the deformation of the core occurs with practically no volume change, i.e., the reduction (DV) of the volume of the core is equal to the ground volume The examined tunnels have a section area A = 0.75 D 2 . The length (L) of the core was calculated by correlating the average radial displacement (U R ) at the tunnel face computed in the numerical analyses with the radial displacement predicted from the face extrusion (U h ) via Eq. (10). The best fit is achieved for L = 0.38 D (Fig. 8).
Thus, Eq. (10) gives (using Eq. 7): and the volume loss (VL) can be expressed as (combining Eq. 9 with 11 and A = 0.75 D 2 ): Figure 9 plots the calculated radial wall convergence (U R ) at the tunnel face (from Eq. 11) versus the face stability parameter (K f ), for typical values of the modulus-to-strength ratio E/p o = 75, 100 and 150. For typical stable faces (K f = 1 -2.5), the calculated radial wall convergence (U R / D) is in the range 0.5 -2.5%. Figure 10 plots the calculated volume loss at the tunnel face (from Eq. 12) versus the face stability parameter (K f ), for typical values of the modulus-tostrength ratio E/p o = 75, 100 and 150. For typical stable faces (K f = 1 -2.5), the calculated volume loss is in the range 0.5 -2.5%.
In 2D (plane strain) numerical analyses, the deconfinement ratio (k) is used to calculate a fictitious radial internal pressure (p i ) which produces the same inward radial convergence (U R ) of the tunnel wall as a corresponding 3D model which, unlike the 2D model, includes the effects of the excavation face (Fig. 7). By definition, the internal pressure (p i ) is related to the deconfinement ratio (k) by the formula: The relationship between (U R ) and (p i ) (or k) is the convergence-confinement relationship, calculated using several methods, such as Duncan Fama (1993), Panet (1995), Kavvadas (1998), Carranza-Torres et al. (2002) and Carranza-Torres (2004).
The deconfinement ratio (k) and the corresponding internal pressure (p i ) vary with the distance (x) from the excavation face. As the radial wall convergence increases along the tunnel axis, and the corresponding internal pressure decreases, the deconfinement ratio varies from k = 0 far ahead of the excavation face (where wall convergence is zero) to k = 1 far behind the excavation face (where wall convergence is stabilized to the maximum value). Several semiempirical formulae have been proposed for the calculation of (k) at various distances (x) from the excavation face. These formulae are produced by equating the radial convergence (U R ) of the tunnel wall from 2D analyses (applying an internal pressure p i ) with the corresponding radial convergence profile along the tunnel axis from 3D finite element analyses (e.g., Panet 1995;Chern et al. 1998;Vlachopoulos & Diederichs 2009).
The methodology developed above can provide an empirical relationship between the deconfinement ratio (k) at the tunnel face and the corresponding face stability parameter (K f ), by correlating the value of (k) at the tunnel face, computed using the above convergence-confinement relationships, with the face stability parameter (K f ), for each of the 87 numerical analyses studied. Figure 11 plots the results of this correlation using four alternative convergence-confinement methods. The best fit curve of the correlation is: Stable tunnel faces (K f [ 1) have deconfinement ratios k = 0.30-0.70, with higher k values for unstable tunnel faces (K f \ 1).

Conclusions
The paper studies the stability of unsupported tunnel faces by analyzing the results of a large set (87 Nos) of 3D numerical analyses of tunnel faces, in various ground conditions and overburden depths. The analyses calculate the average face extrusion (U h ) by averaging the axial displacement over the tunnel face. Limiting face stability occurs when the average face extrusion becomes very large and algorithmic convergence becomes problematic. Using the results of the analyses, a dimensionless ''face stability parameter'' (K f ) is defined (Eq. 3) which depends on a suitable combination of ground strength (r cm ), overburden depth (H) and tunnel width (D). The (K f ) parameter correlates well with critical tunnel face parameters, like the safety factor against face instability (Eq. 6, Fig. 5), the average face extrusion (Eq. 7, Fig. 6), the radial convergence of the tunnel wall at the excavation face (Eq. 11, Fig. 9), the volume loss (Eq. 12, Fig. 10) and the deconfinement ratio at the tunnel face (Eq. 14, Fig. 11). Thus, semi-empirical formulae are proposed for the calculation of these parameters in terms of the face stability parameter. Since the face stability parameter can be easily calculated from basic tunnel and ground parameters, the above critical tunnel parameters can be calculated, and conclusions can be drawn about tunnel face stability, volume loss and the deconfinement ratio at the excavation face which can be useful in preliminary assessments of tunnel behaviour. Furthermore, the calculated volume loss can be used to estimate ground surface settlements in shallow tunnels, while the deconfinement ratio can be used in 2D numerical analyses of tunnel excavation and support.
The numerical analyses used to obtain the proposed correlations have not been verified by actual data from real tunnels, because systematic tunnel face deformation and ground stiffness measurements (required to calibrate the 3D numerical models) are very sparse. Furthermore, although the analyses include both shallow and deep tunnels in a wide range of ground conditions, the produced correlations do not include very shallow (H/D \ 2.5) or very deep (H/D [ 20) tunnels, where face instability mechanisms may be different. Finally, the results do not include cases of face instability controlled by structural discontinuities, because the numerical analyses treat ground as a continuum, and thus are applicable in cases of stiff soils and very weathered or heavily fractured rockmasses. Fig. 11 Correlation of the deconfinement ratio (k) at the tunnel face (computed using four convergence-confinement methods) with the face stability parameter (K f ) for each of the 87 numerical analyses studied. The figure also shows the best-fit curve (Eq. 14)