A Review of Seven Tunnel Face Stability Models

: This paper presents a review of selected tunnel stability models that have been developed and used in calculating the minimum tunnel face pressure as described by original authors. Further-more, this paper provides a comparison of required tunnel face pressure obtained from analytical models, based either on limit equilibrium method or the limit analysis method (upper bound theorem) and numerical models using the ﬁnite element method. The numerical results are presented in charts for the comparative study to discuss the inﬂuence of cover depth to tunnel diameter ratio (C/D), internal friction of the soil ( ϕ ), and cohesion (c) on normalized support pressure ( p u / γ D ) for each model. To verify the accuracy of the selected models, a comparison of the results of seven tunnel stability models with the results of the physical models is carried out. In a ground composed of two layers, a comparison of the required tunnel face pressure is presented. The results show that the wedge–silo models provide higher support pressure than the conical block models. Moreover, the support pressure using the conical block models is only dependent on the friction angle and not on the C/D ratio. Finally, the results of wedge-silo models indicate more signiﬁcant dependence of the required support pressure on the C/D ratio especially for the lower friction angle.


Introduction
During excavation and construction of a tunnel, the stress history of in situ soil surrounding the excavation area is disturbed. Thereafter, a new stress redistribution is induced in the vicinity of the tunnel face due to deformation of the soil mass. The mechanism that leads to this stress redistribution is the so-called arching phenomena (Terzaghi [1]). This mechanism yields a reduction of vertical stress. The arching soil effect has been reported by Terzaghi [2] in trap door experiments and also by Chambon and Corte [3] and Chen et al. [4] conducting a series of physical model tests to investigate the stability of tunnel face.
For preventing the collapse of the tunnel face during the excavation process, adequate support pressure (e.g., compressed air, slurry, or earth pressure support) at the tunnel face is required. Meanwhile, the soil arching effect is associated with a decrease in the required support pressure. To consider the 3D arching soil effect in evaluating the required support pressure, analytical methods are commonly used. The analytical approaches can basically be divided into two groups, namely, Limit Equilibrium Method (LEM) and Limit Analysis Method (LAM). The Limit Equilibrium Method is widely used as a theoretical analysis technique for the stability of the tunnel face. On the other hand, the Limit Equilibrium Method applies the static equilibrium between the forces acting on the soil mass for the critical collapse mechanism. In order to assess the required support pressure, many types of limit equilibrium models have been setup on a 3D wedge-silo mechanism (e.g., Horn [5]; Jancsecz and Steiner [6]; Anagnostou and Kovari [7]; Broere [8]; Kirsch and Kolymbas [9]; Chen et al. [10] ), see Figure 1. However, to incorporate the contribution of the 3D soil arching effect in the calculation of required support pressure, the activation of shear forces on the silo slip surfaces, as well as consideration of shear forces acting on the flanks of the prismatic wedge, are taken into account in the static equilibrium of the forces. The main problem of equilibrium is associated with estimating the shear forces at the slip surfaces. Indeed, the value of shear forces/ stresses depends on the horizontal stresses, which cannot be calculated from equilibrium conditions (Anagnostou [11]). Thus, to overcome this statically indeterminate task, the horizontal stress (σ h ) is assumed to be proportionally dependent on the vertical stress (σ v ). The ratio between the horizontal stress and the corresponding vertical stress represents what is known as the lateral earth pressure coefficient. The value of the lateral earth pressure coefficient must be assumed in advance. Recently, the upper bound theorem of limit analysis became an effective method for studying the stability of the tunnel face. The kinematic approach of this theorem is based on constructing admissible failure mechanisms. To apply the upper bound theorem, some assumptions must be specified (Chen [12]): (1) The strain rates resulting from the velocity field must satisfy the flow rule; (2) For a kinematically admissible failure mechanism, the velocity along a plastic deformed surface must make an angle ϕ with the discontinuity velocity, which fulfils the normality rule; (3) The rate of work of the external forces is equal to the rate of internal energy dissipation. Nevertheless, due to the assumptions in the upper bound theorem, it is difficult to simulate the 3D arching phenomena at the front of the tunnel face. To overcome this difficultly in the analysis, one or more truncated conical sliding wedges is assumed to be the failure mechanism, leading to a reduction of the self weight and thereby to include the 3D soil arching effect in evaluation of the support pressure (e.g., Leca and Dormieux [13]; Mollon et al. [14]; Tang et al. [15]; Ibrahim et al. [16]; Zou et al. [17]; Li and Zhang [18]).
The application of numerical methods with advanced constitutive models has been shown to substantially improve the analysis of the tunnel face stability, using different scenarios for the 3D aspect of tunnel. Finite Element Method (FEM), Discrete Element Method (DEM) and Finite Difference Method (FDM) have proven useful tools in simulating the stability of tunnel face (e.g., Peila [19]; Ohta and Kiya [20]; Vermeer et al. [21]; Kirsch [22]). However, the simulation of the tunnel face failure, using three dimensional DEM, FEM or FDM models, is very time consuming and expensive to model.
More recently, the Kinematical Element Method (KEM) has been utilized to calculate the active pressure of the tunnel face for three-dimensional analysis (e.g., Qarmout et al. [23]; Qarmout et al. [24]).
For investigating the ground soil respond to tunneling, various physical model tests were conducted, using centrifuge model or small scale model tests (e.g., Chambon and Corte [3]; Takano et al. [25]; Kirsch [22]; Idinger et al. [26]; Lü et al. [27]). The physical models provide a beneficial information about ground surface settlement, the required face support and insight view into the 3D arching failure mechanism at the front of the tunnel face. Moreover, the results of the experimental model tests can be used to validate the numerical and analytical models.
This paper presents a review of seven tunnel face stability models that have been used to calculate the minimum support pressure, covering three different approaches, Limit Equilibrium Method (LEM), Limit Analysis Method (LAM) and Finite Element Method (FEM). Using the results for selected tunnel face stability models as described by the original authors, compression results were assembled to enable an assessment of each model. The comparative study includes the cover depth to tunnel diameter ratio (C/D), internal friction of the soil (ϕ) and cohesion (c). Besides, the results given by each model are compared with the results of physical model tests available in the literature. Furthermore, in the case of a ground composed of two soil layers, a comparison of the required support pressure of these models is presented.

Overview of Recent Advances in Tunnel Face Stability Analysis
As mention before, a variety of modeling strategies have been introduced by researchers for simulating the stability of the tunnel face, such as analytical approaches, numerical methods, and experimental model tests. The research activities for each method will be briefly summarized in the following subsections.

Limit Equilibrium Method (LEM)
The first systematic study for 3D face stability of the tunnel was performed by Horn [5]. Horn [5] presented a 3D wedge-silo failure mechanism, replacing the circular shape of the tunnel face with a square shape, as shown in Figure 1a. The horn failure mechanism consists of a prismatic wedge and a vertical silo above the wedge. In the horn failure mechanism a number of proper choices are proposed for adopting the distribution of vertical stress with depth. Latterly, the wedge-silo model has been used by several researchers as a basis for further development. Using the Limit Equilibrium Method, Anagnostou and Kovari [7], Jancsecz and Steiner [6], Broere [8], and Anagnostou [11] proposed a 3D tunnel face stability model based on the Horn failure mechanism, including the effect of soil arching in the silo by using Janssen's silo theory [28].
The silo theory composed the following equation for calculating the vertical stress acting on the base of the silo: where γ is the density of the soil, ϕ is the friction angle, a is the ratio of the area over circumference of horizontal plane of the silo, K silo is the lateral pressure coefficient and z is the soil depth from the ground surface.
Two assumptions are adopted in the wedge-silo model based on Janssen's analysis of soil arching. The first assumption is that the lateral earth pressure coefficient of the soil in the silo is assumed in advance; it is also constant over the tunnel depth. The second assumption is that the vertical stresses are uniformly distributed across any horizontal section of the silo as well as at the base of the silo. In addition to the previous assumptions, the shear force between the silo and the wedge is omitted in the equilibrium of the forces.
A number of authors have suggested a different value for K wedge based on practical experience. Anagnostou and Kovari [7] assumed K wedge = 0.4. Jancsecz and Steiner [6] suggested the lateral earth pressure coefficient for the wedge as the average of K a and K 0 (K 0 = 1 − sin ϕ, Jaky [29], K wedge = (K a + K 0 )/2, see Figure 2. Figure 2 shows the value of K silo and K wedge with different friction angles using various methods. According to Jancsecz and Steiner [6] model, the K silo and K wedge values decrease with the increase in friction angle, this tendency is also revealed by Broere [8]. Whereas, in the Anagnostou and Kovari [7] and Anagnostou [11] models, K silo and K wedge are assumed to be a constant where they are independent of friction angle. As shown in Figure 2a, the value of K silo obtained from the Anagnostou and Kovari [7] model is larger than the other models for the same friction angle.  In order to figure out the distribution of horizontal stress σ h (z) inside the wedge, Anagnostou and Kovari [7] and Jancsecz and Steiner [6] were assumed to be linearly proportional between the horizontal stress σ h (z) and the vertical stress σ v (z). The vertical stress inside the wedge is calculated as follows: where C is the cover depth. The average shear stress (τ wedge ) acting on each triangle slip surfaces of the sliding wedge is obtained by integrating K wedge · σ v (z)· tanϕ over the slip surfaces. The average shear stress over triangle slip surfaces is calculated using the following equation: Anagnostou and Kovari [7] assumed K wedge = 0.4. Jancsecz and Steiner [6] suggested the lateral earth pressure coefficient for the wedge as the average of the summation of K a and K 0 (K 0 = 1 − sin ϕ; Jaky [29]), K wedge = (K a + K 0 )/2, see Figure 2b. Jancsecz and Steiner [6] described the relation between horizontal and vertical stress at the tunnel axis by 3D earth pressure coefficient, the minimum support pressure (p) can be calculated as: Broere [8] and Anagnostou [11] used the infinitesimal slices method in the wedge, described by Walz and Prager [30] for slurry-filled trenches. The wedge is subdivided into smaller horizontal slices bodies i (i is the number of each slices), possibly of different soil conditions, see Figure 1c. The forces equilibrium on each horizontal slices are formulated. The minimum support pressure is calculated based on integrating the horizontal and vertical forces equilibrium over the whole wedge for each horizontal slice. The method of slices assumed linear proportion between the horizontal stress σ h (i) and the vertical stress Broere [8] proposed K (i) [29].Anagnostou [11] assumed K (i) wedge = 0.5. Kirsch and Kolymbas [9] considered the horizontal shear force into the wedge model between the wedge and silo. Kirsch and Kolymbas [9] constructed a simple equation for calculating the support pressure, which follows from the equilibrium condition by assuming a quadratic parabolic distribution of vertical stress above the tunnel crown. Kirsch and Kolymbas [9] assumed K silo = K wedge = K 0 (K 0 = 1 − sin ϕ). The vertical stress distribution and the failure mechanism for different arching implementations shown in Figure 3. Unlike the wedge-silo models, Krause [31] established a 3D shell failure mechanism regarding the internal stability of the soil at front of the tunnel face. The shear and cohesion forces along the sliding surface set up the resistance against the collapse. Based on equilibrium of the forces on the 3D failure body, he deduced a formula for calculating the minimum support pressure.

Limit Analysis Method (LAM)
The limit analysis method is applied to materials which can be idealized as perfectly plastic with the associated plastic flow rule. The concept of limit analysis is based on the theorems of plasticity developed by Drucker et al. [38], namely the lower and upper bound theorem. By using the lower and upper bound theorems, the range in which true collapse load is expected can be found.
The principle of the upper bound theorem is based on the work done by external load in an increment of displacements for a kinematically admissible mechanism being equal to the energy dissipated by internal stresses. This external load is not lower than the true collapse load. For this reason, it represents an upper bound on the actual solution.
The lowest possible upper bound solution is determined with an optimization scheme by trying various possible kinematically admissible failure mechanisms. The lower bound theorem states that if an internal stress field is in equilibrium with external loads without overcoming the yield criterion in the soil mass, the external load is not higher than the true collapse load. The highest possible lower bound solution can be determined by trying different possible statically admissible stress fields.
Many researchers have also used the upper bound method to examine the stability of the tunnel face. A number of authors has assumed different shapes of the failure mechanism to obtain the upper bound solution for calculating the minimum support pressure.
Leca and Dormieux [13] proposed two mechanisms for the failure zone at the front of the tunnel face; one consists of a single conical block and other is composed of two solid conical wedges with elliptic cross sections at the intersection of the tunnel face, see Figure 4. Both failure mechanisms are characterized by only one parameter, namely, the angle α between the axis of the cone and the horizontal tunnel axis. Leca and Dormieux [13] assumed the velocity of each rigid block collinear with the axis of each linked cone. This implies that the plastic energy dissipation along the discontinuities obeys the associated flow rule. The two-blocks mechanism given by Leca and Dormieux [13] is constrained by the normality condition required by plasticity theory. However, this condition does not allow the three-dimensional slip surfaces to develop more freely.
The minimum support pressure at the tunnel axis is expressed by the following equation: where N s and N γ are the non-dimensional coefficients and q s is the surcharge pressure. Leca and Dormieux [13] present two sets of graphs for non-dimensional coefficients N s and N γ with respect to friction angle. An improved failure mechanism composed of several rigid conical blocks was proposed by Mollon et al. [14]. This failure mechanism is an extension of the 3D failure mechanism developed by Leca and Dormieux [13], see Figure 5. Mollon et al. [14] found that a total number of five blocks is sufficient to calculate the minimum support pressure. The improvement of the solution by Mollon et al. [14] is due to the increase in the degree of freedom of the failure mechanism. Moreover, the mechanism is able to account for the whole circular tunnel face.
axis for the upper block is not adequate and leads to nonoptimal collapse pressures.
The failure mechanism presented by Mollon et al. ͑2009͒ and described in more detail in Oberlé ͑1996͒ is an improvement of the two-block collapse mechanism presented by Leca and Dormieux ͑1990͒. This mechanism is a multiblock ͑cf., Fig. 2͒. It is composed of n truncated rigid cones with circular cross sections and with opening angles equal to 2. A mechanism with n=5 is presented in Fig. 2. The geometrical construction of this mechanism is similar to that of Leca and Dormieux ͑1990͒, i.e., each cone is the mirror image of the adjacent cone with respect to the plane that is normal to the contact surface separating these cones. This is a necessary condition to ensure the same elliptical contact area between adjacent cones. In order to make clearer the geometrical construction of the 3D failure mechanism, Fig. 3 shows how the first two truncated conical blocks adjacent to the tunnel face are constructed. The geometrical construction of the remaining truncated conical blocks is straightforward. As for the mechanism by Leca and Dormieux ͑1990͒, Block 1 is a truncated circular cone adjacent to the tunnel face. The intersection of this truncated cone with the tunnel face is an elliptical surface that does not cover the entire circular face of the tunnel. This is a shortcoming not only of the multiblock mechanism by Mollon et al. ͑2009͒ but also of the two-block mechanism by Leca and Dormieux ͑1990͒. On the other hand, Block 1 is truncated with Plane 1 which is inclined at an angle 1 with the vertical direction ͑cf., Fig. 3͒. In order to obtain the same contact area with the adjacent truncated conical block, Block 2 is constructed in such a manner to be the mirror image of Block 1 with respect to the plane that is normal to the surface separating the two blocks ͑i.e., Plane 2 as shown in Fig. 3͒. The mechanism by Mollon et al.
Notice finally that the upper rigid Leca and Dormieux ͑1990͒ and Moll not intersect the ground surface dep values. This phenomenon of no outcr was also pointed out by Chambon and al. ͑2006͒ while they performed expe before, a failure soil mass which has does not necessarily outcrop at the g by these writers.
Both mechanisms by Leca and D et al. ͑2009͒ are translational kinem mechanisms. The different truncate mechanisms move as rigid bodies. T translate with velocities of different ear with the cones axes and make a discontinuity surfaces in order to res required by the limit analysis theory. determined by the condition that the r cones in contact has the direction tha contact surface.
The numerical results obtained by shown that a five-block ͑i.e., n=5͒ m cient since the increase in the number increases ͑i.e., improves͒ the solution provement of the solution by Mollon the one by Leca and Dormieux ͑1990 degree of freedom of the failure m ͑2009͒. Notice however that the solut and those by Leca and Dormieux ͑19 only an inscribed elliptical area to the involved by failure due to the conica the remaining area of the tunnel face and is contrary to what was observe This shortcoming will be removed in nisms developed in this paper.

Kinematical Approach for the of the Tunnel Face Collapse P
The aim of this paper is to compu pressure of a shallow circular tunn shield in a frictional and/or cohesive s based on a three-dimensional multiblo framework of the kinematical appr  The results obtained by Leca and Dormieux [13] and Mollon et al. [14] indicated that for cohesionless or a frictional-cohesive soil with a friction angle greater than or equal to 20 • , the minimum support pressure is independent of the tunnel cover depth.
More recently, the upper bound model used by Leca and Dormieux [13] was evolved for investigating the effect of layered soil on the minimum support pressure. Tang et al. [15] amended the solution of Leca and Dormieux [13] to be applicable in a layered soil. Tang et al. [15] studied the influence of soil properties of the crossed layered soil and the cover layered soil on the minimum support pressure. Their results indicated that the minimum support pressure is highly influenced by the shear strength of crossed soil more than the shear strength of the cover soil.
Zhang et al. [39] proposed a new 3D failure mechanism using the kinematic approach of limit analysis. The new 3D failure mechanism was composed of four conical blocks. The failure mechanisms only consider a portion of the tunnel face (an ellipse on a circular tunnel face). In addition, the numerical simulations using an FLAC3D were performed to analyze the stability of a tunnel with different tunnel cover to diameter ratios. Finally, the minimum support pressure obtained through the FLAC3D were compared to results obtained by the upper bound solution. It was found that the results from the upper bound solution were in very good agreement with the results of numerical simulations. Moreover, the failure zones are compared and the results corresponded with one another.
Using an FLAC3D, Li et al. [40] carried out numerical simulations to investigate the characteristics of the velocity distribution at the front of the tunnel face. Based on the velocity field obtained from the numerical simulations, an improved failure mechanism is constructed using the spatial discretization technique, which takes into account the soil arching effect. Both the minimum support pressure and the failure pattern were compared with the results of the numerical simulations. It was found that the critical face pressure provided by the improved failure mechanism corresponded well to that of the numerical simulation.
Ibrahim et al. [16] improved the 3D failure mechanism of Mollon et al. [14] to compute the minimum support pressure in dry multilayered purely frictional soil. The improved 3D failure mechanism is proposed for two and three soil layers.
Senent and Jimenez [41] extended the solution of Mollon et al. [14] to study the possibility of partial collapse in layered soils. The proposed model by Senent and Jimenez [41] examined the influence of soil properties of the crossed soil and the cover soil on the minimum support pressure.
Khezri et al. [42] investigated the effect of linear variation of cohesion with depth on the minimum support pressure. Their results show that adopting the mean soil cohesion that does not vary with depth would lead to conservative predictions for tunnel face support pressure. However, adopting the cohesive determined for the centreline of the tunnel underestimating the tunnel face support pressure and leads to an unsafe design.
Based on the kinematic approach of limit analysis, Pan and Dias [43] proposed a 3D failure mechanism to study the effect of anisotropy and cohesion on the stability of the tunnel face, including two types of non-homogeneous cohesion-linear variation with depth and layered soil. The results obtained from the presented approach show that both anisotropy and non-homogeneity have a significant effect on the support pressure, especially when the cohesion is relatively large.
Han et al. [44] proposed a 3D multiblocks failure mechanism for multilayered cohesivefrictional soils (see Figure 6). Their failure mechanism combines the silo theory (upper part) with the upper bound solution (lower part). The failure mechanism is composed of five truncated cones in the wedge. The distributed force acting on top of the truncated cone is calculated using silo theory with K silo = K 0 . The minimum support pressure is obtained using the upper bound solution in failure mechanism of the wedge. Fig. 1. Improved failure mechanism. Lee and Nam [45] included the effect of seepage forces emerging from the groundwater flow in the upper bound solution. They found that the minimum support pressure for the face stability is equal to the sum of the effective support pressure obtained from the upper bound solution and the seepage pressure acting on the tunnel face.
Pan and Dias [46] studied the stability of the tunnel face in water-bearing ground under steady seepage flow conditions. The distributions of to pore water pressure resulting from numerical calculation by FLAC3D are adopted to interpolate on the 3D failure mechanism of Mollon et al. [14]. The results indicated that the effective support pressure increases with water table elevation. The influence of anisotropic permeability on tunnel face stability is also discussed, showing that the isotropic model leads to an overestimation of the necessary tunnel face pressure for anisotropic soils.

Comparative Calculations Concerning Minimum Support Pressure
Analogous to the method proposed by Terzaghi [1] for bearing capacity analysis, the normalized support pressure (p u /γD) is represented by the following form (Vermeer et al. [21], Anagnostou [11], Mollon et al. [14]): where the contribution of different loads and soil parameters including self-weight (γ), internal friction angle (ϕ), surface surcharge (q) and cohesion (c) are expressed by the nondimensional bearing capacity coefficients N γ , N q and N c as a function of the friction angle of the soil. Due to the active conditions around the tunnel face the cohesion is reducing the necessary support pressure.

Homogeneous Soil
The comparative studies are performed on seven tunnel stability models, covering three different approaches-Limit Equilibrium Method (LEM), Limit Analysis Method (LAM), and Finite Element Method (FEM) calculations. The effect of friction angle (ϕ), cohesion (c), and depth to diameter ratio (C/D) on normalized support pressure (p u /γD) will be discussed in the following subsections. The soil conditions are assumed to be homogeneous. The dry unit weight is set to γ = 20 kN/m 3 , and the tunnel diameter is set to D = 10 m. For simplicity, the surcharge is neglected (q = 0 kPa).

Influence of the Friction Angle and Cohesion
Using the results calculated from seven tunnel stability models, Figure 7 represents the relationship between the normalized support pressure (p u /γD) and normalized cohesion (c/γD). The analysis of the normalized support pressure is conducted for ϕ = 20 • and 40 • . C/D ratio is assumed to be 2.  Figure 7 shows that for ϕ = 20 • and ϕ = 40 • , the slope of the line represents the solution of Krause [31] has higher value than the other solutions indicating that the value of N c has much more effect on the support pressure compared to the other approaches.

c=.D [-]
In addition, the solutions of Broere [8], Anagnostou and Kovari [7] and Anagnostou [11] have a close slope value. Therefore, the value of N c for these models has a very close effect on the support pressure.
It can be seen from Figure 7 that for ϕ = 20 • and c/γD ≥ 0 (cohesionless and frictionalcohesive soil), the solution of Anagnostou and Kovari [7] using the limit equilibrium method overestimates the value p u /γD with regard to the other approaches. However, at ϕ = 40 • and c/γD ≥ 0, the solution of Broere [8] is above the solution of Anagnostou and Kovari [7]. The alteration of Anagnostou and Kovari [7] results comparing to Broere [8] results can be attributed to two reasons; firstly, the simplified way of considering the vertical stress distribution along the sides of the wedge, i.e. linear vertical stress distribution. Secondly, for ϕ > 35 • Anagnostou and Kovari [7] assumed K wedge = 0.4 which is higher than the assumed value of Broere [8] (K wedge = K 0 ). Consequently, in the Anagnostou and Kovari [7] model, the vertical stress and the frictional resistance acting on the wedge are considerably higher than vertical stress and the frictional resistance acting on the wedge obtained using Broere [8] model. Finally, at ϕ = 40 • , the integration of the two previews reasons during the optimization process leads to reduce distinctly the calculated minimum support pressure for Anagnostou and Kovari [7] model comparing to Broere [8] model.
Moreover, for ϕ = 40 • , the values of p u /γD obtained from Vermeer et al. [21] are much more close to Mollon et al. [14] model, the difference of p u /γD between these two methods is less than 5 %. Whereas, for ϕ = 20 • , the values of p u /γD obtained from Vermeer et al. [21] and Mollon et al. [14] results are almost identical.
According to Figure 7, for cohesionless and frictional-cohesive soil with the two friction angles (ϕ = 20 • and 40 • ), the numerical results of Anagnostou [11] which are based on the infinitesimally thin slices method and the results of Mollon et al. [14] using multi-blocks mechanism are located between the results of Broere [8] and the upper bound solution of Leca and Dormieux [13].
For the two friction angles (ϕ = 20 • and 40 • ) and c/γD ≥ 0, the values of p u /γD obtained from Leca and Dormieux [13] using the upper bound method are clearly below the results given by the others, this is due to the shape of the tunnel face being considered as an elliptic cross section inscribed to the circular face. The two-blocks mechanism given by Leca and Dormieux [13] is constrained by the normality condition, required by plasticity theory, the velocity vector must make an angle ϕ with discontinuity surfaces along the sliding surfaces. However, this condition does not allow the three-dimensional slip surfaces to develop more freely. Figure 8 presents the influence of cover to diameter ratio (C/D) on normalized support pressure (p u /γD). From Figure 8, we can see that with any friction angle, the normalized support pressure (p u /γD) predicted from Leca and Dormieux [13] and Mollon et al. [14] solutions remain constant for C/D ≥ 1. The constant support pressure in Leca and Dormieux [13] and Mollon et al. [14] upper bound solutions explained by both failure mechanisms of Leca and Dormieux [13], and Mollon et al. [14] are translational kinematically admissible failure mechanisms, the 3D rigid conical blocks of these mechanisms move with velocities, which are collinear with the cones axes, making an angle ϕ with the conical discontinuity surfaces, and fulfilling the normality rule required by the limit analysis theory. Due to this condition, the failure mechanism obtained from the maximization process does not change with C/D ≥ 1 [14].

Influence of C/D
In Krause [31] solution, the support pressure is not influenced by the cover to diameter ratio (C/D), which is attributed to the Krause [31] model being based on the internal equilibrium of a half sphere failure mechanism at the front of the tunnel face, where the soil weight of the cover depth is not taken into account in calculating the support pressure. Similarly, in Vermeer et al. [21] numerical results, the normalized support pressure (p u /γD) with C/D ≥ 1 is independent of the cover to diameter ratio.
In limit equilibrium models, the vertical stress tends exponentially to asymptotic value at specific tunnel depth. When the vertical stress reaches its asymptotic value, the minimum support pressure is only dependent on the friction angle and not on the depth of the tunnel, see Equation (1). The rate of approach to the asymptote of the vertical stress is a function of K silo · tanϕ and therefore differs in Anagnostou and Kovari [7], Broere [8], and Anagnostou [11] solutions with different friction angles.
For the case of ϕ = 20 • , the normalized support pressure of Anagnostou and Kovari [7], Broere [8] and Anagnostou [11] solutions increase nonlinearly for increasing C/D and the increasing rate became flatter in deeper tunnel depth.
Moreover, for ϕ = 20 • and ϕ = 40 • , the normalized support pressure calculated by Anagnostou [11] model is slightly impacted by the depth ratio if it is less than 2. Once the cover to diameter ratio is bigger than 2, the support pressure gets constant and is no longer influenced by the cover to diameter ratio. We can also conclude from Figure 8, for ϕ = 40 • , in Anagnostou [11] solution, once the cover to diameter ratio is bigger than 2, the support pressure gets constant and is no longer influenced by the cover to diameter ratio. However, in Broere [8], the effect of C/D starts to vanish for C/D > 2.5.

Verification by Physical Model Tests
Three series of physical model tests from the literature have been chosen to verify the accuracy of selected models for predicting the minimum support pressure. The specifications of the physical model tests are summarized in Table 1.  [22] 1g test Ottendorf-Okrilla sand 32.5 0 0.5, 0.75, 1, 1.5, 2 Chen et al. [47] 1g test Yangtze River sand 37 0-0.5 0.5, 1, 2 The results of normalized support pressure obtained from Chambon and Corte [3], Kirsch [22] and Chen et al. [47] physical models are compared with the results of the wedgesilo models by Anagnostou and Kovari [7], Broere [8], Anagnostou [11], the upper bound solution by Mollon et al. [14], Leca and Dormieux [13], and the finite element method by Vermeer et al. [21]. The predicted results are shown in Figure 9a-f.
From Figure 9a we can see that the normalized support pressure obtained from 1 g tests of Kirsch [22] increases slightly with the increase in C/D ratio. This dependency of C/D ratio was also detected by Anagnostou and Kovari [7] and Broere [8] but is not revealed by Mollon et al. [14], Leca and Dormieux [13], Krause [31], Vermeer et al. [21] and Anagnostou [11]. Besides, the results obtained from Vermeer et al. [21], Mollon et al. [14] and Anagnostou [11] models agree very well with 1g tests of Kirsch [22] results, when the soil is cohesionless soil (c = 0 kPa). According to Figure 9b, we can also see that Leca and Dormieux [13], Mollon et al. [14] and Anagnostou [11] models show a good agreement with the results obtained by the 1g test results of Chen et al. [47] for friction-cohesive soil. As shown in Figure 9c,e, the normalized pressure (p u /γD) obtained from ng tests of Chambon and Corte [3] are less than the results calculated by the seven stability models in cohesionless soil. However, if the cohesion is considered in the calculation of the pressure (c = 5 kPa), see Figure 9d,f, Leca and Dormieux [13], Anagnostou and Kovari [7], Vermeer et al. [21], Mollon et al. [14] and Anagnostou [11] predict the minimum support pressure much more closer to the 1g test results of Chambon and Corte [3]. Moreover, it can be seen that Anagnostou and Kovari [7], Broere [8], and Krause [31] relatively overestimate the (p u /γD) value compared to the results of Chambon and Corte [3], Kirsch [22] and Chen et al. [47].

Layered Soil
In this section, the tunnel face stability in layered soils is investigated using different approaches; the model considers two layers (cover and cross layer) of the soil is shown in Figure 10. The cover and the cross layer are assumed to be located above the groundwater table. The soil strength parameters and geometry of the tunnel are described in Table 2.
The effect of the soil parameters of the cover layer on the minimum support pressure is investigated by varying the friction angle of the cover layer, while the soil properties of the cross layer are kept constant. The variation of the normalized minimum support pressure with the friction angle of the cover layer is shown in Figure 11.   It is clear from Figure 11 that the results obtained from the different approaches are quite different. It can be noticed that the wedge-silo models (Jancsecz and Steiner [6]; Anagnostou and Kovari [7]; Broere [8]) give much larger support pressures than the limit analysis models (Tang et al. [15]; Han et al. [44]). The results of Han et al. [44] are above the results of Tang et al. [15]. This can be attributed to the fact that the 3D failure mechanism of Han et al. [44] is composed of five truncated cones which offer many more degrees of freedom. However, the failure mechanism assumed by Tang et al. [15] is composed of two cones proposed by Leca and Dormieux [13], which is restricted in the number of degrees of freedom.
Furthermore, it can be seen from Figure 11 that the minimum support pressure decreases with increasing friction angle in the case of the upper bound solution of Han et al. [44], Tang et al. [15] and the Anagnostou and Kovari [7] model. The steepest trend is predicted by the model of Anagnostou and Kovari [7]. Surprisingly, applying the solution of Broere [8], the minimum support pressure decreases up to a minimum value at ϕ = 40 • before increasing again with further increase of ϕ. A similar trend is obtained from the approach of Jancsecz and Steiner [6]. In that case, the minimum support pressure decreases slightly to the minimum value at ϕ = 30 • , while it increases subsequently.
The trends of the support pressure predicted by the wedge-silo models (Jancsecz and Steiner [6]; Anagnostou and Kovari [7]; Broere [8]) can be explained as follows. From silo theory (Equation (1)), the value of vertical stress above the tunnel crown is mainly dependent on the value of K silo · tan ϕ (for the same soil strength parameters). For the model of Anagnostou and Kovari [7], this value will be 0.8 · tan ϕ, for the Broere [8] solution it will be (1-sin ϕ) · tan ϕ, and for the Jancsecz and Steiner [6] solution (tan 2 (45 − ϕ/2)) · tan ϕ. As shown in Figure 12, the value of K silo · tan ϕ in the solution of Anagnostou and Kovari [7] increases exponentially. According to the approach of Broere [8], the value of K silo · tan ϕ increases up to a maximum value of ϕ = 40 • before it decreases again. The same trend is obtained for the equations of Jancsecz and Steiner [6], where the maximum value is reached at ϕ = 30 • .
The effect of K silo · tan ϕ on the vertical stress distribution in the different model is shown in Figure 13. From those diagrams one it conclude that in the models of Jancsecz and Steiner [6] and Broere [8], the vertical stress for ϕ = 45 • is close or slightly higher than that for ϕ = 35 • at any depth. In contrast, in the approach of Anagnostou and Kovari [7], the vertical stress for ϕ = 35 • is always higher than that for ϕ = 45 • .  Based on the previous results, it can be concluded that for the set of soil parameters and geometry of the tunnel used in this study (see in Table 2) the approaches of Jancsecz and Steiner [6] and Broere [8] predict a trend of the minimum support pressure with friction angle which contradicts practical experience. The discrepancies between the support pressure predicted by the different approaches are obvious in Figure 11.

Conclusions
This paper presents a review of analytical models based either on the principle of limit equilibrium or the limiting analysis, which have both been used to analyze the stability of the tunnel face for dry frictional and frictional-cohesive soil. The following main conclusions can be drawn:

1.
To use Janssen's solution in the wedge-silo model, it is required to have explicit value for the lateral earth pressure coefficient of the silo. However, the advised values of it vary in a wide range. Due to this, the face support pressures calculated by the existing wedge-silo models are quite different. So, the question about the proper value for the lateral earth pressure K silo remains unanswered.

2.
According to the results of comparative calculations, the wedge-silo models give higher support pressure than conical block models, but the differences in the two methods for the prediction of the support pressure are small for higher friction angles.

3.
Comparing the results of physical model tests to the results of Anagnostou and Kovari [7] and Broere [8] models, indicating that both models are relatively conservative models for estimating the minimum support pressure. The results obtained from Vermeer et al. [21], Mollon et al. [14], and Anagnostou [11] models show a good agreement with those of the experiments. However, comparing to the experiments results, the upper bound solutions proposed by Leca and Dormieux [13] underestimate the minimum support pressure especially when soil is cohesionless. 4.
The solution of the limit equilibrium method revealed that the cover to diameter ratio has significant influence on the minimum support pressure. This effect is also deduced by Chen et al. [4] and Chambon and Corte [3], whereas the solution of the upper bound theorem for the minimum support pressure is only independent, when the failure mechanism of the upper bound theorem does not intersect at the ground surface in a certain soil condition (e.g., C/D ≥ 1), this outcome has also been reported by Vermeer et al. [21]. 5.
In the case of a ground composed of two soil layers, the minimum support pressure obtained from wedge-silo models is higher than that predicted by the upper bound solution. In addition, the discrepancies in the results between the support pressure predicted by the different models are obvious. 6.
Finally, the analytical approaches such as the Limit Equilibrium Method (LEM) and Limiting Analysis Method (LAM) are used to assess the stability of the tunnel face assuming various failure mechanisms. However, the results are quite different. Therefore, there are still considerable potential efforts for calculating the support pressure more accurately. Funding: This research received no external funding.
Data Availability Statement: Not applicable.

Conflicts of Interest:
The authors declare no conflict of interest.