Research Article  Open Access
Numerical Simulation of Bedrock Sagging Sinkholes in StrainSoftening Rock Induced by Embankment Construction
Abstract
A bedrock sagging sinkhole occurred in Jiangxi Province of China when constructing the Changli freeway above shallow karst caves. It was chosen as a case to investigate the failure mechanism and potential evolution. The in situ stress of the study area was measured and numerically reproduced. The Hoek–Brown strength parameters were obtained by laboratory tests. A strainsoftening constitutive model was established according to the strainsoftening behaviour exhibited by the specimens in the triaxial test. The stressstrain curves of the specimens were reproduced by numerical methods. Then, bedrock sagging sinkholes in strainsoftening rock induced by embankment construction were simulated. The occurrence of the strainsoftening zone and its transition to the residual zone were observed and classified into four stages. The stress paths of the four stages were analysed. Interestingly, in this case, the supports at both ends of the bedrock began to yield from the top and extended downward, while the midspan position began to yield from the bottom and extended upward, and the reasons for yielding were related to tension. Further analysis found that the failure mode was basically consistent with the size and direction of the bending moment. In fact, this failure mode was quite similar to a fixed supported beam. Then, the feasibility of calculating the stability of karst caves based on beam assumptions was discussed. Finally, potential evolution of the bedrock sagging sinkhole was discussed.
1. Introduction
In karst areas, sagging of bedrock is generally induced by dissolution of soluble rock or engineering construction disturbance. These subsidence mechanisms produce bedrock sagging sinkholes [1]. Carbonate and evaporite rocks are two main types of soluble rocks. Because of the differences in lithology, the sinkholes in carbonate karst areas are different from those in evaporite karst areas. In evaporite karst areas, sinkholes display a significantly higher activity than in carbonate karst areas largely due to the fact that evaporites dissolve much more rapidly than carbonates [2]. In China, soluble rocks are mainly carbonate rocks (e.g., limestone, dolomite) with relatively low solubility and high mechanical strength [3]. However, with the rapid urbanization, sinkholes induced by engineering construction have become an important issue.
To investigate the subsidence mechanism of the Snowy Mountains Highway in Australia, drilling investigations, electromagnetic method, and ground penetrating radar (GPR) were used by Rumbens to determine the extent of the underground karst cavities [4]. During construction of Solan Tunnel in Kangwon Province, South Korea, water inrushes and support collapses accompanying with sinkholes and subsidence were identified by Song et al. in 2012 [5]. Geophysical methods were used to identify the extent of cavity networks. Fan et al. analysed the karst characteristics and treatment methods of YichangWanzhou railway tunnel in China [6]. Pile foundation, reinforced steel pipe piling structure, fullsection pregrouting reinforcement, anchor net spray protection, and other schemes are introduced. In Guangdong Metro Line 9, leakage of the diaphragm wall was observed when excavating in the karst formation. Drainage, grouting, and other schemes are proposed by Cui et al. to mitigate losses [7]. In Jili Village of Guangxi, China, geographical information system (GIS) was used by Zhou et al. to investigate the susceptibility of secondtime karst sinkholes [8]. Ten major sinkholes related factors were evaluated with logistic regression model. By using upper bound theorem, the formula to determine the minimum safe distance between the karst cave and the tunnel was found by Huang et al. in 2017 [9].
To enrich the study of sinkholes induced by engineering construction, the bedrock sagging sinkholes induced by the construction of Changli highway was selected as a case study. Changli highway is a main road connecting Nanchang city and Shangli city. The karst caves beneath the highway were surveyed before embankment construction. To prevent a highway collapse caused by potential sinkholes, a 0.5 m thick continuous reinforced concrete slab was constructed above the embankment. As a passive engineering measure, this procedure can prevent the sudden collapse of highways.
To analyse the bedrock sagging sinkhole occurring beneath the Changli highway, geological conditions and the in situ stresses of the study area were surveyed. Based on unconfined compressive strength tests, Brazilian splitting tests, and triaxial tests, Hoek–Brown strength parameters were obtained, and a strainsoftening model was established. The strainsoftening behaviour observed by triaxial tests was reproduced numerically. Then, the bedrock sagging sinkhole in strainsoftening rock induced by embankment construction was simulated. According to the softening index, the bedrock was divided into an elastic zone, a strainsoftening zone, and a residual zone. The occurrence of the strainsoftening zone and its transition to the residual zone were analysed. The stress paths of the midspan position and the supports at both ends of the bedrock were analysed. The similarities and differences of the stress path of the strainsoftening model and that of the elastoplastic model were summarized. The bending moments of the strainsoftening numerical model, the simply supported beam model, and the fixed supported beam model were calculated. Finally, further evolution of the bedrock sagging sinkhole was discussed.
2. Geological Conditions
2.1. Lithology
Changli highway is located in the PingxiangLeping EWtrending depression zone, where the Yangtze plate meets the South China Folds Belt. The PingxiangLeping depression zone mainly consists of sedimentary rock (limestone, conglomerate, argillaceous siltstone, and sandstone) and metamorphic rock (phyllite, tuffaceous sandstone, and slate) in late Palaeozoic Era and Permian period. The main topographies are low mountains and hills formed by erosion and denudation.
The study area is approximately 150 km west of Nanchang, China. The outcrop of the study area is mainly Quaternary silty clay sediment (Q_{4}^{el}) and Permian Lower Mouth Group limestone (P_{lm}), as shown in Figure 1. The lithology from top to bottom is the following:(i)Quaternary silty clay sediment: it is composed by silty clay with a small amount of gravel. The predominant colour is yellow or brownish red. The sedimentation thickness varies from 0 to 12 m.(ii)Permian Lower Mouth Group limestone: karst caves are formed in medium to finegrained limestone with a uniform composition. The rock mass is characterized by a blocky structure and moderately weathered surface. The predominant colour is grey.
2.2. Hydrology
Surface water is mainly surface runoff, which is derived from meteoric water and bedrock fissure water. The groundwater is mainly bedrock fissure water, which is hosted in weathered bedrock. It is mainly recharged by pore phreatic water. Surface water and groundwater are scarce. The groundwater level detected in the borehole is below the cave.
In evaporite karst areas, bedrock sagging sinkholes are generally induced by dissolution of evaporites [10–12]. While in carbonate karst areas, the dissolution of carbonate rock is much slower than that of evaporites due to the fact that the solubilities of carbonate rock in water are generally lower than 0.1 g/L [1]. Taking into account the above factors, this paper assumes that the dissolution of limestone could be ignored in the construction and service time of the highway.
2.3. Karst Features
The carbonate formation in Jiangxi Province, China, is mostly buried by Quaternary sediment formation. The buried karst had once been exposed in Permian period and then was buried in Quaternary period. Most karst caves remain stable under natural conditions until they are disturbed by surface engineering construction. There have been few reports on sinkholes in Jiangxi Province in the past. However, they have gradually increased in recent years, mainly focusing on sinkholes induced by engineering construction. For example, a bedrock sagging sinkhole was induced by highway construction in Fengcheng city, with a maximum of 12 cm of settlement [13]. In Lianhua County, a building being constructed induced a bedrock sagging sinkhole with a maximum of 4 mm of settlement [14]. Along the Changjin highway, a bedrock sagging sinkhole with a maximum of 6 cm of settlement was observed at site K940 + 200, while a bedrock collapse sinkhole at a depth of 10 m was observed at site K940 + 750 [15].
Borehole drilling and electrical resistivity tomography (ERT) were used to determine the extent of the karst caves. The ERT surveys were carried out using the Wenner array. The survey line is arranged along the geological profile (Figure 1). The in situ measured apparent resistivity values were inverted to obtain the ERT profile with smoothly varying resistivity values using the inversion procedure Res2Dinv. Borehole drilling results were integrated with the ERT profile to determine the geological profile, as shown in Figure 2.
The geological profile is shown in Figure 3. There are no geological problems, such as landslides and debris flows, in the study area. However, embankment construction may destroy the underground karst caves and induce sinkholes. Caves A and C are empty caves with deep burial and thick bedrock. Cave B is a halffilled cave with shallow burial and thin bedrock. The filling of the cave B is mainly silty clay of Quaternary age. After embankment construction, both cave A and C remain stable, while remarkable ground subsidence is observed above the cave B. This paper mainly focuses on cave B because the bedrock above this cave sagged during embankment construction. Several ground fissures (Figure 1) and a maximum of approximately 14 mm of settlement were observed at this site.
3. Discussion of Hoek–Brown Strength Parameters
3.1. Hoek–Brown Strength Criterion
Many classical rock mass classification systems, such as the rock mass rating (RMR) system, Qsystem, and geological strength index (GSI) system, have been developed for years. More recently, a novel rock mass classification system in karst has been put forward by Andriani and Parise in 2017 [16]. In this paper, the classical GSI system and Hoek–Brown strength criterion is used.
The Hoek–Brown strength criterion can be expressed aswhere σ_{ci} is the uniaxial compressive strength of the intact rock specimen. The parameters s and a are determined by the geological strength index GSI, and the parameter m_{b} is determined by GSI and the material constant m_{i}, as shown in equation (2).
The GSI value is related to weathering conditions and rock structure. The m_{i} value reflects the hardness of rocks and is only related to lithology. Hoek has published tables [17] for determining GSI and m_{i} based on weathering, rock structure, and lithology. However, these tables only provide empirical values and might not conform to the sitespecific situation.where D is the disturbance index. The value of D ranges from 0 (undisturbed rock mass) to 1 (disturbed rock mass). Hoek and CarranzaTorres have proposed several empirical guidelines for the selection of D [18]. If excavation by Tunnel Boring Machine or blasting is controlled in excellent quality and the rock disturbance is minimal, the value of D is suggested to be 0. If blasting is in poor quality and the rocks suffer significant disturbance, the value of D is suggested to be 1. Generally, laboratory tests do not involve disturbances such as blasting or excavation. Therefore, the disturbance factor D of laboratory tests was assumed to be 0. In this paper, the embankment is filled above the ground. The embankment and the carbonate formation are separated by a thick Quaternary sediment formation. Therefore, it is assumed that the disturbance of embankment construction is very small. The disturbance factor D of embankment construction was assumed to be 0.1 [19].
Two types of core samples can be drilled from the site: intact rock samples and rock mass samples with randomly distributed joints and fractures, as shown in Figure 4. Generally, the GSI value of intact rock samples can be assumed to be 100, while that of rock mass samples is not easy to determine. Since m_{i} is only related to lithology, the intact rock samples and rock mass samples have the same m_{i} value.
(a)
(b)
To simplify writing, here, we define a parameter vector Θ = (GSI, m_{i}). A simple way to determine Θ is to test the rock strength under different confining pressures. Then, use equation (1) to fit the strength curve. The solution that makes R^{2} [Θ (GSI, m_{i})] the largest is marked as Θ^{∗} = (GSI^{∗}, ), i.e., Θ^{∗} = arg max {R^{2} [Θ (GSI, m_{i})]}. However, R^{2} [Θ (GSI, m_{i})] tends to have multiple extreme values in most conditions. It is quite common to encounter multiple solutions that make R^{2} [Θ (GSI, m_{i})] reach extreme values during fitting, which makes the parameter vector Θ difficult to determine.
To overcome this problem, an indirect method to determine Θ is presented in this paper, as shown in Figure 5. The assumptions include the following: (1) the disturbance factor D of laboratory tests is 0; (2) the GSI value of intact rock is 100; and (3) the intact rock and rock mass have the same σ_{ci} value and m_{i} value [18].
For intact rock, the GSI value is assumed to be 100. According to equation (2), we obtain m_{b} = m_{i}, s = 1 and a = 0.5. Therefore, only one unknown parameter, m_{i}, needs to be determined. After fitting the failure curve of the intact rock, the value of can be determined by the equation = arg max {R^{2} [Θ (GSI = 100, m_{i})]}. Then, the parameter vector can be obtained by the equation Θ^{∗}(intact rock) = (GSI = 100, ).
The m_{i} value of rock mass is assumed to be equal with that of intact rock. Therefore, only one unknown parameter, GSI, needs to be determined. After fitting the failure curve of the rock mass, the value of GSI^{∗} can be determined by the equation GSI^{∗} = arg max {R^{2} [Θ (GSI, m_{i} = )]}. Then, the parameter vector can be obtained by the equation Θ^{∗}(rock mass) = (GSI^{∗}, ).
3.2. Determination of Hoek–Brown Strength Parameters
Unconfined compressive strength test, Brazilian split test, and triaxial test were conducted on rock samples. Unconsolidated undrained triaxial tests were conducted on embankment filling and silty clay samples. Embankment construction was completed within one month. However, it generally takes several years for the embankment filling and silty clay to be fully consolidated. The permeability of silty clay is very small, and its drainage during construction is neglected. Therefore, unconsolidated undrained triaxial tests were carried out.
To ensure the reliability of material parameters, three samples were tested at a time, and then the test results were averaged to eliminate errors. Three unconfined compressive strength tests, three Brazilian splitting tests, and 18 triaxial tests were conducted on intact rock and rock mass samples, respectively. Cylindrical samples with a diameter of 50 mm and a height of 25 mm were used in the Brazilian splitting tests. Cylindrical samples with a diameter of 50 mm and a height of 100 mm were used in the unconfined compressive strength and triaxial tests. The confining pressures of the triaxial tests were 5, 10, 15, 20, 30, and 40 MPa, as shown in Figure 6. A total of 15 triaxial tests were conducted on embankment filling and silty clay samples, respectively. The confining pressures of the triaxial tests were 100, 150, 200, 250, and 300 kPa.
The test results are summarized in Table 1. In the table, E is Young’s modulus, is Poisson’s ratio, ρ is density, φ is friction angle, c is cohesion, UCS is unconfined compressive strength, and σ_{ti} is splitting strength.

Intact rock samples with slight weathered surface were collected from BH2 and BH3, 3 m beneath the caves B and C. Rock mass samples with moderately weathered surface were collected from BH2 and BH3, 3 m above the caves A and C.
According to the uniaxial tests, σ_{ci} = 120 MPa. After the fitting calculation (Figure 6), we obtain = arg max {R^{2} [Θ (GSI = 100, m_{i})]} = 13; R^{2} [Θ (GSI = 100, )] = 0.98. The fitting curve is the black solid line in Figure 6. According to Hoek’s published tables [17, 20], the m_{i} value of medium fine limestone is 12 + 3, which is quite close to the fitting result in this paper.
According to the assumptions, rock mass and intact rock have the same σ_{ci} value and m_{i} value. Therefore, for rock mass, only one unknown parameter, GSI, needs to be determined. After the fitting calculation, we obtain GSI^{∗} = arg max {R^{2} [Θ (GSI, m_{i} = )]} = 88; R^{2} [Θ (GSI^{∗}, )] = 0.95. The fitting curve is the red solid line in Figure 6.
Substitute Θ^{∗}(intact rock) = (GSI = 100, ) and Θ^{∗}(rock mass) = (GSI^{∗}, ) into equation (2), and the Hoek–Brown strength parameters of both intact rock and rock mass can be determined. The determined Hoek–Brown strength parameters are shown in Table 2.

3.3. Comparison with Hoek’s Published Tables
In the absence of test data, the values of GSI and m_{i} are generally determined by Hoek’s published tables [17, 20]. Using the GSI and m_{i} values determined by Hoek’s published tables, the estimated strength ranges are plotted in Figure 6. In the figure, the grey area represents the estimated strength range of the intact rock, and the blue area represents the estimated strength range of the rock mass. The estimated strength range of the intact rock was basically consistent with the fitting result and the test results. However, the estimated strength range of the rock mass was significantly lower than the fitting result and the test results, showing that the empirical values in Hoek’s published tables might not always conform to sitespecific situations. In this case, the GSI range of rock mass determined by Hoek’s published tables was 50–60, which was significantly smaller than the fitting result (GSI = 88), as illustrated in Table 3. Similarly, Andriani and Parise suggested that classical rock mass classification systems might not be suitable in karst areas, especially when the rocks are affected by karst processes [16].

3.4. Comparison with the Mohr–Coulomb Strength Criterion
The Mohr–Coulomb strength criterion was used to fit the model, and the results were compared with those of the Hoek–Brown strength criterion, as shown in Figure 6. The black dotted line represents the fitting result of the Mohr–Coulomb strength criterion to the intact rock test results. The R^{2} = 0.91 was obviously lower than the fitting results of the Hoek–Brown strength criterion. The red dotted line represents the fitting results of the Mohr–Coulomb strength criterion to the rock mass test results. The R^{2} = 0.86 was also obviously lower than the fitting results of the Hoek–Brown strength criterion.
According to the fitting curves, the UCS and σ_{ti} values were predicted and are listed in Table 4. The result shows that the UCS value predicted by the Mohr–Coulomb strength criterion was approximately 10% less than the test value, and σ_{ti} value predicted by the Mohr–Coulomb strength criterion was twice as high as the test value. The predicted values of the Hoek–Brown strength criterion were generally close to the experimental values except that the predicted UCS value of rock mass was slightly smaller than the test value.

4. Testing and Simulation of the StrainSoftening Behaviour
4.1. StrainSoftening Behaviour Observed in Laboratory Tests
Figure 7 illustrates a stressstrain curve logged in a triaxial test. In the test, the rock experienced an elastic state, a strainsoftening state, and a residual state [21] successively. In the elastic state, the stress increases with strain until the yield limit σ_{1}^{peak} is reached. In the strainsoftening state, the stress gradually decreases to the residual strength σ_{1}^{res}. In the residual state, the stress fluctuates around the residual strength σ_{1}^{res}.
Figure 8 shows the relationship between σ_{1}^{peak}, σ_{1}^{res}, and confining pressure σ_{3}. The GSI^{peak} value of the elastic state and the GSI^{res} value of the residual state are obtained by fitting the strength data with the Hoek–Brown strength criterion. The GSI value and strength curve decreased significantly after strain softening.
(a)
(b)
In the test, Young’s modulus was denoted by E, while the drop modulus was denoted by M, as shown in Figure 7. Figure 9 shows the trend of M. With the increase of confining pressure, M is noted to have decreased gradually, the strainsoftening state lengthened, and the postpeak stress drop velocity slowed down.
4.2. StrainSoftening Model
The rock samples drilled at the site show significant strainsoftening behaviour in laboratory tests. Obviously, this stressstrain relationship cannot simply be represented by an elastoplastic model. Therefore, a strainsoftening model was chosen, as shown in Figure 10.
(a)
(b)
The softening index η can be defined in different ways. A popular method is to define it as the plastic shear strain (), which is obtained by subtracting the minor principal plastic strain () from the major principal plastic strain (), i.e.,
The softening index that marks the transition from the softening state to the residual state is denoted as η^{∗}. As illustrated by Figure 10(a), in the elastic state, the value of the softening index would be maintained at zero. In the strainsoftening state, the value ranges between zero and η^{∗}. In the residual state, the value would be greater than η^{∗}. As suggested by Alonso et al. [22], the parameter η^{∗} could be estimated withwhere K_{ψ} is the dilation coefficient given bywhere ψ denotes dilatancy. The coefficient ψ could be determined by means of the following expression [17]:
The strength parameters m_{b}, s, and a are assumed to be piecewise functions of the softening index η:where ω represents a strength parameter that can be replaced by m_{b}, s, or a. The peak value ω^{peak} is obtained by substituting GSI = GSI^{peak} into equation (2). The residual value ω^{res} is obtained by substituting GSI = GSI^{res} into equation (2).
4.3. Numerical Reproduction of the StrainSoftening Behaviour
It is often helpful to run a simple test of the selected material model before integrating it into a fullscale numerical model. For this aim, the selected strainsoftening model was integrated into a numerical model to simulate the triaxial tests using FLAC^{3D} software. Cylindrical models with a diameter of 50 mm and a height of 100 mm were used in the numerical simulation. The material parameters are shown in Tables 1 and 2. The nodes at the bottom boundary of the model were fixed vertically and horizontally. To simulate confining pressure, stress boundary conditions (σ_{3}) were applied to the nodes of the cylindrical boundary. To simulate the compression process, a constant velocity boundary condition (10^{−8} m/s) was applied to the nodes at the top boundary.
The simulated stressstrain curves were compared with the measured stressstrain curves, as shown in Figure 11. In the figure, the measured stressstrain curves are represented by black solid lines, and the simulated curves are represented by red dotted lines. The results show that the simulated curves under different confining pressures were basically in good agreement with the experimental curves, which indicates that the selected strainsoftening model works well in the numerical model.
The shearstrain increment contour map and the displacement vector map are plotted against the sample loaded to damage, as shown in Figure 12. The result shows that shear failure occurred along the inclined plane, and the region with a large shearstrain increment is basically consistent with the sectors where sample damage occurs. The displacement vector map further shows that after the sample yields, two intersecting yield faces are produced. The fragments generated after yielding have both a tendency to slip along the yield surface and a tendency to bulge outward.
(a)
(b)
(c)
5. Measurement and Simulation of the In Situ Stress
5.1. Measurement of the In Situ Stress
The borehole deformation gauge was used to measure the in situ stress in the study area. The layout of the survey points is shown in Figure 13, and the survey results are shown in Table 5. In Figure 13, the three lines intersecting at one point are the projections of three principal stress vectors on the horizontal plane. The length of the line represents the magnitude of the principal stress. In other words, the longest line represents the projection of major principal stress, and the minimum line segment represents the projection of minor principal stress. Table 5 shows that the major principal stress σ_{h,max} and the intermediate principal stress σ_{h,min} are almost horizontal, and the minimum principal stress σ_{v} is almost vertical. The in situ stress of the study area consisted mainly of tectonic stress, with the horizontal stress exceeding the vertical stress.

5.2. Simulation of the In Situ Stress
The linear regression equation (8) is obtained by fitting the survey results. The in situ stress in the study area was approximately linear with depth, as shown in Figure 14.where H represents the depth.
Equation (8) is used as the initial stress and integrated into the geological model. The simulation result is shown in Figure 15. In the figure, the black triangle indicates the survey results. Due to the disturbance of the cave, the integrated in situ stress inevitably is redistributed. The disturbance is most severe above and below the cave. Fortunately, the redistributed results are still generally consistent with the survey results in most places, despite some biased details. From K50 + 620 to K50 + 800, the simulation results are in good agreement with the survey results; from K50 + 500 to K50 + 580, the simulation results are approximately 1.0 MPa greater than the survey results.
6. Numerical Simulation of the Bedrock Sagging Sinkhole
6.1. Numerical Model
The embankment construction was simulated using FLAC^{3D} software, as shown in Figure 16. The displacement of the bottom and the horizontal displacement of the left and right edges can be neglected. Therefore, the nodes at the bottom boundary of the model were fixed vertically and horizontally. The nodes at the left and right edges were fixed horizontally. To simulate the in situ stress in the rock formation, equation (8) was used as the initial stress. According to the site survey, no controlled structural plane is found, so the rock formation was simulated as a continuous medium. The constitutive model and yield criteria are shown in Section 2. The material parameters are shown in Tables 1 and 2. The simulation process is as follows: (1) establish a freefield model containing only rock and soil strata and integrate the in situ stress, as shown in Figure 15; (2) the EC1–EC9 embankment construction term is simulated in succession, as shown in Figure 17. During each construction term, an approximately 1 m thick embankment is constructed.
6.2. Settlement Computation
To compute the settlement of the karst cave C (Figure 3), the following four schemes are adopted: (1) assuming that the surrounding medium is composed of intact rock and follows an elastoplastic constitutive model; (2) assuming that the surrounding medium is composed of intact rock and follows a strainsoftening constitutive model; (3) assuming that the surrounding medium is composed of rock mass and follows an elastoplastic constitutive model; and (4) assuming that the surrounding medium is composed of rock mass and follows a strainsoftening constitutive model. Figure 18 is a comparison between the computation results of 4 schemes and the measurements at site.
The computation results of schemes 3 and 4 are completely consistent before EC3, but divergence emerges after EC3. The divergence between them rapidly enlarges as the embankment construction proceeds. This is because the surrounding medium was in an elastic state before EC3. Therefore, the elastoplastic model and the strainsoftening model exhibit the same mechanical behaviour. After EC3, the strength of the strainsoftening rock in scheme 4 decreases significantly. Therefore, the deformation of strainsoftening rock is more severe than that of elastoplastic rock after EC3.
Different constitutive models were adopted in schemes 1 and 2, but their computational results are completely consistent from beginning to end. This is because both schemes assume that the surrounding medium is composed of intact rock, which has a high strength. The assumed intact rock remains in an elastic state from beginning to end, and therefore, both schemes exhibit the same mechanical behaviour.
In Figure 18, a diagonal parallel pattern is used to represent the role of strain softening on settlement, and a grid pattern is used to represent the difference between the rock mass assumption and the complete rock assumption. The result shows that deformation will be greatly underestimated if the surrounding medium is assumed to be completely composed of intact rock or the strainsoftening behaviour is neglected. The greater the construction load, the more divergence there is. Scheme 4 assumes that the surrounding medium is composed of rock mass and that the strainsoftening behaviour is integrated as well. The settlement computation is closest to the site measurements. Unless otherwise specified, the numerical simulation was carried out using scheme 4.
7. Results and Discussion
7.1. The Occurrence of a StrainSoftening Zone and Its Transition to a Residual Zone
As shown in Figure 19, the surrounding medium was divided into the elastic zone, strainsoftening zone, and residual zone according to the contour map of softening index η. In the elastic zone, η = 0; in the strainsoftening zone, 0 < η < η^{∗}; in the residual zone, η ≥ η^{∗}.
(a)
(b)
(c)
(d)
The occurrence of the strainsoftening zone and its transition to the residual zone were observed and classified into four stages. During EC1–EC3, the surrounding medium remained elastic, and this period was defined as stage 1 for convenience. During EC3–EC5, the supports at both ends of the bedrock began to yield from the top and extended downward. This period was defined as stage 2. During EC5–EC7, the residual zone emerged in the supports at both ends of the bedrock, and the midspan position began to yield from the bottom and extended upward. This period was defined as stage 3. During EC7–EC9, the strainsoftening and residual zones expanded in the supports at both ends of the bedrock, and the residual zone emerged in the midspan position. This period was defined as stage 4.
In summary, in this case, the time when the supports at both ends of the bedrock began to yield was earlier than that at the midspan position, but the later development was slower than that at the midspan position. The time when the midspan position began to yield was later than that in the supports at both ends of the bedrock, but the yield zone enlarged quickly.
7.2. Stress Path of the Bedrock
The strainsoftening model was used to analyse the stress paths of the midspan and the supports at both ends, as shown with the black solid line in Figure 20. For comparison, the stress paths of the elastoplastic model are also plotted in Figure 20, as the solid grey line. The strainsoftening model and the elastoplastic model have the same yield curve (see the blue dotted line). The red dashed line represents the residual strength of the strainsoftening model.
(a)
(b)
(c)
7.2.1. The Stress Paths of the StrainSoftening Model
The stress paths of the supports at both ends were quite similar. As the bedrock tilted to the left, the left end experienced more load, so the stress on the left end was slightly higher than that on the right end. In stage 1, σ_{1} of the two ends was maintained at approximately 0.5 MPa, and σ_{3} increased along the loading path to approximately −2.6 MPa. The two ends yielded due to tension at the end of stage 1. In stage 2, due to the strainsoftening behaviour, σ_{3} decreased along the unloading path to approximately −0.3 MPa and finally transitioned to the residual zone. In stages 3 and 4, as the embankment construction continued, the stress developed along the residual strength curve.
The stress path of the midspan was obviously different from those of the supports at both ends. In phase 1, σ_{1} increased along the loading path to approximately 6 MPa, and σ_{3} was maintained at approximately 0.2 MPa. In stage 2, σ_{3} increased along the loading path to −2.5 MPa. The midspan yielded due to tension at the end of stage 2. In stage 3, due to the strainsoftening behaviour, σ_{3} decreased along the unloading path to approximately 0 MPa and finally transitioned to the residual zone. In stage 4, as the embankment construction continued, the stress could only develop along the residual strength curve.
7.2.2. The Stress Paths of the Elastoplastic Model
In stage 1, the stress path of the elastoplastic model was exactly the same as that of the strainsoftening model. The supports at both ends were gradually loaded to the tensile strength and eventually yielded at the end of stage 1. After yielding, the stress at both ends continued to develop along the yield curve. The stress of the midspan increased along the loading path, but the loading amplitude was significantly less than that of the strainsoftening model. The midspan still remained in an elastic state until the end of stage 4. Therefore, the stability of caves in elasticplastic strata is obviously higher than that in strainsoftening strata. If strainsoftening behaviour is neglected, the deformation of the cave will be significantly underestimated, and the stability will be significantly overestimated.
7.3. The Feasibility of Calculating the Stability of Karst Caves Based on the Beam Assumption
It is generally believed that the underground cavern may lose stability due to tensile yield in the midspan or shear yield in the support at the two ends [23]. However, in this case, the midspan and the two ends both yielded due to tensile load (Figure 20). Another interesting matter observed from Figure 19 is that the supports at both ends of the bedrock begin to yield from the top and extended downward, while the midspan position begins to yield from the bottom and extended upward. In fact, such a failure mode is quite similar to a fixed supported beam. Calculating the bending moment is helpful to understand the failure mode, especially the reason why the two ends yielded from the top.
The bending moments were calculated with the following three models: (1) the strainsoftening numerical model; (2) the simply supported beam model; and (3) the fixed supported beam model. The results are shown in Figure 21. The black dotted line, grey solid line, and red solid line represent the bending moments of the strainsoftening model, the simply supported beam model, and fixed supported beam model, respectively. Table 6 gives the calculation schemes of the bending moment. In the table, x, y, and z represent the coordinate position, and the coordinate axis has been drawn in Figure 21. σ_{x} represents the normal stress. A(x) indicates the cross section of the bedrock, and it varies with its position x. l represents the cave span, and q represents the external load.
(a)
(b)
(c)
(d)

7.3.1. The Bending Moment of the StrainSoftening Numerical Model
In stage 1, the supports at both ends were in tension on the top with a bending moment of approximately −18000 kN·m and eventually yielded due to tension; the midspan was in tension on the bottom with a bending moment of approximately 8000 kN·m. At the end of stage 1, the midspan was still elastic. In stage 2, the supports at both ends were in a strainsoftening state with a bending moment of approximately −23000 kN·m and eventually entered the residual state; the midspan experienced an increased bending moment of up to 8000 kN·m and eventually yielded due to tension. In stages 3 and 4, the bending moment of the supports at both ends increased to approximately −27000 kN·m after they entered the residual state; the midspan entered the strainsoftening state and the residual state successively with a rapidly increased bending moment of up to 24000 kN·m.
In short, the supports at both ends were in tension on the top. They yielded earlier than the midspan. The midspan was in tension at its bottom. Due to the rapidly increasing bending moment, the yield zone enlarged quickly at the midspan position. Therefore, the failure mode of the bedrock shown in Figure 19 is basically consistent with the size and direction of the bending moment shown in Figure 21.
7.3.2. The Bending Moment of the Simply Supported Beam Model
In the four stages, the supports at both ends did not bear any bending moment; the midspan was in tension at its bottom, with the bending moment increasing from 18000 kN·m to 30000 kN·m. Therefore, if the cave stability is calculated using the simply supported beam model, the stability of the midspan position would be underestimated, but the stability of the supports at both ends would be overestimated. The local regulation, Chinese technical guidance for highway foundation design and construction in a karstified area, recommends using the simply supported beam model to calculate the cave stability; however, the simply supported beam model may result in unrealistic deviations.
7.3.3. The Bending Moment of the Fixed Supported Beam Model
In stage 1, the bending moment of the fixed supported beam model was very close to that of the strainsoftening numerical model, with a difference less than 10%. In stage 2, the bedrock began to yield; thus, the difference gradually increased to 20%. In stages 3 and 4, the strainsoftening state and residual state emerged in the midspan. The bending moment calculated by the fixed supported beam model was far less than the strainsoftening numerical model. Therefore, in the elastic state, it would be acceptable to calculate the cave stability using the fixed supported beam model. However, after yield, the bending moment would be underestimated, and the cave stability would be overestimated.
7.3.4. The Occurrence of Bedrock Sagging Sinkhole and Its Potential Evolution
During embankment construction, the studied sinkhole occurred along with a gradual downward movement of the bedrock, leading to the progressive settlement of the Quaternary sediment. As shown in Figure 22, a failure plane began to develop as the inclined yield band C1 propagated through the bedrock. Shortly after, the bedrock M1 started to slide on the inclined yield band C1. Owing to the slide movement of M1, the support of the overlaying soil was reduced significantly, and the vertical shear band C2 propagated through the soil layer. Then, the soil mass M2 began to subside along the vertical shear band C2. The sliding movement of M2 reduced the support on the remaining embankment, which resulted in the formation of the horizontal tensile band C3 and the embankment mass M3.
(a)
(b)
(c)
(d)
An interesting distinction is that shear band C2 propagated vertically along the damaged mass M2, but the tensile band C3 propagated horizontally along the damaged mass M3, although both damaged masses were formed due to the reduction of support. A possible explanation for this is that the damaged mass M2 continued to be compressed by the pressure of the damaged mass M3 during subsidence, in contrast to the tensile stress state of the damaged mass M3. The constructed embankment was finally destroyed by several layered horizontal tensile bands C3.
According to the main sinkhole classifications [1, 24], this type of sinkhole belongs to the bedrock sagging type. Waltham and Fookes suggested that these sinkholes are more likely to occur when constructing above a cave roof with thickness less than 70% of the cave width [25]. Generally, bedrock sagging sinkholes are more common in evaporite areas than in carbonate karst areas [1, 10] because evaporite rocks are able to show a plastic behaviour to the stresses. However, the case in this paper shows that bedrock sagging sinkholes might happen in carbonate areas as well. Similar cases are poorly documented except the subsidence mentioned by Wigham in the Permian limestones of County Durham, England [26], and that mentioned by Sunwoo over an abandoned underground limestone mine [27]. In the case of this paper, the author believes that the strainsoftening behaviour is one of the main causes of large deformation of carbonate rocks. In the settlement analysis, the settlement of the strainsoftened rock can be several times that of the elastic rock (Figure 17).
In fact, the bedrock sagging sinkhole might be a precursor of potential collapse due to the potential further downward movement of bedrock and the overlaid Quaternary sediment. As suggested by Parise and Lollino, failures of underground caves do not occur without warning [23]. In the early stage of the collapsing development, it may show surface fissures, inner failure planes, and gradual subsidence. Furthermore, after the bedrock is penetrated by the residual zone, downward movement of the bedrock may lead to the generation of a bedrock collapse sinkhole, bringing more disturbances to the karst environment and more serious hazards to engineering construction.
A progressive downward movement of the bedrock sagging sinkhole with overlaid Quaternary sediment was simulated to show a potential evolution to bedrock collapse sinkhole, as illustrated in Figure 23. The evolution of a bedrock sagging sinkhole may be triggered by a weathering process, further disturbance, or construction [28–30]. The damaged mass collapses progressively downwards into the cave (Figure 23(a)). Next, we would see the failure along the margin of the sinkhole and the consequent collapse in the same area, followed by the progressive enlargement of the sinkhole (Figure 23(b)). Eventually, the bedrock sagging sinkhole would evolve into a bedrock collapse sinkhole (Figures 23(c) and 23(d)). The geometry of the simulated bedrock collapse sinkhole seems to be consistent with the model depicted by other authors [1, 23].
(a)
(b)
(c)
(d)
8. Conclusions
In the initial stage of embankment construction, the surrounding medium remains in an elastic state. The elastoplastic rock and the strainsoftening rock exhibit the same mechanical behaviour and deformation. However, in later stages after yield, the strength of strainsoftening rock decreases significantly, and the deformation would be greatly underestimated if strainsoftening behaviour were to be neglected.
In this case, the supports at both ends of the bedrock begin to yield from the top and extended downward, while the midspan position begins to yield from the bottom and extended upward, and the reasons for yielding are related to tension. Further analysis found that the failure mode is basically consistent with the size and direction of the bending moment. In fact, this failure mode is quite similar to a fixed supported beam.
In the elastic state, it would be acceptable to calculate cave stability using the fixed supported beam model. However, after yield, the bending moment would be underestimated, and the cave stability would be overestimated. Otherwise, if the cave stability is calculated using the simply supported beam model, the stability of the midspan position would be underestimated, but the stability of the supports at both ends would be overestimated.
According to the main sinkhole classifications [1, 24, 25], sinkholes of this type could be classified as the bedrock sagging type. Furthermore, after the bedrock is penetrated by the residual zone, the bedrock sagging sinkhole may even evolve to a more hazardous bedrock collapse sinkhole.
Data Availability
All data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
Yaoting Zhu, Wenhua Hu, Huquan Wu, and Chaoqun Liu are the funders of the Technology Program of Jiangxi Transportation Hall (Grant no. 2015C0022). The authors would like to acknowledge their supports. This research was funded by the Technology Program of Jiangxi Transportation Hall (Grant no. 2015C0022) and the Fundamental Research Foundation of the Central Universities (Grant no. 310821175026).
References
 F. Gutiérrez, M. Parise, J. De Waele, and H. Jourde, “A review on natural and humaninduced geohazards and impacts in karst,” EarthScience Reviews, vol. 138, pp. 61–88, 2014. View at: Publisher Site  Google Scholar
 J. P. Galve, F. Gutiérrez, P. Lucha et al., “Sinkholes in the saltbearing evaporite karst of the Ebro River valley upstream of Zaragoza city (NE Spain),” Geomorphology, vol. 108, no. 34, pp. 145–158, 2009. View at: Publisher Site  Google Scholar
 Q.L. Cui, H.N. Wu, S.L. Shen, Y.S. Xu, and G.L. Ye, “Chinese karst geology and measures to prevent geohazards during shield tunnelling in karst region with caves,” Natural Hazards, vol. 77, no. 1, pp. 129–152, 2015. View at: Publisher Site  Google Scholar
 A. J. Rumbens, “Detection of cavities in karstic terrain. Road subsidencesnowy mountains highway near yarrangobilly, state of new South walesAustralia,” Exploration Geophysics, vol. 21, no. 12, pp. 121–124, 1990. View at: Publisher Site  Google Scholar
 K.I. Song, G.C. Cho, and S.B. Chang, “Identification, remediation, and analysis of karst sinkholes in the longest railroad tunnel in South Korea,” Engineering Geology, vol. 135136, no. 7, pp. 92–105, 2012. View at: Publisher Site  Google Scholar
 H. Fan, Y. Zhang, S. He, K. Wang, X. Wang, and H. Wang, “Hazards and treatment of karst tunneling in QinlingDaba mountainous area: overview and lessons learnt from Yichang–Wanzhou railway system,” Environmental Earth Sciences, vol. 77, no. 19, p. 679, 2018. View at: Publisher Site  Google Scholar
 Q.L. Cui, S.L. Shen, Y.S. Xu, H.N. Wu, and Z.Y. Yin, “Mitigation of geohazards during deep excavations in karst regions with caverns: a case study,” Engineering Geology, vol. 195, pp. 16–27, 2015. View at: Publisher Site  Google Scholar
 G. Zhou, H. Yan, K. Chen, and R. Zhang, “Spatial analysis for susceptibility of secondtime karst sinkholes: a case study of Jili Village in Guangxi, China,” Computers & Geosciences, vol. 89, pp. 144–160, 2016. View at: Publisher Site  Google Scholar
 F. Huang, L. Zhao, T. Ling, and X. Yang, “Rock mass collapse mechanism of concealed karst cave beneath deep tunnel,” International Journal of Rock Mechanics and Mining Sciences, vol. 91, pp. 133–138, 2017. View at: Publisher Site  Google Scholar
 F. Gutiérrez, J. P. Galve, P. Lucha, J. Bonachea, L. Jordá, and R. Jordá, “Investigation of a large collapse sinkhole affecting a multistorey building by means of geophysics and the trenching technique (Zaragoza City, NE Spain),” Environmental Geology, vol. 58, no. 5, pp. 1107–1122, 2008. View at: Publisher Site  Google Scholar
 M. Parise, L. Pisano, and C. Vennari, Sinkhole Occurrence in Consequence of Heavy Rainstorms, EGU, Munich, Germany, 2016.
 S. Margiotta, S. Negri, M. Parise, and T. A. M. Quarta, “Karst geosites at risk of collapse: the sinkholes at Nociglia (Apulia, SE Italy),” Environmental Earth Sciences, vol. 75, no. 1, pp. 1–10, 2016. View at: Publisher Site  Google Scholar
 Z. X. Han, W. P. Hua, Y. Zhang, and B. Xia, “Formation mechanism of subsidence and cracks on a highway in fengcheng, Jiangxi,” Journal of Geological Hazards & Environment Preservation, vol. 24, no. 4, pp. 31–6, 2013. View at: Google Scholar
 J. Liu, Y. Huang, M. Yang, H. Wei, and X. Chu, “Karst collapse developmental conditions and genetic analysis in the urban area of Lianhua County,Jiangxi,” Chinese Journal of Geological Hazard & Control, vol. 28, no. 3, pp. 58–65, 2017. View at: Google Scholar
 Q. Liu and M. Zheng, “The causes analysis and treatment of karst subgrade diseases—luxi section of ChangboJinyushi expressway,” Journal of East China Jiaotong University, vol. 31, no. 1, pp. 135–40, 2014. View at: Google Scholar
 G. F. Andriani and M. Parise, “Applying rock mass classifications to carbonate rocks for engineering purposes with a new approach using the rock engineering system,” Journal of Rock Mechanics and Geotechnical Engineering, vol. 9, no. 2, pp. 364–369, 2017. View at: Publisher Site  Google Scholar
 E. Hoek and E. Brown, “Practical estimates of rock mass strength,” International Journal of Rock Mechanics and Mining Science & Geomechanics Abstracts, vol. 34, no. 8, pp. 1165–1186, 1997. View at: Publisher Site  Google Scholar
 E. Hoek and C. CarranzaTorres, “HoekBrown failure criterion—2002 edition,” in Proceedings of Fifth North American Rock Mechanics Symposium, pp. 18–22, Toronto, Canada, July 2002. View at: Google Scholar
 F. Wei, Z. D. Chen, Z. F. Chen et al., “Numerical simulation of the mechanical characteristic and failure mode of karst subgrade,” China Journal of Highway & Transport, vol. 31, no. 6, pp. 195–206, 2018. View at: Google Scholar
 P. Marinos and E. Hoek, “Estimating the geotechnical properties of heterogeneous rock masses such as flysch,” Bulletin of Engineering Geology and the Environment, vol. 60, no. 2, pp. 85–92, 2001. View at: Publisher Site  Google Scholar
 X. Wu, Y. Jiang, and Z. Guan, “A modified strainsoftening model with multipostpeak behaviours and its application in circular tunnel,” Engineering Geology, vol. 240, pp. 21–33, 2018. View at: Publisher Site  Google Scholar
 E. Alonso, L. R. Alejano, F. Varas, G. FdezManin, and C. CarranzaTorres, “Ground response curves for rock masses exhibiting strain‐softening behaviour,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 27, no. 13, pp. 1153–1185, 2003. View at: Publisher Site  Google Scholar
 M. Parise and P. Lollino, “A preliminary analysis of failure mechanisms in karst and manmade underground caves in Southern Italy,” Geomorphology, vol. 134, no. 12, pp. 132–143, 2011. View at: Publisher Site  Google Scholar
 F. Gutiérrez, J. Guerrero, and P. Lucha, “A genetic classification of sinkholes illustrated from evaporite paleokarst exposures in Spain,” Environmental Geology, vol. 53, no. 5, pp. 993–1006, 2007. View at: Publisher Site  Google Scholar
 A. C. Waltham and P. G. Fookes, “Engineering classification of karst ground conditions,” Quarterly Journal of Engineering Geology and Hydrogeology, vol. 36, no. 2, pp. 101–118, 2003. View at: Publisher Site  Google Scholar
 D. Wigham, “Occurrence of mininginduced open fissures and shear walls in the Permian limestones of County Durham,” Mining Technology, vol. 109, no. 3, pp. 172–178, 2013. View at: Publisher Site  Google Scholar
 C. Sunwoo, W.K. Song, and D.W. Ryu, “A case study of subsidence over an abandoned underground limestone mine,” Geosystem Engineering, vol. 13, no. 4, pp. 147–152, 2010. View at: Publisher Site  Google Scholar
 P. Van Beynen and K. Townsend, “A disturbance index for karst environments,” Environmental Management, vol. 36, no. 1, pp. 101–116, 2005. View at: Publisher Site  Google Scholar
 M. Parise and J. Gunn, Natural and Anthropogenic Hazards in Karst Areas: Recognition, Analysis and Mitigation, Geological Society, London, UK, 2007.
 L. A. North, P. E. van Beynen, and M. Parise, “Interregional comparison of karst disturbance: westcentral Florida and southeast Italy,” Journal of Environmental Management, vol. 90, no. 5, pp. 1770–1781, 2009. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2019 Zhen Zhang and Zhongda Chen. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.