Undrained stability of active and passive trapdoors

The recent growth in the number of sinkhole occurrences due to human activities has highlighted the need for better understanding and prediction of the problem. This paper investigates the use of Broms and Bennermark’s original stability number for trapdoor problems in cohesive soil. The shear-strength-reduction method built in a finite-difference method software program (FLAC) is used to obtain the factor of safety (FOS) under different combinations of pressures for collapse and blowout. Unlike previous research on the use of critical pressure ratios, the FOS results are now functions of the original stability number and depth ratio. The obtained numerical results are compared and validated by using rigorous upperand lower-bound finite-element limit analysis, as well as other existing solutions available in the literature. Surface failure extents are also examined in the paper. The dimensionless ratios employed in this study are useful for preparing design charts with a broad range of trapdoor geometries and soil parameters.

The recent growth in the number of sinkhole occurrences due to human activities has highlighted the need for better understanding and prediction of the problem. This paper investigates the use of Broms and Bennermark's original stability number for trapdoor problems in cohesive soil. The shear-strength-reduction method built in a finite-difference method software program (FLAC) is used to obtain the factor of safety (FOS) under different combinations of pressures for collapse and blowout. Unlike previous research on the use of critical pressure ratios, the FOS results are now functions of the original stability number and depth ratio. The obtained numerical results are compared and validated by using rigorous upper-and lower-bound finite-element limit analysis, as well as other existing solutions available in the literature. Surface failure extents are also examined in the paper. The dimensionless ratios employed in this study are useful for preparing design charts with a broad range of trapdoor geometries and soil parameters.

Introduction
Sinkholes present environmental risk through subsidence or sudden ground collapse, leading to loss of life and infrastructure. The recent growth in the number of sinkhole occurrences due to human activities, such as urbanisation, mining and agricultural development, has highlighted the need for better understanding and prediction of the problem (Drumm et al., 2009). Sowers (1996) outlined the sub-profile of karst soil and described the process of forming a sinkhole. It was suggested that in limestone areas, the gradual erosion of rock at a depth caused by the passing of underground water leads to subsidence of overburden deposited soil, resulting in a saucer-shaped depression. Field investigation studies (Newton, 1976;Sowers, 1996) also suggested that the underground voids, created either naturally or by humans, initiated in cracks between the underground rocks. As indicated by Tharp (2003), the initial size of the cavity does not reflect the actual size of the trapdoor at collapse because the size of the initial cavity will grow further due to internal erosion and will create a reverse-funnel shape.
A considerable number of studies have been published on the stability of trapdoors. A study was initiated by Terzaghi (1936), who experimentally investigated the effect of distributed stress in sand. The study categorised the failure as either an active or a passive mode and described active failure as occuring due to overburden pressure and passive failure occurring as an uplifting force such as an anchor.
Through laboratory experiments and field data collection, Broms and Bennermark (1967) stated that the support pressure required to maintain the stability of an opening on a vertical wall should equate to overburden pressures (surcharge and self-weight) and the undrained shear strength of the soil multiplied by a 'factor'. The stability number (N) was therefore defined in the following equation where s s is the surface surcharge pressure; s t is the support pressure; H is the depth of the opening; g represents the soil unit weight; and S u is the undrained shear strength of the soil. Broms and Bennermark (1967) also concluded that the opening would be unstable when the overburden pressure (s s + gH) is six times greater than the undrained shear strength of the soil.
One of the most influential studies of underground stability comes from Davis et al. (1980), who used an analytical approach to study tunnel heading stability. Davis et al. (1980) used the limit analysis theorem to determine upper-bound and lower-bound solutions to the problem. In their approach, unlike the original Broms and Bennermark (1967) stability number approach, the problem was approached differently by using a critical pressure ratio (s s − s t )/S u , which is a function of the strength ratio (gD/S u ) and depth ratio (H/W), as indicated in Equation 2. Numerous studies have since been performed using this approach (Augarde et al., 2003;Drumm et al., 2009).
Koutsabeloulis and Griffiths (1989) used the displacement finiteelement method to investigate soil displacement in active and passive modes of the trapdoor problem. Martin (2009) introduced a new slipline solution for a shallow trapdoor. Craig (1990) utilised a centrifuge model to investigate the critical stability of a circular cavity. Using the centrifugal approach, many experimental investigations of trapdoor stability have been carried out by researchers, such as Abdulla and Goodings (1996) and Jacobsz (2016). Recently, Keawsawasvong and Ukritchon (2017) studied active trapdoor problems with a linear increase in undrained shear strength with depth using finite-element limit analysis (FELA).
Although extensive research has been carried out on the stability of trapdoors in the past, most of the studies have predominately focused on the use of the critical pressure ratio ((s s − s t )/S u ). Very few studies have used the original stability number approach (Broms and Bennermark, 1967) to study soil stability and further explore the relationship between the stability number and the factor of safety (FOS). In this paper, the shear-strength-reduction method (SSRM) is used with the finite-difference method (FDM) to obtain the FOS for a wide range of stability numbers (N) and depth ratios (H/W). This study also investigates the extent of sinkhole collapse on the ground surface. The results are validated by using FELA and other previous published studies. These FOS results are used to produce comprehensive design charts for the problem of trapdoor stability.

Statement of the problem and modelling technique
The development of cover-collapse sinkholes is a complex procedure due to the continuous expansion of the cavity size over the time. To simplify the problem, it is assumed that the cavity is in the critical stage where the failure is imminent. Figure 1 shows the problem definition of an idealised horizontal trapdoor underlying a homogeneous layer of cohesive soil. The undrained soil is modelled as uniform Mohr-Coulomb material with a zero soil internal friction angle (F = 0). S u is the undrained shear strength, and g is the soil unit weight. The trapdoor opening width is W, and the depth from the surface to the trapdoor opening is notated as H.
Note that the combination of surcharge pressure (s s ), overburden pressure (gH) and support pressure (s t ) can produce failure in either a collapse or a blowout. For undrained clay without volume loss during plastic shearing, stability results are independent of loading directions, and Broms and Bennermark's original undrained stability number is a suitable design parameter (Shiau and Al-Asadi, 2018). For drained soils, the original stability number is not applicable.
A broad range of stability numbers (N = −15-15) and depth ratios (H/W = 1-10) have been chosen to cover all possible investigations of collapse and blowout. Note that the actual values of s s , s t , g, S u , H, and W used in the analyses are insignificant and are not to be reported here due to the nature of the dimensionless definition. The SSRM is adopted to solve the FOS, which is a function of the stability number N and the depth ratio H/W, as shown in the following equation Due to the lack of previous literature on implementing the FOS approach as well as correlating the stability number N to the FOS, it is important to use two techniques for the comparison of results. For this reason, both the FDM and FELA were used to analyse the problem in this paper.
The FDM is one of the oldest techniques used in numerical studies and is a powerful method for analysing complex geotechnical stability problems involving non-linear solutions (Itasca, 2003). A typical grid for simulating the trapdoor stability problem in the FDM is shown in Figure 2. To improve the computational efficiency, a symmetrical condition is considered. This symmetrical condition is particularly important for deep cases that normally require more central processing unit time. An effective domain should be such that it is large enough to present the entire velocity field. Both the left and the right boundaries of the mesh were fixed in the x-direction, allowing the soil to move in a vertical direction. The lower boundary of the mesh was restrained in both the x-and y-directions except where the trapdoor opening is positioned. The trapdoor opening was not restrained so the soil body can freely move downwards into the cavity. A Fish script was developed to assist in auto mesh generation and problem solving of the trapdoor problem using SSRM (Shiau and Sams, 2019;Shiau et al., 2018 particularly important tool in this study, as it allows parametric study to be conducted efficiently. The numerical process of the FELA and SSRM is based on the limit theorems of classical plasticity, which were described by Sloan (2002a, 2002b) and Sloan (2013). The details of the formulation will not be repeated here. It is worth noting that the new technique utilises the finite-element discretisation to deal with complicated geometry and loading conditions and the plastic bounding theorems to bracket the true limit load using upper-and lower-bound solutions. Figure 3 shows the mesh used in the paper. The domain sizes of the models were carefully chosen by observing the non-zero velocity fields, so as to minimise boundary effects. Note that the symmetrical faces are fixed only in the normal direction (i.e. x-direction) to allow vertical movement, so as the outer face boundary. Similar to the FDM mesh, the lower face is restrained in the x-and y-directions and the top face is free to displace in all directions. The solutions are then triggered by involving the shear-strength-reduction technique stated in the paper by Krabbenhoft and Lyamin (2015).
The SSRM was utilised as early as 1975 by Zienkiewicz et al. (1975) and many other researchers to investigate various geotechnical engineering problems (Griffiths and Lane, 1999;Krabbenhoft and Lyamin, 2015;Matsui and San, 1992;Michalowski, 2002;Ugai and Leshchinsky, 1995;Yang and Drumm, 2002). Although the SSRM provides a straightforward solution to many geotechnical problems, such as slopes and retaining walls, this method has seldom been used in the analysis of underground stability problems (Shiau et al., 2017(Shiau et al., , 2018. In this study, the SSRM and FOS approach is adopted to analyse the stability of the trapdoor in collapse and blowout conditions. A total of 690 trapdoor cases are studied using the FDM, upperbound FELA (FELA UB) and lower-bound FELA (FELA LB). Numerical results of the extensive investigation are presented in the form of design charts and equations.

Results and discussions
Comprehensive findings of this study are presented in Tables 1-3 for a broad range of stability numbers (N = −15-15) and depth ratios (H/W = 1-10). Using the data in the tables, Figure 4 plots the FOS results of collapse and blowout from the analyses of lower bound, upper bound and finite differences for the depth ratio of H/W = 3. The results show that the curves are in hyperbolic form where FOS and N are the vertical and horizontal asymptotes, respectively. The general equation of the curve is presented as follows Equation 4 suggests that for a given depth ratio (H/W = 3), any combination of FOS and N on the curve yields a unique value. This unique value is the critical stability number (N c ), which corresponds to an FOS of 1. By drawing an FOS = 1 horizontal line in Figure 4, the two intersection points give an N C value of 4·925 for the collapse and −4·930 for the blowout.
When the supporting pressure ratio (SPR, s t /S u ) is greater than the overburden pressure ratio (OPR ((s s + gH)/S u )), the negative value of N represents a blowout movement. Contrary to this, a positive value of N indicates that the soil moves in the collapse condition. This occurs when the overburden pressure ratio ((s s + gH)/S u )) is greater than the supporting pressure ratio (s t /S u ). As N further increases, an incipient collapse is reached where FOS = 1 and the corresponding N is the critical N c . When the supporting pressure ratio (s t /S u ) is equal to the overburden pressure ratio ((s s + gH)/ S u )), N is equal to zero and FOS is at a maximum (infinite) where a 'stressless' scenario exists on the asymptote line.
The finite-difference results of N c were chosen for the regression analysis. These are presented in Equations 6 and 7 for collapse and blowout, respectively, with a correlation coefficient r 2 = 0·998.
Using Equations 1 and 4, the FOS can be determined with the following equation Substituting Equations 6 and 7 into Equation 8, Equation 9 can be used to determine the FOS for known design parameters (s s , s t , g, H, W and S u ).
Equation 9 is also presented graphically in Figure 6. The design contour map of FOS was constructed based on the FDM numerical solutions.
By rearranging Equation 8, one can determine the required support pressure s t for a given FOS using the following equation Comparison Figure 7 and Table 4 compare the N c values obtained in this paper with those in the published literature. The comparison shows that the FDM results of N c are considerably lower than the analytical upper-bound solutions of Davis (1968) for large H/W. Although Davis's investigation had been improved by Gunn (1980) by using the three rigid block parameters, the analytical upper-bound solutions of Gunn (1980) are still 4·0-12·9% larger than the current FDM.
The comparison in Figure 7 also shows some 0·65-15·00% difference between the FDM results and the results by Sloan et al. (1990) for the range of depth ratios (H/W = 1-10) studied.
Although Sloan et al. (1990) carried out extensive research on trapdoor stability, the use of linear programming was a major drawback on the solution accuracy. In particular, the lower bound shows a large variation in comparison with FDM. The slip-line solutions of Martin (2009)

Failure extent
Results of the failure extent investigation are shown in Figure 8. The distance of failure extent was determined by inspection of the velocity vector plots produced in the program. Figure 8 suggests that the failure extent ratio (E/W) is linearly proportional to the depth ratio (H/W). The linear relationship is presented in the following equation

11.
A practical conclusion can be drawn from ■ Since there is no internal pressure, only the collapse failure should be considered. ■ The stability number is N = gH/S u = 4·21. ■ Using H/W = 6 and N = 4·21, Equation 9 gives an FOS of 1·58 for the collapse. ■ Using H/W = 6 and N = 4·21, Figure 6 gives an approximate FOS of 1·6. ■ An actual computer analysis of this case gives an FOS of 1·63.
Design of a supported cavity (s t ) An FOS of 4 is required for the design of an underground military bunker where the surcharge pressure is given as s s = 50 kPa. The following parameters are known: S u = 25 kPa, g = 18 kN/m 3 , H = 40 m and W = 30 m.
■ Using Equation 6, the critical stability number is N c = 2·89. Note that Figure 5 can also be used to find the N c value. ■ Substitute the N c value into Equation 10, s t = s s + gH − (N c × S u /FOS) = 50 + (18 × 40) − (2·9 × 25/4) = 752 kPa. ■ The required pressure to support the cavity for an FOS of 4 is 752 kPa.

Conclusion
This study successfully investigated the stability of trapdoor problems using Broms and Bennermark's original stability number. Numerical results of this study were obtained by utilising the SSRM. Three numerical techniques were usednamely, the FDM, FELA UB and FELA LB.
The FOS was found to be a function of the depth ratio (H/W) and stability number (N). The numerical results suggest that the FOS increases when the stability number (N) decreases. The FOS becomes very large when the stability N is very small. Further investigation on failure mechanisms indicates a linear relationship between the failure extent ratio (E/W) and the depth ratio (H/W). The failure angle (q) measured from the centre of the opening (W) to the outer boundary of the failure surface is approximately equal to 55°for all depth ratios (H/W).

Geotechnical Research
Undrained stability of active and passive trapdoors Shiau and Hassan