Determination of Safety Factor for Rock Mass Surrounding Tunnel by Sudden Change of Equivalent Plastic Strain in Strength Reduction Method

The safety factor of rock mass surrounding tunnel can be determined using the strength reduction method (SRM), furthermore, it is very important to solve the criterion of critical state. For the stability estimation of rock mass surrounding tunnel, there is need to discuss that it is preferable to use the same criteria for the slope stability, such as non-convergence of finite element calculation, penetration of plastic strain and sudden change of horizontal displacement. The safety factor can be determined by sudden change of equivalent plastic strain in relationship between a reduction coefficient of strength parameter and equivalent plastic strain. This method is based on the elasto-plastic finite element method (FEM) and the SRM by ABAQUS and Mohr–Coulomb yield criterion. Simulation results using this method show how a safety factor varies with geometries, friction angles and cohesions for circle and square tunnels. Simulation results also show the safety factor varying with quality change of rock mass, pore water pressure and tunnel depth.


Introduction
One of the most important problems for structural calculation of tunnel is how to estimate the stability of rock mass surrounding the tunnel.
The stability of circular tunnels has been extensively studied at the Cambridge since the 1970s: see, for example, the work reported by Atkinson and Cairncross (1973), Davis et al. (1980), Mair (1979) and Seneviratne (1979). Before 1990s, most researches on the tunnel stability focused on the undrained stability of circular tunnels in clay. Davis et al. (1980) derived upper and lower bound stability solutions for the collapse under undrained conditions in consideration of three different shapes of shallow underground opening relevant to tunneling. Assadi and Sloan (1991) investigated the undrained stability of a square tunnel in a cohesive soil under plane Abstract The safety factor of rock mass surrounding tunnel can be determined using the strength reduction method (SRM), furthermore, it is very important to solve the criterion of critical state. For the stability estimation of rock mass surrounding tunnel, there is need to discuss that it is preferable to use the same criteria for the slope stability, such as non-convergence of finite element calculation, penetration of plastic strain and sudden change of horizontal displacement. The safety factor can be determined by sudden change of equivalent plastic strain in relationship between a reduction coefficient of strength parameter and equivalent plastic strain. This method is based on the elasto-plastic finite element method (FEM) and the SRM by ABAQUS and strain. Yamamoto et al. (2011Yamamoto et al. ( , 2013 studied the stability of one and dual circular tunnels in cohesivefrictional soils subjected to surcharge loading. They described the ultimate surcharge load of cohesivefrictional soils with inclusion of shallow square tunnels in terms of the dimensionless load parameter σ s /c′ as a function of three variables (φ′ γB/c′ and H/B). But there are so many factors affecting the stability of a tunnel, which implies that there is no unique criterion of the stability estimation (Lee et al. 2003).
To estimate the stability of rock mass surrounding tunnel, researchers have studied how to apply the strength reduction method (SRM). The SRM was first proposed by Zienkiewicz et al. (1975), and now, it has been widely used in the estimation of the slope stability. The definition of the safety factor for the slope using the SRM is considered to be a ratio of actual shear strength to the lowest shear strength of geomaterials (soil and rock) required to maintain the equilibrium of a slope (Dawson et al. 1999;Ugai 1989;Griffiths and Lane 1999;Liu et al. 2015).
Many researchers have compared the various criteria of the critical state of a slope, which is the most important in the estimation of slope stability using the SRM (Fu and Liao 2010;Khosravi and Khabbazian 2012;Ri et al. 2022;Yuan et al. 2020;Zheng et al. 2009a, b). Various criteria were developed, however, there were mainly three methods by the penetration of equivalent plastic strain, non-convergence of numerical calculation and sudden change of horizontal displacement at a specific point, respectively. Several researchers have sought for alternatives by establishing correlation between the assessments of tunnel and slope stability. Huang et al. (2012) determined the safety factor for a shallow tunnel in combination with the upper bound theorem of limit analysis and SRM. Jiang et al. (2009) judged critical failure state of tunnel by the penetration of equivalent plastic strain. Goh and Zhang (2012) and Yin et al. (2007) considered the non-convergence to be the criterion of the critical state, which produces the infinite plastic strain on the slide surface, and the finite element calculation is not converged, when tunnel lost the stability and collapsed. The earthquake stability for rock and soil masses surrounding tunnel was estimated by ANSYS program and SRM (Cheng et al. 2014). An et al. (2011) and Han et al (2021) proposed that vault settlement rate or ground subsidence should be a major criterion for assessment of tunnel stability. Zhang (2009) suggested that it was difficult to confirm the instability for rock mass surrounding tunnel by the distribution of plastic zone or the non-convergence of the finite element calculation. And he proposed the sign of collapse in rock mass, indicating that the slide move infinitely and the displacement of the slide surface change and develop infinitely. Liang (2019) analyzed the failure mechanism around deep excavations using rock failure process analysis (RFPA) with SRM.
There are several advantages for the application of the SRM in geotechnical stability analysis, and the SRM is very attractive and commonly accepted approach among the geotechnical researchers and engineers (Huang et al. 2012;Shen and Karakus 2013). Presently, the safety factors used in the geotechnical stability analysis are mostly calculated by the finite element method (FEM) with shear strength reduction technique. However, the key for the estimation of tunnel stability using the SRM is how to judge the instability and how to determine the safety factor. Therefore, there is need to consider whether the criteria for slope can be also used in tunnel and need to develop a more suitable criterion for determining the safety factor of rock mass surrounding tunnel by the SRM.
In this paper, various criteria such as the nonconvergence of numerical calculation, sudden changes of horizontal displacement and plastic strain are compared and analyzed, respectively. Then, we present a new criterion for determining the safety factor of rock mass surrounding tunnel by a sudden change of equivalent plastic strain on the relationship between the shear strength reduction coefficient (SRC) and the equivalent plastic strain (PEEQ) at a specific point. Simulation results using the proposed method show how the safety factor varies with geometries, friction angles and cohesions for circle and square tunnels, respectively, and how the safety factor varies with quality change of rock mass, pore water pressure and tunnel depth, respectively.

Determination of the Safety Factor
There are mainly two well-accepted definitions of the safety factor for rock or soil mass (Zheng et al. 2006). One is the safety factor as strength reserve, i.e., when the strength of rock or soil mass is continuously reduced until it fails, a reduced ratio of the strength is just the safety factor, while the other is the safety factor as load reserve. This method increases continuously the gravity of rock or soil mass until it fails, i.e., an increased ratio of gravity is just the safety factor. For tunnel, a failure typically occurs due to the strength reduction of rock mass caused by drill and blast excavation method, therefore, the safety factor as strength reduction is more suitable for tunnel. For shear failure mode, failure mode of a tunnel is similar to that of a slope, and the the safety factor is a ratio of actual strength to yielding strength of rock or soil mass on shear failure surface. If shear strength parameters are continuously reduced using finite-difference or finite-element SRM until rock or soil mass reaches critical failure state, the reduction coefficient corresponding to critical sate is just the whole safety factor of rock mass surrounding tunnel.
Stability estimation proposed in this section means a judgment of critical state and a determination of the safety factor, which is based on the Mohr-Coulomb strength criterion. The Mohr-Coulomb strength criterion is commonly used in a calculation of rock or soil mass. The safety factor of rock mass surrounding a tunnel can be defined as a ratio of actual shear strength parameters to critical failure shear strength parameters. Here, the stability estimation is performed with analysis results of elasto-plastic FEM and the SRM using ABAQUS.

Simulation of Strength Reduction Process in ABAQUS
The shear strength reduction technique can be adopted to solve for a safety factor of tunnel stability. In this technique, a series of trial safety factors are used to adjust the cohesion c and the friction angle φ of geomaterials, respectively, as Eq. (1). Surrounding rock mass can be idealized as an elastic, perfectly plastic material with cohesion and frictional angle. Shear strength parameters of rock mass surrounding tunnel are continuously reduced, and so, the rock mass turns into critical state and the calculation stops. From the simulation results, the plastic zone contribution and the factor of strength reserve can be obtained, and the development process of equivalent plastic strain can be monitored. In ABAQUS, the shear strength parameters of geomaterial can be reduced by using field variable F v , which equals to an SRC. The simulation of the strength reduction process in ABAQUS is as follows: In the strength reduction calculation by ABAQUS, range of the SRC (F v ) must be defined so that the shear strength parameters can be reduced. Commonly, this range is equal to one of the SRCs. If the shear strength parameters of material are too small, non-convergence may occur. Therefore, a small value of F v should be chosen to avoid non-convergence at the beginning of calculation, so that the strength parameters can be changed into the large value. According to our calculation experiences, F v may have the range of approximately 0.25 ~ 4.0 or 1.0 ~ 10.0, that is, the minimum and maximum values (F v min and F v max ) of the SRC may be chosen approximately as 0.25-1.0 and 4.0-10.0, respectively. If rock mass surrounding tunnel is soft, generally shear strength parameters are small. Consequently, at the begining of the strength reduction, the rock mass can lose the stability, and the calculation of elastoplastic FEM for shear strength reduction can be also diverged. Therefore, it is more reasonable to choose the minimum value (F v min ) of the SRC as an enough small one than 1.0 and start in calculation. For the hard rock mass surrounding tunnel, the minimum value (F v min ) of 1.0 may be chosen, while the maximum value (F v max ) may be chosen as large value (e.g., 4.0-10.0). Only under such conditions, numerical calculation can be sufficiently performed without divergence, and then, total process of plastic strain development can be sufficiently monitored.
According to the SRC (F v ), change of the shear strength parameters for geomaterial such as cohesion c, inner frictional angle φ and shear dilatation angle ψ can be separately defined. These parameters are related to the SRC (F v ) by the following equation. (1) where c′, φ′ and ψ′ is reduced cohesion, frictional angle and shear dilatation angle for geomaterial, respectively. The gravity (or volume load) must be applied, and the equilibrium stress state must be created for the model. The SRC (F v ) should linearly increase up to F v max in the later analysis. When the elasto-plasticity finite element calculation is ceased (i.e., the divergence in this calculation), the results with total process of the calculation should be monitored so that a critical state can be judged, and a safety factor can be determined.

Determination of Critical Failure State and Safety Factor
The stability estimation of rock mass surrounding tunnel is a process which judges a failure state of rock mass, and at the same time, determines safety factor of strength reserve. Critical failure state can be judged from the observation of development process of equivalent plastic strain (PEEQ). In other words, the moment when plastic strain starts to generate, is just the critical failure state. Next, it is also necessary to search for a sudden change point on the relative curve between F v and PEEQ so that the safety factor can be determined. This phenomenon of sudden change conforms to the appearance moment of equivalent plastic strain in the excavation part of tunnel. Finally, it is suggested that the SRC (F v ) corresponding to the sudden change point can be determined as the safety factor ( Fig. 1). When the sudden change point on the curve is considered, attention is to confirm a relation between F v and PEEQ corresponding to the sudden change point through numerical value. That is, it is necessary to make sure that the beginning point of plastic strain value occurring from zero coincides with the sudden change point. The above process is described in detail through numerical simulation for examples of Sect. 3.

Numerical Simulation for Examples
Examples of tunnel models using ABAQUS are simulated to calculate the safety factor based on the above-described methodology. The tunnels are modeled as unlined circular tunnel and square tunnel, respectively, which have no internal pressure. The deformation is assumed to be the plane strain state. Besides, the surrounding rock mass is modeled as a uniform Mohr-Coulomb material with a unit weight γ, cohesion c, inner friction angle φ, shear dilatation angle ψ, elastic modulus E and Poisson's ratio μ, respectively (Figs. 2 and 3). The surrounding rock mass is modeled using the base feature of 2D Planar shell. In case of no considering the pore water pressure, the element type of CPE6 (i.e., 6-node modified quadratic plane strain triangle) is used, while, for considering the pore water pressure, the element type of CPE6M (i.e., 6-node modified quadratic plane strain triangle, pore pressure, hourglass control) is used. The range of field variable F v is defined as from 0.25 to 4.0.  Figure 3 shows the finite element mesh along with the assumed boundary conditions. Considering vertical symmetry, only half tunnel is used for numerical analysis. The edge along the center axis of tunnel is fixed in the horizontal direction, i.e., horizontal displacement u is zero (u = 0), while the outside boundary part of circle arc is fixed both in the horizontal and vertical directions, i.e., horizontal displacement u and vertical one v are all zero (u = v = 0). And as for the load condition, it is assumed that the weight load of the ground is applied, while both outer and inner pore pressure are neglected. The boundary of tunnel circumference is free. Physical and mechanical parameters of material for the circle tunnel geometries of H/D = 2.0 and D = 5.0 m are shown as Table 1. Analysis procedures are divided into 3 steps, consisting of creation step of initial stress, removal one of excavation elements and reduction one of the shear strength parameters, respectivley. The calculation results in the visualization module are considered after finishing the finite element calculation.
This numerical example converged up to the increment step 100 of analysis step 3 (i.e., reduction step of the shear strength parameters) and terminated the calculation of the plastic strain increments at this time. The equivalent plastic strain isn't occurred in the increment steps 0 to 43, however, it begins to occur from the increment step 44. The equivalent plastic strain contours are shown at the increment step 54, the increment step 64 and the increment step 74, respectively. The distribution region of equivalent plastic strain is gradually extended, and Fig. 4 shows the equivalent plastic strain at the last increment step 100. As a result, at the start, the equivalent plastic strain occurs from the side wall of excavation part, while it increases gradually with the increase of the reduction coefficient. Because the development state of equivalent plastic strain varies with the position on the surface of the excavated tunnel, the specific 5 points can be chosen at the top, the bottom and the sides of excavation part of tunnel, respectively, as shown in Fig. 5. For these 5 points, there are different relationships between the SRC (F v ) and the equivalent plastic strain (PEEQ) or between the SRC (F v ) and the horizontal displacement (U 1 ). These curves are plotted as shown in Fig. 6a-e.
There are several remarkable results by considering above-calculated data of finite element SRM for circle tunnel.
(1) It is difficult to determine the safety factor using non-convergence of finite element calculation, because all the calculations finished at the SRC F v of 4, and at this time, all the safety factors are determined to be 4 depending on the non-convergence of finite element calculation. Therefore, this can cause a very big error for determining the safety factor. (2) It is difficult to determine the safety factor using the penetration of plastic strain. In other words, plastic strain value is very small, and potential  failure surface (i.e., plastic strain penetration zone) happens in soil slope but that doesn't happen in rock mass surrounding tunnel.
(3) It is impossible to determine the safety factor using the sudden change of horizontal displacement (U 1 ), because the sudden change point can't be found in all the relationships between F v and U 1 for these 5 points. Therefore, the above 3 reasons indicate that it is difficult to determine using the preceding criteria for the slope such as nonconvergence of finite element calculation, penetration of plastic strain and sudden change of horizontal displacement, respectively. (4) The sudden change point can be found in the relationships between F v and PEEQ, therefore, the safety factor can be determined by sudden change of equivalent plastic strain in relationship between a reduction coefficient of strength parameter and equivalent plastic strain.
However, the position of the sudden change point varied with the position of specific point as shown in Fig. 6, while there is the most obvious sudden change phenomenon in the equivalent plastic strain curve for the point 3 as shown in Fig. 6c. While the SRC F v increases more than 2.005 and the increment step increases more than 43, the value of equivalent plastic strain keeps almost linear increase till the calculation stops, and the distribution region of equivalent plastic strain also expands. In this case, F v is 2.005 as the minimum among the reduction coefficients corresponding the sudden change points. Finally, the safety factor can be determined as about 2.0.  Vol.: (0123456789) Fig. 6 Relation curves between the SRC (F v ) and the equivalent plastic strain (PEEQ), between the SRC (F v ) and the horizontal displacement (U 1 ) on the specific 5 points of a circle tunnel Then, on the basis of the study above, a compact set of safety factor charts can be calculated and presented for circle tunnels. Parameters considered in this section are for circle tunnel geometries of H/D = 1.0-5.0, D = 5.0 m, friction angles of φ = 5-35° and γD/c = 0.5-3.0. Table 2 shows the comparison of the safety factor solutions obtained from the sudden change of equivalent plastic strain for the calculation arrange of the above. As shown in Table 2, the change of the safety factor can be obtained with geometries H/D, friction angles φ and cohesions c, respectively, for circle tunnels. Also, it can be seen that the larger H/D is, the smaller safety factor becomes, while, for the same H/D, the smaller friction angle and cohesion are, the smaller the safety factors become.

Determination of the Safety Factor for Square Tunnel
Then the safety factor for square tunnel can be determined. Figures 7 and 8 show the finite element mesh along with the assumed boundary conditions for half of square tunnel. The geometric model and boundary conditions are equal to those of circle tunnel. Likewise, physical and mechanical parameters of material and analysis procedures for square tunnel of H/B = 2.0, B = 5.0 m are defined as circle tunnel. The calculation results in the visualization module are considered after finishing the finite element calculation. The calculation of the plastic strain increments in this example terminated when calculation step converged up to the increment step 45 of analysis step 3 (i.e., reduction step of the shear strength parameters). The equivalent plastic strain isn't occurred in the increment step 0-25. However, the equivalent plastic strain begins to occur from the increment step 26. The equivalent plastic strain contours are shown at the increment step 36, the increment step 40 and the increment step 43, respectively. The distribution region of the equivalent plastic strain is gradually extended, and it is shown that at the last increment step 45, the equivalent plastic strain is contoured as in Fig. 9. As a result, at the start, the equivalent plastic strain occurs from the corners of top and bottom of excavation part while increases gradually with the increase of the reduction coefficient. In the same way as described for circle tunnel, the specific 5 points can be chosen at the top, the bottom and the sides of excavation part of square tunnel, respectively, as shown in Fig. 10. For these 5 points, there are different relationships between the SRC (F v ) and the equivalent plastic strain (PEEQ) or between the SRC (F v ) and the horizontal displacement (U 1 ). These curves are plotted as shown in Fig. 11a-e.
Like for circle tunnel, there are several remarkable results through the considering above calculation data of finite-element SRM for square tunnel.
(1) It is difficult to determine the safety factor using non-convergence of finite element calculation. The calculations for point 1 and point 5 finished at the reduction coefficient F v of 4, and at this time, the safety factors are determined to be 4 depending on the non-convergence of finite element calculation, while the calculations for point 2, point 3 and point 4 finished at the reduction coefficient F v of 1.80, and at this time, the safety factors are determined to be 1.80 depending on the non-convergence of finite element calculation. But this can cause a very big error during the determination of the safety factor. (2) It is difficult to determine the safety factor using the penetration of plastic strain. In other words, plastic strain value is very small, and potential failure surface (i.e., plastic strain penetration zone) happens in soil slope but that doesn't happen in surrounding rock mass. (3) It can't determine the safety factor using the sudden of horizontal displacement (U 1 ) because the sudden change point can't be found in all the rela-tionships between F v and U 1 for these 5 points. Therefore, the above 3 reasons indicate that it is difficult to determine using the preceding criteria for the slope such as non-convergence of finite element calculation, penetration of plastic strain and sudden change of horizontal displacement, respectively. (4) The sudden change point can be found in the relationships between F v and PEEQ, therefore the safety factor can be determined by sudden change of equivalent plastic strain in relationship between a reduction coefficient of strength parameter and equivalent plastic strain.
However, the position of the sudden change point varied with the position of specific point, as shown in Fig. 11, while there is the obvious sudden change phenomenon in the equivalent plastic strain curve for the point 2, point 3 and point 4, respectively, as shown    Fig. 11b-d. Then, the F v of the sudden change point for point 2 is the minimum as 1.1875 among the reduction coefficients corresponding the sudden change points for points 2, 3 and 4, respectively, and the increment step of this F v is 25. While the reduction coefficient F v increases more than 1.1875 and the increment step increases more than 25, the value of equivalent plastic strain continues to increase almost linearly till the calculation stops, and the distribution region of equivalent plastic strain also expands. Finally, the safety factor can be defined as about 1.19.
Then, on the basis of the study above, a compact set of safety factor charts can be calculated and presented for square tunnel. Parameters considered in this section are for square tunnel geometries of H/B = 1.0-5.0, B = 5.0 m, friction angles of φ = 5-35°, and γB/c = 0.5-3.0, respectively. Table 3 shows the comparison of the safety factors obtained using the sudden change of equivalent plastic strain within the above parameter ranges. It can be found from Table 3 shows that the larger H/B is, the smaller safety factor becomes, while, for the same H/B, the smaller friction angle and cohesion are, the smaller the safety factors become.

Change of the Safety Factor with Various Conditions of Rock Mass
There are many factors affecting stability of surrounding rock mass. However, the most important factors Fig. 10 Selection of specific 5 points in excavation part of circle tunnel Fig. 11 Relation curves between the SRC (F v ) and the equivalent plastic strain (PEEQ), between the SRC (F v ) and the horizontal displacement (U1) on the specific 5 points of a square tunnel are geostructure and property within surrounding rock mass, groundwater and in-situ stress, respectively. The safety factor of surrounding rock mass is the changeable indices depending on the above conditions of surrounding rock mass. If the safety factor is obtained, it will be simple to define excavation limit and lining thickness of tunnel. The larger safety factor is, the larger excavation space will be and the thinner lining thickness will be. That is because the safety factor relates to geological conditions, tunnel depth and seepage action and so on. Namely, this means that the safety factor of surrounding rock mass is an important index reflecting the condition of surrounding rock mass. Therefore, in this section, it is useful to simulate how the safety factor varies with the quality class of the surrounding rock mass, pore water pressure and tunnel depth, respectively. The quality of surrounding rock mass is divided into degree I to degree V, and the calculation cases are divided into nine. The surrounding rock mass for circle tunnel of diameter D = 6.6 m is modeled as homogeneous Mohr-Coulomb material. The property parameters of surrounding rock mass are shown in Table 4, and action mode of load is shown in Fig. 12. As shown in Fig. 12, p0 increases linearly with the water depth (H w ) as pore water pressure activating surrounding rock mass, while σ H increases linearly with the ground depth (H) as in-situ stress of surrounding rock mass. First of all, it is useful to simulate how the safety factors vary with the quality class of surrounding rock mass for the tunnel with depth of H = 30, 65, 100, 200, 300 and 400 m, respectively. The values of the safety factor in each case are shown in Fig. 13, respectively.
As shown in Fig. 13, the safety factor has a close relation with the quality class of surrounding rock mass and the tunnel depth. The lower the quality class of rock mass is and the higher the tunnel depth is, smaller it becomes. That is, the softer rock mass is and the deeper tunnel depth is, the smaller the safety factor becomes. As shown in Fig. 13, the safety factor is reduced from above 1.0 to below 1.0 with the decrease of the quality class of rock mass. If the tunnel is shallow, the safety factor is rapidly reduced, however, if the tunnel is deep, the safety factor is slowly reduced. To put it concretely, if the tunnel depth is 30 m, the safety factor is rapidly reduced from 7.07 to 0.77, while if the tunnel depth is 65 m, the safety factor is rapidly reduced from 4.01 to 0.51. If the tunnel depth is 100 m, the safety factor is reduced from 2.97 to 0.4, while if the tunnel depth is 200 m, the safety factor is reduced from 1.98 to 0.31. If the tunnel depth is 300 m, the safety factor is reduced from 1.62 to 0.25, while if the tunnel depth is 400 m, the safety factor is reduced from 1.40 to 0.25.
With the increase of the tunnel depth, the safety factor becomes smaller and smaller. If the quality class of rock mass is high (i.e., hard rock mass), the safety factor is rapidly reduced, however, if the quality class of rock mass is low (i.e., soft rock mass), the safety factor is slowly reduced. To put it concretely, when the tunnel depth is increased from 30 to 400 m in the calculation case 1 (i.e., very hard rock mass), the safety factor is rapidly reduced from 7.07 to 1.4, and in the calculation case 2 (i.e., hard rock mass), the safety factor is rapidly reduced from 6.29 to 1.31. Also, in the calculation case 3, the safety factor is rapidly reduced from 5.57 to 1.19, while in the calculation case 4, the safety factor is rapidly reduced from 5.45 to 1.12. However, in the calculation case 7, all the safety factors are below 3.0, and the safety factor is slowly reduced from 2.21 to 0.54, while in the calculation case 8, all the safety factors are below 2.0, and the safety factor is slowly reduced from 1.91 to 0.42. And in the calculation case 9 (i.e., very soft rock mass), all the safety factors are smaller than 1.0, and the safety factor is slowly reduced from 0.77 to 0.25. In all the calculation cases, when the tunnel depth is smaller than 100 m, the safety factors show a marked trend toward reduction.
Then, it is useful to simulate how the safety factors vary with the outside pore water pressure and the class of quality of surrounding rock mass for the tunnel of depth with H = 30, 65 and 100 m, respectively. For action heads of H w = 0, 60, 80 and 100 m, respectively, the pore water pressure of 0, 588, 784 and 980 kPa are activated on the bottom of model, respectively. Action head is the value by which pore water pressure activated on the bottom is converted into water head. At that time, the pore water pressure activities as triangle distribution load. Figure 14 shows how the safety factors vary with the action head of H w = 0, 60, 80 and 100 m, respectively, for the calculation case of 1-9 and the tunnel depth of H = 30 m. Figure 15 shows how the safety factors vary with the action head of H w = 0, 60, 80 and 100 m, respectively, for the calculation case of 1-9 and the tunnel depth of H = 65 m. As in Figs. 14 and 15, 16 illustrates the curves by which the safety factors are changed for the tunnel depth with H = 100 m, and in general, the safety factors are smaller than the cases of H = 30 and 65 m, respectively. As shown in Figs. 14, 15, 16, the safety factors relate not only to the quality class of rock mass and tunnel depth but also to the pore water pressure. Therefore, the higher pore water pressure is, the smaller safety factor becomes. When the pore water pressure increases, the safety factor reduces, while the higher the quality class of rock mass is and the shallower tunnel depth is, the more obvious this trend is. For example, when the tunnel depth is 30 m, it can be seen from Fig. 14 that, higher the quality class of rock mass is, the trend of the safety factor to being reduced is obvious with the increase of pore water pressure, however, if the quality class of rock mass is lower, this trend is vague. On the other hand, as shown in Figs. 15 and 16, the higher pore water pressure is, the vaguer the reducing trend

Discussion and conclusion
The major focus of this paper was to present a new method for determining the safety factor of rock mass surrounding tunnel. First, various criteria used in slope stability analysis, such as the non-convergence of finite element calculation, penetration of plastic strain and sudden change of horizontal displacement, were compared and analyzed, respectively. Thus, it was found that the methods for estimating slope stability couldn't be used in tunnel stability.
Second, the sudden change of equivalent plastic strain in the SRM is used to estimate the stability of rock mass surrounding tunnel. For circle tunnel, the sudden change phenomenon of equivalent plastic strain was very obvious on side wall of excavation part for any circle tunnels, and the minimum value of the reduction coefficient was found at the sudden change point. Namely, the best solution for the safety factor can be obtained from the plastic strain curve for the specific point on side wall of circle tunnel. While for square tunnel, the sudden change phenomenon of equivalent plastic strain was quite obvious on the corners of top and bottom of excavation part for any square tunnels, and then, the minimum value of the reduction coefficient was found at the sudden change point as well. This method is relatively effective and simple for judging the critical state and for determining the safety factor of the rock mass surrounding tunnel and can overcome the weakness of previous criteria.
Third, on the basis of this method, a compact set of the safety factor charts is proposed for circle and square tunnels. Tables 2 and 3 show that the safety factors are varied with geometries (H/D or H/B), friction angles (φ) and cohesions (c) for the circle and square tunnels. In the other word, Tables 2 and 3, respectively, show that the larger H/D or H/B is, the smaller safety factor becomes while for the same H/D or H/B the smaller friction angle and cohesion are, the smaller the safety factors become.
Fourth, it is important to understand the implication of the safety factors shown in Figs. 14, 15, 16. These figures imply that the safety factor has a close relativity with the quality class of surrounding rock mass, the tunnel depth and pore pressure. As shown in these figures, the deeper the tunnel, the softer the surrounding rock mass and the higher pore pressure is, the smaller safety factors become. The safety factors under every calculation condition can be classified smaller than 1.0, 1.5 or larger than 2.0 using the above-described calculation results. It is interesting to note that the safety factors can be obtained depending on the quality class of surrounding rock mass, the pore water pressure and tunnel depth, respectively. This study will be exploited in any rock mass surrounding tunnel where the stability estimation is needed, and it seems that this is a new prospecting standard for judging the critical failure state. It is important to carefully monitor plastic strain behavior of the surrounding rock mass in different parts of the tunnel. On the basis of the results, it is necessary to estimate not only the local instability but the total instability induced by the local instability. It appears that this is However, since compared calculation models are not sufficient, these results were not examined for all cases and may be influenced by some external factors. Our research may have a limitation. The stability estimation is limited for only surrounding rock mass of homogeneous material. Further studies are needed in order to judge more rationally critical failure state for the stability estimation of tunnel or surrounding rock mass with any structure type and any rock mass condition. The research for estimation of tunnel stability has direct influence on the investment of tunnel construction and is of great significance in protecting environment or decreasing environmental pollution.