Numerical Simulation of Bedrock Sagging Sinkholes in Strain-Softening Rock Induced by Embankment Construction

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. ,e in situ stress of the study area was measured and numerically reproduced. ,e Hoek–Brown strength parameters were obtained by laboratory tests. A strainsoftening constitutive model was established according to the strain-softening behaviour exhibited by the specimens in the triaxial test. ,e stress-strain curves of the specimens were reproduced by numerical methods. ,en, bedrock sagging sinkholes in strainsoftening rock induced by embankment construction were simulated. ,e occurrence of the strain-softening zone and its transition to the residual zone were observed and classified into four stages. ,e 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. ,en, 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.


Introduction
In karst areas, sagging of bedrock is generally induced by dissolution of soluble rock or engineering construction disturbance.
ese 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 Yichang-Wanzhou railway tunnel in China [6].Pile foundation, reinforced steel pipe piling structure, full-section 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 second-time 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.e 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 strain-softening model was established.
e strain-softening behaviour observed by triaxial tests was reproduced numerically.en, the bedrock sagging sinkhole in strain-softening 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.
e occurrence of the strain-softening zone and its transition to the residual zone were analysed.e stress paths of the midspan position and the supports at both ends of the bedrock were analysed.
e similarities and differences of the stress path of the strain-softening model and that of the elastoplastic model were summarized.
e 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.In evaporite karst areas, bedrock sagging sinkholes are generally induced by dissolution of evaporites [10][11][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.

Karst Features.
e carbonate formation in Jiangxi Province, China, is mostly buried by Quaternary sediment formation.
e 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.
ere 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.
e ERT surveys were carried out using the Wenner array.e survey line is arranged along the geological profile (Figure 1).e 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.
e geological profile is shown in Figure 3. ere 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.1) and a maximum of approximately 14 mm of settlement were observed at this site.

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, Q-system, 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.e Hoek-Brown strength criterion can be expressed as where σ ci is the uniaxial compressive strength of the intact rock specimen.e 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). e GSI value is related to weathering conditions and rock structure.e 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 site-specific situation.
where D is the disturbance index.e value of D ranges from 0 (undisturbed rock mass) to 1 (disturbed rock mass).Hoek and Carranza-Torres 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.erefore, the disturbance factor D of laboratory tests was assumed to be 0. In this paper, the embankment is filled above the ground.e embankment and the carbonate formation are separated by a thick Quaternary sediment formation.erefore, it is assumed that the disturbance of embankment construction is very small.e 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.
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.
en, use equation ( 1) to fit the strength curve.e solution that makes R 2 [Θ (GSI, m i )] the largest is marked as Θ * � (GSI * , m * i ), 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.
e 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.erefore, only one unknown parameter, m i , needs to be determined.After fitting the failure curve of the intact rock, the value of m * i can be determined by the equation m * i � arg max {R 2 [Θ (GSI � 100, m i )]}.en, the parameter vector can be obtained by the equation Θ * (intact rock) � (GSI � 100, m * i ). e m i value of rock mass is assumed to be equal with that of intact rock.erefore, 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 � m * i )]}.en, the parameter vector can be obtained by the equation Θ * (rock mass) � (GSI * , m * i ).

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.e permeability of silty clay is very small, and its drainage during construction is neglected.erefore, unconsolidated undrained triaxial tests were carried out.Advances in Civil Engineering To ensure the reliability of material parameters, three samples were tested at a time, and then the test results were averaged to eliminate errors.ree uncon ned 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 uncon ned compressive strength and triaxial tests.e con ning 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 lling and silty clay samples, respectively.e con ning pressures of the triaxial tests were 100, 150, 200, 250, and 300 kPa.
e test results are summarized in Table 1.In the table, E is Young's modulus, v is Poisson's ratio, ρ is density, φ is friction angle, c is cohesion, UCS is uncon ned 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 tting calculation (Figure 6), we obtain 98. e tting curve is the black solid line in Figure 6.According to Hoek's published tables [17,20], the m i value of medium ne limestone is 12 + 3, which is quite close to the tting result in this paper.
According to the assumptions, rock mass and intact rock have the same σ ci value and m i value.erefore, for rock mass, only one unknown parameter, GSI, needs to be determined.After the tting calculation, we obtain GSI * arg max {R 2 e tting curve is the red solid line in Figure 6.
Substitute Θ * (intact rock) (GSI 100, m * i ) and 2), and the Hoek-Brown strength parameters of both intact rock and rock mass can be determined.e determined Hoek-Brown strength parameters are shown in Table 2.

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 gure, 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.
e estimated strength range of the intact rock was basically consistent with the tting result and the test results.However, the estimated strength range of the rock mass was signi cantly lower than the tting result and the test results, showing that the empirical values in Hoek's published tables might Fitting the failure curve of intact rock Fitting the failure curve of rock mass Advances in Civil Engineering not always conform to site-speci c situations.In this case, the GSI range of rock mass determined by Hoek's published tables was 50-60, which was signi cantly smaller than the tting result (GSI 88), as illustrated in Table 3.Similarly, Andriani and Parise suggested that classical rock mass classi cation systems might not be suitable in karst areas, especially when the rocks are a ected by karst processes [16].

Comparison with the Mohr-Coulomb Strength Criterion.
e Mohr-Coulomb strength criterion was used to t the model, and the results were compared with those of the Hoek-Brown strength criterion, as shown in Figure 6. e black dotted line represents the tting result of the Mohr-Coulomb strength criterion to the intact rock test results.e R 2 0.91 was obviously lower than the tting results of the Hoek-Brown strength criterion.e red dotted line represents the tting results of the Mohr-Coulomb strength criterion to the rock mass test results.e R 2 0.86 was also obviously lower than the tting results of the Hoek-Brown strength criterion.
According to the tting curves, the UCS and σ ti values were predicted and are listed in Table 4. e 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.e 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.Figure 7 illustrates a stress-strain 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 strain-softening state, the stress gradually decreases to the residual strength σ 1 res .In the residual state, the stress uctuates around the residual strength σ 1 res .Figure 8 shows the relationship between σ 1 peak , σ 1 res , and con ning pressure σ 3 .e GSI peak value of the elastic state and the GSI res value of the residual state are obtained by tting the strength data with the Hoek-Brown strength criterion.e GSI value and strength curve decreased signi cantly after strain softening.
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 conning pressure, M is noted to have decreased gradually, the strain-softening state lengthened, and the postpeak stress drop velocity slowed down.

Strain-Softening Model.
e rock samples drilled at the site show signi cant strain-softening behaviour in laboratory tests.Obviously, this stress-strain relationship cannot simply be represented by an elastoplastic model.erefore, a strain-softening model was chosen, as shown in Figure 10.
e softening index η can be de ned in di erent ways.A popular method is to de ne it as the plastic shear strain (c p ), which is obtained by subtracting the minor principal plastic strain (ε p 3 ) from the major principal plastic strain (ε p 1 ), i.e., ( e 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 with where K ψ is the dilation coe cient given by where ψ denotes dilatancy.e coe cient ψ could be determined by means of the following expression [17]: e 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. e peak value ω peak is obtained by substituting GSI GSI peak into equation ( 2). e residual value ω res is obtained by substituting GSI GSI res into equation (2).

Numerical Reproduction of the Strain-Softening
Behaviour.It is often helpful to run a simple test of the selected material model before integrating it into a full-scale numerical model.For this aim, the selected strain-softening 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.e material parameters are shown in Tables 1 and 2. e nodes at the bottom boundary of the model were xed vertically and horizontally.To simulate con ning pressure, stress boundary conditions (σ 3 )  Advances in Civil Engineering 7 Elastic state 8 Advances in Civil Engineering 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.e simulated stress-strain curves were compared with the measured stress-strain curves, as shown in Figure 11.In the gure, the measured stress-strain curves are represented by black solid lines, and the simulated curves are represented by red dotted lines.
e results show that the simulated curves under di erent con ning pressures were basically in good agreement with the experimental curves, which indicates that the selected strain-softening model works well in the numerical model.
e shear-strain increment contour map and the displacement vector map are plotted against the sample loaded to damage, as shown in Figure 12. e result shows that shear failure occurred along the inclined plane, and the region with a large shear-strain increment is basically consistent with the sectors where sample damage occurs.e displacement vector map further shows that after the sample yields, two intersecting yield faces are produced.e fragments generated after yielding have both a tendency to slip along the yield surface and a tendency to bulge outward.

Measurement of the In Situ Stress.
e borehole deformation gauge was used to measure the in situ stress in the study area.
e 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.e 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.e in situ stress of the study area consisted mainly of tectonic stress, with the horizontal stress exceeding the vertical stress.

Simulation of the In Situ Stress.
e linear regression equation ( 8) is obtained by tting the survey results.e in situ stress in the study area was approximately linear with depth, as shown in Figure 14.
σ h,max 0.0276H + 2.4594, σ h,min 0.0298H + 0.4475, where H represents the depth.Equation ( 8) is used as the initial stress and integrated into the geological model.e simulation result is shown in Figure 15.In the gure, the black triangle indicates the survey results.Due to the disturbance of the cave, the integrated in situ stress inevitably is redistributed.e 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.

Numerical Simulation of the Bedrock Sagging Sinkhole
6.1.Numerical Model.e embankment construction was simulated using FLAC 3D software, as shown in Figure 16.
e displacement of the bottom and the horizontal displacement of the left and right edges can be neglected.
erefore, the nodes at the bottom boundary of the model were xed vertically and horizontally.e nodes at the left and right edges were xed 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.e constitutive model and yield criteria are shown in Section 2. e material parameters are shown in Tables 1 and 2. e simulation process is as follows: (1) establish a free-eld 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.

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 strain-softening constitutive model; (3) assuming that the surrounding Advances in Civil Engineering   4) assuming that the surrounding medium is composed of rock mass and follows a strain-softening constitutive model.Figure 18 is a comparison between the computation results of 4 schemes and the measurements at site.e computation results of schemes 3 and 4 are completely consistent before EC3, but divergence emerges after EC3. e divergence between them rapidly enlarges as the embankment construction proceeds.
is is because the surrounding medium was in an elastic state before EC3.
erefore, the elastoplastic model and the strain-softening model exhibit the same mechanical behaviour.After EC3, the strength of the strain-softening rock in scheme 4 decreases signi cantly.erefore, the deformation of strainsoftening rock is more severe than that of elastoplastic rock after EC3.
Di erent constitutive models were adopted in schemes 1 and 2, but their computational results are completely consistent from beginning to end. is is because both schemes assume that the surrounding medium is composed of intact rock, which has a high strength.e 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 di erence between the rock mass assumption and the complete rock assumption.
e result shows that deformation will be greatly underestimated if the surrounding medium is assumed to be completely composed of intact rock or the strain-softening behaviour is neglected.e greater the construction load, the more divergence there is.Scheme 4 assumes that the surrounding medium is composed of rock mass and that the strain-softening behaviour is integrated as well.e settlement computation is closest to the site measurements.
Unless otherwise speci ed, the numerical simulation was carried out using scheme 4.

e Occurrence of a Strain-Softening 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 strain-softening zone, 0 < η < η * ; in the residual zone, η ≥ η * .
e occurrence of the strain-softening zone and its transition to the residual zone were observed and classi ed into four stages.During EC1-EC3, the surrounding medium remained elastic, and this period was de ned 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.
is period was de ned 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.
is period was de ned as stage 3.During EC7-EC9, the strain-softening and residual zones expanded in the supports at both ends of the bedrock, and the residual zone emerged in the midspan position.is period was de ned 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.e 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.

Stress Path of the Bedrock.
e strain-softening 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.e strain-softening model and the elastoplastic model have the same yield curve (see the blue dotted line).e red dashed line represents the residual strength of the strainsoftening model.

e Stress Paths of the Strain-Softening Model.
e 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. e two ends yielded due to tension at the end of stage 1.In stage 2, due to the strain-softening behaviour, σ 3 decreased along the unloading path to approximately −0.3 MPa and nally transitioned to the residual zone.In stages 3 and 4, as the embankment construction continued, the stress developed along the residual strength curve.12 Advances in Civil Engineering 3, due to the strain-softening behaviour, σ 3 decreased along the unloading path to approximately 0 MPa and nally transitioned to the residual zone.In stage 4, as the embankment construction continued, the stress could only develop along the residual strength curve.

e Stress Paths of the Elastoplastic Model.
In stage 1, the stress path of the elastoplastic model was exactly the same as that of the strain-softening model.e 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. e stress of the midspan increased along the loading path, but the loading amplitude was signi cantly less than that of the strain-softening model.e midspan still remained in an elastic state until the end of stage 4.
erefore, the stability of caves in elastic-plastic strata is obviously higher than that in strain-softening strata.If strain-softening behaviour is neglected, the deformation of the cave will be signi cantly underestimated, and the stability will be signi cantly overestimated.

e 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 xed supported beam.Calculating the bending moment is helpful to understand the failure mode, especially the reason why the two ends yielded from the top.
e bending moments were calculated with the following three models: (1) the strain-softening numerical model; (2) the simply supported beam model; and (3) the xed supported beam model.e results are shown in Figure 21.e black dotted line, grey solid line, and red solid line represent the bending moments of the strain-softening model, the simply supported beam model, and xed 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.

e Bending Moment of the Strain-Softening 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 strain-softening 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 strain-softening 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.ey yielded earlier than the midspan.e midspan was in tension at its bottom.Due to the rapidly increasing bending moment, the yield zone enlarged quickly at the midspan position.erefore, 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.

e 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    Advances in Civil Engineering 18000 kN•m to 30000 kN•m.erefore, 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.e local regulation, Chinese technical guidance for highway foundation design and construction in a karsti ed area, recommends using the simply supported beam model to calculate the cave stability; however, the simply supported beam model may result in unrealistic deviations.

e Bending Moment of the Fixed Supported Beam
Model.In stage 1, the bending moment of the xed supported beam model was very close to that of the strainsoftening numerical model, with a di erence less than 10%.In stage 2, the bedrock began to yield; thus, the di erence gradually increased to 20%.In stages 3 and 4, the strainsoftening state and residual state emerged in the midspan.
e bending moment calculated by the xed supported beam model was far less than the strain-softening numerical model.erefore, in the elastic state, it would be acceptable to calculate the cave stability using the xed supported beam model.However, after yield, the bending moment would be underestimated, and the cave stability would be overestimated.

e 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 signi cantly, and the vertical shear band C2 propagated through the soil layer.en, the soil mass M2 began to subside along the vertical shear band C2. e 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.
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. e constructed embankment was nally destroyed by several layered horizontal tensile bands C3.
According to the main sinkhole classi cations [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 strain-softening behaviour is one of the main causes of large deformation of Advances in Civil Engineering carbonate rocks.In the settlement analysis, the settlement of the strain-softened 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 ssures, 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.e evolution of a bedrock sagging sinkhole may be triggered by a weathering process, further disturbance, or construction [28][29][30].
e damaged mass collapses progressively  e geometry of the simulated bedrock collapse sinkhole seems to be consistent with the model depicted by other authors [1,23].

Conclusions
In the initial stage of embankment construction, the surrounding medium remains in an elastic state.e elastoplastic rock and the strain-softening rock exhibit the same mechanical behaviour and deformation.However, in later stages after yield, the strength of strain-softening rock decreases signi cantly, and the deformation would be greatly underestimated if strain-softening 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 xed supported beam.
In the elastic state, it would be acceptable to calculate cave stability using the xed 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 classi cations [1,24,25], sinkholes of this type could be classi ed 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.

MFigure 9 :
Figure 9: e trend of the drop modulus M.

Figure 14 :Figure 15 :Figure 16 :
Figure 14: Trend of in situ stress along the depth.

Figure 20 :
Figure 20: Stress path of the bedrock.(a) e support at the left end.(b) e support at the right end.(c) Midspan.
2.1.Lithology.Changli highway is located in the Pingxiang-Leping EW-trending depression zone, where the Yangtze plate meets the South China Folds Belt.
lm ), as shown in Figure1.e 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.e predominant

Table 3 :
Comparison between the tting method and the published estimating method.

Table 4 :
Comparison between Hoek-Brown criterion and Mohr-Coulomb criterion.

Table 5 :
Survey results of in situ stress.