Post-failure behavior of tunnel heading collapse by MPM simulation

Tunnel collapse presents a serious threat to the safety of urban construction. The traditional approach adopted to assess this risk is to evaluate the factor of safety against failure. However, this analysis only determines on whether the tunnel will collapse or not, and does not provide information on the magnitude of the post-failure behavior (for example, catastrophic or progressive) if the tunnel collapse occurs. In this study, a meshless method based on the material point method (MPM) was used to investigate the post-failure behavior of tunnel heading collapse in two-dimensional plane-strain conditions. The capability and accuracy of MPM were verified by comparing the elicited results to centrifuge test data and to analytical solutions obtained from limit state methods. MPM simulations were conducted at different soil conditions (clay or sand) and profiles (homogenous or linear increasing strength) as well as at different tunnel geometries (i.e. tunnel depth and unlined length). The differences in the post-failure behavior and mechanisms are examined and reported.


Introduction
An increasing number of tunnels are being built all over the world. Although substantial technical progress has been made during the past decades to avoid any catastrophic failures during tunnel construction, there are still major incidents being reported. For instance, in Seattle, USA, a tunnel constructed with a tunnel boring machine (TBM) collapsed on the 8 March 2009 due to over-excavation, and created a sinkhole that was 4.5 m in depth and 9 m in width. In Prague, Czech Republic, a tunnel constructed with the new Austria tunneling method (NATM) collapsed on the 6 July 2010. A crater was created that was 20 m wide and 30 m long. In Guangzhou, China, the construction of a tunnel using mine methods caused a very large sinkhole on 28 January 2013. It was approximately 10 m deep and 690 m 2 in area on the surface, and two neighboring buildings were totally damaged.
Some of the past tunnel failures are well documented. Seidenfuß [1] collected data from 110 tunnel collapses and relevant incidents throughout the world during a 70-year period before 2006. The Health and Safety Executive of the UK (HSE) [2] established a list of 39 NATM tunnel incidents from 1973 to 1994, and the Civil Engineering and Development Department of Hong Kong (CEDD) [3] summarized a catalogue with 42 cases of notable tunnel failures. They all provide useful guidance on developing conceptual models for tunnel collapse mechanisms and for estimating the associated potential degrees of damage in the cases at which collapse occurs. For example, HSE [2] identified a number of general collapse mechanisms of NATM tunnels in London clay.
In conventional geotechnical engineering, defining a failure mechanism is the basis for establishing the stability analysis. Past work mainly focused on the determination of various failure mechanisms of tunnels and on the evaluation of the associated minimum support pressures to avoid failure occurrence. Apart from field observations, they are often determined from the results of physical model tests and numerical simulations. Centrifuge models in particular, have led to useful insights into how a tunnel fails. For example, Mair [4] conducted a series of centrifuge tests to study the tunnel heading collapse mechanism in homogenous clay. He then proposed a design chart to demonstrate the effect of heading geometry (tunnel depth and unlined length) on the tunnel stability ratio. Potts [5] and Atkinson and Potts [6] carried out several lined and unlined plain-strain circular tunnel centrifuge tests using sand, and showed that the measured collapse pressures were bracketed closely by the upper and lower bound mechanisms they proposed. Chambon and Corté [7] showed from their centrifuge test results that the collapse mechanism of tunnel heading in sandy soil could be represented by a bulb shape, and that the shape was influenced by the tunnel depth.
The use of computational methods is also popular in evaluating tunnel failures. The finite element method (FEM) has been used to find a critical mechanism with regard to the onset of tunnel collapse and the required support pressure. Many researchers contributed to this since the 1970s [4,5,[8][9][10][11][12][13][14][15]. Sloan et al. [8,12,13,[15][16][17][18][19][20][21] used the upper-and lower bound theorems of classical plasticity in conjunction with FEM to solve the stability problem of the plane-strain tunnel section and the plane-strain tunnel heading. Various geometries (circular, square, dual circular, dual square, rectangular, etc.) and different soil conditions (clay with undrained shear strength uniformly distributed or increasing linearly with depth, and cohesive-frictional soil) were considered by them, and corresponding design charts were proposed for use by practicing engineers for the determination of the tunnel stability. More recently, discrete element methods (DEM), such as the universal distinct element code (UDEC) and the particle flow code (PFC) have been used to simulate tunnel collapse in sandy soils [22][23][24].
The determination of support pressure in the working chamber is a crucial parameter for tunnel safety in shield tunneling, such as the slurry shield and the earth pressure balanced (EPB) shield [25]. Insufficient support pressure can induce large settlements in the ground or even cause collapse of the tunnel face. If a tunnel collapse was to occur, it could be catastrophic or progressive. Although such information could be useful in the design and upon the assessment of the risk of tunnel collapse, most of the past studies give no or limited information on the post-failure mechanism of tunnel collapse. This is partly due to the dif-ferent aims or, in the case of numerical simulations, due to their limited ability to simulate large deformation behavior. In this study, a meshless method, called the material point method (MPM) was used to evaluate post-failure mechanisms of tunnel heading collapse. MPM is a powerful numerical technique for large deformation simulations [26] and can adopt well-developed continuum-based constitutive models.
In this study, the performance of this method was verified first by comparing the results to centrifuge test data and analytical solutions derived from limit state methods. The MPM simulations were then conducted by varying several important factors, such as soil conditions (clay or sand) and profiles (homogenous or linear increasing strength) as well as different tunnel geometries (i.e., tunnel depth and unlined length) in order to examine in detail how these factors would affect the post-failure mechanism and the extent of damage.

Material point method
The MPM is a type of a meshless method combining Lagrangian and Eulerian computational schemes [26,27]. In MPM, objects are discretized into a collection of material points (or particles) using the Lagrangian description, as shown in Figure 1. State variables (e.g., mass, velocity, stress, strain, and history variables) are all traced at these material points throughout the entire simulation process. The procedures of the MPM algorithm at each time step are shown in Figure 1. At the beginning of each step, some variables, such as mass, velocities, and forces, are transferred from the material points to the grid nodes of a predefined Eulerian background mesh, as shown in Figure 1(a). The convection of information between the material points and the background mesh nodes is achieved by various mapping techniques, such as linear interpolation methods used in standard MPM, and higher-order interpolation methods used in the generalized interpolation material-point method (GIMP) [28]. The momentum equation of motion is then solved and the velocities at the mesh nodes are updated as shown in Figure 1(b). From the computed velocity field of the background mesh, strain increments can be computed at material points using the mapping techniques. Subsequently, the stresses of the material points are computed based on the constitutive model. At the same time, the velocities of the material points are computed from the velocities estimated at the mesh nodes. These procedures are shown in Figure  1(c). The material points are finally moved to new positions according to their velocities, as shown in Figure 1(d).
The major advantage of MPM is that the mesh distortion problem in a purely Lagrangian method and the usual convection term in the acceleration associated with the Eulerian formulation are avoided [27,29]. Due to the large deformation calculation capacity, MPM has been utilized to solve a wide range of problems, such as anchor pulling [30], landslide [31,32], retaining wall [26,32], fluid membrane interaction [33,34], silo discharging and filling, the extrusion and cutting problems [32], and others. This study utilizes MPM to examine the post-failure mechanisms of tunnel heading collapses in both clay and sandy soils. The Cambridge MPM code was used for the investigation. Further mathematical details of the method can be found in [35].

Verification of MPM simulation
In order to verify the feasibility and accuracy of MPM in simulating tunnel-heading collapse, the results from MPM were compared to those from centrifuge tests and analytical solutions derived from limit state methods. A centrifuge test for tunnel face collapse was selected from the literature as a reference to check on whether the MPM can obtain a relatively realistic collapse mechanism, especially when simulating large-deformation problems. Furthermore, the MPM was used to calculate the deformation of the tunnel face supported by different internal pressures. The results were compared to those derived from classical analytical methods and limit state-based computational methods.

Simulation of centrifuge test of tunnel collapse
Kamata and Mashimo [22] conducted a series of centrifuge tests in sand to investigate the influence of bolting on face stability. The collapse pattern of the tunnel without any reinforcement (the photo and their corresponding DEM simulation results) is shown in Figure 2. The tunnel section shown in Figure 2(a) was simulated in plane-strain conditions by MPM, whereas the actual tunnel model was a circular opening.
The plane geometry of the MPM model was the same as the cross-section at the window boundary of the centrifuge test, as shown in Figure 3. The diameter of the tunnel was 80 mm and the burial depth was also 80 mm. The background mesh element was 5 mm × 5 mm, and each element consisted of four particles. There were 14592 particles in total in the model. The Mohr-Coulomb model parameters of the sand were determined based on the information given in [22]: c = 0 kPa, φ′ = 34.5°, ψ = 0°, E′= 5 MPa, ν = 0.3, and  = 1510 kg/m 3 . The horizontal velocities of the grid nodes at the two lateral boundaries were set to zero, and the vertical edge surfaces were assigned to a friction coefficient equal to 0.55, which is based on the friction angle of the sand (0.8 × tanφ′). For the bottom boundary, the vertical direction of the grid nodes was fixed, and the horizontal direction was frictional with the same friction coefficient, equal to 0.55. The linings of the tunnel were also modeled by fixing the vertical direction and applying the same horizontal friction as other boundaries on the grid nodes. A vertical boundary was placed at the tunnel face, and the movement in the horizontal direction was controlled to simulate the aluminum plate, as shown in Figure 2(a). When the acceleration due to gravity reached 30 g, the aluminum plate was pulled out to trigger the collapse of the tunnel face. The initial stresses of the MPM model were therefore set to be the geostatic stresses of the soil under 30 g with a lateral pressure coefficient of K 0 = 1-sinφ′ = 0.43.
As the tunnel boundary was moved, the soil began to flow into the tunnel. Figure 4 illustrates the total displacement contour results after the collapse using MPM. The collapse configuration in the centrifuge test is denoted by the red lines. In the MPM simulation, soil near the excavation face flowed into the tunnel and caused a settlement trough at the ground surface. Two slip surfaces of the failure body started from the top and bottom of the excavation face separately. The surface that started from the crown extended almost vertically to the ground surface. The other one, spreading from the invert, extended obliquely upwards to the ground surface. The collapse pattern of the MPM simulation was close to that from the centrifuge test. That is, the two slip surfaces and the failure body of the MPM results were similar to those represented by the red lines. However, the chimney shape of the failure body was wider than that in the centrifuge test. This is considered to be due to the geometry effect (plane-strain tunnel MPM model versus circular opening centrifuge model). An extension of the MPM to 3D is needed to model the centrifuge test more accurately. Nevertheless, the MPM was able to simulate the chimney failure mechanism observed in the centrifuge test.

Comparison with analytical solutions for tunnel stability
In this section, a series of MPM simulations were conducted to examine the post-failure behavior of a tunnel supported by various pressures, and the results were compared to the solutions of the limit state stability methods. The geometry of the model used in the simulations is shown in Figure 5, in which C is the tunnel depth, and D is the diameter of the tunnel (equal to 6 m). The background mesh element was 1 m × 1 m, initially consisting of two particles in each dimension. Overall, there were 22560 particles in total. For the clay cases, an undrained total stress analysis was conducted, whereas for the sand cases, a drained effective stress analysis was conducted. For the left and right boundaries, horizontal velocities of grid nodes were fixed, and vertical velocities were free. For the bottom boundary, both horizontal  and vertical velocities were fixed. A geostatic stress field was applied as the initial condition. As an example, the initial vertical stress contour is shown in Figure 5 for the C/D = 3 sand case. The friction angle of sand used in the following simulations was 25°, and in order to ensure that sand and clay models had relatively similar conditions, a lateral pressure coefficient of K 0 = 1-sinφ′ = 0.58 was adopted for both models. The simulation of the linings was the same as previous models, and the friction coefficient between the soil and tunnel lining was set to be 0.37 (0.8 × tanφ′).

Clay
Davis et al. [36] studied the stability of the plain-strain tunnel heading using the lower and upper bound limit analysis. Sloan and Assadi [8] did similar work based on limit theorems but in conjunction with the FEM. Their upper and lower bound solutions were adopted in this paper to examine the MPM results. The clay was homogenous with c u = 30 kPa, G = 3 MPa, ν = 0.495, and  = 1800 kg/m 3 . The Poisson's ratio ν was set to be 0.495 to achieve the constant volume condition, and to avoid the numerical oscillations at the same time. Four C/D ratios were considered, ranging from 1 to 4. Figure 6 shows the deformation patterns of the tunnel supported by different pressures when C/D = 2. When the support pressure increased, the collapse was mitigated. The maximal horizontal movement of the excavation face and the maximum settlement of the ground surface are selected as measures of the degree of collapse, and the results are shown in Figure 7. It can be seen from Figure 7 that the displacement of the excavation face is always much larger than the settlement of the surface, and the difference between them becomes larger quickly when the support pressure decreases. As the support pressure decreases, the deformation of the tunnel increases, and the rate becomes larger. When the support pressure is larger than a certain trigger value (turning points of the curves were approximately 50 kPa for C/D = 1, 130 kPa for C/D = 2, 220 kPa for C/D = 3, and 310 kPa for C/D = 4), the deformations of the face and ground decrease almost linearly as the support pressure increases. Therefore, these trigger pressures were considered as a critical pressure from the MPM simulations. The turning point of the curve can be determined as the calculation point where the slope of the curve obviously decreases and starts to change slightly.  The trigger pressures, σ T , normalized by the undrained strength, c u , obtained from the MPM simulations, are compared to the limit pressures computed by the two limit analysis methods in Figure 8. The critical pressures derived by the different methods all increase linearly when the C/D ratio increases. The MPM results are located between the upper and lower bound solutions derived from ref. [36]. The limit analysis solutions derived by Sloan and Assadi [8], which were based on relatively denser finite element meshes, are more accurate and have a smaller bracket than the solutions by Davis et al. [36]. The MPM results are close to their lower bound solution. In addition, the collapse patterns of the MPM results are similar to the failure mechanism adopted by Davis et al., depicted in Figure 6 by the red lines.

Figure 8
Critical support pressure obtained from the MPM and two limit analysis methods in clay. Figure 7, for different C/D ratios, the curves have similar shapes and are nearly parallel on both sides of the turning point, i.e., the critical pressure. No matter how deep the tunnel is, when the support pressure approaches the critical pressure, the magnitude of the deformation of the excavation face is at the same level, and so is the settlement of the ground surface. However, the displacement of the excavation face at the critical pressure increases slightly and almost linearly when the C/D ratio increases (as denoted by the dashed line), while the settlement of the surface is maintained nearly constant at different C/D ratio values(as denoted by the dash-dotted line).

Sand
Krause [37] presented some formulas for calculating the minimal support pressure of a 2D tunnel heading in sandy material, based on the limit equilibrium mechanisms proposed by Broms [38], as shown in Figure 9. He pointed out that the quarter-circle mechanism might not be a realistic representation of the failure body, and it could lead to an unreasonably high minimal support pressure.
The parameters for the Mohr-Coulomb model were: c = 0 kPa, φ′ = 25°, ψ = 0°, E′= 20 MPa, ν = 0.3, and  = 1800 kg/m 3 . A small friction angle was selected to observe the collapsed mechanism more clearly. Four C/D ratios between 1 and 4 were considered. The support pressures obtained by the limit equilibrium method for the two different collapse mechanisms were 18.6 and 77.2 kPa for the half-circle and the quarter-circle, respectively. Both solutions showed that the required supporting pressures were independent of the tunnel depth. However, the difference between these two results is significant. According to the conclusion given by Krause [37], 77.2 kPa might be an unreasonably high value and 18.6 kPa is closer to the accurate value. Figure 10 shows the deformation patterns of the tunnel supported by different pressures when C/D = 4, and the two collapse mechanisms, i.e., the half-circle and the quarter-circle mechanisms are indicated using white line and purple dashed line, respectively. However, neither the half-circle nor the quarter-circle mechanisms are similar to the numerical results. The mechanism indicated by the displacement contour has an elliptical form. Figure 11 shows the curves of the maximal horizontal displacement of the tunnel face, and the settlement of the ground surface versus the support pressure calculated by   MPM. The deformation of the excavation face is much larger than the settlement of the surface. From these curves, the turning points can be considered as the critical pressures and the turning points are especially noticeable at small C/D ratios. The curve of the excavation face has a nearly identical turning point with that of the ground surface for each C/D ratio. The critical pressures can be estimated to be respectively 20, 25, 30, and 40 kPa, for C/D ratio values from 1 to 4. Compared with the results derived by the limit equilibrium methods, MPM results are much smaller in value than the corresponding values elicited by the quarter-circle mechanism (77.2 kPa), and are larger in value than the values elicited by the half-circle mechanism (18.6 kPa). Additionally, the critical pressure calculated by MPM increases with the C/D ratio. This is different from the concept of limit equilibrium methods, at least for shallow tunnels, as shown in Figure 12. The centrifuge test results derived by Chambon and Corté also indicated that the critical pressure increases slightly with the C/D ratio value [7,39], which is consistent with the MPM results. Considering that the slip mechanisms used in the limit equilibrium analysis are simplified models and cannot reflect the influence factors, such as the arching effect, thickness of overburden soil, and oth-ers, the MPM results should be more accurate than the limit equilibrium results.

Differences in post-tunnel collapse mechanism between sand and clay
The general post-failure behavior observed in the MPM simulations is shown in Figure 13(a) and (b) for sand and clay, respectively. In sand, soil mass initially moves uniformly perpendicular to the tunnel face, but in time, more mass moves at the bottom half of the tunnel. A sloping surface then gradually forms inside the tunnel under the self-weight of the mass. As the sloping surface is established by reaching its equilibrium state, more soil flows near the crown location and a chimney-shaped failure zone is created.
In clay, the soil mass tends to move uniformly in the horizontal direction from the tunnel face into the tunnel as an integrated block. As more mass flows into the tunnel, the zone of influence extends laterally into the soil formation. In turn, this creates a much wider surface settlement compared to what is observed in sand.
The post-tunnel collapse behavior depends on various factors, such as the tunnel geometry and soil properties. The next section describes the results of the parametric study that examined these effects. Mair [4], Chambon and Corté [7] showed, based on their centrifuge test results that the sizes of the failed chimney  and the settlement trough vary significantly at increasing tunnel depths. The post-failure behavior of the tunnel heading collapse is expected to be influenced by the tunnel depth. To obtain a better understanding of this aspect, a series of collapse simulations were carried out by MPM for tunnels with C/D ranging from 1 to 4.5. The geometries and boundary conditions of the MPM models in this section are the same as those described in Section 3.2. The only difference is that, in order to initiate the tunnel collapse, no support pressure was imposed on the excavation face. Different soil conditions are considered herein.  Figure 14 shows the total displacement contours after collapse in four cases (C/D = 1, 2, 3, and 4). It can be seen that in all the C/D cases, the excavation faces could not remain stable and soil mass flowed into the tunnel. The flowing soil body above the tunnel crown has a chimney-like shape. The displacement contours close to the excavation face are nearly horizontal and with a semi-elliptic shape. Hence, the width of the failed chimney and the settlement trough expands dramatically when the tunnel is deeper, as illustrated in Figures 14 and 15, and the distance between the lowest point of the surface crater and the excavation face becomes  longer.

Clay
The depth and width of the surface crater and flow-in volume are adopted in this paper as indices that reflect the severity of the post-failure behavior. The flow-in volume was defined as the volume of the soil flowing into the tunnel, and it can be obtained approximately by counting the particle numbers inside the tunnel. The width of the surface crater was defined as the length of the ground surface where the settlement was larger than 10 cm. The plots that show the variation of these indices with the C/D ratio values for different soil conditions are shown in Figure 16.
As shown in Figure 16, the flow-in volume, the sizes of the surface crater in width and depth, are very small when the tunnel is shallow, and become larger when the tunnel is deeper. The flow-in volume and the width of the crater both increase almost linearly as the C/D ratio increases, while the depth of the crater increases linearly at small C/D ratio values (i.e., 1 to 2) but gradually becomes constant at large C/D ratio values (i.e., 3 to 4). The results imply that the stability of the tunnel decreases significantly when the burial depth increases.

Clayey ground profile with undrained shear strength c u linearly increasing with depth
In soft clay conditions, the undrained shear strength typically increases with depth. Herein, MPM simulations were conducted in ground conditions, which had an undrained strength profile defined by (5 + 1.6z) kPa, where z is the depth. This ground condition will be called linear clay profile in this paper. Ladd [40] discussed that c u of normally consolidated clay derived from triaxial compression is approximately 0.32 times of the vertical consolidation effective stress, but the undrained shear strength in triaxial extension conditions is 0.5-0.75 times of the compression conditions. Taking this strength anisotropy effect into account, the c u profile adopted here is approximately 0.2 times the vertical consolidation effective stress. The shear modulus was assigned to be 100 times the undrained strength, and hence it increased linearly with depth.
As shown in Figure 17, the total displacement vectors that are very close to the tunnel face are horizontal, which is similar to the homogenous clay case. However, the vectors gradually become inclined as the movement propagates up to the surface. This is in contrast to the homogenous clay cases, in which the movement propagated more laterally. Because of this reason, the widths of failed chimney do not vary significantly at different burial depths.
As shown in Figures 17 and 18, due to the low strength near the ground surface, the collapse is severe for a shallow tunnel in linear clay. It can be seen from Figure 16 that  when the tunnel is relative shallow (C/D < 2.5), the depth and width of the surface crater and flow-in volume of the linear clay case are all much larger than that of the homogenous clay case. Additionally, the relationships of the indices that vary with burial depth in the linear clay case are significantly different to those in the homogenous clay case. The flow-in volume is nearly independent of the tunnel depth. The width of the surface crater increases with the C/D ratio, but the rate of increase is much smaller than that in the homogenous clay case. In contrast to the homogenous clay case, the depth of the crater decreases when the C/D ratio increases.
Based on the differences of the results between the homogenous clay case and the linear clay case, it can be concluded that the soil strength distribution, especially the overburden soil layer, has a significant influence on the tunnel face collapse. When homogenous clay and linear clay have the same strength at the tunnel level, the linear clay exhibits a collapse behavior that is much more severe. In reality, the undrained strength of clay usually increases with depth within the soft soil area, so it might be unsafe to predict the post-failure behavior and stability of the tunnel in such a soil distribution using the conclusions derived from the homogenous clay case. In addition, a large portion of research studies on tunnel stability in clay are currently based on a homogenous distribution, such as the stability ratio charts from the centrifuge tests in [4]. It is advised that particular caution must be exercised when the results derived from the homogenous clay case are applied to the linear clay case.
In a lot of soft soil areas, soil conditions are close to the linear clay profile used in this research. In order to get a better understanding of the potential risk of tunnel collapse in such areas, additional cases were calculated at C/D ratio values ranging from 5 to 9. Additionally, the change rule of the measurement indices was proposed (width δ w , and depth δ d of the surface crater normalized by the tunnel diameter) for different C/D ratio values, as shown in Figure 19. The relationship between the normalized δ w and the C/D ratio is linear, as shown in eq. (1). The decreasing rate of the normalized δ d gradually becomes smaller when the C/D ratio becomes larger, and the relationship can be expressed by the exponential function, as shown in eq. (2). Based on the above analysis, the soil profile influences the degree of the tunnel collapse significantly. It can be argued that, for a given tunnel depth, the homogenous clay case with c u = 30 kPa will be different since the undrained shear strength at the tunnel depth will be different in the linear clay case. Hence, another series of simulations using homogenous clay conditions was conducted by assigning the undrained strength of the homogenous clay model to be equal to the strength at the central line of the tunnel in the linear clay case, as shown in Table 1. This series of simulations is denoted as homogenous clay 2 in this paper. Figures 20 and 21 show the total displacement contours and ground settlements for this series of simulations. Comparing Figure 20 with Figure 17, the difference of the total displacement vectors for the two soil profiles, i.e., the linear clay and the homogenous clay 2, can be observed. For C/D = 1, the contours are similar, but in the linear clay case, the failed chimney is much wider than that of the homogenous case. For larger C/D ratio values, the shapes of the contours in the two cases are quite different. In the linear clay case, the total displacement contours near the tunnel face (such as the area in red color) stretch far towards the ground surface. This indicates that more soil exhibits very large movements Table 1 Undrained strengths (c u ) for homogenous clay 2 models  and flows into the tunnel. On the contrary, in the homogenous clay 2 case, the total displacement contours from the tunnel face extend horizontally for a relatively short distance, and hence less soil flows into the tunnel. This difference is mainly because the homogenous clay cases possess a stronger soil above the tunnel axis. The plots shown in Figure 16 confirm the analysis listed above. The trends of the flow-in volume, and the width and depth of the settlement trough observed in the homogenous clay 2 cases, are similar with those in the linear clay cases. However, all these values are approximately 30% less, indicating that the degrees of damage in the homogenous clay 2 cases are smaller than those elicited in linear clay cases. The results further indicate that the soil strength distribution has a significant influence on the tunnel heading collapse, and the degree of damage of the collapse will be underestimated if the linear clay profile is simplified to a homogenous clay profile.

Sand
For the sand cases, soil parameters are exactly the same as the sand models listed in Section 3.2. In all the C/D cases, the excavation faces would fail. As shown in Figure 22, a sloping surface is present inside the tunnel after the collapse. The contours inside the failed soil body are bulb-shaped. The shape change begins from the excavation face and expands almost vertically towards the ground surface. The width of the failed chimney becomes wider as the tunnel becomes deeper. When the tunnel is extremely shallow, like the C/D = 1 case in Figure 22(a), the failure reaches the ground surface, and creates a deep hollow. As the tunnel depth increases, the surface settlements decrease and the ground surface becomes flat, as shown in Figure 23. Figure 16 also summarizes how the flow-in volume and the width and depth of the surface crater change with C/D ratio values in the sand cases. The flow-in volume is almost independent of the C/D value, exhibiting a trend that is similar to that observed in the linear clay case. The width of the surface crater increases almost linearly proportional to C/D, as shown in Figure 16(b). Simultaneously, the depth of the surface crater decreases rapidly, as shown in Figure 16(c). These opposite trends in the width and the depth are related to the fact that the flow-in volume is independent of the C/D. Figure 24 shows the maximum principal stress contours. It also shows the possible location of the force arching developed for different C/D cases. For C/D=1, no arching can be observed, and the collapse is exposed at the ground surface, as shown in Figure 24(a). When the tunnel depth increases, arching begins to form over the tunnel crown, as shown in Figure 24(b)-(d), and the soil weight of the upper part of the failure body is transferred to two chimney sides, i.e., the tunnel liner on the left hand side, and the stable soil on the right hand side. The arching prevents failure from propagating to the ground, and ensures that the ground   surface remains almost flat.
In an urban environment, the effects of potential failure risk are especially important for buildings, municipal pipes, and citizens. The deflection ratio (defined as the depth over the length of the settlement trough) of the ground surface is directly related to the damage degree of the environment. As shown in Figure 23, when the tunnel is shallow, the deflection ratio is quite large, even though the influenced area is small. Therefore, the damage to the ground, adjacent buildings, and pipelines, could be very severe and irreparable due to the steep and deep sinking. For a relatively deep tunnel, settlements of the surface ground and the deep soil layer (besides the part that is very close to the excavation face) are much smoother, so the damage is limited and may be recovered easily. However, it should be noted that the analysis was based on a sandy soil layer without groundwa-ter, and did not consider the influences of the groundwater flow. For the tunnel located at the soil layer under the groundwater level, the tunnel heading collapse behavior should be analyzed further.
Dilation is an important mechanical feature of sandy soils and the development of arching depends very much on this. In this study, a series of models at different dilation angles was used. Two types of dilatancy models were considered and compared. The first is a constant dilation angle case, in which the dilation angle is kept constant through the entire shearing process. The other is a variable dilation angle case in which the dilation angle varies with the plastic shear strain after yielding, and becomes zero at large strains (i.e., at the critical state). It is assumed that the dilation angle varies linearly with plastic shear strain after yielding, and that when the plastic shear strain is larger than 0.15, sand attains the critical state (see Figure 25). This type of dilation model is called the critical state dilation case in this paper. For both types of dilatancy, dilation angles (or maximum dilation angle for the critical state dilation) of 5°, 10°, 15°, and 20°, were simulated for the cases where C/D = 1 and 3. Figure 26 illustrates some collapse states of different dilation model cases when C/D = 1. Results show that the dilation characteristics have significant influence on the degree of the collapse. In both dilation models, when the dilation angle increases, the surface crater becomes smaller and shallower, and the failed chimney becomes narrower, which means that the impact of the collapse to the ground is mitigated. For the same peak dilation angle, the degree of collapse with the constant dilation case is much smaller than that of the critical state dilation case, as shown in Figure  26(b) and (c).
As shown in Figure 27, the arching mechanism discussed   When the dilation angle is larger, arching develops at a smaller soil displacement and grows slightly in size. Based on the above discussion, it can be concluded that with the aid of the dilation effect, the force arching above the tunnel heading will develop at smaller soil movements. The arching can then effectively prevent the soil flow from moving from the excavation face to the ground, and the collapse is therefore localized near the tunnel face. Figure 28 shows the flow-in volumes and the depths of the surface crater for different dilation cases. Results show that the higher the dilation angle is, the smaller the degree of collapse becomes. The constant dilation angle cases have much greater influence on the degree of tunnel collapse than the critical state dilation case. In addition, when compared to Figures 28(a) and (b), the effect of dilation on the ground settlement (e.g. depth of crater) is larger than that on the tunnel face collapse (e.g., the flow-in volume). This phenomenon confirms the assumption that the dilation effect influences the tunnel collapse mainly through the development of a force arching mechanism. The assumption of a constant dilation angle can lead to an overestimation of the role of the dilation effect on the tunnel collapse. Hence, when the risk of tunnel collapse is assessed, accurate modeling of dilation should be carefully considered.

Effect of tunnel geometry
In conventional tunneling (i.e., the NATM or sprayed concrete tunneling), the excavation face is only a few meters beyond the lined tunnel. This method takes advantage of the standstill ability and arching effect of the soil at the unlined section of the tunnel. In this section, the effect of the unlined length on the post-failure mechanism is investigated The amount of soil moving from the crown then rapidly increases, as shown in Figure 29(b)②. Finally, the soil mass falling downwards from the crown fills the unlined section, and starts to flow horizontally towards the inner parts of the tunnel, as shown in Figure 29(b)③. Figure 30 shows the final post-failure states for four different P/D values in sand, and the unlined section of the tunnel is depicted by white lines. The shape of the soil inside the lined tunnel (i.e., the sloping surface) is very similar. The unlined section is completely filled with the collapsed soil. As shown in Figure 31, the total flow-in volume represents the soil mass volume that flowed into the excavated part, including both the lined and the unlined sections,   and the volume in the lined tunnel only represents the former one. The total flow-in volume increases linearly with the P/D value. However, the volume in the lined tunnel changes slightly. The main differences of the total flow-in volume at different P/D values originate from the length of the unlined section. The width of the surface crater also increases linearly at different P/D values. On the other hand, the depth of the surface crater increases at small P/D values, but the rate of increase then decreases at large P/D values. As shown by the displacement contours, this is because the majority of the collapsed soil moves from the zone within a distance of one tunnel diameter from the crown of the lining, and there is less soil movement beyond that zone. Figure 32 shows the post-failure patterns at different P/D values in homogenous clay. For the cases of small P/D values, the lateral movement of the soil from the face is dominant. For the cases of large P/D values, the vertically downward movements of soil from the crown become dominant. However, it is also important to note that the upheave movement at the invert is much greater than what was observed in the sand cases. The red lines in the figures are the locations of the collapsed zone of the semi-circular tunnels from the centrifuge tests by Mair [4]. In general, the MPM results are similar to the centrifuge test results. However, the widths of the failed chimneys obtained in the MPM simulations are slightly wider because simulations were conducted under 2D plane-strain conditions, which do not allow the development of arching in the circumferential direction of the tunnel. Hence, the difference between the simulation and centrifuge model data becomes greater as P/D increases because the effect of arching in the longitudinal direction diminishes with P/D increases. The linear clay model cases show similar behavior to the homogenous clay cases.
As shown in Figure 31, the total flow-in volume and the width of the surface crater increase with P/D for both homogenous clay and linear clay. For homogenous clay, the flow-in volume in the lined section decreases as the P/D ratio increases. However, for linear clay the flow-in volume in the lined section decreases when the P/D ratio is less than 1, and increases for larger P/D values. The depth of the surface crater was more or less independent of P/D. From the displacement contours in Figure 32, soil flows into the tunnel mainly from the excavation face in the cases where P/D values were small. However, in large P/D value cases, soil flows from the top of the unlined section, and the overburden soil moves nearly downwards. Although the post-failure mechanisms are different, the magnitude of the maximum settlements turns out to be similar.

Conclusions
In this paper, MPM, which is capable of simulating large deformation behavior, was adopted to investigate post-failure mechanisms of tunnel collapse in clay and sand. In particular, the collapse processes and the final state of the collapsed tunnel were examined in detail. The pattern of collapse simulated by MPM was similar to what was observed in a centrifuge test. In addition, the collapse pressures computed by MPM were similar to those computed by various limit state analyses reported in the literature. These results provided confidence in using MPM to pursue the aims of this study. MPM simulations were conducted at different soil conditions (clay or sand) and profiles (homogenous or linear increasing strength) as well as different tunnel geometries (i.e., at different C/D and P/D values). The main findings from the simulation results are summarized below.
(1) In homogenous sand, the width of the surface crater becomes larger and the depth becomes shallower when the tunnel is deeper. The collapse is always localized near the tunnel face, and flow-in volume does not vary with the tunnel depth. When the tunnel depth is shallow, the collapse can outcrop on the ground surface. However, when the tunnel depth increases, a smooth settlement trough is observed due to the arching that develops above the collapsed zone. The development of arching is governed by the dilatancy characteristics of the sand, and the stability of a tunnel may be overestimated if a constant dilation angle is used in the analysis. Therefore, it is recommended to model the dilation characteristics of the sand carefully when the post-failure behavior of tunnel heading collapse needs to be assessed.
(2) When a tunnel is excavated in a homogenous clay profile, the magnitude of collapse increases dramatically when the tunnel becomes deeper. The size of the surface crater and the flow-in volume both increase at increasing tunnel depths. Soil mass tends to move uniformly in the horizontal direction from the tunnel face into the tunnel. As more soil mass flows into the tunnel, the zone of influence extends laterally into the soil formation. This in turn creates a much wider surface settlement compared to what is observed in sand. Therefore, tunneling-induced collapse in clayey soils can lead to a more catastrophic failure compared to that in sandy soil.
(3) When a tunnel is excavated in a linear clay profile, the tunnel collapse is relatively severe for the cases of small C/D ratio values. The relationship between the magnitude of the collapse and the tunnel depth becomes similar to what was observed in sand because both soil profiles yield an increasing strength with depth. The flow-in volume was approximately independent of C/D. As C/D increased, the width of the surface crater increased but the depth decreased because the flow-in volumes are similar in all C/D cases. Because the linear clay profile was close to the soil condition in soft soil areas in practical engineering terms, expressions of the relationship between the depth and width of the surface crater and the C/D ratio are proposed for reference for similar engineering applications. Additionally, results suggest that the soil strength distribution has a significant influence on the post-failure behavior of a collapsed tunnel, and it might be unsafe to evaluate the stability of a tunnel in a linear clay profile based on the results using the equivalent homogenous clay profile.
(4) Unlined length in conventional tunneling can dramatically change the post-failure mechanism. For a short unlined length, soil flows from the tunnel face, but when the unlined length becomes larger, soil mainly flows from the tunnel crown. However, the extent of soil mass moving into the lined section of the tunnel is similar for all P/D cases. In sand, the majority of soil moves in the zone within a distance of one tunnel diameter from the crown of the lining, and there is less soil movement beyond that zone. In clay, the upheave movement at the invert is much greater than what is observed in the sand cases. In turn, this makes the extent of the soil moving into the tunnel larger and leads to wider surface settlements. The flow-in volume and the width of the surface crater increase at increasing unlined lengths. In clay, the depth of the surface crater is approximately independent of P/D. In sand, the depth of the surface crater increases with P/D at small P/D values. However, the rate of increase then decreases at larger P/D values. This is because the main soil movements occur in the zone within a distance of one tunnel diameter from the crown of the lining.
It should be noted that the MPM simulations were conducted in 2D plane-strain geometry, which is obviously different from the actual 3D tunnel geometry. Although the findings reported in this paper may provide insights into the post-failure behavior of tunnel heading collapse, further investigation is needed to quantify more precisely the magnitude of damages. Three-dimensional MPM simulations are currently underway and will be reported in future publications.