After the excavation of the tunnel, the maximum principal stress of surrounding rock is approximately symmetrical, and the maximum value appears at the vault, inflected arch and side wall. Figure 6A illustrates the distribution of the maximum principal stress of the model after geostress equilibrium. As a whole, the tunnel model is under pressure, and the maximum principal stress is approximately linear along the depth direction. After tunnel excavation, the surrounding rock is affected by excavation unloading effect, and the value of the maximum principal stress around the excavated section is greatly reduced, as shown in Figure 6B . The far end of the excavated section is less affected, and the maximum principal stress does not change. With the continuous action of the excavation disturbance, the stress disturbance area radiates annular from the excavation section. And the stress above and below the tunnel excavation section is symmetrical, as is the stress on the left and right sides, as shown in Figure 6C . At the 20,000th time step after tunnel excavation, the surrounding rock reaches the strength limit, and a crack germinates from the arch foot of the tunnel section. Meanwhile, the crack propagation in turn affects the distribution of the maximum principal stress, resulting in the stress concentration near the crack tip. The stress concentration phenomenon will induce the crack to continue to develop, so that the crack continues to expand to the deep of the rock mass ( Figure 6D ), and the location of the maximum principal stress concentration also moves to the deep of the rock mass simultaneously. Finally, the model is calculated to be convergent, cracks fully develop and expand, large deformation occurs on both sides of the spandrel, collapse happens at the vault, slight uplift occurs in the inflected arch, and tension areas are formed in the vault and inflected arch ( Figure 6E ). At this time, the maximum principal stress distribution tends to be stable, and the maximum principal stress is still symmetrically distributed about the center line of tunnel with the phenomenon of stress concentration.

After excavation, the stress of rock mass in the tunnel is released, and the radial compressive stress provided by the rock mass in the cavern disappears, resulting in the surrounding rock to transition from a state of two-direction compressive stress to unidirectional compressive stress. And the stress borne by the rock mass in the cavern is transferred to the surrounding rock, causing the increase of tangential stress in the surrounding rock. In essence, tunnel excavation is a process in which tangential pressure of surrounding rock increases and radial pressure decreases. The surrounding rock breaks when the tangential stress exceeds the compressive strength corresponding to the radial stress. During the failure process after excavation, cracks and spalling occur on both sides of the side wall first, the fracture zone develops along the radial direction, and cracks extend and develop on both sides of the cavern ( Figure 7A ). Subsequently, the cracks continued to grow at the spandrel and intersected at the vault, resulting in settlement deformation of the vault with significant collapse ( Figure 7B ). Eventually, the surrounding rock breaks and expands, and the rock mass moves toward the center of the cavern, causing the deformation of the tunnel section and the reduction of the area.

When the tunnel face was excavated, stress concentration occurred at the arch foot of the tunnel, and the crack was first generated at this position ( Figure 8A ). With the advance of tunnel excavation time, the stress concentration near both sides of the side wall intensifies, and the crack continues to expand to the deeper rock mass until the concentrated stress at the crack tip reaches the ultimate equilibrium state. As the whole surrounding rock is in a state of compression, the vertical geostress is much greater than the horizontal geostress, which makes the cracks mainly expand at the side wall, and the cracks are mainly tensile-shear mixed crack ( Figure 9 ), showing a conjugate cross distribution ( Figure 8B ). After the shear crack has fully extended, the rock mass on both sides of the side wall gradually loosens and falls off, converging and deforming into the cavern. The broken rock blocks extrude each other, causing them to overturn into the cavern. This results in a decrease in stress in the fracture zone, forming a local tensile stress zone. Subsequently, cracks on both sides of the side wall gradually extended to the vault and intersected, causing the vault to collapse ( Figure 8C ). After tunnel excavation, the failure zone is mainly distributed in the spandrel, taking on an inverted V shape. Tensile cracks are the primary cracks in the shallow part of the surrounding rock, while shear cracks are the primary cracks in the deep. The conjugate shear angle gradually decreases and annihilates in the deep part.

The finite element model with the same mechanical parameters and geometric dimensions as the FDEM model is established to analyze the tunnel stability, in which the surrounding rock adopts the Mohr-Coulomb constitutive model. The plastic zone distribution after tunnel excavation was calculated, and the specific simulation results are presented in Figure 10 . By comparing the plastic zone and fracture zone calculated by both methods, it can be observed that the failure location and degree are essentially identical, and the contours of the failure area are in good agreement. This verifies the validity and rationality of FDEM for tunnel stability analysis. Moreover, FDEM successfully simulates the gradual transformation process of the surrounding rock from continuous to discontinuous under the influence of excavation unloading effect. It also reproduces the process of crack germination, propagation, and collapse in the process of tunnel surrounding rock failure after excavation, which is more consistent with the failure mode of tunnel surrounding rock as a discontinuous medium. The simulation results are more realistic, intuitive, and consistent with actual tunnel excavation failure mode.

The model of a high-speed railway tunnel with a design speed of 350km/h in IV-level surrounding rock is utilized. The ratio of the uniaxial saturated compressive strength ( R c ) of the rock to the maximum initial geostress ( σ max ) in the vertical axis direction is used to reflect the initial geostress state, in accordance with the evaluation standard of the initial geostress state in the Code for Design of Railway Tunnel (TB10003-2016). Other calculation parameters are held constant, and the R c / σ max value is varied to 1.5, 3, 4.5, 6, 7.5, and 9, resulting in a total of six working conditions. This is used to simulate three stress states: general stress, high stress, and very high stress.

Based on the engineering geological conditions of the Wufeng Tunnel, Silurian shale has been chosen as the representative rock for the tunnel surrounding rock. The uniaxial saturated compressive strength of the shale has been determined to be 18MPa based on engineering data and serves as the benchmark for calculations. And then the simulated working conditions under various initial geostress are confirmed ( Table 3 ). Since the deep buried tunnel is mainly affected by geostress, gravity is not considered in the model. The application method of initial geostress is shown in Figure 11 .

Three different conditions of geostress were chosen for analysis: general R c / σ max = 9 , high R c / σ max = 6 , and very high R c / σ max = 3 . The simulation results are presented in Figure 12 . When the geostress is low ( Figure 12A ), the tunnel’s surrounding rock can maintain a stable state following excavation, and the tension area is limited to within 2 meters near the tunnel section. Meanwhile, the excavation’s disturbance does not propagate deep into the surrounding rock, resulting in a small stress disturbance area. When R c / σ max = 6 ( Figure 12B ), the surrounding rock collapses, causing the tunnel’s critical stable section to shift deeper. The distribution range of the tension zone increases to twice the tunnel diameter, and the stress disturbance zone caused by excavation forms an obvious evenly-distributed ring. When the geostress reaches a very high state ( R c / σ max = 3 , Figure 12C ), the collapse range of the surrounding rock is further expanded, and the distribution range of the tension zone extends to three times the diameter of the tunnel. The stress disturbance caused by excavation also spreads to the boundary of the model. In conclusion, under a geostress field with a lateral pressure coefficient of 1, the stress value of the surrounding rock near the excavation section will sharply drop after excavation. The tension zone will appear near the excavation section, and the stress disturbance zone will be distributed annularly around the contour surface of the cavern. Then the surrounding rock stress gradually transition from tensile stress in shallow part to the original stress in deep part of tunnel face. Compared with the simulation results of gravity stress field, the maximum principal stress distribution in geostress field does not exhibit an obvious stress concentration phenomenon.

The working condition with high geostress ( R c / σ max = 6 ) was selected for analysis, and the simulation results are presented in Figure 13 . Following tunnel excavation, failure occurred around the tunnel face, and the overall failure form was circular. Initially, cracking and loosening occurred in the vault, side wall and inflected arch, and the fracture area gradually developed to the depth of surrounding rock ( Figure 13A ). The volume of rock block produced by shear crack cutting in deep surrounding rock increases gradually. Additionally, the deep crushed rock compresses the shallow rock, making the shallow rock turn over to the cavern ( Figure 13B ). The large overturning movement of rock blocks and the contact non-occlusion between the blocks make the surrounding rock produce the large deformation of crushing expansion in meter scale. Compared to the failure mode under the gravity stress field, the failure of the surrounding rock under the geostress field is influenced by both horizontal and vertical compressive stresses. There was significant uplift in the inflected arch and extensive collapse in the vault. Moreover, the failure zone of the surrounding rock exhibits no apparent directional preference. These observations suggest that the horizontal geostress plays a critical role in determining both the failure mode and the extent of surrounding rock.

High geostress ( R c / σ max = 6 ) is selected for analysis. The simulation results demonstrate that the cracks produced due to tunnel excavation are primarily tensile-shear mixed crack, followed by shear cracks and a small proportion of tensile cracks. The overall distribution pattern of cracks appears to be annular and conjugate, with a higher density of cracks observed near the excavation section and a lower density in the deep surrounding rock. After the tunnel excavation, cracks appear in the shallow part of the surrounding rock due to the action of both vertical and horizontal geostress. These cracks are mainly tensile-shear mixed crack ( Figure 14A ). Subsequently, the cracks progressively extend to the depth of surrounding rock, with the maximum propagation depth being 30.4 m, approximately twice the tunnel’s equivalent diameter. The cracks observed in the deeper surrounding rock are mainly shear cracks, exhibiting a distinct conjugate distribution pattern ( Figure 14B ). Eventually, the shear angle of the crack decreases gradually and becomes annihilated along the depth direction, which leads to the cessation of failure of the surrounding rock.

Figure 15 illustrates the failure of the surrounding rock after excavation under six different geostress conditions. It is observed that only at a geostress of 2 MPa, the surrounding rock can maintain the stability. Under all other conditions, the tunnel experiences deformations at the meter-level after excavation, and large deformations occur under very high geostress conditions. This emphasizes the critical role of the stress state of the tunnel in maintaining the stability of the surrounding rock. When the geostress is 2 MPa ( Figure 15A ), the displacement of the surrounding rock after excavation is limited to the millimeter level, indicating robust tunnel stability. It is noteworthy that the maximum displacement of the surrounding rock is concentrated at the vault and inflected arch, whereas the maximum displacement of the tunnel after excavation under the gravity stress field is generated in both sides of the side wall. This is due to the presence of horizontal geostress that alter the stress state of the cavern. Therefore, when the tunnel is located in an engineering rock mass with a significant geostress field, special attention should be paid to the stress on the vault and inflected arch during the design and construction process. As the geostress increases to 4 MPa ( Figure 15B ), the surrounding rock begins to collapse in the shallow section, with the failure area distributed evenly around the cavern contour, and the maximum failure depth reaching 7.4 m. When the geostress gradually increases to 4.8MPa and 6MPa (high geostress), the failure zone begins to propagate to the depth of the surrounding rock. The deep broken rock began to extrude the shallow rock, forcing the shallow rock to overturn to the cavern. The maximum depth of the corresponding failure zone reaches 16m and 25m respectively ( Figures 15C, D ). Under very high geostress ( Figures 15E, F ), the failure zone of surrounding rock further radiates outward, and the maximum depth of the failure zone under 8MPa and 12MPa conditions reaches 31 meters and 42 meters respectively. Meanwhile, the surrounding rock at the excavated section is more broken, and loose and collapsed rock blocks jam the entire tunnel face.

In general, the IV-level surrounding rock can maintain a certain stability after excavation under general geostress conditions. However, as geostress increases, a failure zone appears around the excavation section and gradually develops deeper into the surrounding rock. At very high geostress conditions, the failure zone can extend more than three times of the tunnel diameter. At this point, the rock generated from collapse completely blocks the cavern, resulting in the overall collapse and destabilization of the tunnel.

Based on the strength reduction method, the FDEM is utilized to determine the required minimum supporting force for maintaining the equilibrium of the failure zone of the surrounding rock at the tunnel vault. Then 1.4 times of the minimum supporting force is taken as the design value of surrounding rock pressure. Figure 16 shows the schematic diagram of the calculation principle of the design value of surrounding rock pressure. The specific steps are as follows:

(1) A numerical model of FDEM is established. Considering the weakening effect of groundwater, joint surface and other factors on the strength of surrounding rock, the strength of mechanical parameter of surrounding rock is reduced, and the reduction coefficient is generally 1.15 (Xiao, 2020). After the tunnel excavation, the model is calculated to converge, and the crack zone and failure zone of the surrounding rock are computed separately.

(2) The gravity of surrounding rock failure zone within the span of the cavern is calculated, and it is equivalent to the uniform load q 0 along the span direction.

(3) Following the tunnel model described in (1) an excavation simulation is again conducted. Normal supporting force of initial value q 0 is applied to the tunnel section, and the supporting force is gradually adjusted until there is no failure zone around the section. The supporting force is recorded at this point as the minimum supporting force q min . Then, the horizontal load is calculated by multiplying the vertical load q min with the lateral pressure coefficient.

(4) The design value q of the surrounding rock pressure is determined by multiplying the minimum supporting force with the safety factor k. k is usually greater than 1.4 and can be adjusted based on factors such as the importance of the engineering project and the degree of deformation control required.

In the above steps, the failure zone is identified based on displacement. When the displacement of the surrounding rock near the excavation section reaches a specific value, it is considered locally unstable and the failure zone is formed. The displacement criterion refers to the allowable relative convergence value around the cavern ( Table 4 ) in the Technical Code for Engineering of Ground Anchorages and Shotcrete Support (GB50086-2015). The relative convergence value around the cavern represents the ratio of the measured displacement between two measuring points and the distance between them, or the ratio of the measured displacement of the vault and the span of the tunnel.

Taking a high-speed railway tunnel with a burial depth of 200m and a design speed of 350km/h in the IV-level surrounding rock as an example, the design value of surrounding rock pressure is calculated. After the tunnel excavation, the equivalent uniform load of gravity of the failure zone in the tunnel vault is calculated as 16kPa. Subsequently, excavation simulation is conducted again, followed by the application of supporting force, and the model is calculated until it converges. Finally, the supporting force is gradually adjusted until the surrounding rock no longer produces a failure zone after excavation.

The simulation results under different supporting forces are shown in Figure 17 . When the supporting force is 16kPa ( Figure 17A ), the simulation results show that when the supporting force is small, the cracks will intersect at the tunnel vault, and the side walls on both sides will converge and deform towards the cavern. The surrounding rock around the cavern will collapse, and the excavation disturbance area will be large. With 25kPa supporting force applied ( Figure 17B ), the excavation disturbance area decreases significantly. Although cracks still appear at the vault and intersect at the right spandrel, no rock falls off from the rock mass. When the supporting force increases to 54kPa ( Figure 17C ), the excavation disturbance area is concentrated in a smaller region around the cavern. The cracks in the surrounding rock decreased obviously and mainly concentrated on both sides of the side wall. The cracks do not intersect at the vault, which remains stable, and no rock spalling occurs on both sides of the side wall. The relative convergence value of the tunnel is 0.99%, and the tunnel remains stable. Subsequently, the supporting force continues to increase, but the relative convergence value of the tunnel remains stable at about 1% ( Figure 18 ), meeting the requirements of surrounding rock stability. Therefore, 54kPa is the minimum supporting force required for the high-speed railway tunnel.

When the buried depth of the tunnel is no less than 10 ~ 15 times of the tunnel diameter, the formula for calculating the surrounding rock pressure is as follows for surrounding rock that conforms to the Mohr-Coulomb yielding criteria (Cai, 2002).

Vertical uniform load:

Horizontal uniform load:

where, γ is the volumetric weight of surrounding rock; λ is the lateral pressure coefficient of surrounding rock; α and β are the adjustment coefficients of surrounding rock pressure at the arch and side of the tunnel respectively; and R p d is the radius of the tunnel plastic zone at 45° when the supporting force P i =0. p 0 , c , and φ represent the initial stress, cohesive force, and internal friction angle of the surrounding rock respectively. θ is the included angle with the horizontal axis of the tunnel; R 0 is the radius of tunnel. If the section is not circular, the radius of the outer circle of the tunnel section is taken. a is the distance from the center of the outer circle of the tunnel section to the tunnel excavation face in the direction of 45°.

When applying the aforementioned theoretical formula to resolve the surrounding rock pressure, the outer circle of tunnel section should be taken as the equivalent circular tunnel first. Subsequently, the theoretical solution for the plastic zone of the equivalent circular tunnel is computed, along with the determination of the depth of the plastic zone in the 45° direction. The gravitational force exerted by the surrounding rock within this depth range serves as the fundamental value of the surrounding rock pressure at the top of the tunnel. Ultimately, the approximate value of the surrounding rock pressure is obtained by multiplying the fundamental value by the adjustment coefficient.

The relevant parameters of the tunnel model in Section 5.1.2 are applied to equations (10) and (12), resulting in the calculation of a surrounding rock pressure of 62 kPa at the top of the tunnel. Although slightly higher than the minimum supporting force determined by the FDEM method in Section 5.1.2, the two values are close. The differences in the calculation results between the two methods are primarily attributed to the parameter related to the importance of the tunnel, which is the adjustment coefficient in the theoretical formula, as well as the allowable relative convergence value around the cavern in the FDEM-based calculation method. When the allowable relative convergence value is set to 0.75%, the calculation results from both methods can be equal ( Figure 18 ). Therefore, it is important to select an appropriate allowable relative convergence value based on the significance of the tunnel. This finding further validates the accuracy of the FDEM-based calculation method for determining the design value of surrounding rock pressure. Moreover, it should be noted that the theoretical formula can solely provide the surrounding rock pressure for a circular tunnel in a homogeneous stratum, whereas the FDEM-based approach allows for the determination of the design value of surrounding rock pressure under complex section and surrounding rock conditions.

The calculation method, based on FDEM as proposed in Section 5.1, was utilized for the calculations. Table 5 presents the calculation results for the design values of surrounding rock pressure of the high-speed railway tunnel with a design speed of 350km/h in the IV-level surrounding rock, considering various geostress. It should be noted that the safety factor “k” in the table is maintained at 1.4.

Based on the calculation results of the design value of surrounding rock pressure, it is evident that engineering support measures need to be implemented in all five working conditions, except when the geostress is 2MPa, in order to meet the design requirements. As depicted in Figure 19 , it can be observed that the equivalent uniform load, the minimum supporting force, and the design value of surrounding rock pressure exhibit a significant increase with the augmentation of geostress. Specifically, when the geostress varies from 2MPa to 12MPa, the corresponding equivalent uniform load ranges from 18KPa to 275KPa, the minimum supporting force escalates from 35KPa to 920KPa, and the design value of surrounding rock pressure rises from 49KPa to 1288KPa. The curve that represents the design value of the surrounding rock pressure in response to geostress approximately conforms to the quadratic function relation. By fitting the formula, the equation q = 14.0637 σ 2 − 68.5260 σ + 85.8181 is derived. Notably, all three indicators steadily increase in the general and high geostress stages but sharply increase in the very high geostress stage. Moreover, the curve that illustrates the design value of the surrounding rock pressure shows obvious nonlinearity. This is because under extreme geostress conditions, the tunnel’s surrounding rock generally approaches the limit state. This results in a severe stress-release phenomenon after tunnel excavation, leading to a broad distribution of surrounding rock failure areas. To control the destabilization and collapse of the large-scale failure area, significant supporting force is necessary.

Additionally, it should be noted that there is an approximate linear relationship between the minimum supporting force and the equivalent uniform load ( Figure 20 ), and the corresponding fitting formula is q = 3.1473 q 0 − 12.5888 . This supporting force can restrict the deformation of the surrounding rock of the tunnel after excavation to a small scale and ensure the integrity of the tunnel section. This calculation result also proves the correctness and reasonableness of computing the design value of the surrounding rock pressure based on the equivalent uniform load of the surrounding rock gravity in the failure zone within the cavern span.

In summary, geostress directly affects the design value of the surrounding rock pressure of tunnels, particularly under very high geostress. An obvious correlation exists between the minimum supporting force and the equivalent uniform load of the surrounding rock gravity in the failure zone within the cavern span, so the uniform load can be used as the reference value to calculate and verify the final design value of surrounding rock pressure.