Stability analysis of shield tunnel face under complex ground conditions

Current underground space development features ever-increasing scales. More and more planned tunnel lines have alignments that are characterized by varying passing depths. Both conditions make it more frequent that the shield tunnel face is exposed to heterogeneous ground consisting of distinct soil layers and hydraulic loads, which makes it remarkably challenging to maintain a stable shield face. This work presents a limit equilibrium model to analyze the face stability of shield tunnel under complex ground conditions. The model is consistent with the existing method of slices for uniform soil profile in that the effect of horizontal arching is considered without the need of an a priori assumption on the distribution of horizontal stresses, but extends it to the stratified ground by allowing the variation of failure mechanism across soil layers. Three-dimensional coupled pore fluid flow and stress finite element analyses are performed to study the tunnel face stability in heterogeneous ground consisting of horizontal layers and subjected to hydraulic loads. The results show that the multilayer wedge model can compute limit support pressures that are reasonably close to finite element solutions. Lastly, the finite element method and the multilayer wedge model are employed to analyze the stability of tunnel face in aquifer-aquitard alternate layers with confined water.


Introduction
With the increase of underground space excavation scale, shield tunneling in a coastal area that is largely composed of clay deposits often encounters complex geological conditions exemplified by confined aquifer layers. The encountered aquifers can bear confined water with high pressure (i.e., artesian aquifer) that can trigger seepage towards shield face if a sufficient piezometric head is not maintained in the excavator chamber. In particular, how to determine a proper face support pressure for shield tunneling under such geological conditions is significant, since insufficient and excessive support pressure can result in the tunnel face collapse and induce large ground movements and  Perazzelli et al., 2014).
Theoretical models based on limit equilibrium and limit analysis have been widely employed to assess the stability conditions of shield tunnel (Anagnostou & Kovari, 1996; Anagnostou,2012; Mollon et al., 2010;Zhang et al., 2019). For instance, by assuming a three-dimensional wedge-shaped failure mechanism, limit equilibrium solutions were proposed by Anagnostou and Kovari (1996) and Anagnostou (2012) to compute the limit face support pressure. Wilson et al. (2014) used finite element limit analysis and semi-analytical rigid body method to study the face stability of dual circular tunnel in undrained clays with heterogeneous strength. Based on the upper-bound theorem of limit analysis, Tang et al. (2014) obtained an analytical model for estimating the minimum support pressure for shield tunneling in layered soils. These general analytical techniques, however, are not constructed for complex ground conditions like confined aquifer ground.
Numerical modeling represents a powerful alternative technique for investigating the face stability of shield tunnels due to its advantages in handling complex geological and hydrological conditions (Chen et al., 2011;Vermeer et al., 2002). Perazzelli and Anagnostou (2013) analyzed the stability of shield tunnel face reinforced by bolts using the finite difference software FLAC3D. Perazzelli et al. (2014) incorporated the effects of seepage into analyzing tunnel face stability based on the distribution of hydraulic head around tunnel face in homogenous ground computed from numerical simulations. Lü et al. (2020) considered the effects of seepage on the tunnel face stability through the strength reduction finite element analyses. To the authors' knowledge, there are, however, few numerical studies aimed to analyze the face stability for tunneling in layered ground with confined aquifer. Most of the research is based on homogeneous soil or layered sand ground to study the stability of the tunnel face. But for the layered ground with confined aquifer, there is difference in the hydraulic head distribution and the method that estimate limit support pressure.
This work presents a limit equilibrium model to estimate the limit support pressure of the tunnel face. The model is consistent with the existing method of slices for uniform soil profile in that the effect of horizontal arching is considered without the need of an a priori assumption on the distribution of horizontal stresses, but extends it to the stratified ground by allowing the variation of failure mechanism across soil layers. Three-dimensional coupled pore fluid flow and stress finite element analyses are performed to study the tunnel face stability and validate the proposed model. Lastly, the finite element method and the multilayer wedge model are employed to analyze the stability of tunnel face in layered ground with confined water.

Multilayer wedge model
Anagnostou (2012) based on the equilibrium conditions of an infinitesimally thin slice proposed the wedge model to estimate the limit support pressure of tunnel face in homogenous ground. We extend the method to the analysis in layered ground. Figure 1 sketches the basic concepts related to the multilayer wedge model. The most important feature of the proposed model is that it allows the change of collapse mechanism across soil strata. The wedge body in front of the tunnel face, therefore, is discretized into multiple blocks in accordance with soil stratification, while the individual block can be characterized by different wedge inclinations (i.e., angle ωj). The limit equilibrium of each wedge block is analyzed by using the method of slices such that horizontal arching can be accounted for without the need for an a priori assumption on stress distribution.  Figure 1. The multilayer wedge model: (a) longitudinal and cross section of failure mechanism in layered ground; (b) force equilibrium of infinitesimally thin prism slice and wedge slice.
When the seepage is considered, the stress of an infinitesimal slice in prism is shown in Fig. 1(b), the corresponding equilibrium conditions can be expressed as： where B and t are the width and length, respectively (see Fig. 1(a)). The force components are: (5) where the coefficient of the earth pressure at rest is λi = 1-sinφ'i, c'i, φ'i are the cohesive strength and friction angle, γ'i is the effective gravity. The Eqs.
(2) to (5) are substituted into the Eq. (1) and simplified, and the vertical equilibrium equation for an infinitesimal slice as follows： tan 0 In Eq. (6),  indicates the area-to-circumference ratio of the cross-section of prism, and the seepage force of unit volume fzi is: where hp denotes the approximated distribution of the hydraulic head in prism. It can be expressed as: (1 ) Where h is the difference between the far-field hydraulic head and the piezometric head applied to the tunnel face from inside the tunnel, D is the height of the tunnel face, while x and z are the horizontal and vertical coordinates as depicted in Fig. 1(a). Lü et al. (2017) suggest that the value of a is fixed at 3.5, while the value of b varies in accordance with H0/D. The integration of Eq. (6) leads to the effective vertical pressure in prism ( ) vi z  and the prism load on the wedge silo V . As shown in Fig. 1(b), in order to calculating the vertical force of the wedge and the limit support pressure, the equilibrium of an infinitesimal thin wedge slice in layer soils is considered. The tangential and normal equilibrium equations of the sliding surface are as follows: The expressions for the force components in Eqs. (9) and (10) and more detailed discussion on their meanings can be found in Anagnostou (2012) or Perazzelli et al. (2014). The vertical stress vj  can be expressed as: Where the variables uj and pj indicate the z coordinate and the largest x coordinate at the bottom of the j th wedge block.
The support pressure j s can be obtained by solve the equilibrium equation Eq. (9) and Eq. (10) with considering the boundary conditions of V(0) = 0 (at the bottom of tunnel face) and V(D)= Vsilo (at the top boundary of the wedge). The support pressure s  is:

Finite element analysis
The finite element analyses are conducted using the commercial software ABAQUS/STANDARD (Abaqus, 2014). The finite element model as Fig. 2(a) is adopted in this study. The selected element has 8-node with linear interpolation for displacement and pore pressure (C3D8P). All displacement components at the model bottom and tunnel wall are fixed, whereas the four sides of the model were only laterally constrained. In Fig. 2, the variable C and D denote the cover depth and tunnel diameter, respectively. Accordingly, to reduce the size of the finite element model and also constrain the number of control variables, this work is focused on the scenario of C/D=1 and D=10m. The H is water head at the crown of the tunnel. The pore pressure boundary conditions at the ground surface are specifically setting to modelled the different water table. To consider the influence of the seepage at the tunnel face, a fully coupled soil deformation-water flow analysis is performed according to Biot's consolidation theory (Biot, 1941). Seepage flow analysis in the FEM includes two computation stages. During the first step, initial pore pressure field prior to excavation is computed by allowing equilibrium under the aforementioned pore pressure boundary conditions. The second step involves prescribing zero-valued pore pressure (open drainage boundary) along tunnel face and computes the hydraulic head distribution due to shield tunneling.
The cases that the face in homogenous ground and layered ground as in Figs. 2(c) to (e) are considered in this study. In the simulations presented later, the soil layer is modeled by an elastoperfectly-plastic model with a Mohr-Coulomb yield surface. The non-associated flow rule with zero dilatation angle is adopted for FEM analysis in this study.
The excavation step involves removing displacement constraints along the tunnel face, applying support pressure equal to the earth pressure at rest (calculated according to the K0 value), and subsequently gradually decreasing the support until reaching a limit state. Following previous studies (e.g. Ibrahim et al. 2015; Alagha & Chapman 2019), the minimal face pressure, as marked in Fig. 2(b), is defined as the critical face support pressure.  Figure 2. Finite element modeling: (a) mesh and boundary conditions; (b) example of computed face support pressure-displacement curve and the way to define the limit support pressure; (c) schematics illustrating the shield tunnel face passes the homogenous ground; (d) schematics illustrating the shield tunnel face in layered ground; (e) schematics illustrating the shield tunnel face in layered ground with confined aquifer.

Model validation
To assess the performance of the proposed multilayer wedge model, the comparison between the limit face support pressures computed by the multilayer wedge model and the finite element solutions in homogenous ground and layered ground is shown in this section. The saturated unit weight of γsat=20 kN/m 3 is adopted in this section.

homogeneous ground
We first use the proposed multilayer wedge model to analyze the limit support pressures for shield tunneling in homogeneous ground. This relatively simple case has been investigated by previous studies (e.g., Lü et al. 2017), hence can be used to benchmark the proposed model.
It is seen from Fig. 3(a) that the proposed model consistently gives results closer to FEM solutions than its predecessor. The improvement tends to grow when the effective friction angle of soils is relatively small and/or tunnel excavation at relatively great depths below the water table. We can attribute such enhanced performance to the better consideration of the horizontal arching effects under groundwater flow, as depicted in Fig. 3(b). The classical postulation of linearly varying vertical stresses is independent of the underground water flow conditions and can potentially underestimate the horizontal arching when strong seepage flow is present. The required face support pressures accordingly are overestimated. Moreover, Fig. 3(b) shows that the differences between the stress distributions computed by using the method of slices and the conventional assumption become greater with the increase in the piezometric head at the tunnel face. Table 1 shows the comparison between the results of the proposed method and physical tests. The computed support pressures also match reasonably with the tests data.  Table 1. Comparison between limit support pressures in homogenous ground computed by proposed LEM and those measured in model tests (note that the reported tunnel diameter D in centrifuge tests is in prototype; cover ratio C/D=1. in the cases where dry sand is employed in model tests, the dry unit weight of sand is used to normalize the critical support pressures, while when seepage flow is present, the effective unit weight of soils is included in the normalization).  Fig. 4(a) and (b) compare the limit support pressures computed by using the multilayer wedge model and FEM in layered ground (see Fig. 2(d)). The permeability ratio k1:k2=1 is adopted in this case. It is seen that the results of the proposed model are reasonably close to the results of FEM (again the multilayer wedge model predicts larger support pressures) and noticeably improve the solution of the reference wedge model. As discussed before, the better performance can be attributed to using the method of slices to consider the horizontal arching effects. Fig. 4(a) and (b) respectively indicate that the aforementioned improvement becomes more significant with the decrease in the overall strength of the two-layered soil system (i.e., Δφ' becomes smaller) and with the increase in the piezometric head at the tunnel face. The results indicate that the proposed multilayer wedge model can reasonably predict the critical face support pressures for shield tunneling under complex geological conditions.

The stability analysis in confined aquifer ground
In this section, the hydraulic head distribution in confined aquifer ground is discussed and considered in the computation of the limit support pressure. The soil parameters as Table 2 is adopted in this section.

Hydraulic head distribution
Confined aquifer strata typically contain upper and lower aquitards and interlayer confined aquifer. The confined water within aquifer can bear high pressure and such an artesian head can range from -5m to 6.1m (the datum is set at the ground surface), as reported by case histories worldwide  Fig. 2(e)) is the water head in aquitard and aquifer respectively at the crown of the tunnel. This study mainly discusses the influence of the water head in aquifer H2. The phreatic water table is set at the ground surface (i.e., the H1=10m is adopted for the whole study), while the confined water is modelled by specifying pore pressure boundary at the bottom of the aquifer layer according to given artesian heads, and allows for pore pressure equilibrium before tunnel excavation.
The far-field head distribution is different in complex ground with confined aquifer, and it is significant for approximate hydraulic head adjacent to the tunnel face (Perazzelli et al., 2014). The farfield hydraulic head of the cases as Fig. 2(e) obtain by FEM as shown in Fig. 5(a). According to the results of FEM, the far field hydraulic head in aquitard can be represented by a linear distribution. The far-field hydraulic head in aquifer is uniform as Fig. 5(a). The hydraulic head distribution can be expressed as z is the vertical coordinate of the ground surface. The 2 z is the vertical coordinate of the interface of aqutard and aquifer.
The comparison of the distribution of the hydraulic head along the vertical axis above the tunnel face (x=0) compute by FEM, approximate distribution as Eq. (8) with constant far-field head (Perazzelli et al., 2014) and approximate distribution with variable far-field head is shown in Fig. 5(b). The results by the method of this study is are closer to those of FEM. When the value of the z/D increases, the result of the method of Perazzelli et al. (2014) is larger. As the distance from the tunnel face increases (i.e., the increase of z/D), the approximated hydraulic head is close to the far field hydraulic head according to FEM results. The constant far-field hydraulic head h=H2 is considered in the approximate distribution, and the approximate head is greater than H1 when the H2/H1>1. This can lead to overestimation of seepage force when estimating limit support pressure.  Fig. 6 presents the variation of limit support pressures with the water head in aquifer for the case that the face passes the bottom interface of the aquitard layer (see Fig. 2(e)). In this case, the constant far-field hydraulic head and the linear distribution of far-field hydraulic head is considered in Eq. (8) to approximate hydraulic head adjacent to the tunnel face. For consistency of comparison, the γ'=10 kN/m 3 is used for normalized limit support pressure. As Fig. 6, the limit support pressure compute by FEM is positively correlated with the confined water head. Compared to the multilayer wedge model with the constant far-field head distribution, the method considers the linear distribution of far-field hydraulic head can more reasonably approximate the finite element results.  Figure 6. Comparison between the limit support pressures computed by multilayer wedge model and FEM for the cases that the shield tunnel face in layered ground with confined aquifer (see Fig. 2(e)).

Conclusion
This work presents a limit equilibrium (LE) model for evaluating the critical face support pressures for shield tunneling in layered ground under the presence of groundwater flow. Three-dimensional finite element analyses are conducted to study the tunnel face stability and validate the proposed model. The results show that the multilayer wedge model can compute limit support pressures that are reasonably close to finite element solutions. Lastly, the finite element method and the multilayer wedge model are employed to analyse the stability of tunnel face in layered ground with confined water.
The multilayer wedge model which does not need a priori assumption of stress distribution can obtain more reasonable stress distribution accounts seepage flow, and then can reasonably estimate the limit support pressure of tunnel face in homogeneous ground. Meanwhile, the model can consider the difference of wedge inclination angle caused by the mechanical properties of different soil layers. Through the optimization calculation of the inclination angle, the reasonable limit support pressure of tunnel face in layered ground can be obtained.
Different from in homogeneous ground or in layered ground with the same permeability, the farfield hydraulic head is not constant in layered ground with confined aquifer, and it is linear distribution in aquitard, which leads to the difference of near-field hydraulic head distribution. Considering the characteristics of the far-field hydraulic head distribution in layered ground with confined aquifer, the multilayer wedge model can obtain the limit support pressure which is closer to the FEM results.