Three-dimensional sinkhole stability of spherical cavity

Sinkhole occurrences due to underground water mains operations have piqued people’s curiosity. Most research works were in relation to geophysical practices to discover the subsurface cavity. Very few works can be found in relation to the investigation of soil stability due to underground cavity shapes. The actual shape of an underground cavity and its transformation is difficult to predict, though the sinkhole failures are mostly circular in shape on the ground surface. This study explores the three-dimensional collapse stability of three distinct idealized cavity geometries, namely the circular, semi-spherical, and spherical cavities. For an active failure, dimensionless parameters are used to investigate the combined impacts of soil cover, surcharge pressure, soil weight, and internal pressure using advanced finite element limit analysis. Numerical results are compared with the two-dimensional axisymmetric results, and design charts presented to cover a wide range of design parameters for practical applications.


Introduction
Sinkholes are frequent seen in karst terrain. The presence of a subterranean soil void causes a near surface sinkhole collapse. It is a geological condition where bedrock may be dissolved by slightly acidic groundwater [18]. Sinkholes can have various sizes and depths. The formation can be relatively slow and yet progressive. It can also be instantaneous and catastrophic. A comprehensive review on natural and human-induced geohazards and impacts in karst can be found in Gutiérrez et al. [18]. Previous researchers have examined geological case studies of sinkhole stabilities using various numerical techniques [1,4,13,[24][25][26][27][28][29][30]. It is not uncommon to see global news reports on road-related sinkholes. Such sinkhole events have been on the rise in recent decades, and it is proving difficult for concerned authorities to address the issue since the fundamental parameters influencing the magnitude of soil erosion caused by pipe faults are unclear [16].
Subterranean mining, tunnelling, and underground pipes are the primary sources of human-induced sinkhole formation in metropolitan regions [5]. Sinkhole may be caused by several influential factors such as a lack of maintenance, increasing water pressure, differential settling, root damage, and corrosion. Kwak et al. [21] reported that pipeline leaks have caused a shortage of domestic water supply in many developed nations in recent years. The after-effects of these pipeline leaks are soil erosion, subsidence and sinkholes, and they have detrimental impact to the existing infrastructure [15]. Water main damage can have severe humanitarian, social, economic, and environmental implications, and in certain circumstances, they can put people's lives in jeopardy [14]. Figure 1a presents a schematic diagram of a pipeline defectrelated sinkhole problem, while Fig. 1b shows a recent sinkhole event in Brisbane, Australia [49].
There have been numerous studies related to sinkhole failures, even though many of them are confined to very fundamental analyses. For example, flat planar trapdoor stability under active plane strain conditions was investigated by Sloan et al. [46], Martin [22], Wang et al. [47], Keawsawasvong and Ukritchon [20] and Shiau and Hassan [39]. In addition, a circular cavity or tunnel in plain strain condition was studied by Drumm et al. [12], Wilson et al. [48], Yamamoto et al. [50], and Carranza-Torres et al. [8,9] to derive the active stability solutions.
The problems of flat rectangular and circular trapdoors under three-dimensional conditions were considered by Shiau et al. [43]. Shiau and Al-Asadi [31-33, 36, 37] have also studied both the collapse and the blowout problem for two-and three-dimensional tunnels using the latest nonlinear finite element limit analysis technique. Very recently, Shiau et al. [38] studied the passive blowout stability for three different opening shapes of trapdoors under plane strain conditions. Design charts were provided to evaluate soil stability, and it was concluded that the effect of opening shapes on the blowout stability is significant. Other numerical studies using finite element, finite difference, and distinct methods can be found in [2,3,19,34,35,37,[40][41][42]44].
Very few past studies were related to soil stability due to the changes in opening shapes for an active failure condition in three dimensions. Therefore, this paper aims to present advanced three-dimensional upper-and lowerbound active sinkhole failure solutions using finite elements and limit analysis theorems. Soil stability of three cavity shapes, namely the flat circular, semi-spherical, and spherical cavities, with a wide range of dimensionless parameters such as depth ratios and shear strength ratios are investigated in this paper. Two-dimensional axisymmetric models were also established in the study to compare numerical results with the actual 3D analyses. The produced stability design charts and tables can be used in the preliminary design stage of cavity formations by practitioners.

Problem statement and FELA modelling
Three-dimensional stability of three idealized cavity shapes, namely the flat circular, the semi-sphere, and the sphere are studied in the paper using finite element limit analysis. Figures 2, 3, and 4 present the problem definitions, finite element meshes, and failure contour plots for the three idealized shapes representing three stages of cavity formation, respectively. Figure 2 is taken as an example to discuss the detailed 3D modelling. A similar discussion can be applied to the semi-sphere cavity and the full-sphere cavity in Figs the idealized flat circular opening in 3D, representing the initial stage of a pipeline burst. It is to be noted that C is the soil cover from the ground surface to the top of the cavity, and D is the diameter of the opening. The opening face is subjected to a normal internal pressure (r t ), while the ground surface is subjected to an external surcharge pressure (r s ). Tresca soil model with the undrained shear strength (S u ) and the unit weight (c) is considered as the failure criteria in this study. The soil internal friction angle is assumed to be zero throughout the study. Figure 2b shows a typical adaptive mesh and absolute velocity (|u|) contour plot of the idealized flat circular opening. Due to the model's symmetrical state, only a quarter of the 3D mesh is shown along with the bottom view of the circular aperture. The problem's boundary requirements are that the four sides be fixed in the normal direction, whereas the bottom side be fixed in all directions. The top surface is a free surface that is not constrained in any way. On the other hand, the symmetrical plane has all nodes fixed in the normal direction, allowing for tangential movement. The model is established with a large domain to ensure that the provided boundary conditions do not affect the development of total velocity fields [6]. This has been tested by the authors extensively before a final model is determined. In general, the domain extent has been chosen as 5 times the diameter of the trapdoor. The latest adaptive mesh refinement technique was used to determine the 3D limit load [10]. The same boundary conditions are applicable to the semi-sphere cavity and the full-sphere cavity in Figs. 3 and 4.
Broms and Bennermark [7] introduced a stability number (N) in their vertical trapdoor problem in which the effect of surface surcharge, soil self-weight, and the trapdoor pressure are combined into a single dimensionless stability number. A new way to represent the stability number results is to use a pressure ratio {PR = (r s -r t )/ S u } that is a function of soil strength ratio (SR = cD/S u ) and depth-diameter ratio (DR = C/D), as proposed by Davis et al. [11]. The current study adopts the critical pressure ratio to present upper-and lower-bound stability solutions throughout the paper. This is shown in Eq. (1).
In Eq. (1), a standard site investigation such as SPT test will provide values for the soil properties (c and S u ). The values of (C) and (D) can be obtained using geophysics techniques [17]. The surface surcharge (r s ) is a fixed pressure to be specified and the supporting pressure (r t ) is the required parameter to be optimized using upper and lower bound limit analysis techniques. The combined group {(r s -r t )/S u } is characterized as the pressure ratio (PR) that is a function of the depth ratio (DR) and the strength ratio (SR). It delivers a critical information of the pressure difference between the surface surcharge (r s ) and the cavity pressure (r t ), as shown in Eq. (1).
A series of depth ratios (C/D = 0.5-4) and shear strength ratios (cD/S u = 0-2) are used to determine the lower and upper bound limits of the pressure ratio (PR). A total of seventy cases were investigated for the flat circular opening problem where the objective function to be optimized is the critical supporting pressure (r t ). The value of pressure ratio (PR) was computed directly by substituting the obtained critical supporting pressure (r t ) into the (PR) equation, i.e. Eq. (1), which is a function of the other design parameters (C, D, r s , c and S u ).
Limit analysis is most useful when both upper bound (UB) and lower bound (LB) estimations can be used to represent the actual collapse load [45]. The underlying bound theorems assume a rigid-perfectly plastic material and have been effectively utilized to study a variety of drained and undrained stability problems in geotechnical engineering. The latest development of the technique can be found in the OptumG3 FELA software [23]. In the LB method, four-node tetrahedron elements are used in the analysis. Each tetrahedron element has the six nodal stresses that are set to be the basic unknown variables. The statically admissible stress discontinuities are allowed for producing the continuity of normal and shear stresses along with the interfaces of all the elements. The conditions of stress equilibrium, stress boundary condition, and the Tresca failure criterion are all constraints in a typical LB analysis, in which the objective function is to maximize the collapse load of problems.
The upper bound theorem requires a kinematically admissible velocity field where the external work is greater or equal to the plastic shear dissipation. In UB method, tennode tetrahedron elements are used in the formulation. At each node of the element, there are three velocities defined as the basic unknown variables. The setting of kinematically admissible velocity discontinuities is applied at the interfaces of all the elements. The material is set to obey the associated flow rule which is satisfied along any velocity discontinuity. The constraints involved in this procedure are nonlinear and non-smooth but remain convex and amenable to analysis. These two theorems are perfectly fitted to the nonlinear programming optimization problems using the second-order cone programming (SOCP).

Results and discussion
The stability of soils due to water main bursts is studied in this section using the latest three-dimensional FELA technique with upper and lower theorems. Numerical results of the soil stability above the three different cavity shapes are described by the pressure ratio (PR = {(r s -r t )/S u }).

The three idealized shapes
For the phase 1 problem with a flat circular trapdoor (see Fig. 2), the relationships between the pressure ratio (PR) and the depth-diameter ratio (C/D) for the five various shear strength ratios (SR = 0-2) are shown in Fig. 8. Results have shown that the pressure ratio (PR) increases nonlinearly with increasing depth-diameter ratios (C/D). An increase in the critical value of PR means that a larger surcharge pressure r s is required to fail the system, according to the PR = {(r s -r t )/S u } definition. Note that the nonlinear increase occurs for most shear strength ratios (SR = cD/S u ), except for large values of SR = cD/S u (SR [ = 1.5), which may indicate a heavier (c) or a wider (D) system. In this case, the soil stability (PR) eventually decreases as the depth-diameter ratio C/D is greater than 2.5 (see Fig. 5a for cD/S u = 2.0) * a much less geometrical arching support was developed. In general, the upper bounds (r t ) are larger than those of the lower bounds (r t ), and they can always bracket the true solution within a few percentages.
A design contour chart for the flat circular opening problem is shown in Fig. 5b. In this stability design chart, the x-axis represents the depth ratio (C/D), whereas the strength ratio (SR, cD/S u ) is shown in the y-axis. The pressure ratio (PR) is represented by the contour values in the chart. One would enter the known design parameters such as (C/D) and (cD/S u ) to obtain a critical PR value or vice versus. While comparing numerical PR results with the stage one problem (flat circular opening), similar trends are found for stage two problem (see Fig. 6 for the semi-sphere opening) as well as for stage three (see Fig. 7 for the fullsphere opening). It is to be noted that the PR value decreases as the shape transforms from flat circular opening to full-sphere opening. Besides, for larger values of cD/ S u (greater unsupported length D), the reduction in geometrical arching support is significant as the shape transforms from stage 1 to stage 3, resulting in a decrease in PR value as the depth ratio increases.

Comparison of the three idealized shapes
Even though the three idealized shapes are used in the study to represent the internal soil erosion process, it is essential to compare the corresponding effects concerning their stability against collapse failure. Figure 8 presents the pressure ratio (PR) comparison of the three different cavity shapes. Three different strength ratios (SR = 0, 1 and 2) are also shown in the figure.
For shallow depths (C/D \ 1.5, all values of cD/S u ), the both the semi and full-sphere openings give very close PR values. It is to be noted that the PR value of the two shapes (semi-and full-spheres) openings is greater than that of the flat circular opening for the shallow case (C/D \ 1.5). Interestingly, a trend reversal is seen as (C/D) increases, and the flat circular opening has the largest PR among the three opening shapes. As discussed before, owning to the decrease in soil arching for large SR (or D value (e.g. cD/ S u = 2.0), the pressure ratio (PR) decreases as C/D increases (see Fig. 8). It is evident that careful consideration of soil stability for the various opening shapes is essential in the preliminary stage of cavity formations. Figure 9 presents the absolute velocity (|u|) contour plot for a full-sphere cavity. Note that the real values of the velocity fields are not important. Instead, it is more relevant to discuss the overall nonzero velocity fields as they represent the associated possible failure mechanism. The plot highlights the various collapse failure mechanisms of a full-sphere cavity for various depth ratios (C/D = 1-4) with the strength ratio (SR = cD/S u = 1). Interestingly, the results have shown that the mechanism tends to transform from a global failure to a local one as the depth ratio (C/ D) increases. Figure 9a, C/D = 1, shows a global failure that extends from the cavity to the ground surface. While in Fig. 9d, where C/D = 4, most of absolute velocity (|u|) are confined near the cavity indicating a possible local failure mechanism. This is mostly due to the development of geometrical arching in the system. The geometric arching is significant in a 3D analysis, although 2D analysis oversimplifies the solution that is mostly conservative in a geotechnical stability analysis. More detailed study on the associated 3D failure mechanisms as the depth-diameter ratio C/D increases and the 3D soil arching is recommended for a future study.  Table 1. The comparison is made for a flat circular opening. The comparison has shown that the 3D flat circular trapdoor results reported by Shiau et al. [43] agrees well with our current 3D and AX results. Shiau and Al-Asadi [31] studied the tunnel heading blowout and collapse stability using the critical stability number of Broms and Bennermark [7]. Their study is literally a vertical trapdoor problem which is different from the current study. Nevertheless, it is still valuable to present here for the sake of completeness. It was noted that the vertical trapdoor problem (3D tunnel heading) in Shiau and Al-Asadi [31] produces results that are consistently greater than our flat circular trapdoor ones. The larger the C/ D value is, the smaller the difference between the two problems. This is mostly due to the strong geometrical arching in the 3D tunnel heading problem studied in Shiau and Al-Asadi [31].

Comparison with two-dimensional axisymmetric models
Axisymmetric study reduces the computing time by reducing 3D problems into 2D ones. To improve the confidence in the current 3D results, axisymmetric models (AX) for the three stages of cavity development are established.
A typical axisymmetric (AX) adaptive mesh for the flat circular trapdoor is shown in Fig. 11a. The boundary condition for the axial line is such that while allowing vertical movements, the nodes along the line of axial symmetry are not permitted to move in the normal direction. Very accurate results are obtained as ten thousand elements are used for the axisymmetric model. Numerical results of PR in Fig. 11b have shown that very tight bound can be achieved using the two-dimensional AX (UB) and (LB), in comparison with the three-dimensional results. Similar studies have also been implemented for semispherical cavity (Fig. 12) and full-spherical cavity (Fig. 13). The excellent agreement between the AX (2D) and 3D results has greatly enhanced the confidence in producing the design charts in the paper.

Application examples
A simple example is proposed to demonstrate how to use the produced results to evaluate soil stability above various cavity shapes. Figure 14 shows a schematic representation of the problem. The design parameters are given as unit weight c = 20 kPa and undrained shear strength of soil S u = 50 kPa. The recommended tyre pressure of a standard car is 30 psi, i.e. about 207 kPa. With r s = 207 kPa, evaluate the soil stability for the following cavity shapes.  Using Fig. 9, the critical pressure ratio PR = (r s -r t )/ S u = 9.3. 4. Given r s = 207 kPa, S u = 50 kPa, the supporting pressure is calculated as r t = -258 kPa. 5. The negative value of r t = -258 kPa indicates that a downward pressure is needed to cause the system to a failure. 6. In other word, the system is stable even without this ''pulling'' pressure (i.e. when r t = 0).  Using Fig. 11, the critical pressure ratio PR = (r s--r t )/S u = 6.1 4. Given r s = 207 kPa, S u = 50 kPa, the supporting pressure is calculated as r t = -98 kPa 5. The negative sign of r t = -98 kPa indicates that a downward pressure is needed to bring the system to a failure. 6. In other word, the system is stable without this ''pulling'' pressure (i.e. when r t = 0)  Using Fig. 13, the critical pressure ratio PR = (r s--r t )/S u = 3.6 4. Given r s = 207 kPa, S u = 50 kPa, the supporting pressure is calculated as r t = 27 kPa 5. The positive value of r t = 27 kPa indicates that an upward pressure is needed to prevent the system from failure. 6. In other word, the system is unstable without this ''supporting'' pressure. When r t = 0, soil collapse would occur.
As the cavity shape changes, the soil cover C decreases and the cavity diameter or width D increases. Therefore, the depth ratio C/D decreases. Consequently, the soil stability is weakened. The above examples have successfully demonstrated the process of evaluation as well as the usefulness of the design charts.

Conclusion
This paper has successfully used three-dimensional limit analysis with the lower and upper bound theorems to study soil stability above various cavity shapes. The effect of cover depth ratio (C/D) and shear strength ratio (SR) on the pressure ratio (PR) was investigated. The three-dimensional stability solutions were compared with those using axisymmetric models and the comparisons are in excellent agreement. Several design charts were developed for practical uses and an example is used to demonstrate the usefulness of the charts. This research has made significant contributions to our knowledge of soil stability above various cavity shapes and provides a novel sinkhole modelling approach for future research and scientific reference for actual engineering projects. The current work shall be expanded to study layered soil condition over various cavity shapes, in which the finite element limit analysis under axisymmetric conditions can be adopted. In addition, the blowout collapses of cavities are also the immediate future works. Note that this study employs the finite element limit analysis which is commonly based on the continuum mechanics. The possible future works can be the use of the discrete element Funding Open Access funding enabled and organized by CAUL and its Member Institutions.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.