Analysis of the Limit Support Pressure of a Shallow Shield Tunnel in Sandy Soil Considering the Inﬂuence of Seepage

: The objective was to optimize the existing solution for the limit support pressure of a tunnel face. Firstly, based on the numerical simulation results, the existing three-dimensional analytical solution for pore water pressure distribution is expanded to a three-dimensional solution considering the pore water pressure distribution in the upper formation behind the tunnel face. Then, according to the results of physical model tests, a failure model considering the failure range in the upper formation behind the tunnel face is established, and the newly established three-dimensional solution for pore water pressure is introduced into the model, and then the limit e ﬀ ective support pressure of the tunnel face considering seepage is obtained by the method of soil–water joint calculation. Finally, the calculation results in this paper are compared with the experimental results, numerical simulation results and existing theoretical solutions. The major ﬁndings are as follows. The distribution of pore water pressure in the front and back strata above the tunnel face is basically symmetrical. The limit e ﬀ ective support pressure of the tunnel face will increase linearly with an increase in the hydraulic head di ﬀ erence between the tunnel face and the ground surface. The calculated results of the new limit equilibrium theory are obviously larger than those of the existing theory and numerical simulation and closer to the results of the physical model tests. Therefore, the new limit equilibrium model can better predict the limit e ﬀ ective support pressure of the tunnel face considering seepage and provide a reference for actual projects.


Introduction
It is easy to form a seepage field near a tunnel face when the shield is excavated in a saturated sandy soil layer, which is not conducive to the stability of the tunnel face. The determination of the limit support pressure of the tunnel face under seepage conditions is very important to maintain the stability of the tunnel face in saturated sandy soil. Physical model tests [1][2][3] and numerical simulations [4][5][6] can be used to study the limit support pressure of the tunnel face under seepage conditions, but these methods are often limited to specific seepage fields and soil physical and mechanical parameters. However, theoretical calculation can often predict the limit support pressure of the tunnel face under different seepage conditions and soil physical and mechanical parameters. Theoretical calculation of the limit support pressure of the tunnel face under seepage conditions mainly includes limit analysis [7][8][9] and limit equilibrium analysis, and the limit equilibrium method is simple and easy to understand and more convenient to apply to actual projects.

Three-Dimensional Analytical Solution for Pore Water Pressure Distribution
The distribution of pore water pressure is key to study of the limit support pressure of the tunnel face considering seepage by the limit equilibrium method. Figure 1 shows the seepage model of the tunnel face, where, C is the cover thickness of the tunnel; D is the tunnel diameter; h t is the hydraulic head of the tunnel face; h 0 is the hydraulic head of the ground surface; and ∆ is the hydraulic head difference between the tunnel face and the ground surface. Taking the bottom of the tunnel as the origin of the coordinates, the tunnel excavation direction is the positive x-axis direction, and the vertical upward direction is the positive z-axis direction.
Cao et al. [17] obtained a three-dimensional solution for pore water pressure distribution in front of the tunnel face through Equations (1) and (2), but did not consider the distribution of pore water pressure behind the tunnel face. Here, u is the pore water pressure; γ w is the unit weight of water; H is the equivalent height of the tunnel in this paper, H = D.
u(x, y, z) = γ w ( ) In order to expand the above three-dimensional solution without considering the pore water pressure behind the tunnel face to the three-dimensional solution considering the pore water pressure behind the tunnel face, a numerical model, as shown in Figure 2, is established by FLAC 3D software, considering the symmetry of the tunnel. Here, the diameter of the tunnel D is 6m, the thickness of the covering soil C is 12m, the top surface is free and the hydraulic head pressure is fixed at 1m, the bottom surface and the side wall of the tunnel are fixed and impermeable, the surrounding normal direction is fixed and impermeable, the tunnel face is fixed in the normal direction, the hydraulic head difference between the tunnel face and the ground surface is fixed at 0.5, 1.0 or 2.0 D , and the friction angle of the soil ϕ is 37°, 31° or 25°. Then, the distribution of pore water pressure is obtained by three-dimensional seepage calculation. In order to expand the above three-dimensional solution without considering the pore water pressure behind the tunnel face to the three-dimensional solution considering the pore water pressure behind the tunnel face, a numerical model, as shown in Figure 2, is established by FLAC 3D software, considering the symmetry of the tunnel. Here, the diameter of the tunnel D is 6 m, the thickness of the covering soil C is 12 m, the top surface is free and the hydraulic head pressure is fixed at 1 m, the bottom surface and the side wall of the tunnel are fixed and impermeable, the surrounding normal direction is fixed and impermeable, the tunnel face is fixed in the normal direction, the hydraulic head difference between the tunnel face and the ground surface is fixed at 0.5, 1.0 or 2.0 D, and the friction angle of the soil ϕ is 37 • , 31 • or 25 • . Then, the distribution of pore water pressure is obtained by three-dimensional seepage calculation.   Figure 3 shows the cloud image of the pore water pressure distribution obtained by numerical simulation. According to the numerical simulation results, it is found that the pore water pressure in the front and back of the tunnel face is symmetrical. Therefore, this paper expands the three-dimensional solution for the pore water pressure distribution of Equations (1) and (2), which did not consider the pore water pressure distribution behind the tunnel face, into Equation (3), which comprehensively considers the pore water pressure distribution directly in front of the tunnel face, the upper formation in front of the tunnel face and the upper formation behind the tunnel face.  Figure 3 shows the cloud image of the pore water pressure distribution obtained by numerical simulation. According to the numerical simulation results, it is found that the pore water pressure in the front and back of the tunnel face is symmetrical. Therefore, this paper expands the three-dimensional solution for the pore water pressure distribution of Equations (1) and (2), which did not consider the pore water pressure distribution behind the tunnel face, into Equation (3), which comprehensively considers the pore water pressure distribution directly in front of the tunnel face, the upper formation in front of the tunnel face and the upper formation behind the tunnel face.  Figure 3 shows the cloud image of the pore water pressure distribution obtained by numerical simulation. According to the numerical simulation results, it is found that the pore water pressure in the front and back of the tunnel face is symmetrical. Therefore, this paper expands the three-dimensional solution for the pore water pressure distribution of Equations (1) and (2), which did not consider the pore water pressure distribution behind the tunnel face, into Equation (3), which comprehensively considers the pore water pressure distribution directly in front of the tunnel face, the upper formation in front of the tunnel face and the upper formation behind the tunnel face.
(3) Figure 4 shows the distribution of pore water pressure on the vertical symmetrical plane of the tunnel obtained by numerical simulation and theoretical calculation in this paper. It can be seen from Figure 4 that the pore water pressure distribution law above the tunnel face obtained by the two methods is the same (the maximum relative error is less than 15%), both of which take the tunnel face as the symmetrical plane. Directly above the tunnel face, the pore water pressure gradient calculated by theory is obviously larger than that of the numerical simulation, which will make the seepage force in the theoretical calculation larger and make the limit support pressure of the tunnel face larger, and give the calculation result more safety. As the ground directly in front of the tunnel face is mainly affected by the horizontal seepage force [18], the pore water pressure on the central axis of the tunnel represents the pore water pressure on the ground directly in front of the tunnel face. According to the distribution of pore water pressure on the central axis of the tunnel in Figure 4, the theoretical calculation is basically consistent with the numerical simulation results. Based on the above analysis, it is safe to apply the pore water pressure calculated according to Equation (3) to calculation of the limit support pressure of the tunnel face.

Limit Support Pressure of the Tunnel Face Considering the Influence of Seepage
In this section, firstly, according to the physical model tests of Mi and Xiang [3], a failure model considering the failure range behind the tunnel face is established, then the newly established three-dimensional solution for the pore water pressure is introduced into the model, and, finally, the limit support pressure of the tunnel face is solved by the limit equilibrium method.

Limit Support Pressure of the Tunnel Face Considering the Influence of Seepage
In this section, firstly, according to the physical model tests of Mi and Xiang [3], a failure model considering the failure range behind the tunnel face is established, then the newly established Symmetry 2020, 12, 1023 6 of 13 three-dimensional solution for the pore water pressure is introduced into the model, and, finally, the limit support pressure of the tunnel face is solved by the limit equilibrium method. Figure 5 shows the test model of excavation-seepage instability of the shield tunnel established by Mi and Xiang [3]. Different seepage fields are generated by regulating the pump switch. By rotating the runner of the tunnel model, the support plate moves backward gradually until the ground is unstable and damaged. In the process of the support plate moving backward, the ground settlement and the support pressure of the tunnel face are monitored. Based on analysis of the monitoring data, the ground collapse range and the limit support pressure of the tunnel face are obtained.    Figure 6 shows the ground collapse range of the vertical symmetry plane in the center of the tunnel and the plane in front of the tunnel face under different seepage conditions based on the physical model tests of Mi and Xiang [3]. Here, ∆ is the hydraulic head difference between the tunnel face and the ground surface and D is the tunnel diameter. It can be seen from Figure 6 that when the support pressure of the tunnel face is not enough, the ground in front and behind the tunnel face will lose stability within a certain range. The collapse range boundary in front of the tunnel face extends from the bottom of the tunnel to the upper part of the tunnel at a certain angle and reaches a certain height, and then extends vertically to the ground surface. The collapse range boundary behind the tunnel face extends from the top of the tunnel to the rear and up at a certain angle and reaches a certain height, then extends vertically upward to the ground surface. The horizontal collapse range boundary extends from the bottom of the tunnel to the upper side at a certain angle and reaches a certain height, then extends to the ground surface vertically.

Establishment of the Failure Model
According to the ground collapse range in Figure 6, considering the symmetry of the tunnel, a failure model considering the influence of seepage, as shown in Figure 7, is established. This model consists of three parts: a pyramid, a prismoid and a prism from bottom to top. Here, H is the equivalent height of the tunnel (in this paper, H = D); C is the cover thickness of the tunnel; h 0 is the hydraulic head of the ground surface; h t is the hydraulic head of the tunnel face; q is the additional pressure on Symmetry 2020, 12, 1023 7 of 13 the ground surface; ω 1 is the break angle in front of the tunnel face; ω 2 is the horizontal break angle; ω 3 is the break angle behind the tunnel face.  According to the ground collapse range in Figure 6, considering the symmetry of the tunnel, a failure model considering the influence of seepage, as shown in Figure 7, is established. This model consists of three parts: a pyramid, a prismoid and a prism from bottom to top. Here, H is the equivalent height of the tunnel (in this paper, H D = ); C is the cover thickness of the tunnel; h 0 is the hydraulic head of the ground surface; t h is the hydraulic head of the tunnel face; q is the additional pressure on the ground surface; 1 ω is the break angle in front of the tunnel face; 2 ω is the horizontal break angle; 3 ω is the break angle behind the tunnel face.

Calculation of the Limit Support Pressure of the Tunnel Face
It is assumed that the soil in front of and behind the tunnel face is uniform and obeys the Mohr-Coulomb failure criterion. Considering that the divided calculation of water and soil needs to solve the seepage force, the process of solution is rather complicated, so this paper uses the method of soil-water joint calculation to analyze the stress of the failure model. Figure 8 shows the stress state of the prism micro-element, where, S dp 1 and S dp 2 are the

Calculation of the Limit Support Pressure of the Tunnel Face
It is assumed that the soil in front of and behind the tunnel face is uniform and obeys the Mohr-Coulomb failure criterion. Considering that the divided calculation of water and soil needs to solve the seepage force, the process of solution is rather complicated, so this paper uses the method of soil-water joint calculation to analyze the stress of the failure model. Figure 8 shows the stress state of the prism micro-element, where, dp S1 and dp S2 are the water pressure on the two sides of the prism micro-element, respectively; dN 1 and dN 2 are the normal effective pressure on the two sides of the prism micro-element, respectively. From the static equilibrium conditions of the prism micro-elements in the vertical direction, it can be obtained that dV + dp V + dG = 2dS 1 + 2dS 2 (4) where dV is the vertical effective pressure increment of the prism micro-element; where σ z is the average vertical effective compressive stress on the z plane of the prism; dp V is the vertical water pressure increment of the prism micro-element; where u(z) is the average water pressure on the z plane of the prism; dG is the gravity of the prism micro-element, where γ sat is the saturated weight of the soil; dS 1 is the friction resistance on the front and back sides of the prism micro-element, dS 1 = (c + kσ z tan ϕ)(4H tan ω 2 )dz (8) where c is the cohesion of the soil; k is the pressure coefficient (in this paper, k = 1 − sin ϕ); ϕ is the friction angle of soil; dS 2 is the friction resistance on the left and right sides of the prism micro-element: Substituting Equations (5)-(9) into Equation (4) and combining the surface boundary conditions σ z | z=H+C = q and u(H + C) = γ w (h 0 − H − C), the vertical effective compressive stress of the prism σ z can be obtained, u(x, y, z) dxdy (10) where R is the area-perimeter ratio of the cross-section of the prism, where R is the area-perimeter ratio of the cross-section of the prism,  In order to simplify the calculation, it is assumed that the vertical effective compressive stress in the prismoid continues the change rule in the prism; that is, within the range H z H C ≤ ≤ + , the In order to simplify the calculation, it is assumed that the vertical effective compressive stress in the prismoid continues the change rule in the prism; that is, within the range H ≤ z ≤ H + C, the vertical effective compressive stress σ z (z) can be obtained by Equation (10). When z = H, the vertical effective compressive stress of the interface between the prismoid and the pyramid can be obtained, that is, σ p = σ z (H). Figure 9 shows the stress state of the whole pyramid. According to the static equilibrium conditions in the x and z directions of the pyramid, it can be obtained that T + p T + 2S s sin α + S f sin ω 1 = N f cos ω 1 + p f cos ω 1 (12) V + p V + G = 2S s cos α cos ω 2 + S f cos ω 1 + 2N s sin ω 2 + 2p s sin ω 2 + N f sin ω 1 + p f sin ω 1 (13) where N f is the normal effective pressure on the slope in front of the pyramid; S f is the friction resistance on the slope in front of the pyramid. The relationship between them is as follows: Substitute Equation (14) into Equations (12) and (13), respectively, and eliminate N f simultaneously to obtain the effective support pressure of the tunnel face T , where α is the angle between two side edges of the pyramid, α = arctan(tan ω 1 cos ω 2 ) V is the vertical effective pressure on the top surface of the pyramid, It can be seen in Figure 3 that the break angle behind the tunnel face, 3 ω , does not change much and can be approximately taken as σ is the limit effective support pressure stress of the tunnel face considering the influence of seepage. Figure 9. The stress state of the whole pyramid (soil-water joint calculation). Figure 9. The stress state of the whole pyramid (soil-water joint calculation).

Comparative Analysis of the Limit Effective Support Pressure of the Tunnel Face Obtained by Different Methods
In order to verify the correctness of the theoretical solution for the limit equilibrium in this paper, the physical model test parameters of Mi and Xiang [3] are used for calculation, and the calculation results are compared with the physical model test results [3], the existing theoretical calculation results of the limit equilibrium [11] and the numerical simulation results [3], respectively. The diameter of the tunnel D is 0.3m; the covering thickness C is 0.6m; the hydraulic head of ground surface h 0 is 1.0m; the hydraulic head of the tunnel face t h is 1.0, 0.9, 0.8, 0.7 or 0.5m; the saturated weight of soil sat γ is 21.2 kN m 3 / ; the friction angle of the soil ϕ is 37°; the cohesion of the soil c is 0; the additional pressure on the ground surface q is 0. Figure 11 shows the limit effective support pressure of the tunnel face under different seepage conditions obtained by different methods, where t =h h 0 Δ − is the hydraulic head difference between the tunnel face and the ground surface. It can be seen that there is a linear relationship between the limit effective support pressure T lim ' σ and the hydraulic head difference Δ , that is, T lim ' σ will increase linearly with the increase in Δ . The theoretical calculation results of the limit equilibrium in this paper are obviously larger than those of Perazzelli et al. [11] and the numerical simulation results of Mi and Xiang [3], which are closer to the physical model test results of Mi and Xiang [3]. The larger the hydraulic head difference Δ is, the closer the prediction result calculated in this paper is to the actual test result. The theoretical model of limit equilibrium established in this paper can better predict the limit effective support pressure of the tunnel face considering seepage and can provide a reference for practical engineering.

Comparative Analysis of the Limit Effective Support Pressure of the Tunnel Face Obtained by Different Methods
In order to verify the correctness of the theoretical solution for the limit equilibrium in this paper, the physical model test parameters of Mi and Xiang [3] are used for calculation, and the calculation results are compared with the physical model test results [3], the existing theoretical calculation results of the limit equilibrium [11] and the numerical simulation results [3], respectively. The diameter of the tunnel D is 0.3 m; the covering thickness C is 0.6 m; the hydraulic head of ground surface h 0 is 1.0 m; the hydraulic head of the tunnel face h t is 1.0, 0.9, 0.8, 0.7 or 0.5 m; the saturated weight of soil γ sat is 21.2 kN/m 3 ; the friction angle of the soil ϕ is 37 • ; the cohesion of the soil c is 0; the additional pressure on the ground surface q is 0. Figure 11 shows the limit effective support pressure of the tunnel face under different seepage conditions obtained by different methods, where ∆= h 0 − h t is the hydraulic head difference between the tunnel face and the ground surface. It can be seen that there is a linear relationship between the limit effective support pressure σ Tlim and the hydraulic head difference ∆, that is, σ Tlim will increase linearly with the increase in ∆. The theoretical calculation results of the limit equilibrium in this paper are obviously larger than those of Perazzelli et al. [11] and the numerical simulation results of Mi and Xiang [3], which are closer to the physical model test results of Mi and Xiang [3]. The larger the hydraulic head difference ∆ is, the closer the prediction result calculated in this paper is to the actual test result. The theoretical model of limit equilibrium established in this paper can better predict the limit effective support pressure of the tunnel face considering seepage and can provide a reference for practical engineering.

Conclusions
The objective was to study the limit support pressure of a tunnel face in saturated sandy soil under seepage conditions. Firstly, according to the numerical simulation results, the existing three-dimensional analytical solution for the pore water pressure distribution not considering the distribution of pore water pressure behind the tunnel face is extended to the three-dimensional solution considering the distribution behind the tunnel face. Secondly, based on the results of the physical model tests, a failure model considering the failure range behind the tunnel face was established, and the newly established three-dimensional pore water pressure solution was introduced into the model to obtain the limit support pressure of the tunnel face with the method of soil-water joint calculation. Finally, the calculation results were compared with the experimental results, numerical simulation results and existing theoretical solutions. The main conclusions are as follows.
(1) The distribution of pore water pressure in the front and back strata above the tunnel face is basically symmetrical.
(2) When the support pressure of the tunnel face is not enough, both the ground in front and that behind the tunnel face will lose stability within a certain range. Therefore, when establishing the failure model, the failure range of the ground directly in front of the tunnel face, the upper formation in front of the tunnel face, and the upper formation behind the tunnel face should be considered comprehensively.
(3) The limit effective support pressure of the tunnel face will increase linearly with the increase in the hydraulic head difference between the tunnel face and the ground surface.
(4) The calculated results in the new limit equilibrium theory are obviously larger than those in the existing theory and numerical simulation and are closer to the results of the physical model tests. Therefore, the new limit equilibrium model can better predict the limit effective support pressure of the tunnel face considering seepage and provide a reference for practical engineering.
It is worth noting that the ground studied in this paper is sandy soil, and the limit effective support pressure of the tunnel face for cohesive soil needs further study.
Author Contributions: B.M. and Y.X. contributed equally to this work. All authors have read and agreed to the published version of the manuscript.

Conclusions
The objective was to study the limit support pressure of a tunnel face in saturated sandy soil under seepage conditions. Firstly, according to the numerical simulation results, the existing three-dimensional analytical solution for the pore water pressure distribution not considering the distribution of pore water pressure behind the tunnel face is extended to the three-dimensional solution considering the distribution behind the tunnel face. Secondly, based on the results of the physical model tests, a failure model considering the failure range behind the tunnel face was established, and the newly established three-dimensional pore water pressure solution was introduced into the model to obtain the limit support pressure of the tunnel face with the method of soil-water joint calculation. Finally, the calculation results were compared with the experimental results, numerical simulation results and existing theoretical solutions. The main conclusions are as follows.
(1) The distribution of pore water pressure in the front and back strata above the tunnel face is basically symmetrical.
(2) When the support pressure of the tunnel face is not enough, both the ground in front and that behind the tunnel face will lose stability within a certain range. Therefore, when establishing the failure model, the failure range of the ground directly in front of the tunnel face, the upper formation in front of the tunnel face, and the upper formation behind the tunnel face should be considered comprehensively.
(3) The limit effective support pressure of the tunnel face will increase linearly with the increase in the hydraulic head difference between the tunnel face and the ground surface.
(4) The calculated results in the new limit equilibrium theory are obviously larger than those in the existing theory and numerical simulation and are closer to the results of the physical model tests. Therefore, the new limit equilibrium model can better predict the limit effective support pressure of the tunnel face considering seepage and provide a reference for practical engineering.
It is worth noting that the ground studied in this paper is sandy soil, and the limit effective support pressure of the tunnel face for cohesive soil needs further study.