Improved Protodyakonov ’ s Method of the Tunnel Surrounding Rock Pressure under the Seepage Condition of Weak Interlayer

,


Introduction
From 2009 to 2021, the construction of China's road tunnel has shown a rapid growth in speed and scale [1][2][3] (Figure 1). With the complex geological conditions in China, the construction of large-scale tunnels would inevitably lead to disasters. Construction causes technological disasters such as landslides and flash floods, which disrupt constructions, delay operations, and place human lives and properties under serious threat [4]. Among them, groundwater and weak interlayer are typical representatives. Therefore, because of this engineering geological situation, it is essential to put forward a reasonable calculation method for calculating the surrounding rock pressure of tunnels considering their seepage and weakness.
To calculate the pressure on the surrounding rock of tunnels, it was initially thought that the pressure on the sup-port system mainly came from the gravity of the overlying soil, and the formula γH was used for the calculation. However, the computed results differ significantly from the actual situation. With the continuous development of subsequent results, the collapse vault theory plus Protodyakonov's theory and K. T. Erzahi's idea were used together [5]. Austrian scholar L.V. Rabcewicz revised the New Austrian tunneling method, a more practical tunnel construction method, by summarizing the method used to calculate the surrounding rock pressure [6]. Furthermore, the energy support theory was formed through scholars' systematic results and discussion [7,8]. According to this theory, the sum of the energy released by the surrounding rock and the energy absorbed by the support system are fixed. Therefore, the support system should be adjusted to minimize the energy absorbed by the surrounding rock for it to be stable. To calculate a large section tunnel with a lower flattening rate, Tricot et al. [9] used the method for calculating the surrounding rock pressure; it is a more suitable method used for practical engineering. Based on the actual monitoring engineering data, Chen [10] obtained the law of stress change in the circumferential and radial directions of a rock mass and studied the self-supporting state of the pressure arch of the surrounding rock. Paolo and InggLorenzo [11] considered that if the support scheme is not reasonable, the residual stress of some parts of a tunnel will be released during the excavation process. This makes the relatively weak surrounding rock easily gets to the compressive or compressive ultimate yield state, resulting in the instability of the tunnel [12][13][14]. Wang studied the method used for calculating the small clearances of tunnels with large cross-sections. Based on the calculation result of the surrounding rock pressure obtained by the code method, the pressure of the surrounding rock is amplified, or the elastic resistance of the middle sandwich wall is reduced according to the influence law of the smalldiameter distance [15][16][17][18]. As the method used for calculating the surrounding rock loses pressure in a shallow tunnel, a new calculation model used for the surrounding rock pressure in the folded zone is proposed [19]. Sarfarazi et. al., Zhang et al., and Zhang and Qui reported that with the multifactor surrounding rock pressure calculation theory, the surrounding rock pressure of arbitrary shape and span can be easily calculated and analyzed, as well as its stability and safety [20][21][22]. Wen et al. [23] introduced the weight coefficients of tunnel depth H, tunnel span B, internal friction angle φ, and rock weight γ and deduced the calculation formula of the surrounding rock pressure for extensive excava-tion of section tunnels. The shear strength of surrounding rock is always affected by joints. Sun et al. and Jiang et al. [24][25][26] visually observed the failure process through a discrete element numerical direct shear test, and the shear strength is related to the failure form and mechanism. Numerical simulation combined with engineering experiments verify the effectiveness of the formula. The finite element analysis is given to determine the formula of the type 1 fracture toughness using the compact tension (CT) test. The hollow center-cracked Brazillian sample was used for validation, simulation, and numerical simulation. The cross-validation test confirms the results have high reliability.
To sum up, the related previous scholars have done indepth and extensive studies on calculating the surrounding rock pressure of tunnels and have achieved specific results. Tesaghi's theory, Protodyakonov's theory, and code methods are commonly used in the computation of surrounding rock pressure in underground engineering design in China [27][28][29][30][31]. However, little result has been obtained from the calculation method used for the surrounding rock pressure due to the simultaneous effect of weak interlayer and seepage. Compared with the above three calculation schemes of the surrounding rock pressure of tunnels, it is found that the calculation method used for the surrounding rock pressure of tunnels based on Tesaghi's theory is more suitable for shallow buried tunnels in terms of principle and hypothesis, and the calculation results are relatively conservative. Protodyakonov's theory and the code method can address the characteristics of deeply buried tunnels, while

Geofluids
Protodyakonov's theory can more precisely reflect the excavation process. This work developed a modified Protodyakonov's method (MPM) for calculating the surrounding rock pressure of a tunnel under the coupled condition of seepage and weak interlayer. Based on Protodyakonov's theory, the finite difference method is adopted to determine the modified method of calculating the surrounding rock pressure by taking the angle and position of the interlayer as the entry point. The case project of this is the surrounding rock of a tunnel in Chongqing affected by the seepage of its weak interlayer. The findings of three alternative calculation methods were compared with the monitoring results obtained from Chongqing Mountain tunnel calculated with the MPM to confirm its efficacy. The results show that the calculation results of the modified method are closer to those of the field monitoring, and their values are more significant than those of the field monitoring. They meet the requirements of the support strength. Furthermore, this work discussed the research results of the modified coefficients considering the angle and distance of the weak interlayer to the general significance of the engineering application.

Developing MPM for Calculating Surrounding Rock Pressure
This paper mainly used the modified coefficient to measure the adverse influence of the corrected weak interlayer and its seepage on the tunnel support. The assumption conditions of the modified formula obey the assumption conditions of Protodyakonov's theory and meet the requirement that the application range contains a single main weak sandwich, where the dip angle of the sandwich is x 1 and the content is 0~90°. Furthermore, the distance from the tunnel vault is x 2 . Within 0~0.5 times the span, the interlayer is in a saturated state. k 1 is introduced considering sthe influence of the dip angle of the interlayer on the surrounding rock pressure of the tunnel, and k 2 considers assessing the impact of the interplay at different positions on the surrounding rock pressure of the tunnel; the modified simplified formula is shown as where k 1 is the modified coefficient of the dip angle of the interlayer on the surrounding rock pressure, and k 2 is the modified coefficient of the interlayer position on the rock pressure. The numerical values of k 1 and k 2 are fitted and determined using the finite difference method's numerical simulation of monitoring data.

Permeability Coefficient of Dynamic Change under
Seepage Condition. The contact surface relationship between the soft and weak interlayer and the surrounding rock is treated by the interface, and the seepage effect during the excavation of the surrounding rock of the tunnel is simulated using the flow-solid coupling medium. Finally, the dynamic change of the permeability coefficient under the seepage state of the surrounding rock during the excavation is realized by Fish language. The permeability coefficient in weak FLAC3D differs from the soil mechanics concept. The international unit of permeability coefficient K is m 2 /Pa/sec, and the conversion relation between the permeability coefficient K (cm/s) in soil mechanics is expressed as Weak FLAC3D also has the advantage of implementing some postprocessing functions in the Fish language. The permeability coefficient in all numerical simulations in this paper is obtained from the model test multiplied by 1.02e-6 to perform the simulation [32]. The Kozeny-Carman equation is widely used, and the Kozeny-Carman equation for the variation of the permeability coefficient with porosity is shown as where K is the permeability coefficient of rock masses; K z is constant, approximately equal to 5; Φ is porosity, %; Σ is body surface area; S p is pore surface area. It can be seen from Equation (3) that a proportional relationship exists between the porosity and permeability coefficient of rock. The relationship between volumetric strain and porosity is shown as where Φ 0 is original porosity %, and ε v is under volumetric strain.
The relationship between the permeability coefficient and volumetric strain is shown as Through the correlation between the stress-strain state of the surrounding rock and the permeability coefficient, the equation of permeability coefficient and volumetric strain is written into the software as a command stream in a Fish language for calculation. The real-time permeability coefficient is obtained by monitoring the volumetric strain of the model unit during the simulation calculation in FLAC3D [33][34][35].

Establishment of the Numerical Simulation Model and
Setting of Boundary Conditions. The numerical simulation is based on a mountain tunnel project in Chongqing, and FLAC3D finite difference software is used for the numerical 3 Geofluids simulation. For the numerical simulation, hexahedral, wedge, and column grids are used to form the essential grid cells of the model, with a total number of 39324 cells and 43654 nodes. The physical parameters of the surrounding rock and the weak interlayer are selected as shown in Table 1, and the contact between the weak interlayer and the adjacent surrounding rock is simulated by the contact surface cell of the Coulomb shear intrinsic model (in Figure 2). The dimension of the numerical simulation is 54 m × 54 m × 18 m. The tunnel is a single horseshoe-shaped tunnel with a span diameter of 13.5 m, a clear height of 9.45 m, and an interlayer thickness of 0.675 m ( Figure 3). Taking the center of the elevation arch as the origin, the tunnel extends by 27 m along the x-axis in both positive and negative directions, and 18 m along the y-axis, which is the longitudinal length of the tunnel. More-Coulomb elastoplastic constitutive model is used for the simulation of the surrounding rock and its weak interlayer. The method adopted for the constitutive model is the "mixed discrete method" which is used to simulate the plastic characteristics of the surrounding rock and its interlayer. This method is more accurate and applicable than the "discrete integrated method" commonly used in finite elements.
Stress boundary conditions. The numerical simulation models the entire submerged layer, i.e., 61.5 m to 115.5 m below the ground. The model above the submerged layer is   Seepage boundary and stress boundary generation. In the control boundary conditions, first the calculation of stress field is turned off, and then the seepage field is opened, followed by seepage field balance calculation, and so on. Each node and the unit of pore water pressure field are generated after turning off the seepage mechanics equilibrium calculation based on the seepage boundary individual balance with mechanics. Then, the corresponding stress field is developed. In other words, the stress limit on the water is added to the seepage limit after the initial seepage equilibrium.
After the first stress field is formed on the ground and the seepage field, excavation occurs. The mechanical equilibrium is first performed separately after completing the excavation. This is done to simulate the process of stress redistribution to form a new stress field after excavation unloading. Then, the seepage field and stress field are opened simultaneously to simulate the transient seepage calculation using a fluid-solid coupling, and the excavation cycle is repeated eight times. Each working condition is consistent with the model boundary in the calculation process, and the simulation calculation of constant water level and constant ground stress field is considered. Monitoring points at the left arch waist, right arch waist, and arch top are set up to monitor the surrounding rock pressure value during the whole excavation process.

Influence of Dip Angle of Interlayer under Seepage
Condition. To further analyze the influence of the interlayer dip angle, the interlayer angle was altered to 0~90°, and the surrounding rock pressure was observed at the corresponding position.
(1) The influence rule of angle change on the pressure of left hance, right hance, and the surrounding rock of the vault is shown in Table 2   Table 2 shows the change of the surrounding rock pressure of different angle interlayers in an anhydrous state. It is possible to see that when the infill angle is less than 45°, the pressure of the surrounding rock at an angle of inclination increases. The surrounding rock pressure decreases when the interlayer's angle is more significant above 45°. The interlayer's angle has an apparent influence on the surrounding rock pressure. The maximum value is reached in the interlayer's angle range of 40°to 60°. The surrounding rock pressure of the interlayer with different dip angles presents a Gaussian distribution. Therefore, Gaussian distribution is intended to be used to fit the surrounding rock (2) Gaussian distribution is used to fit the surrounding rock pressure of the left hance of the interlayer at different dip angles (Figure 4), where y 0 = 0:413, A = 32:3, ω = 32:5, and R 2 = 0:961 > 0:95. FLAC3D simulation results show that the pressure of left hance surrounding rock without interlayer is 0.095 Mpa, so the modified coefficient k l1 of the interlayer with different dip angles is obtained. The calculation formula is shown in Equation (7). Gaussian distribution is used to fit the surrounding rock pressure of the right hance of the interlayer at different dip angles ( Figure 5), where y 0 = 0:094, A = 37:2, ω = 45:9, and R 2 = 0:974 > 0:95. FLAC3D simulation results show that the right hance surrounding rock pressure without interlayer is 0.095 MPa; the modified coefficient k r1 of interlayer with different dip angles is obtained.
The calculation formula is shown in Equation (8).
Gaussian distribution was used to fit the surrounding rock pressure of the sandwich vault at different dip angles (in Figure 6) To further explore the influence of the interlayer position, the distance between the interlayer and vault is changed by 0~0.5 m, and the surrounding rock pressure is monitored at the relevant position.
(1) The influence rule of position change on the pressure of the left hance, right hance, and surrounding rock of the vault is shown in Table 3   Table 3 shows the change of the surrounding rock pressure of the interlayer with different distances from the vault in an anhydrous state. We can observe that with the increase in the spacing of the vault, the pressure around the rock   7 Geofluids changing with the angle of the interlayer is obtained. The exponential function is adopted for this fitting, is shown as (2) Exponential distribution is used to fit the pressure of the surrounding rock on the left hance of the interlayer at different angles (in Figure 7), where y 0 = 0:434, A 1 = 0:704, t 1 = 0:171, and R 2 = 0:982 > 0:95.
According to the FLAC3D simulation results, the surrounding rock pressure of the left hance of the weak interlayer tunnel with zero span spacing to the tunnel wall is 1.086 Mpa, and the modified coefficient k l2 of the position of different weak interlayers is obtained. The calculation formula is shown in Equation (11). Exponential distribution is used to fit the surrounding rock pressure of the right hance of the interlayer at different angles (Figure 8), where y 0 = 0:086, A 1 = 0:513, t 1 = 0:161, and R 2 = 0:976 > 0:95. According to the FLAC3D simulation results, the surrounding rock pressure of the right hance of the weak interlayer tunnel with zero span spacing to the tunnel wall is 0.571 MPa, and the modified coefficient k r2 of the different interlayer positions is obtained. The calculation formula is shown in Equation (12). The exponential distribution is used to fit the surrounding rock pressure of the interlayer vault at different angles (in Figure 9), where y 0 = 1:58, A 1 = 4:01, t 1 = 0:421, and R 2 = 0:968 > 0:95. According to the FLAC3D simulation results, the surrounding rock pressure of the tunnel vault without a weak interlayer is 5.371 Mpa, and the modified coefficient k t2 of the interlayer with different dip angles is obtained. The calculation formula is shown as

Case Application of the MPM
The reasonableness and effectiveness of the MPM used for the mountainous tunnel in Chongqing region is validated. The surrounding rock pressure of the project's field is monitored. The monitoring results are compared with the results of the surrounding rock pressure calculated with the commonly used code methods, Protodyakonov's method, and MPM.

Overview of the Case Project.
The case project is a mountain tunnel with single-hole double-track in Chongqing. The tunnel length is 853 m, its buried depth is 110 m, span length is 23.60 m, and height is 20.83 m. The tunnel is a curved wall-round vault section. The survey area is within the structural belt of the southeast arc of Sichuan, and the site is in the western wing of the Nanwenquan anticline. The strata dip is about 290°, and the dip angle is about 60°. The interface plane of sandstone and sandy mudstone is poorly combined, and the phenomenon of thin-bedded marginalization appears. According to the engineering characteristics of the proposed tunnel, the diving layer is 40 m below the ground, and there is mainly bedrock vein-like fissure water and open layer pore water (upper stagnant water) related to the proposed project. A typical AK29+834 section was selected to study the influence of the weak interlayer on the surrounding rock pressure of the tunnel structure. The upper step section of the surrounding rock in the layered excavation is shown in Figure 10. A weak interlayer with a thickness of 1.2 m~1.6 m and a dip angle of 70°crosses the vault between the inclined rock layers. According to geological prospecting data and field observations, the lithology of an interlayer is mostly strongly weathered mudstone, which has low strength and large compression modulus, and the water strength of an interlayer soil is significantly reduced or even lost.
To study the surrounding rock pressure of the initial support, the surrounding rock pressure of the initial support at the left and right hance is monitored from the excavation process to the excavation stability period. The earth pressure box is buried in the surrounding rock in contact with the ini-tial support. The monitoring points are distributed in the same section, and the monitoring points arranged are firm and reliable to ensure that the monitoring data are accurate and effective. The detailed distribution section of the earth pressure box is shown in Figure 11.
The data obtained during the excavation are shown in Table 4. The surrounding rock pressure monitoring values of the left and right hance and vault reached a fixed deal   after October 23 and tended to be stable. The curve of the surrounding rock pressure change is shown in Figure 12.

Analysis of the Results of Three Calculation Methods.
The monitoring results and the values of the surrounding rock pressure calculated by the commonly used code method, Protodyakonov's method, and the MPM derived are listed in Table 5. The comparative analysis in Table 5 shows that the monitoring values in the field are generally larger than those calculated in the Code for Design of Road Tunnel, indicating that the surrounding rock pressure calculated by the code method is smaller than the monitoring values in the case of weak interlayers. The surrounding rock pressure values of the vaulted roof and the left and right hance of the tunnel with weak interlayer calculated using the Presbyterian theory are smaller than those monitored in the field. Therefore, the surrounding rock pressure calculated by the code and Protodyakonov's theory method cannot reach the support strength of the surrounding rock of the tunnel with a weak interlayer, and the relative safety risk is significant. At the same time, the pressure value of the surrounding rock calculated by the modified formula is greater than the field monitoring value, and the two values are close. The errors of the surrounding rock pressure calculated by the MPM with the area monitoring results are 38.8%, 51.4%, and 166%, respectively. Therefore, to satisfy the support strength of the vault and the right hance, it is more appropriate to conservatively consider the calculation scheme of the support strength of the left hance. Compared with the usual calculation method used for the surrounding rock pressure, the result obtained from the surrounding rock pressure calculated with the MPM has more guiding significance for the support design of similar engineering with the weak interlayer tunnel.

Discussion
The calculation result of the surrounding rock pressure of a mountain tunnel with single-hole double-track in Chongqing obtained by MPM is greater than the monitoring values (Table 5); while the calculation results of the other two methods are less than the monitoring values. To ensure that the supporting measures meet the requirements of strength, MPM could be used as a recommended method to calculate the surrounding rock pressure. It is safer than the other two methods and relatively suitable. Hence, the calculation results of the surrounding rock are partially conservative. As the pressure value of the surrounding rock calculated by the modified general theory is 1.66 times larger than the field monitoring value, the following points are discussed plus the limitations and shortcomings of the resulting work.

Application Effect of Examples
(1) The calculation results of surrounding rock pressure obtained by normative method and Protodyakonov's theory are generally lower than 50% ( (1) From the aspects of geology, rock mass correction can be done separately, considering the boundary conditions (different head and different states of in situ stress, deformation, and bearing load), and more abundant weak intercalation characteristics (such as the thickness of the interlayer, clip layer upon layer, mezzanine type, quantity, distribution, and extended range), the surrounding rock pressure calculation method of correction can be perfected for solving a more comprehensive range of practical engineering problems (2) From the engineering aspect of the correction, the form and stiffness of supporting structures and the sequence of controlled blasting and excavation in construction can be investigated based on the cave factors (such as the shape of a three-center vault, semicircular vault, and cut round vault; depth, scale, thickness of overburden layer, and axis direction of the tunnel). These are the breakthrough points for future results

Applicability of the MPM.
According to the survey report, the method applies only to the rock mass class, which is IV level, weak formation period of V level, and single weak intercalation in situ stress. The phreatic aquifer change relative to the base of deep-buried tunnel mode is not large. Thus, the study determines the water head and in situ stress state affected by the action of seepage under the weak intercalations of horseshoe-shaped deep-buried tunnels.
As the results of the numerical simulation of the fitting modified coefficient are more influenced by engineering and geological effects (such as buried tunnel depth, span, sectional shape, the trend of the tunnel, and lithology of the surrounding rock combined with the shape), this paper, considering the weak interlayer's angle and distance modified coefficient, mainly provides a train of thought. However, the application layer of universal significance needs further discussion.

Conclusions
For the surrounding rock of the tunnel under the seepage condition of weak interlayer, this work proporsed MPM, which is a novel calculation method of tunnel surrounding rock pressure considering the influence of interlayer angle and position.
(1) In the process of determining the modified coefficient, the corresponding surrounding rock pressure under the influence of the interlayer's angle presents the distribution rule of the Gaussian state, and the change of the surrounding rock pressure under the influence of the interlayer position presents the distribution rule of the exponential function. The maximum value is reached in the interlayer's angle range of 40°to 60°. With the increase in the distance 12 Geofluids between the sandwich vault, the surrounding rock pressure decreases, and the decrease is more obvious when the distance is less than 0.2 D. Therefore, the variation of the vertical surrounding rock pressure of sandwich vault with different spacing presents the distribution form of the exponential function (2) The surrounding rock pressure monitoring value in engineering is compared with the surrounding rock pressure calculated by three schemes, namely, the code method, Protodyakonov's theory method, and the MPM of Protodyakonov's theory. It is found that the result obtained from the surrounding rock pressure calculated with the MPM is the most suitable (3) The analysis results show that the surrounding rock pressure calculated by the code method and the Protodyakonov's theory method is too small to reach the support strength of the surrounding rock of the tunnel with a weak interlayer, and the relative safety risk is significant. The pressure values calculated by the modified formula are greater than those monitored in the field, and the values calculated by the two schemes are close. On the premise of satisfying the support strength of the right vault and the left vault, it is more appropriate to conservatively consider the calculation of the support strength of the excellent vault

Data Availability
The data presented in this study can be available on request from the first author.

Conflicts of Interest
The authors declare that they have no conflicts of interest regarding the publication of this paper.