Working Face Stability of Box Shield Tunneling under Non-Uniform Support Pressure

: This paper proposes a theoretical model for the stability analysis of a box tunnel face in non-cohesive soils considering the uneven distribution of support pressure caused by multiple cutter heads and screw conveyors. The support pressure distribution on the tunnel face is concave. Accordingly, the failure mechanism is composed of a prism and a wedge, both including three blocks. The relatively smaller support pressure acting on the middle blocks lead to the tendency of slide. Assuming that the support pressure acting on the side blocks is obtained using the active earth pressure coefﬁcient, the support pressure acting on block II can be achieved by limit equilibrium analysis considering the interactions between the blocks. The inﬂuences of strength parameters and geometric parameters on the tunnel face stability are discussed in the parametric analysis. For comparison, numerical analysis is conducted in the commercial software OptumG3. It is found that the results given by the proposed model agree well with those from the numerical model. Therefore, the rationality of the proposed model in predicting the collapse geometry is veriﬁed.


Introduction
Rectangular tunnels are widely used in the construction of urban underground projects such as subways, which have advantages such as high space utilization and little impacts on ground traffic, etc.As a trenchless technique, box jacking is often adopted in urban construction.The stability of the tunnel face is one of the most important factors in tunneling.It is necessary to ensure face stability during the excavation of a tunnel.The cross section of a rectangular tunnel tends to be large, and the support pressure is unevenly distributed, both of which greatly increase the risk of tunnel face instability.
The methods to study the stability of a tunnel boring face mainly include theoretical analysis , model testing, and numerical simulation [31][32][33][34][35][36][37][38][39][40][41][42].Theoretical analysis presents a clear force model, and the force analysis is scientific and accurate, which is ideal for engineering applications.Commonly used theoretical methods mainly include upper limit analysis and limit equilibrium analysis.
Limit analysis methods [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15] construct a damage maneuver mechanism according to the plastic flow law for soil in the critical state and use an energy balance to derive the judgment criteria for the destabilization of the excavation surface.Davis et al. [1] studied the active failure mode of shield tunneling in undrained clay soil under 2D plane strain conditions and proposed the calculation method of the ultimate support pressure and stability coefficient.Leca & Dormieux [2] proposed a 3D cone model, constructed a translational maneuver field that satisfies the motion permission, and derived the upper limit solution of the limit support pressure.According to Leca & Dormieux [2], Soubra et al. [3] improved the cone model, and a 3D failure mechanism composed of a rigid cone and a logarithmic spiral shear body was proposed.Mollon et al. [4][5][6] focused on the elliptical failure area and the circle of the cone model and constructed the 3D logarithmic spiral failure mechanism for the first time using spatial discretization technology to achieve the full-section instability of the failure mechanism.
Based on the spatial discretization technology, Pan et al. [7][8][9] extended an advanced 3D failure mechanism into the face stability analysis in anisotropic and nonhomogeneous soils, achieved to generate a failure mechanism for a non-circular tunnel face to estimate the safety factor of a tunnel face.Senent et al. [10] used both planar triangular blocks and spiral blocks to construct an active instability maneuvering model and studied the influence of the tunnel excavation step on the stability of the driving face.Chen et al. [11] constructed the passive instability failure mechanism of a shallow tunnel driving face and considered the pore water, the influence of pressure, and an uneven overload on the stability of the driving face.Liu et al. [12,13] developed 2D and 3D kinematically admissible mechanisms, and the composite mechanism comprises rotational and translational blocks for analyzing the passive face stability during shield tunneling.Liu et al. [14] considered rectangular cross-sectional shield tunneling and derived the upper-bound solution of support pressure applied to maintain the face stability.Di et al. [15] proposed an analytical model of face stability for a tunnel under anisotropic seepage based on the limit analysis upper-bound method.
However, the maneuver model construction rule in limit analysis is influenced by the boundary conditions, and the construction model is complicated, which reduces the applications of limit analysis in engineering.In contrast, limit equilibrium analysis is more adaptable to the boundary conditions and easier to analyze by constructing a rigid-body limit equilibrium model.Jancsecz et al. [16] established a 3D wedge model based on the limit equilibrium method using the silo theory and proposed a calculation method for the limit support pressure.Broere [17] extended the calculation model to layered soil conditions and analyzed the earth pressure and mud-water balance shield construction limit support pressure.Anagnostou [18] studied the influence of the soil arching effect on the stability of the tunnel excavation surface, optimized the calculation model, and put forward the expression of support pressure.Xu et al. [19] considered the face stability of rectangular pipe jacking under the condition of gravel formation and used the trapezoidal wedge model to derive the expression of the active failure limit support pressure of the excavation face.Cheng et al. [20] assumed that the failure area is an arc-shaped wedge and derived the active limit support force calculation formula for the lower excavation face under dry sand conditions based on the silo theory.Lin et al. [21], based on the limit equilibrium method, proposed a logarithmic spiral prism model considering the unsupported length to analyze the stability of the shallow tunnel driving face in sand.Vu et al. [22] extended the calculation model to stratified soil conditions and derived the ultimate support pressure for the earth pressure and mud-water balance shield construction.Wei et al. [23,24] modified the 3D wedge model by assuming the sliding block shape to be a trapezoidal wedge when the excavation surface is destabilized and obtained better calculation results for the support force.Sun et al. [25][26][27] investigated the analytical solutions for the responses of tunnels excavated in rock masses and verified the proposed solutions by numerical simulations.The existing theoretical studies mainly focus on the stability of tunnel boring faces in circular cross-sections, and there is less research on rectangular cross-sections, and there are especially few theoretical studies on the uneven distribution of support pressure in the boring face.Liu et al. [28] focused on the situation in which the working face chamber is exposed to free air due to obstruction removal.In this situation, the support pressure distribution on the working face is characterized by zero in the middle.At this point, the soil in the middle of the working face has a great risk of destabilization because, without support pressure, the soil is likely to have an insufficient shear strength to maintain stability.
This paper investigates the stability of the box tunnel working face.An improved wedge-prism model is established for the uneven distribution of support pressure at the tunnel face due to the large cross-section and the characteristic behavior of the cutter head.In the proposed model, the wedge and prim blocks are each divided into three blocks according to the characteristics of the support pressure distribution.The support pressure distribution on the working face is characterized as concave.Therefore, the soil in the middle with a small support pressure has a downward sliding tendency relative to the soil on both sides because the support pressure is not sufficient to maintain stability, and even active damage may occur.Assuming that the support pressure on both sides of the tunnel face is known and taking into account the interactions between the blocks, the smaller support pressure in the middle can be calculated by force analysis.The rationality of the proposed model is investigated, and the results given by the proposed method are compared with those obtained by numerical simulation.

Distribution Characteristics of Support Pressure in Box Tunnels
As shown in Figure 1, due to its large section, the boring machine often uses multiple cutter heads and double-screw conveyors to improve the efficiency of excavation and soil removal.The gates of the conveyor are marked with yellow circles in Figure 1a-e.When the tunnel is driven with cutter heads rotating, the cut-down soil enters the chamber through the conveyor gates.A double-screw conveyor mounted behind the bulkhead then discharges the cut soil.(Figure 1f).Since the opening rate of the cutter disc is large, earth pressure in the chamber is highly possible to be fluctuating and unevenly distributed.tunnel face due to the large cross-section and the characteristic behavior of the cutter head.
In the proposed model, the wedge and prim blocks are each divided into three blocks according to the characteristics of the support pressure distribution.The support pressure distribution on the working face is characterized as concave.Therefore, the soil in the middle with a small support pressure has a downward sliding tendency relative to the soil on both sides because the support pressure is not sufficient to maintain stability, and even active damage may occur.Assuming that the support pressure on both sides of the tunnel face is known and taking into account the interactions between the blocks, the smaller support pressure in the middle can be calculated by force analysis.The rationality of the proposed model is investigated, and the results given by the proposed method are compared with those obtained by numerical simulation.

Distribution Characteristics of Support Pressure in Box Tunnels
As shown in Figure 1, due to its large section, the boring machine often uses multiple cutter heads and double-screw conveyors to improve the efficiency of excavation and soil removal.The gates of the conveyor are marked with yellow circles in Figure 1a-e.When the tunnel is driven with cutter heads rotating, the cut-down soil enters the chamber through the conveyor gates.A double-screw conveyor mounted behind the bulkhead then discharges the cut soil.(Figure 1f).Since the opening rate of the cutter disc is large, earth pressure in the chamber is highly possible to be fluctuating and unevenly distributed.Based on the characteristics of the equipment, such as cutter and screw conveyors, the cut soil is conveyed out from the screw conveyor gates, which are at the bottom of the working face.As shown in Figure 2a, with the discharge of soil in the area near the outlet of the double-screw conveyors, the soil in the upper parts and the surrounding areas collapses.The discharge of the screw conveyors makes the air pressure in the central area Based on the characteristics of the equipment, such as cutter and screw conveyors, the cut soil is conveyed out from the screw conveyor gates, which are at the bottom of the working face.As shown in Figure 2a, with the discharge of soil in the area near the outlet of the double-screw conveyors, the soil in the upper parts and the surrounding areas collapses.The discharge of the screw conveyors makes the air pressure in the central area small and unstable.Therefore, the central area is significantly impacted, while the neighboring two parts are almost unaffected.The residual soil on both sides serves as a stable support for the working face, and the support force is defined as S 1 .The air pressure in the middle area is low and unstable, forming a pressure fluctuation zone.The support force is defined as S 2 , which is smaller than S 1 .
Based on the characteristics of the equipment, such as cutter and screw conveyors, the cut soil is conveyed out from the screw conveyor gates, which are at the bottom of the working face.As shown in Figure 2a, with the discharge of soil in the area near the outlet of the double-screw conveyors, the soil in the upper parts and the surrounding areas collapses.The discharge of the screw conveyors makes the air pressure in the central area small and unstable.Therefore, the central area is significantly impacted, while the neighboring two parts are almost unaffected.The residual soil on both sides serves as a stable support for the working face, and the support force is defined as 1 S .The air pressure in the middle area is low and unstable, forming a pressure fluctuation zone.The support force is defined as 2 S , which is smaller than 1 S .
Therefore, the working face is non-uniformly supported, and the support pressure is large on both sides and small in the middle.As shown in Figure 2b, divide the working face into three parts, and simplify the central irregular area of the unstable support pressure into a rectangle (the width of which is B , and the support force is denoted as 2 S ).
The two adjacent areas of stable support pressure are also simplified to rectangles with support force 1 S .In Figure 3, a tunnel excavation of height D and width L under a cover depth C in homogenous soil is considered.The soil obeys the Mohr-Coulomb failure condi- tion, defined by its effective cohesion ( c ′ ), effective angle of internal friction ( ϕ ′ ), and the total unit weight (γ ).At the working face, the middle part with a width of B is the area of unstable support pressure.In other words, this part (II) is not subject to the full support force.However, the two parts of the working face on both sides (I and III) are entirely supported.The ground surcharge s σ is not taken into account; thus, s σ is as- sumed to be 0 kPa.Therefore, the working face is non-uniformly supported, and the support pressure is large on both sides and small in the middle.As shown in Figure 2b, divide the working face into three parts, and simplify the central irregular area of the unstable support pressure into a rectangle (the width of which is B, and the support force is denoted as S 2 ).The two adjacent areas of stable support pressure are also simplified to rectangles with support force S 1 .
In Figure 3, a tunnel excavation of height D and width L under a cover depth C in homogenous soil is considered.The soil obeys the Mohr-Coulomb failure condition, defined by its effective cohesion (c ), effective angle of internal friction (ϕ ), and the total unit weight (γ).At the working face, the middle part with a width of B is the area of unstable support pressure.In other words, this part (II) is not subject to the full support force.However, the two parts of the working face on both sides (I and III) are entirely supported.The ground surcharge σ s is not taken into account; thus, σ s is assumed to be 0 kPa.

Modified Model
Face stability in homogeneous ground can be assessed by considering the mechanism which consists of a wedge and a prism that extends from the tunnel crown to the surface.Figure 4 shows the profile of the proposed model.The width, height, and depth of the

Modified Model
Face stability in homogeneous ground can be assessed by considering the mechanism which consists of a wedge and a prism that extends from the tunnel crown to the surface.Figure 4 shows the profile of the proposed model.The width, height, and depth of the tunnel, respectively, are L, D, and C. In the proposed model, the width of the partially supported area is B. Due to symmetry, only block I, block II, and their interaction are necessary to analyze.There are interaction frictions, T 12 and T 21 , and normal forces, N 12 and N 21 , on the contact surface of block I and block II.When block II has a downward trend due to the reduction in the support force, it will produce downward friction T 12 on block I which will cause the risk of instability.Therefore, the analysis of the limit equilibrium state of the two blocks and the interaction forces are used so that the minimum support force S 2 to keep block I stable can be deduced when support force S 1 is known.

Modified Model
Face stability in homogeneous ground can be assessed by considering the mechanism which consists of a wedge and a prism that extends from the tunnel crown to the surface.Figure 4 shows the profile of the proposed model.The width, height, and depth of the tunnel, respectively, are L, D , and C .In the proposed model, the width of the par- tially supported area is B .Due to symmetry, only block I, block II, and their interaction are necessary to analyze.There are interaction frictions, 12 T and 21 T , and normal forces,    The resultant vertical force of the prism 1 P on the top surface o a e is given by: The resultant vertical force of the prism P 1 on the top surface oae is given by: where θ denotes the wedge angle between slanted failure plane and the horizontal.The resultant normal force N 1 and the shear force T 1 along the failure surface bde are given by: The resultant normal force N 3 and the shear force T 3 along the failure surface abe are given by: where α denotes the bottom corner of isosceles trapezoid on the top of the wedge and can be solved as α = arctan 2D cot θ L−B ; K 0 is a lateral stress coefficient and can be solved as K 0 = 1 − sin ϕ ; and the vertical stress σ z along abe changes linearly with depth and can be solved as The support force S 1 is defined as: where s 1 denotes the support pressure acting on surface oabd.
where K a is the active earth pressure coefficient and can be solved as K a = tan 2 π 4 − ϕ 2 .The weight of block I is: where V 1 represents the volume of block I.The equilibrium conditions in the y and z directions lead to: Its solution reads as follows: As is shown in Figure 6, block II is analyzed by the same method.Compared with block I, block II is acted upon by the support force S 2 acting on surface oghd.
( ) As is shown in Figure 6, block II is analyzed by the same method.Compared with block I, block II is acted upon by the support force 2 S acting on surface oghd .The vertical force at the wedge-prism interface oefg reads as follows: where v σ denotes the vertical stress.Considering the contribution of arching, v σ can be calculated as: The vertical force at the wedge-prism interface oe f g reads as follows: where σ v denotes the vertical stress.Considering the contribution of arching, σ v can be calculated as: where σ s denotes the surface load and is set equal to zero; R equals the ratio of the volume of the prism to its side area: R = BD cot θ 2(B+D cot θ) .The resultant normal force N 2 and the shear force T 2 along the failure surface de f h are given by: The support force S 2 is defined as: where s 2 denotes the support pressure acting on surface oghd.
The self-weight of block II is: where V 2 represents the volume of block II.
The resultant normal force and the shear force of the two lateral slip surfaces, oed and g f h, are equal: Taking into account Equation ( 14) and the equilibrium conditions in the y and z directions leads to: Its solution reads as follows: Due to interaction frictions T 12 = T 21 , combining Equation ( 9) and Equation ( 18) makes it possible to eliminate T 12 and T 21 ; the required support force S 2 can now be found from: Due to Equations ( 1)-( 7) and ( 10)-( 15), we obtain from Equation (20): This simplifies to: where Here, the shorthand notation is introduced: As is shown in Equation ( 22), the formula consists of two elements.The first element is related to the weight of the soil.N γ is a coefficient related to γ.The second element is related to the cohesion of the soil.N c is a coefficient related to c .These influencing factors are selected for discussion in the following parametric analysis.

Comparisons of the Support Pressure
In this section, the finite element limit analysis code OptumG3 is used for numerical validation.OptumG3 is favorable for 3D analysis and has the advantage of automatically searching for optimal faults through mathematical optimization [43].The failure mechanisms of tunnels with different strength parameters and geometric parameters have been investigated.The results of the proposed model-the value of the collapse pressure and the failure geometry-are compared to those of a numerical simulation.
Due to symmetry, half of the tunnel is selected for computation.Figure 7 illustrates the numerical model constructed to validate the analytical solution.The dimensions of the model are 30 m × 24 m × 40 m.The tunnel is excavated to half of its length (20 m).The number of grids is 3000.Standard boundary conditions are used.The support pressure s 1 acting on block I is calculated by Equation ( 6).Accordingly, the support pressure s 2 acting on block II is calculated.
ing on block II is calculated.
The cover depth of the tunnel is 6.0 m, i.e.,     The cover depth of the tunnel is 6.0 m, i.e., D = 6.0 m.In addition, the surface surcharge is zero.The partially supported ratio B/L varies from 0.4 to 0.7.For shallowly buried tunnels, the depth ratio C/D varies from 1 to 2. For box tunnels, the width-height ratio L/D varies from 1.5 to 2.0.The soil is set as an elastic-plastic material following the Mohr-Coulomb failure criterion and associated flow rule.Typical non-cohesive soil is considered.The mechanical parameters of the soil are listed in Table 1. Figure 8 indicates the variation in the normalized support pressure with increasing ϕ .Both the proposed model and the numerical model predict that s 2 /γD decreases nonlinearly as ϕ increases.At constant ϕ , s 2 /γD associated with a higher B/L has a larger result.This is because the increase in B/L makes the fully supported area on the working face smaller and the partially supported area larger, thus leading to an increase in the risk of instability.At constant ϕ and B/L, the proposed model predicts a higher s 2 /γD than the numerical model does.The results of the proposed model are more secure for the engineering design.Figure 8a,b shows that s 2 /γD with different depth ratios decreases in a similar pattern with an increasing ϕ .A higher C/D makes a higher s 2 /γD, which can be observed by comparison.
Table 2 shows the variation in the ratio of the two support pressures s 2 /s 1 with an increasing ϕ .It can be seen that s 2 /s 1 decreases in a similar pattern as ϕ changes with s 2 /γD in Figure 8.However, the changes in s 2 /γD in Figure 8 are more pronounced than the corresponding changes in s 2 /s 1 in Table 2.It is known that s 1 is related to ϕ by Equation (6).s 1 decreases with the increase in ϕ , while γD is constant.Under the same conditions, s 2 /s 1 is larger than s 2 /γD, especially when ϕ is greater.Table 2 shows that s 2 /s 1 with different depth ratios decreases in a similar pattern with an increasing ϕ .A higher C/D makes a lower s 2 /s 1 , as can be seen by comparison.This is because s 1 is related to C by Equation (6).A higher C/D makes a higher s 1 .The results of the proposed model are greater than the results of the numerical model.The reliability of the proposed model can be verified.Figure 9 shows the variation in S 2 depending on B/L for a typical ϕ = 25 • and L/D = 2. S 2 grows rapidly with an increasing B/L; the change curve is steep, and the range of change is large.Figure 10 shows the variation in s 2 /γD depending on B/L for a typical ϕ = 25 • .In Figure 10a (C/D = 1), both models predict that s 2 /γD increases with an increasing B/L.The curves are relatively flat with less variation.Comparing Figures 9 and 10, it can be found that under the same conditions, S 2 increases by almost 100%, while s 2 /γD increases by approximately 10%.Due to s 2 = S 2 /BD (Equation ( 15)), the change in s 2 /γD is related to B when γD is constant.When the increase in S 2 is apportioned to the increase in B, the change in s 2 /γD will be relatively small.This is the reason why s 2 /γD is less sensitive to B/L compared to S 2 .In Figures 10 and 11, both s 2 /γD and s 2 /s 1 increase in a similar pattern with an increasing B/L.The curves indicating the results of the proposed model are mostly above the curves of the numerical model results.The comparison indicates that the proposed model agrees well with the numerical model and predicts higher and safer pressures.The error mainly comes from the assumption that the support pressure acting on block I is known and can be obtained by the active earth pressure coefficient K a , which is shown in Equation ( 6).This assumption makes s 1 smaller and s 2 larger.The findings from the proposed model are larger than those from the numerical model.It can also be found that at an identical B/L, the results increase with a higher ratio of L/D.This probably occurs because the increase in the width of the working surface weakens the soil arching effect, making the overlying soil pressure increase, and the area partially supported increases accordingly, and it will be more difficult to maintain soil stability.

Comparisons of the Failure Geometry
Figure 12 compares the collapse geometries obtained with the analytical solution and with the numerical model.It can be seen that θ increases with an increasing ϕ at a constant B/L, while it decreases with an increasing B/L at a constant ϕ .As shown in Figure 12, in the case of a small ϕ value (ϕ= 15 • ), the results of the proposed model deviate slightly from numerical results.When ϕ is beyond this range, the proposed model gives almost the same value as the numerical simulation.The results of the numerical method show that the failure mechanism in front of the tunnel face is wedge-shaped.The shape of the collapse zone obtained from the proposed model agrees well with the result of the numerical model.Therefore, the proposed model is capable of predicting a reasonable collapse zone.

Parametric Analysis
In this part, a typical situation of D = 6 m and L/D = 2 is considered.The partially supported ratio B/L varies from 0.4 to 0.7.For shallowly buried tunnels, the depth ratio C/D varies from 1 to 2, and the total stress concept is accepted for calculating the soil stress.The soil weight was chosen as γ = 19.8kN/m 3 .The frictional angle ϕ varies from 15 • to 45 • .

Comparisons of the Failure Geometry
Figure 12 compares the collapse geometries obtained with the analytical solution and with the numerical model.It can be seen that θ increases with an increasing ϕ′ at a constant / B L, while it decreases with an increasing / B L at a constant ϕ′ .As shown in Figure 12, in the case of a small ϕ′ value ( =15 ϕ °), the results of the proposed model deviate slightly from numerical results.When ϕ′ is beyond this range, the proposed model gives almost the same value as the numerical simulation.The results of the numerical method show that the failure mechanism in front of the tunnel face is wedge-shaped.The shape of the collapse zone obtained from the proposed model agrees well with the result of the numerical model.Therefore, the proposed model is capable of predicting a reasonable collapse zone.Figure 13 shows the variation in N γ depending on ϕ .For C/D= 1 (in Figure 13a), N γ decreases nonlinearly with an increasing ϕ .When ϕ <30 • , N γ increases with an increasing B/L at a constant ϕ , and they become close to each other when ϕ is between 30 • and 35 • .N γ decreases with an increasing B/L at a constant ϕ , which is less than 30 • .In Figure 13b, similar variations in N γ at C/D= 2 can be observed.A comparison between Figure 13a and Figure 13b shows the influence of C/D on N γ .A higher C/D makes a higher coefficient N γ .Figure 14 shows the variations in N c depending on ϕ .For C/D= 1 (in Figure 14a), N c decreases nonlinearly with an increasing ϕ .As B/L becomes smaller, the change in N c becomes more pronounced.When ϕ is relatively small, the discrepancy between N c values associated with different B/L values becomes more significant.Figure 14b demonstrates a similar tendency of change in N c .A comparison between Figure 14a and Figure 14b shows that a higher C/D makes a higher a coefficient N c .

Parametric Analysis
In this part, a typical situation of Figure 13 shows the variation in N γ depending on ϕ′ .For =1 C D (in Figure 13a), N γ decreases nonlinearly with an increasing ϕ′ .When 30 ϕ′ °＜ , N γ increases with an increasing B L at a constant ϕ′ , and they become close to each other when ϕ′ is between 30° and 35°.N γ decreases with an increasing B L at a constant ϕ′ , which is less than 30°.In Figure 13b, similar variations in N γ at =2 C D can be observed.A comparison between Figure 13a and Figure 13b shows the influence of C D on N γ .A higher C D makes a higher coefficient N γ .Figure 14 shows the variations in c N de- pending on ϕ′ .For =1 C D (in Figure 14a), c N decreases nonlinearly with an increas- ing ϕ′ .As B L becomes smaller, the change in c N becomes more pronounced.When ϕ′ is relatively small, the discrepancy between c N values associated with different B L values becomes more significant.Figure 14b demonstrates a similar tendency of change in c N .A comparison between Figure 14a and Figure 14b shows that a higher C D makes a higher a coefficient c N .3 shows the variation in θ , α , and V depending on ϕ′ .V denotes the total volume of the wedge model, that is, the failure zone.V can be calculated by adding up the volumes of block I, II, and III, that is, It can be seen that with an increasing ϕ′, θ increases, while α and V decrease.In particular ( / 1, 0.4 ), as ϕ′ in- creases from 15° to 45°, θ increases from 55.9° to 71.0°, while α decreases from 48.5° to 29.8°, and V decreases from 117.1 m 3 to 59.4 m 3 .The increase in θ indicates that the wedge slope becomes larger.The decrease in α indicates that the failure zone becomes smaller.Changes in θ and α both cause the volume of the failure zone to decrease.
The decrease in V directly confirms the volume of failure zone decreases with an in- creasing ϕ′ .The reason is that the increase in ϕ′ means that the soil has a greater shear strength, and the greater the self-stabilizing capacity of the soil, the better it is for the tunnel face to sustain stability by itself.At a constant ϕ′ and C D , the increasing B L makes θ decrease, but α and V increase.In particular ( 25 ϕ′ = °), as B L increases from 0.4 to 0.7, θ decreases from 61.7° to 59.4°, while α increases from 41.9° to 63.1°, and V increases from 92.9 m 3 to 114.9 m 3 .In addition, at a constant ϕ′ and B L , θ ,  In this part, typical non-cohesive soil (c /(γD) = 0) is considered for the most adverse scenario.As a result, Equation ( 22) is simplified as s 2 = γDN γ .Table 3 shows the variation in θ, α, and V depending on ϕ .V denotes the total volume of the wedge model, that is, the failure zone.V can be calculated by adding up the volumes of block I, II, and III, that is, It can be seen that with an increasing ϕ , θ increases, while α and V decrease.In particular (C/D = 1, B/L = 0.4), as ϕ increases from 15 • to 45 • , θ increases from 55.9 • to 71.0 • , while α decreases from 48.5 • to 29.8 • , and V decreases from 117.1 m 3 to 59.4 m 3 .The increase in θ indicates that the wedge slope becomes larger.The decrease in α indicates that the failure zone becomes smaller.Changes in θ and α both cause the volume of the failure zone to decrease.The decrease in V directly confirms the volume of failure zone decreases with an increasing ϕ .The reason is that the increase in ϕ means that the soil has a greater shear strength, and the greater the self-stabilizing capacity of the soil, the better it is for the tunnel face to sustain stability by itself.At a constant ϕ and C/D, the increasing B/L makes θ decrease, but α and V increase.In particular (ϕ = 25 • ), as B/L increases from 0.4 to 0.7, θ decreases from 61.7 • to 59.4 • , while α increases from 41.9 • to 63.1 • , and V increases from 92.9 m 3 to 114.9 m 3 .In addition, at a constant ϕ and B/L, θ, α, and V remain virtually constant when C/D increases from 1 to 2. The failure zone is hardly affected by C/D.In this part, typical non-cohesive soil (c /(γD) = 0) is considered.Figure 15 shows the variation in support pressure s 2 depending on the angle of internal friction ϕ .In Figure 15a (C/D = 1), the proposed model predicts that s 2 decreases nonlinearly as ϕ increases, and changes gradually flatten out.The change in s 2 becomes more apparent with a higher B/L.The discrepancy between s 2 values associated with different B/L values becomes more significant at a lower ϕ .In Figure 15b (C/D = 2), a similar variation in s 2 depending on ϕ can be seen.A comparison between Figure 15a and Figure 15b shows that a higher C/D makes a higher s 2 .This is because the higher the C/D, the greater the pressure of the overburden, and therefore, the more likely it is to slide unstably.
In Figure 16a (C/D = 2), the proposed model predicts that s 2 /γD decreases nonlinearly with an increasing ϕ .s 2 /γD associated with a higher B/L has a larger variation.Compared to Figure 15b (C/D = 2), due to s 2 = S 2 /BD (Equation ( 15)), the change in s 2 /γD is related to B/L when γD and L are constant.Thus, the variation in s 2 /γD is different from that of S 2 depending on ϕ .Assuming that the S 2 associated with a lower B/L (assumed to be B 1 /L) is named as S 2 | B 1 /L , the S 2 associated with B 2 /L(B 2 >B 1 ) is named as S 2 | B 2 /L , named in the same way as s 2 /γD| B 1 /L and s 2 /γD| B 2 /L .A simple calcu- lation (s 2 = S 2 /BD) leads to the conclusion that at a constant ϕ , when with a higher B L .The discrepancy between 2 s values associated with different B L values becomes more significant at a lower ϕ′ .In Figure 15b ( 2 C D = ), a similar varia- tion in 2 s depending on ϕ′ can be seen.A comparison between Figure 15a and Figure 15b shows that a higher C D makes a higher 2 s .This is because the higher the C D , the greater the pressure of the overburden, and therefore, the more likely it is to slide unstably.In Figure 16a ( 2 C D = ), the proposed model predicts that 2 s D γ decreases non- linearly with an increasing ϕ′ .) leads to the conclusion that at a constant ϕ′ , when

Conclusions
The large cross-section of the tunnel and the features of soil cutting and discharging caused by multiple cutter heads and screw conveyors tend to cause an uneven distribution of support pressure, which has the risk of instability.This study develops a theoretical model based on the limit equilibrium method for the stability analysis of a box tunnel face in non-cohesive soil.A parametric analysis for the solution is performed.The proposed model is verified by the numerical model.The main conclusions can be summarized as follows:

Conclusions
The large cross-section of the tunnel and the features of soil cutting and discharging caused by multiple cutter heads and screw conveyors tend to cause an uneven distribution of support pressure, which has the risk of instability.This study develops a theoretical model based on the limit equilibrium method for the stability analysis of a box tunnel face in non-cohesive soil.A parametric analysis for the proposed solution is performed.The proposed model is verified by the numerical model.The main conclusions can be summarized as follows: (1) The obtained support pressure distribution on the working face was uneven, with a large support pressure on both sides and small one in the middle.Accordingly, the failure mechanism was composed of a prism (upper part) and a wedge (lower part), which were both divided into three blocks.The relatively small support force acting on the middle block (II) led to its tendency to slide and the risk of destabilization.The interaction forces between the blocks were obtained by force analysis.Assuming that the tunnel face was in the active failure limit state, the support pressure acting on block I was obtained using the active earth pressure coefficient; then, the support pressure acting on block II could be calculated by limit equilibrium analysis; (2) The formula of the limit support pressure was obtained.The influences of the strength parameters and geometric parameters on the tunnel face stability were investigated in the parametric analysis.The limit support pressure rose with the increase in the partially supported ratio and the width-height ratio, and the pressure dropped greatly with the increasing internal friction angle; (3) The applicability of the proposed model was verified by comparing it with the numerical model built in OptumG3.The results of the proposed model have a similar pattern of variation with that of the numerical model, showing that the proposed model agrees well with the numerical model.For projects using similar tunnel machines, the theoretical solution presented in this paper can predict the collapse geometry and provide valuable guidance for risk control.

Figure 1 .
Figure 1.Typical box tunneling machine with multiple cutter heads and double-screw conveyors.

Figure 1 .
Figure 1.Typical box tunneling machine with multiple cutter heads and double-screw conveyors.

Figure 2 .
Figure 2. Distribution characteristics of support pressure.

Figure 2 .
Figure 2. Distribution characteristics of support pressure.

Figure 4 .
Figure 4. Modified model.Figure 4. Modified model.As shown in Figure 5, block I is acted upon by (a) the resultant vertical force of the prism on the top surface oae; (b) the resultant normal forces and shear forces along the failure surfaces ode, abe, and bde; (c) the support force S 1 acting on surface oabd; and (d) its weight.

1 T
(a) front view (b) back view

Figure 5 .
Figure 5.The forces applied on block I.

Figure 5 .
Figure 5.The forces applied on block I.

Figure 6 .
Figure 6.Forces applied on block II.

Figure 6 .
Figure 6.Forces applied on block II.

.Figure 7 .
Figure 7. Geometry of the numerical model built in OptumG3.

Figure 8
Figure 8 indicates the variation in the normalized support pressure with increasing ϕ′ .Both the proposed model and the numerical model predict that 2 s D γ decreases nonlinearly as ϕ′ increases.At constant ϕ′ , 2 s D γ associated with a higher / B L

Figure 7 .
Figure 7. Geometry of the numerical model built in OptumG3.

higher 1 s
. The results of the proposed model are greater than the results of the numerical model.The reliability of the proposed model can be verified.

Figure 12 .
Figure 12.Comparison of failure geometries obtained with the proposed mechanism and with the numerical model for different strength parameters.

6 D
= m and / 2 L D = is considered.The par- tially supported ratio / B L varies from 0.4 to 0.7.For shallowly buried tunnels, the depth ratio / C D varies from 1 to 2, and the total stress concept is accepted for calculating the soil stress.The soil weight was chosen as 3 19.8kN/mγ = .The frictional angle ϕ′ varies from 15° to 45°.

Figure 12 .= 2 Figure 13 . 7 Figure 13 .= 2 Figure 13 .
Figure 12.Comparison of failure geometries obtained with the proposed mechanism and with the numerical model for different strength parameters.Appl.Sci.2022, 12, x FOR PEER REVIEW 16 of 23

= 2 Figure 14 .
Figure 14.Variations in c N depending on ϕ ′ .In this part, typical non-cohesive soil ( ( )0 c D γ ′ =) is considered for the most ad- verse scenario.As a result, Equation (22) is simplified as 2 s DN γ γ = .Table3shows the

α
, and V remain virtually constant when C D increases from 1 to 2. The failure zone is hardly affected by C D .

Figure 14 .
Figure 14.Variations in N c depending on ϕ .

γ
is related to / B L when D γ and L are constant.Thus, the varia- tion in 2 s Dγ is different from that of 2 S depending on ϕ′ .Assuming that the 2

. 2 S shows that when 35
For example, two scenarios are shown in Figure16b.The calculation about

Table 1 .
Parameters for numerical modeling.

Table 2
shows the variation in the ratio of the two support pressures 2 1

Table 1 .
Parameters for numerical modeling.

Table 3 .
Variations in the optimal θ, α, and V depending on ϕ .