Installation Time of an Initial Support for Tunnel Excavation upon the Safety Factors of Surrounding Rock

: In engineering practice, the initial support system is commonly installed in the vicinity of the tunnel face after excavation, whereas the self-capacity of rock mass will fail to be utilized and the cost of the initial support system will be expensive. In this study, a methodology is proposed to determine the appropriate timing of initial support installation to ﬁnd out the balance of tunnel safety and construction cost. Firstly, the global safety factor is introduced as the critical indicator to evaluate tunnel stability. Then, the comprehensive graphic relationship between the global safety factor and the distance to the tunnel face is established. Once the global safety factor decreases to an admissible value, the stability of the surrounding rock is in a critical state and the corresponding distance is the recommended location for installing the initial support. In these procedures, the installation time of the initial support at the typical tunnel section can be quickly designed and fed back by a direct indicator during construction. Meanwhile, several cases with di ﬀ erent conditions have been carried out to discuss the regularity of the method.


Introduction
As the important components in tunneling, the support system serves to stabilize the tunnel and provide safe working conditions. The support system can be classified as the initial support structure and the secondary lining by the sequence of installation. The initial supports, such as the steel arch and cable that are installed firstly after excavation, will share most of the loads to control the deformation of the surrounding rock [1]. Therefore, it is valuable to optimize the design of the initial support system.
There are four different approaches for the design of the initial support: (a) empirical method; (b) analytical method; (c) observational method; (d) numerical method. The ground support is commonly selected by the empirical method, such as the rock mass classification systems that are firstly formulated by Terzaghi [2]. Subsequently, numerous rock classification systems are proposed, including the Q-system and Rock Mass Rating (RMR) system [3,4]. These methods have been proven effective in the preliminary stage and widely implemented in practical projects. However, these methods rely on a large number of historic cases and fail to reflect the mechanical behavior and deformation of rock and support structures [5,6]. Then, in the early 1960s, the new Austrian tunneling method (NATM) by Rabecewicz firstly considers the self-capacity of the rock mass as a part of support structure [7].
In summary, the majority of the initial support design methods are aimed to choose appropriate support types and parameters, while the installation time of the initial support is commonly designed by semi-empirical methods. According to the modern support design theory, the early installation of support structures will lead to excessive load on the support structures and fail to make full use of the rock capacity, which is uneconomical. Contrarily, the delayed installation of Figure 1. Schematic representation of the convergence-confinement method [8]. Reproduced with permission from Carranza-Torres C., Tunnelling and Underground Space Technology; published by Elsevier, 2000. In summary, the majority of the initial support design methods are aimed to choose appropriate support types and parameters, while the installation time of the initial support is commonly designed by semi-empirical methods. According to the modern support design theory, the early installation of support structures will lead to excessive load on the support structures and fail to make full use of the rock capacity, which is uneconomical. Contrarily, the delayed installation of support structures will lead to excessive deformation and threaten the safety of the tunnel. Consequently, it is necessary to determine the appropriate installation time of the initial support for the reliability and economy. Theoretically, when the self-bearing capacity of rock mass has been exhausted after excavation, the stability arrives at a critical state and the support structures are supposed to carry on the tunnel. In the engineering, the estimation of the critical state normally depends on the experience with the deformation of surrounding rock or the depth of plastic zone. For instance, Zhang [23] combined the displacement of the monitoring point with the principal stress direction by analyzing three-dimensional stress rotation after excavation as a standard to select the installation time of liner. Su et al. [24] proposed a method to choose optimal installation time of initial ground support by the plastic ratio upon the evolution of the displacement completion rate. Moreover, combining the observational method and the inverse analysis, a methodology for the real-time adaptation of the support systems during tunneling is presented by Miranda et al. [25] by the limit criteria for displacements and stresses. However, in the field monitoring, the limited monitoring points cannot always be representative for the whole section. From the perspective of mechanics, after excavation, one point in rock mass will lead to yielded state when the stress at this point exceeds the strength of rock mass. And restrained by the unyielding part, the intact rock will not be totally damaged. Consequently, it is acknowledged that the instability of rock mass is due to the development of the yield region and continuous slip surface. Among the various criteria for the stability of rock mass, the safety factor method based on the elastic-plastic theory can reflect the ratio between the stress state and the yield strength, which has been widely implemented in slope engineering [26]. Moreover, numerous works have considered the safety factor method in tunnel construction. Li et al. [27] calculated and compared the safety factor based on the Mohr-Coulomb and Drucker-Prager rule after tunnel excavation to judge the stable range of rock mass for designing the reasonable supporting parameter. By small-scale model tests and numerical analysis, Lee [28] utilized the safety factor of the tunnel face stability as an indicator to determine the round length that was defined as the length of unsupported span after excavation.
In this paper, the global safety factor is defined to estimate the stability of the tunnel and a methodology to suggest the installation location of the initial support. Meanwhile, several schemes with diverse lateral pressure coefficients and geological conditions are carried out to discuss the regularity of this method.

Point Safety Factor Based on Mohr-Coulomb Failure Criterion
When nonlinear deformation begins, the stress at this point is defined as yield strength. Yield criterion is a hypothesis concerning the limit of elasticity under any combination of stresses. The Mohr-Coulomb failure criterion is widely used to describe the response of geotechnical materials to shear stress as well as normal stress [29]. The Mohr-Coulomb yield surface divides the principal stress space into two individual areas, elastic zone and plastic zone. And the Mohr-Coulomb yield function can be expressed by principal stresses as Equation (1): where f is the function of the yield criterion; σ 1 and σ 3 are the maximum principal stress and the minimum principal stress, respectively; c is the cohesion and φ is the angle of internal friction. As mentioned above, the point will yield when the stress state at this point exceeds the material strength. To indicate the stability of rock mass, the point safety factor is introduced by an engineer on the basis of elastic theory as Equation (2).
where H is the parameters' function for the internal variable, χ, and F is the point safety factor. F > 1 indicates that the zone is in the non-yielded state; F = 1 indicates that the zone is in critical state; F < 1 indicates that the zone is under the yielded state. Based on the Mohr-Coulomb failure criterion, the distance BC between the stress state and the strength envelope in Figure 2 can be considered as the strength margin. Accordingly, the point safety factor based on Mohr-Coulomb criterion can be defined as Equation (3): Appl. Sci. 2020, 10, x FOR PEER REVIEW 4 of 19 where H is the parameters' function for the internal variable, χ, and F is the point safety factor. F > 1 indicates that the zone is in the non-yielded state; F = 1 indicates that the zone is in critical state; F < 1 indicates that the zone is under the yielded state. Based on the Mohr-Coulomb failure criterion, the distance BC between the stress state and the strength envelope in Figure 2 can be considered as the strength margin. Accordingly, the point safety factor based on Mohr-Coulomb criterion can be defined as Equation (3)

Longitudinal Displacement Profile
The longitudinal displacement profile (LDP) describes the displacement at the monitoring section versus distance away from tunnel face and provides an analytical solution to establish a distance-convergence relationship. As shown in Figure 3a, the origin of the x-axis is located at the monitoring section and the x-axis represents the relative location from monitoring section to tunnel face; the y-axis represents the displacement of tunnel boundary at the monitoring section. When tunnel face proceeds along the x-axis to the coordinate value of x, the displacement at the section is denoted as ux. Before the arrival of tunnel face, a portion of displacement has occurred. When tunnel face arrives at the monitoring section, initial displacement u0, approximately 30% of the ultimate displacement u∞, has been achieved, after which the convergence of surrounding rock will increase dramatically to the peak value and keep stable. In other words, the excavation disturbance effect by the tunnel face will disappear until the tunnel face is far away enough, as shown in Figure 3b. It is commonly believed that the distance between the section and tunnel face is no less than three times the opening diameter. In tunnel design, the appropriate location for the installation of the initial support can be determined by the LDP [8].

Longitudinal Displacement Profile
The longitudinal displacement profile (LDP) describes the displacement at the monitoring section versus distance away from tunnel face and provides an analytical solution to establish a distance-convergence relationship. As shown in Figure 3a, the origin of the x-axis is located at the monitoring section and the x-axis represents the relative location from monitoring section to tunnel face; the y-axis represents the displacement of tunnel boundary at the monitoring section. When tunnel face proceeds along the x-axis to the coordinate value of x, the displacement at the section is denoted as u x . Before the arrival of tunnel face, a portion of displacement has occurred. When tunnel face arrives at the monitoring section, initial displacement u 0 , approximately 30% of the ultimate displacement u ∞ , has been achieved, after which the convergence of surrounding rock will increase dramatically to the peak value and keep stable. In other words, the excavation disturbance effect by the tunnel face will disappear until the tunnel face is far away enough, as shown in Figure 3b. It is commonly believed that the distance between the section and tunnel face is no less than three times the opening diameter. In tunnel design, the appropriate location for the installation of the initial support can be determined by the LDP [8].
To reflect the propagation of the convergence during the tunnel face advancing, the dimensionless displacement, λ, is defined as the ratio of u x to u ∞ to normalize the displacement of surrounding rock in Equation (4):  To reflect the propagation of the convergence during the tunnel face advancing, the dimensionless displacement, λ, is defined as the ratio of ux to u∞ to normalize the displacement of surrounding rock in Equation (4): The longitudinal displacement profile can be obtained by the monitoring data. However, the field data usually just contain the excavated part and fail to monitor the unexcavated rock mass. To establish the intact relationship between the dimensionless displacement and distance to tunnel face, numerous authors have proposed alternative expressions. M. Panet [31] derived the relationship for the normalized LDP from elastic analysis as Equation (5): (5) where R is the tunnel radius; x is the monitoring location away from the tunnel face and x > 0 is behind the tunnel face whereas x < 0 is ahead of the face.
Based on the measured data, Hoek [8] suggested the empirical best-fit formula as Equation (6): The elastic theoretical expression was deduced by Zhao et al. [32] as Equation (7): (7) where X and X1 are constants, and, moreover, Additionally, the expression based on Equation (8) was improved by the numerical analysis as [33]: The longitudinal displacement profile can be obtained by the monitoring data. However, the field data usually just contain the excavated part and fail to monitor the unexcavated rock mass. To establish the intact relationship between the dimensionless displacement and distance to tunnel face, numerous authors have proposed alternative expressions. M. Panet [31] derived the relationship for the normalized LDP from elastic analysis as Equation (5): where R is the tunnel radius; x is the monitoring location away from the tunnel face and x > 0 is behind the tunnel face whereas x < 0 is ahead of the face.
Based on the measured data, Hoek [8] suggested the empirical best-fit formula as Equation (6): The elastic theoretical expression was deduced by Zhao et al. [32] as Equation (7): where X and X 1 are constants, and, moreover, X 1 = λ 0 1−λ 0 X. Additionally, the expression based on Equation (8) was improved by the numerical analysis as [33]: Based on the work of C. Carranza-Torres and C. Fairhurst [8], in which they compared Equations (5) and (6) with the measured data of a tunnel in the Mingtam Power Cavern project, Equation (5) could illustrate the excavated part of the LDP curve and that this elastic approximation still had deviation with the field data; Equation (6) could meet the requirement of the description in the excavated part, but the error in the unexcavated part was obvious. Based on the theoretical analysis of the elastic medium, Equation (7) had distinction with the practical rock material. In our previous work, we improved Equation (7) by the elastoplastic numerical simulation and proposed that Equation (8) that had a better effect [33]. Accordingly, we will adopt Equation (8) to describe the LDP curve.

Proposed Methodology
To fulfill the processes of releasing load on the excavation boundary and advancing the tunnel face, the full-face tunnel excavation will be simulated in two separate ways, the one step and the step by step. The procedures for determining the appropriate installation location of the initial support will be realized by combining these two processes.

The Simulation of Releasing Load on Excavation Boundary
The simulation of one-step excavation is realized by releasing excavation load on the excavation boundary step by step. After one-step excavation, as shown in Figure 4, firstly, the equivalent nodal load f (σ) on the excavation boundary, initially equal to the in-situ stress f (σ 0 ), will be extracted and gradually released with the deformation of rock mass increasing. To achieve the simulation of excavation unloading, f (σ 0 ) will be equally divided into n releasing steps. Therefore, the increment of releasing load ∆f (σ) is f (σ 0 )/n in each step. Moreover, the excavation releasing ratio r, defined as the proportion of released load in the total load, can be calculated as the ratio of i to n, where i is the count of releasing step and ranges from 1 to n. Then, the opposite force f (σ r ) = −r·f (σ r ) will be applied on the excavation boundary as the released load to realize the process of unloading. After cycling to equilibrium, this releasing step is accomplished. The point safety factors of each monitoring objective are calculated by Equation (3) and denoted as F i 1 , F i 2 , . . . , F i k . . . . To evaluate the overall surrounding rock and avoid the specificity from the selection of several monitoring objectives, the geometric mean of all point safety factors are denoted as the global safety factor, F i in Equation (9): Similarly, the dimensionless displacements of each monitoring objectives can be calculated and denoted as λ i 1 , λ i 2 , . . . , λ i k . . . . Then, the geometric mean of all dimensionless displacements is calculated as the global dimensionless displacement, λ in Equation (10): In this way, the curve of the global safety factor with the excavation releasing ratio, r − F, and the curve of global dimensionless displacement with the excavation releasing ratio, r − λ, can be plotted. The internal stress on the excavation boundary is f ( it represents that the unloading should be continued until the stress on the excavation boundary decreases to zero.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 7 of 19 zero, it represents that the unloading should be continued until the stress on the excavation boundary decreases to zero.

The Simulation of an Advancing Tunnel Face
To obtain the deformation during tunneling, the advancement of the tunnel face will be simulated. The tunnel excavation footage is determined and the dimensionless displacements of all monitoring objectives will be recorded to calculate the global dimensionless displacement, λ , at each round of excavation. The distance from the monitoring section to the tunnel face is denoted as x. The curve of global dimensionless displacement with the distance away from the tunnel face, λ x  , can be plotted. Fitted by the LDP equation, the mathematical expression between x and λ can be built. In engineering practice, the longitudinal displacement profile can be substituted by the data from in-situ measurement.

The Determination of Installation Location for the Initial Support
From the above steps, combining the curve of

The Simulation of an Advancing Tunnel Face
To obtain the deformation during tunneling, the advancement of the tunnel face will be simulated. The tunnel excavation footage is determined and the dimensionless displacements of all monitoring objectives will be recorded to calculate the global dimensionless displacement, λ, at each round of excavation. The distance from the monitoring section to the tunnel face is denoted as x. The curve of global dimensionless displacement with the distance away from the tunnel face, x − λ, can be plotted. Fitted by the LDP equation, the mathematical expression between x and λ can be built. In engineering practice, the longitudinal displacement profile can be substituted by the data from in-situ measurement.

The Determination of Installation Location for the Initial Support
From the above steps, combining the curve of r − F, r − λ, and x − λ, the relationship between F and x can be established by the bridge of λ, as the comprehensive graph F − r − λ − x shows in Figure 5a. Suppose that the rock mass will enter the critical state of instability when the global safety factor decreases to a certain admissible value, [F]. Consequently, [F] is served as the indicator of the appropriate installation time for the initial support. As the procedure in Figure 5b could be calculated as the relative distance to the tunnel face, x s . Moreover, according to the given advance velocity, v, the timing of installation for the initial support can be deduced as Equation (11): Appl. Sci. 2020, 10, x FOR PEER REVIEW 8 of 19 procedure in Figure 5b, once the global safety factor declines to [ F ], interpolating ] [F F  into the graph, the recommended installation location could be calculated as the relative distance to the tunnel face, xs. Moreover, according to the given advance velocity, v, the timing of installation for the initial support can be deduced as Equation (11):

Numerical Model
The finite difference software package FLAC 3D (version 6.0, Itasca Consulting Group Inc., Minneapolis, MN, USA, 2017) is employed, and a hydraulic tunnel in China is chosen as the practice example to evaluate the method proposed in this study. As shown in Figure 6, altogether 32,349 grid points and 29,920 zones are meshed; the opening diameter D is 10 m. The distance from the right boundary to the tunnel center is 7 times that of D, as well as the left and top boundary. Additionally, the bottom boundary is 60 m away from the tunnel center. The model length is 80 m along the tunnel axis, 8 times that of D. The extent of the numerical model in this study can meet the requirement of computational accuracy [34]. The boundary conditions are shown in Figure 6b: The nodes on the all boundaries except the top boundary are fixed along the normal direction; the 300 m overlying rock is fulfilled by a distributed pressure on the top boundary. As shown in Figure 6c,d, all nodes and zones on the excavation boundary are served as monitoring objects and numbered in consecutive order, including 40 zones and 40 grid points.
Based on the engineering report, the rock mass of this section is composed of granite with the parameters that density of rock mass is 2700 kg•m −3 , and the deformation modulus and Poisson ratio is 8.5 GPa and 0.27, respectively. The friction angle of 40°, the cohesive strength of 0.8 MPa are

Numerical Model
The finite difference software package FLAC 3D (version 6.0, Itasca Consulting Group Inc., Minneapolis, MN, USA, 2017) is employed, and a hydraulic tunnel in China is chosen as the practice example to evaluate the method proposed in this study. As shown in Figure 6, altogether 32,349 grid points and 29,920 zones are meshed; the opening diameter D is 10 m. The distance from the right boundary to the tunnel center is 7 times that of D, as well as the left and top boundary. Additionally, the bottom boundary is 60 m away from the tunnel center. The model length is 80 m along the tunnel axis, 8 times that of D. The extent of the numerical model in this study can meet the requirement of computational accuracy [34]. The boundary conditions are shown in Figure 6b: The nodes on the all boundaries except the top boundary are fixed along the normal direction; the 300 m overlying rock is fulfilled by a distributed pressure on the top boundary. As shown in Figure 6c,d, all nodes and zones on the excavation boundary are served as monitoring objects and numbered in consecutive order, including 40 zones and 40 grid points.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 9 of 19 adopted in the Mohr-Coulomb failure criterion. The vertical in-situ stress is simulated by the rock self-weight, and the horizontal in-situ stress is calculated by lateral pressure coefficient Kx = Ky = 0.4.

The Evolution of the Global Safety Factor with the Excavation Load Releasing Ratio
After each step of the excavation load releasing, the F of 40 monitoring objects are recorded in Figure 7. With the completion of excavation, the point safety factors around the waist continue to decline, while the factors around the crown and the bottom show a dramatic growth before the continuous decrease. When r reaches 60%, the safety factors at tunnel waist are the first to be less than 1, and half of the monitoring objects are in the state of instability.  Based on the engineering report, the rock mass of this section is composed of granite with the parameters that density of rock mass is 2700 kg·m −3 , and the deformation modulus and Poisson ratio is 8.5 GPa and 0.27, respectively. The friction angle of 40 • , the cohesive strength of 0.8 MPa are adopted in the Mohr-Coulomb failure criterion. The vertical in-situ stress is simulated by the rock self-weight, and the horizontal in-situ stress is calculated by lateral pressure coefficient K x = K y = 0.4.

The Evolution of the Global Safety Factor with the Excavation Load Releasing Ratio
After each step of the excavation load releasing, the F of 40 monitoring objects are recorded in Figure 7. With the completion of excavation, the point safety factors around the waist continue to decline, while the factors around the crown and the bottom show a dramatic growth before the continuous decrease. When r reaches 60%, the safety factors at tunnel waist are the first to be less than 1, and half of the monitoring objects are in the state of instability.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 9 of 19 adopted in the Mohr-Coulomb failure criterion. The vertical in-situ stress is simulated by the rock self-weight, and the horizontal in-situ stress is calculated by lateral pressure coefficient Kx = Ky = 0.4.

The Evolution of the Global Safety Factor with the Excavation Load Releasing Ratio
After each step of the excavation load releasing, the F of 40 monitoring objects are recorded in Figure 7. With the completion of excavation, the point safety factors around the waist continue to decline, while the factors around the crown and the bottom show a dramatic growth before the continuous decrease. When r reaches 60%, the safety factors at tunnel waist are the first to be less than 1, and half of the monitoring objects are in the state of instability.  Considering that the structure and the stress state are symmetrical, the point safety factors of monitoring points in Figure 6d from 1 to 20 after each step of excavation load releasing, which are denoted as F i 1 , F i 2 , . . . , F i 20 , are recorded in Figure 8. All these data are contained to a certain range and the degree of dispersion reaches the maximum when the excavation load releases around 70%. Then, at each step, the global safety factor, F, can be calculated by Equation (9) and plotted as a red line in Figure 8. As shown in Figure 8, there are two y-axes. The primary y-axis represents the point safety factor, F, and the secondary y-axis represents the global safety factor, F. For instance, in Figure 8, when r is 10%, the point safety factors of 20 monitoring points around the boundary are calculated by Equation (3) and recorded as a scatter plot, of which the maximum value is 1.842 and the minimum value is 1.516. Then, the geometric average upon these 20 points, namely global safety factor, F, is calculated as 1.705 by Equation (9) and plotted as a red line. In this way, the global safety factors at each calculation step are calculated to obtain the r − F curve with the downward tendency in Figure 8. 70%. Then, at each step, the global safety factor, F , can be calculated by Equation (9) and plotted as a red line in Figure 8. As shown in Figure 8, there are two y-axes. The primary y-axis represents the point safety factor, F, and the secondary y-axis represents the global safety factor, F . For instance, in Figure 8, when r is 10%, the point safety factors of 20 monitoring points around the boundary are calculated by Equation (3) and recorded as a scatter plot, of which the maximum value is 1.842 and the minimum value is 1.516. Then, the geometric average upon these 20 points, namely global safety factor, F , is calculated as 1.705 by Equation (9) and plotted as a red line. In this way, the global safety factors at each calculation step are calculated to obtain the r-F curve with the downward tendency in Figure 8.

The Evolution of Global Dimensionless Displacement with the Excavation Load Releasing Ratio
After each step of releasing load, the dimensionless displacements of monitoring objectives from 1 to 20 in Figure 6d, 20 , are recorded as a scatter plot in Figure 9, corresponding to the primary y-axis. When r is equal to 0%, λ of all points are 0 by default; when r is equal to 100%, all points are considered to reach ultimate displacement and λ of all points are 1. Then, the global dimensionless displacements are calculated by Equation (10) after each calculation step to plot the curve of r-λ as a red line in Figure 9. Before r reaches 60%, the global dimensionless displacement shows the trend of uniform linear growth. After that, there is a rapid rise.

The Evolution of Global Dimensionless Displacement with the Excavation Load Releasing Ratio
After each step of releasing load, the dimensionless displacements of monitoring objectives from 1 to 20 in Figure 6d, λ i 1 , λ i 2 , . . . , λ i 20 , are recorded as a scatter plot in Figure 9, corresponding to the primary y-axis. When r is equal to 0%, λ of all points are 0 by default; when r is equal to 100%, all points are considered to reach ultimate displacement and λ of all points are 1. Then, the global dimensionless displacements are calculated by Equation (10) after each calculation step to plot the curve of r − λ as a red line in Figure 9. Before r reaches 60%, the global dimensionless displacement shows the trend of uniform linear growth. After that, there is a rapid rise.

The Evolution of Global Dimensionless Displacement with the Tunnel Face Advancing
Due to the lack of the relevant monitoring data, the LDP curve will be obtained by the numerical simulation of step-by-step excavation method. The excavation footage is set as 4 m per day. After each excavation step, the λ of all monitoring objectives are recorded to calculate λ in Figure 10a. Then, the mathematical relationship between λ and x is fitted by Equation (8), which was proposed by Zhang et al. [33]. The parameters in Equation (8)

The Evolution of Global Dimensionless Displacement with the Tunnel Face Advancing
Due to the lack of the relevant monitoring data, the LDP curve will be obtained by the numerical simulation of step-by-step excavation method. The excavation footage is set as 4 m per day. After each excavation step, the λ of all monitoring objectives are recorded to calculate λ in Figure 10a. Then, the mathematical relationship between λ and x is fitted by Equation (8), which was proposed by Zhang et al. [33]. The parameters in Equation (8) are λ 0 = 0.273 and X = 3.291, respectively. As shown in Figure 10b, the fitting result is good, and the function of λ for x could be converted using Equation (12). Once λ is certain, the distance away from the tunnel face will be calculated.

The Evolution of Global Dimensionless Displacement with the Tunnel Face Advancing
Due to the lack of the relevant monitoring data, the LDP curve will be obtained by the numerical simulation of step-by-step excavation method. The excavation footage is set as 4 m per day. After each excavation step, the λ of all monitoring objectives are recorded to calculate λ in Figure 10a. Then, the mathematical relationship between λ and x is fitted by Equation (8), which was proposed by Zhang et al. [33]. The parameters in Equation (8) are λ0 = 0.273 and X = 3.291, respectively. As shown in Figure 10b, the fitting result is good, and the function of λ for x could be converted using Equation (12). Once λ is certain, the distance away from the tunnel face will be calculated.

The Determination of Installation Position for the Initial Support
As an example, the admissible value of global safety factor is set to [F] = 1. As shown in Figure 11, interpolating F= [F] into the r − F curve, the corresponding excavation load releasing ratio r s is 96.5%. Moreover, the corresponding global dimensionless displacement λ s is 0.918 by interpolating r s into the r − λ curve. Lastly, substituting the specific λ s into Equation (12), the installation location x s is calculated as 9.88, namely the initial ground support should be installed behind the tunnel face at a distance of no more than 9.88 m. Besides, according to the excavation speed, the timing of installation for the initial support can be determined by Equation (11)  Similarly, if the admissible value of the global safety factor is set to [F] = 1.2, the corresponding installation location will be behind the tunnel face no more than 2.79 m away, namely in 0.7 day after excavation, as shown the above process.
behind the tunnel face at a distance of no more than 9.88 m. Besides, according to the excavation speed, the timing of installation for the initial support can be determined by Equation (11)  Similarly, if the admissible value of the global safety factor is set to [ F ] = 1.2, the corresponding installation location will be behind the tunnel face no more than 2.79 m away, namely in 0.7 day after excavation, as shown the above process.

Parametric Study
To validate the proposed methodology, the cases with different lithological characters and lateral pressure coefficients in Table 1 are implemented to validate the regularity of the proposed method. The case 0 has been analyzed in Section 4 Numerical Analysis of Case Study. The material parameters in case 1 and case 2 remain unchanged with case 0, while the lateral pressure coefficients in case 3 and case 4 are same as that of case 0. Differing from case 0, the lateral pressure coefficient is 1.0 in case 1 and 1.5 in case 2; the rock mass in case 3 is softer than that of case 0 under the geological condition of cretaceous clay rock, while the rock mass in case 4 is harder quartz diorite than that of case 0.  Figure 11. Installation location of initial ground support.

Parametric Study
To validate the proposed methodology, the cases with different lithological characters and lateral pressure coefficients in Table 1 are implemented to validate the regularity of the proposed method. The case 0 has been analyzed in Section 4 Numerical Analysis of Case Study. The material parameters in case 1 and case 2 remain unchanged with case 0, while the lateral pressure coefficients in case 3 and case 4 are same as that of case 0. Differing from case 0, the lateral pressure coefficient is 1.0 in case 1 and 1.5 in case 2; the rock mass in case 3 is softer than that of case 0 under the geological condition of cretaceous clay rock, while the rock mass in case 4 is harder quartz diorite than that of case 0. According to the above procedures, the monitoring data in four separate cases from the numerical simulations are recorded in Figures 12-14, including the evolution of F with r, the evolution of λ with r, and the evolution of λ with x. Then, the geometric mean at each calculation step is calculated to obtain the curve of r − F, r − λ and x − λ and establish the graphic correspondence between F and x in Figure 15. The curve of x − λ can be fitted by Equation (8) for the explicit mathematical expression and the fitting results are listed in Table 2.
According to the above procedures, the monitoring data in four separate cases from the numerical simulations are recorded in Figures 12-14, including the evolution of F with r, the evolution of λ with r, and the evolution of λ with x. Then, the geometric mean at each calculation step is calculated to obtain the curve of r-F , r-λ and x-λ and establish the graphic correspondence between F and x in Figure 15. The curve of x-λ can be fitted by Equation (8) for the explicit mathematical expression and the fitting results are listed in Table 2.  Figure 15. In case 1, with the lateral pressure coefficient of 1.0, when F decreases to [ F ] = 1.0, the corresponding excavation load releasing ratio r = rs is 90.9% and the global dimensionless displacement s λ λ  is 0.698. Substituting into the fitting formula can calculate that the installation location of the initial support, xs, is equal to 4.09, which represents that the initial support should be installed after excavation behind the tunnel face at a distance of no more than 4.09 m. Similarly, under the condition that the lateral pressure coefficient is 1.5 in case 2, when F decreases to 1.0, the corresponding excavation load releasing ratio r = rs is 89.6%, the recommended installation location, xs, is 3.49. Comparing the installation location in case 0, case 1 and case 2, we can find that the initial support should be installed earlier with the increase in the lateral pressure coefficient in Figure 16a.  can be installed behind tunnel face to infinity. As shown in Figure 16b, if the geological conditions are relatively good, the installation of the initial support can be delayed or even unnecessary.       Figure 15. In case 1, with the lateral pressure coefficient of 1.0, when F decreases to [F] = 1.0, the corresponding excavation load releasing ratio r = r s is 90.9% and the global dimensionless displacement λ = λ s is 0.698. Substituting into the fitting formula can calculate that the installation location of the initial support, x s , is equal to 4.09, which represents that the initial support should be installed after excavation behind the tunnel face at a distance of no more than 4.09 m. Similarly, under the condition that the lateral pressure coefficient is 1.5 in case 2, when F decreases to 1.0, the corresponding excavation load releasing ratio r = r s is 89.6%, the recommended installation location, x s , is 3.49. Comparing the installation location in case 0, case 1 and case 2, we can find that the initial support should be installed earlier with the increase in the lateral pressure coefficient in Figure 16a. Moreover, the material parameters of rock mass in case 4, case 0 and case 3 range from hard rock to soft rock and have been compared to the recommended timing for the installation of the initial support determined by the above procedures. In case 3, when [F] is set as 1.0, the appropriate installation location, x s , is 0.73 m; when [F] is set as 1.2, the installation location, x s , is calculated as −2.64 m, namely that the advance support system is suggested to be installed ahead of the tunnel face by 2.64 m at least. Contrarily, in case 4, under the condition of hard and intact rock mass, the surrounding rock keeps stable during excavation and the global safety factor, F, is greater than [F] all the time, which means that initial support system is unnecessary and the support structure can be installed behind tunnel face to infinity. As shown in Figure 16b, if the geological conditions are relatively good, the installation of the initial support can be delayed or even unnecessary.

Conclusions
In this study, a novel method for the installation time of the initial support is proposed. The point safety factor is introduced to represent the stability of each rock zone on the excavation boundary. To describe the overall surrounding rock and avoid selecting the special monitoring points, the global safety factor is defined as the geometric mean of all point safety factors to estimate the stability of surrounding rock at a certain tunnel section. With the releasing of excavation load on the excavation boundary, the evolution of the global safety factor can be obtained by two-dimensional analysis. Based on the in-situ measurement, the complete longitudinal displacement profile can be fitted by the existing LDP equation. Combing the evolution of global safety factor with the LDP, the comprehensive graphic relationship between the global safety factor and the distance from the monitoring tunnel section to tunnel face can be established. When the global safety factor decreases to an admissible value, the surrounding rock is considered in a critical state and the corresponding distance in the above relationship is the recommended location to install the initial support. Meanwhile, several cases with diverse lateral pressure coefficients and geological conditions have been presented to validate the proposed method with good effect.
Although the proposed method has some restrictions, it also provides quick feedback for the initial support design at the typical section in practical engineering.

Conflicts of Interest:
The authors declare that they have no conflict of interest.