Advances in Civil Engineering
Volume 2019 (2019), Article ID 9426029, 19 pages
https://doi.org/10.1155/2019/9426029

Research Article

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

Zhen Zhang and Zhongda Chen

Key Laboratory for Special Area Highway Engineering of Ministry of Education, Chang’an University, Xi’an 710064, China

Correspondence should be addressed to Zhen Zhang; 470604245@qq.com

Received 1 October 2018; Accepted 28 January 2019; Published 14 February 2019

Academic Editor: Giovanni Biondi

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.

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 strain-softening constitutive model was established according to the strain-softening behaviour exhibited by the specimens in the triaxial test. The stress-strain curves of the specimens were reproduced by numerical methods. Then, bedrock sagging sinkholes in strain-softening rock induced by embankment construction were simulated. The occurrence of the strain-softening 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 Sol-an 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. 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 strain-softening model was established. The strain-softening behaviour observed by triaxial tests was reproduced numerically. Then, 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 strain-softening zone, and a residual zone. The occurrence of the strain-softening 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 strain-softening model and that of the elastoplastic model were summarized. The bending moments of the strain-softening 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 Pingxiang-Leping EW-trending depression zone, where the Yangtze plate meets the South China Folds Belt. The Pingxiang-Leping 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₄el) and Permian Lower Mouth Group limestone (Plm), 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 fine-grained limestone with a uniform composition. The rock mass is characterized by a blocky structure and moderately weathered surface. The predominant colour is grey.

Figure 1: Geological map.

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.

Figure 2: ERT profile.

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 half-filled 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.

Figure 3: Geological profile.

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, 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.

The Hoek–Brown strength criterion can be expressed as

where σ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 mb is determined by GSI and the material constant mi, as shown in equation (2).

The GSI value is related to weathering conditions and rock structure. The mi value reflects the hardness of rocks and is only related to lithology. Hoek has published tables [17] for determining GSI and mi 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. The 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. 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 mi is only related to lithology, the intact rock samples and rock mass samples have the same mi value.

(a)

(b)

Figure 4: Intact rock sample (a) and rock mass sample (b).

To simplify writing, here, we define a parameter vector Θ = (GSI, mi). 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² [Θ (GSI, mi)] the largest is marked as Θ∗ = (GSI∗, ), i.e., Θ∗ = arg max {R² [Θ (GSI, mi)]}. However, R² [Θ (GSI, mi)] tends to have multiple extreme values in most conditions. It is quite common to encounter multiple solutions that make R² [Θ (GSI, mi)] 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 mi value [18].

Figure 5: Determination procedure of Hoek–Brown strength parameters.

For intact rock, the GSI value is assumed to be 100. According to equation (2), we obtain mb = mi, s = 1 and a = 0.5. Therefore, only one unknown parameter, mi, 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² [Θ (GSI = 100, mi)]}. Then, the parameter vector can be obtained by the equation Θ∗(intact rock) = (GSI = 100, ).

The mi 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² [Θ (GSI, mi = )]}. 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.

Figure 6: Fitting curves of triaxial tests.

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.

Table 1: Basic material properties.

Materials
E (MPa)

ρ (kg·m⁻³)
φ (°)
c (kPa)
UCS (MPa)
σti (MPa)

Embankment filling
50
0.35
2100
24
10
–
–

Silty clay
30
0.35
1900
20
50
–
–

Intact rock
60000
0.25
2700
44
23000
120
11

Rock mass
44000
0.25
2700
40
13000
88
3.5

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² [Θ (GSI = 100, mi)]} = 13; R² [Θ (GSI = 100, )] = 0.98. The fitting curve is the black solid line in Figure 6. According to Hoek’s published tables [17, 20], the mi 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 mi value. Therefore, for rock mass, only one unknown parameter, GSI, needs to be determined. After the fitting calculation, we obtain GSI∗ = arg max {R² [Θ (GSI, mi = )]} = 88; R² [Θ (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.

Table 2: Hoek–Brown strength parameters.

Materials
GSI
mi
mb
s
a

Intact rock
100
13
13
1
0.5

Rock mass
88
13
8.5
0.3
0.5

3.3. Comparison with Hoek’s Published Tables

In the absence of test data, the values of GSI and mi are generally determined by Hoek’s published tables [17, 20]. Using the GSI and mi 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 site-specific 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].

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

Method
Θ
Intact rock
Rock mass

Determined by fitting the failure curves
GSI
100
88

mi
13
13

Estimated by Hoek’s published tables
GSI
100
50–60

mi
9–15
9–15

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² = 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² = 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.

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

Samples
Intact rock
Rock mass

UCS (MPa)
Laboratory measurements
120
88

Estimated by Hoek–Brown
120
62

Estimated by Mohr–Coulomb
111
72

σti (MPa)
Laboratory measurements
11
3.5

Estimated by Hoek–Brown
10
4.0

Estimated by Mohr–Coulomb
19
16

4. Testing and Simulation of the Strain-Softening Behaviour

4.1. Strain-Softening Behaviour Observed in Laboratory Tests

Figure 7 illustrates a stress-strain curve logged in a triaxial test. In the test, the rock experienced an elastic state, a strain-softening state, and a residual state [21] successively. In the elastic state, the stress increases with strain until the yield limit σ₁peak is reached. In the strain-softening state, the stress gradually decreases to the residual strength σ₁res. In the residual state, the stress fluctuates around the residual strength σ₁res.

Figure 7: Strain-softening behaviour observed in laboratory tests.

Figure 8 shows the relationship between σ₁peak, σ₁res, and confining pressure σ₃. The GSIpeak value of the elastic state and the GSIres 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)

Figure 8: Envelope of the yield limit and the residual strength. (a) Intact rock. (b) Rock mass.

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 strain-softening state lengthened, and the postpeak stress drop velocity slowed down.

Figure 9: The trend of the drop modulus M.

4.2. Strain-Softening Model

The rock samples drilled at the site show significant strain-softening behaviour in laboratory tests. Obviously, this stress-strain relationship cannot simply be represented by an elastoplastic model. Therefore, a strain-softening model was chosen, as shown in Figure 10.

(a)

(b)

Figure 10: Strain-softening constitutive model. (a) Stress-strain relationship. (b) Strength model.

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 strain-softening 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 coefficient given by

where ψ denotes dilatancy. The coefficient ψ could be determined by means of the following expression [17]:

The strength parameters mb, s, and a are assumed to be piecewise functions of the softening index η:

where ω represents a strength parameter that can be replaced by mb, s, or a. The peak value ωpeak is obtained by substituting GSI = GSIpeak into equation (2). The residual value ωres is obtained by substituting GSI = GSIres into equation (2).

4.3. 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 FLAC3D 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 (σ₃) were applied to the nodes of the cylindrical boundary. To simulate the compression process, a constant velocity boundary condition (10⁻⁸ m/s) was applied to the nodes at the top boundary.

The simulated stress-strain curves were compared with the measured stress-strain curves, as shown in Figure 11. In the figure, the measured stress-strain 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 strain-softening model works well in the numerical model.

Figure 11: Strain-softening behaviour and numerical reproduction (rock mass).

The shear-strain 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 shear-strain 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)

Figure 12: Failure forms in triaxial tests and simulations (rock mass, σ₃ = 10 MPa). (a) Laboratory test. (b) Shear-strain increment. (c) Displacement vector.

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.

Figure 13: Projection of principal stresses.

Table 5: Survey results of in situ stress.

No.
Depth (m)
σh,max
σh,min
σv

Value (MPa)
Direction
Dip angle (°)
Value (MPa)
Direction
Dip angle (°)
Value (MPa)
Direction
Dip angle (°)

S1
37
3.48
N84°E
8
1.43
S7°E
–5
1.25
S50°W
–100

S2
42
3.98
N69°E
13
1.67
S21°E
–2
1.60
S62°W
–103

S3
47
3.76
N65°E
4
1.79
S26°E
–6
1.55
S6°W
–97

S4
57
4.27
N77°E
1
1.91
S13°E
–3
1.83
S6°W
–93

S5
45
3.54
N84°E
5
2.01
S7°E
–13
1.74
S15°W
–104

S6
52
3.77
N82°E
12
2.3
S10°E
–11
2.09
S38°W
–106

S7
48
3.46
N82°E
7
1.8
S9°E
–9
1.65
S27°W
–102

S8
56
4.01
N87°E
9
2.12
S3°E
–1
1.89
S83°W
–99

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.

Figure 14: Trend of in situ stress along 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.

Figure 15: Simulation of the in situ stress.

6. Numerical Simulation of the Bedrock Sagging Sinkhole

6.1. Numerical Model

The embankment construction was simulated using FLAC3D 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 free-field 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.

Figure 16: Numerical model.

Figure 17: Embankment construction term.

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 strain-softening 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 strain-softening constitutive model. Figure 18 is a comparison between the computation results of 4 schemes and the measurements at site.

Figure 18: Settlement computation.

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 strain-softening model exhibit the same mechanical behaviour. After EC3, the strength of the strain-softening rock in scheme 4 decreases significantly. Therefore, the deformation of strain-softening 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 strain-softening 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 strain-softening 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 Strain-Softening Zone and Its Transition to a Residual Zone

As shown in Figure 19, the surrounding medium was divided into the elastic zone, strain-softening 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, η ≥ η∗.

(a)

(b)

(c)

(d)

Figure 19: The occurrence of the strain-softening zone and its transition to the residual zone.

The occurrence of the strain-softening 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 strain-softening 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 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. The strain-softening 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 strain-softening model.

(a)

(b)

(c)

Figure 20: Stress path of the bedrock. (a) The support at the left end. (b) The support at the right end. (c) Midspan.

7.2.1. The Stress Paths of the Strain-Softening 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, σ₁ of the two ends was maintained at approximately 0.5 MPa, and σ₃ 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 strain-softening behaviour, σ₃ 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, σ₁ increased along the loading path to approximately 6 MPa, and σ₃ was maintained at approximately 0.2 MPa. In stage 2, σ₃ 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 strain-softening behaviour, σ₃ 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 strain-softening 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 strain-softening model. The midspan still remained in an elastic state until the end of stage 4. Therefore, 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 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 strain-softening 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 strain-softening 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)

Figure 21: Bending moment of bedrock. (a) Stage 1. (b) Stage 2. (c) Stage 3. (d) Stage 4.

Table 6: Calculating schemes of the bending moment.

Scheme
Bending moment solution

Numerical simulation

Simply supported beam solution

Fixed supported beam solution

7.3.1. The 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. 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 strain-softening 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 strain-softening state and residual state emerged in the midspan. The bending moment calculated by the fixed supported beam model was far less than the strain-softening 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)

Figure 22: Evolution of the bedrock sagging sinkhole during the embankment construction. (a) Stage 1. (b) Stage 2. (c) Stage 3. (d) Stage 4.

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 strain-softening behaviour is one of the main causes of large deformation of 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 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)

Figure 23: Potential evolution from bedrock sagging sinkhole to bedrock collapse sinkhole. (a) Stage 1. (b) Stage 2. (c) Stage 3. (d) Stage 4.

8. Conclusions

In the initial stage of embankment construction, the surrounding medium remains in an elastic state. The 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 significantly, 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 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 human-induced geohazards and impacts in karst,” Earth-Science Reviews, vol. 138, pp. 61–88, 2014.

J. P. Galve, F. Gutiérrez, P. Lucha et al., “Sinkholes in the salt-bearing evaporite karst of the Ebro River valley upstream of Zaragoza city (NE Spain),” Geomorphology, vol. 108, no. 3-4, pp. 145–158, 2009.

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.

A. J. Rumbens, “Detection of cavities in karstic terrain. Road subsidence-snowy mountains highway near yarrangobilly, state of new South wales-Australia,” Exploration Geophysics, vol. 21, no. 1-2, pp. 121–124, 1990.

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. 135-136, no. 7, pp. 92–105, 2012.

H. Fan, Y. Zhang, S. He, K. Wang, X. Wang, and H. Wang, “Hazards and treatment of karst tunneling in Qinling-Daba mountainous area: overview and lessons learnt from Yichang–Wanzhou railway system,” Environmental Earth Sciences, vol. 77, no. 19, p. 679, 2018.

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.

G. Zhou, H. Yan, K. Chen, and R. Zhang, “Spatial analysis for susceptibility of second-time karst sinkholes: a case study of Jili Village in Guangxi, China,” Computers & Geosciences, vol. 89, pp. 144–160, 2016.

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.

F. Gutiérrez, J. P. Galve, P. Lucha, J. Bonachea, L. Jordá, and R. Jordá, “Investigation of a large collapse sinkhole affecting a multi-storey building by means of geophysics and the trenching technique (Zaragoza City, NE Spain),” Environmental Geology, vol. 58, no. 5, pp. 1107–1122, 2008.

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.

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.

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.

Q. Liu and M. Zheng, “The causes analysis and treatment of karst subgrade diseases—luxi section of Changbo-Jinyushi expressway,” Journal of East China Jiaotong University, vol. 31, no. 1, pp. 135–40, 2014.

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.

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.

E. Hoek and C. Carranza-Torres, “Hoek-Brown failure criterion—2002 edition,” in Proceedings of Fifth North American Rock Mechanics Symposium, pp. 18–22, Toronto, Canada, July 2002.

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.

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.

X. Wu, Y. Jiang, and Z. Guan, “A modified strain-softening model with multi-post-peak behaviours and its application in circular tunnel,” Engineering Geology, vol. 240, pp. 21–33, 2018.

E. Alonso, L. R. Alejano, F. Varas, G. Fdez-Manin, and C. Carranza-Torres, “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.

M. Parise and P. Lollino, “A preliminary analysis of failure mechanisms in karst and man-made underground caves in Southern Italy,” Geomorphology, vol. 134, no. 1-2, pp. 132–143, 2011.

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.

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.

D. Wigham, “Occurrence of mining-induced open fissures and shear walls in the Permian limestones of County Durham,” Mining Technology, vol. 109, no. 3, pp. 172–178, 2013.

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.

P. Van Beynen and K. Townsend, “A disturbance index for karst environments,” Environmental Management, vol. 36, no. 1, pp. 101–116, 2005.

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: west-central Florida and southeast Italy,” Journal of Environmental Management, vol. 90, no. 5, pp. 1770–1781, 2009.