# Constraint Embankment Construction to Prevent the Collapse of Underground Caves

Authors: Zhen Zhang; Zhongda Chen
Journal: Advances in Civil Engineering
 (2019)

---

## Abstract

Dozens of underground karst caves were found before constructing the Changli highway. The thickness-to-span ratio of nearly half of the caves is less than 0.05, and the greatest ratio is only 0.35, far less than the value demanded by local construction specifications (0.8). The caves located at K50 + 700 and K178 + 800 are by far the only two caves that have become unstable. Only one passive measure was taken when constructing the highway, i.e., building 0.5 m thick continuous reinforced concrete slabs above the embankment; this measure did not contribute to the improvement of the stability of the underground caves. Numerical solutions based on strength reduction and analytical solutions based on the beam hypothesis are used to assess the stability of underground caves. The capacity of an underground cave to bear embankment construction is observed to be proportional to the tensile strength of the rock mass and the square of the thickness-to-span ratio of the cave roof. The tensile strength of the rock mass isψ times lower than that of the intact rock. The value of ψ is mainly determined by the geological strength index (GSI). To prevent instability of underground caves, the embankment height should be reasonably controlled. However, local construction specifications requiring that the thickness-to-span ratio of underground cave be greater than 0.8 are conservative.


---

## Body

## 1. Introduction

Embankment construction may pose potential hazards to underground caves, especially those caves with a wide span and thin roof. Much concern has been paid to prevent the collapse of underground caves and failure of surface civil engineering works.In practice, the thickness-to-span ratio (i.e., the ratio of the roof thickness to the cave span) is a simple indicator to determine whether the underground cave is stable. Chinese specifications for the design of highway embankments (JTG D30-2015) demand that engineering actions be taken in an underground cave when its thickness-to-span ratio is less than 0.8; this excludes disintegrated rock mass. An informal guideline in the United Kingdom suggests that the underground cave would be stable if the roof thickness is equal to or greater than the cave span (i.e., thickness-to-span ratio ≥ 1.0). However, many civil engineering works do stand on caves with much thinner roofs than recommended by the guidelines [1]. After full-scale loading tests in Nottingham and physical modelling tests in the laboratory, Waltham and Swift [1] suggested that the underground cave with a thickness-to-span ratio greater than 0.7 would be adequate for bearing most surface civil engineering works. The thickness-to-span ratio indicator has allowed for simple estimation of the stability of underground caves that is easy to implement in engineering practice. The drawback of this approach is that it makes no distinction regarding the rock strength and the load size, each of which has a great influence on the stability [2–4].The antibend estimation method based on the assumption of a simply supported beam is one of the recommended methods of the Chinese technical guidance for highway foundation design and construction in a karstified area. The antibend estimation method can be used to judge the stability of the underground cave with different amounts of rock strength and load size; however, it provides a worst case with least bearing capacity because it ignores the fact that the rock at the margins can provide powerful constraints on the bending moment [5]. The analytical solution based on the simply supported beam hypothesis is therefore conservative.In this paper, the stability of underground caves is assessed by numerical solutions based on the strength reduction and analytical solutions based on the beam hypothesis. Several procedures to perform the strength reduction for the Hoek–Brown criterion are discussed, two of which are selected and used for case studies. An approximate relationship between the tensile strength of the rock mass and the unconfined compressive strengthσci, the rock type mi, the disturbance index D, and the geological strength index (GSI) is derived and used as the basic parameter in analytical solutions. The analytical solutions based on beam hypothesis are deduced and compared with the numerical solutions. The feasibility to derive the thickness-to-span ratio from the analytical solution is also discussed.


## 2. Hoek–Brown Criterion

The Hoek–Brown criterion is frequently used to depict the rock strength [6, 7] and can be expressed as follows:(1)fσ1,σ3,GSI=σ1−σ3−σcimbσ3σci+sa,where σ1 and σ3 are the major and minor principal stresses, respectively, and compression is taken as positive. The strength parameters mb, s, and a are functions of the geological strength index GSI:(2)mb=expGSI−10028−14Dmi,s=expGSI−1009−3D,a=12+16exp−GSI15−exp−203,where mi represents the intact rock material constant and D represents the disturbance factor.


## 3. Strength Reduction

### 3.1. A Brief Review

Strength reduction is a convenient tool for calculating stability. Currently, there are four different procedures available to perform the strength reduction for the Hoek–Brown criterion.Method 1 (reduction of equivalent c-φ values).
First, equivalentc-φ values were determined for a set of Hoek–Brown parameters using the equations described by Hoek and Carranza-Torres [8]. In this approach, the area below the Mohr–Coulomb plot (line A in Figure 1) and the area below the Hoek–Brown plot (line B in Figure 1) are balanced by using the equivalent c-φ values:(3)φ=sin−13ambs+mbσ3na−11+a2+a+3ambs+mbσ3na−1,(4)c=σci1+2as+1−ambσ3ns+mbσ3na−11+a2+a1+a2+a+6ambs+mbσ3na−1,where σ3n=σ3max/σci and σ3max is the upper limit of the confining stress. The value of σ3max can be determined by searching the maximum σ3 in the case.
Subsequently, the equivalentc-φ values were used for executing global c-φ-reduction by the classical approach [7, 9–11]. The drawback of this method is that the equivalent c-φ values ignore the local stress state in a specified element.Figure 1
Hoek–Brown criteria and the equivalent Mohr–Coulomb criteria.Method 2 (reduction of local c-φ values).
If the stress states (σ1, σ2, and σ3), particularly σ3, of an element is known, then the corresponding instantaneous c and φ values can be calculated. These instantaneous c and φ values are used for executing local c-φ-reduction (e.g., [12–14]). The key problem of this method is how to determine the instantaneous c and φ values. Fu and Liao derived the instantaneous c and φ values by solving the Newton iterative formula. Despite the complexity of the technique, it can be applied to any nonlinear failure criteria. A recommended simpler method to obtain the instantaneous c and φ values is to use formulas (3) and (4) by setting σ3n=σci−1⋅σ3, as done by Shen and Karakus [14].Method 3 (lowering of equivalent shear-normal envelope).
First, based on relationships developed by previous research studies [8], a shear-normal stress envelope (line A in Figure 2) equivalent to the Hoek–Brown principal stress envelope can be determined. Next, the equivalent shear-normal stress envelope is lowered by a factor of safety (line B in Figure 2), as done by Hammah et al. [15]. More recently, a modified method to conduct the shear strength reduction has been proposed by Clausen and Damkilde [16].Figure 2
Lowering of equivalent shear-normal envelope.Method 4 (reduction of Hoek–Brown parameters).
This method selects one or several Hoek–Brown parameters (e.g., GSI,σci, and mi) to perform strength reduction. The key problem of this method is that the selection of Hoek–Brown parameters influences the factor of safety (FOS). Song et al. [17] found that the FOS value obtained by reduction of GSI and σci converges with that obtained by the limit equilibrium method.
In methods 1 and 3, the equivalent envelope is lowered by a FOS. Then, to determine the reduced Hoek–Brown parametersGSI′ and mi′, the lowered envelope should be fitted to a Hoek–Brown curve using a nonlinear regression algorithm. Therefore, the reduced Hoek–Brown shear curve is an approximate fit to the lowered envelope. It may result in some errors in determining the reduced Hoek–Brown parameters in the fitting process [18]. In method 2, instantaneous c and φ values have to be deduced in each element. Therefore, the safety factor result determined by method 2 is believed to be more accurate than that determined by method 1 [18]. In method 4, neither equivalent nor fitting is needed. Therefore, method 4 provides a more convenient way in practical use than methods 1–3.


### 3.2. Discussion on the Parameter Reduction Schemes

Two different parameter reduction schemes are compared. Both schemes belong to aforementioned method 4.Scheme 1 (reduction of σci and GSI).
The parametersσci and GSI are reduced by FOS, and the reduced parameters σci′ and GSI′ are obtained as follows:(5)σci′=σciFOS,GSI′=GSIFOS.
SubstituteGSI′ into equation (2) to obtain the reduced parameters mb′, s′, and a′:(6)mb′=expGSI−100FOS28−14DFOSmi,s′=expGSI−100FOS9−3DFOS,a′=12+16exp−GSI15FOS−exp−203.Scheme 2 (reduction of σci and expGSI/K).
The parametersσci and expGSI/K are reduced by FOS, and the reduced parameters σci′ and GSI′ are obtained as follows:(7)σci′=σciFOS,GSI′=GSI−KlnFOS,where K=28−14D+9−3D+15/3=52−17D/3.
SubstituteGSI′ into equation (2) to obtain the reduced parameters mb′, s′, and a′:(8)mb′=1FOSK/28−14DexpGSI−10028−14Dmi,s′=1FOSK/9−3DexpGSI−1009−3D,a′=12+161FOS−K/15exp−GSI15−exp−203.
The relationships of instantaneousGSI′, mb′, s′, and a′ with FOS can be used to highlight the reduction process, as shown in Figure 3. The following parameters were used for the calculation: GSI = 60 and D = 0.1; the values of FOS range from 0.8 to 5.0.
Using scheme 1,GSI′, mb′, and s′ descend sharply when FOS ranges from 0.8 to approximately 1.5. Using scheme 2, the sharp curves are significantly eased. In Figure 3(d), the nonlinear growth trend of the parameter a′ becomes approximately linear using scheme 2.Figure 3
The relationship between (a)GSI′ and FOS; (b) mb′ and FOS; (c) s′ and FOS; (d) a′ and FOS.
(a)
(b)
(c)
(d)


### 3.3. Comparison of the Strength Reduction Methods

To validate the strength reduction method used in this work, safety factor analysis was carried out for some examples described in literatures.Example 1 (rock slope).
The factors of safety for 10 m high shale rock slopes [18] with 45° slope angle were calculated with the aforementioned strength reduction methods, as shown in Figure 4(a). The results are illustrated in Figure 4(b). The parameters are as follows: Young’s modulus = 5,000 MPa, Poisson’s ratio = 0.3, unit weight = 0.025 MN/m3, σci = 40 MPa, mi = 2, and D = 0. GSI values range from 5 to 65. In the numerical model, the bottom boundary is fixed in the vertical and horizontal directions and the left and right boundaries are fixed in the horizontal direction.Figure 4
The safety factors of the rock slope. (a) The illustration of the rock slope and (b) the relationship between FOS and GSI.
(a)
(b)Example 2 (rectangular cave).
The factors of safety for 3 m wide and 4 m high rectangular caves were calculated with the aforementioned strength reduction methods, as shown in Figure5(a). The results are illustrated in Figure 5(b). The parameters are as follows: Young’s modulus = 12,000 MPa, Poisson’s ratio = 0.3, unit weight = 0.023 MN/m3, σci = 18 MPa, mi = 8, and D = 0.5. GSI values range from 40 to 100. In the numerical model, the bottom boundary was fixed in the vertical and horizontal directions and the left and right boundaries were fixed in the horizontal direction.
The results show that the FOS values derived from scheme 2 of method 4 are very close to those from the methods 1–3. However, the FOS values derived from scheme 1 of method 4 are significantly smaller than those from other methods. If method 4 is used, the author recommends that the FOS values should be determined by scheme 2 therein.Figure 5
The safety factors of the rectangular cave. (a) The illustration of the rectangular cave and (b) the relationship between FOS and GSI.
(a)
(b)


## 3.1. A Brief Review

Strength reduction is a convenient tool for calculating stability. Currently, there are four different procedures available to perform the strength reduction for the Hoek–Brown criterion.Method 1 (reduction of equivalent c-φ values).
First, equivalentc-φ values were determined for a set of Hoek–Brown parameters using the equations described by Hoek and Carranza-Torres [8]. In this approach, the area below the Mohr–Coulomb plot (line A in Figure 1) and the area below the Hoek–Brown plot (line B in Figure 1) are balanced by using the equivalent c-φ values:(3)φ=sin−13ambs+mbσ3na−11+a2+a+3ambs+mbσ3na−1,(4)c=σci1+2as+1−ambσ3ns+mbσ3na−11+a2+a1+a2+a+6ambs+mbσ3na−1,where σ3n=σ3max/σci and σ3max is the upper limit of the confining stress. The value of σ3max can be determined by searching the maximum σ3 in the case.
Subsequently, the equivalentc-φ values were used for executing global c-φ-reduction by the classical approach [7, 9–11]. The drawback of this method is that the equivalent c-φ values ignore the local stress state in a specified element.Figure 1
Hoek–Brown criteria and the equivalent Mohr–Coulomb criteria.Method 2 (reduction of local c-φ values).
If the stress states (σ1, σ2, and σ3), particularly σ3, of an element is known, then the corresponding instantaneous c and φ values can be calculated. These instantaneous c and φ values are used for executing local c-φ-reduction (e.g., [12–14]). The key problem of this method is how to determine the instantaneous c and φ values. Fu and Liao derived the instantaneous c and φ values by solving the Newton iterative formula. Despite the complexity of the technique, it can be applied to any nonlinear failure criteria. A recommended simpler method to obtain the instantaneous c and φ values is to use formulas (3) and (4) by setting σ3n=σci−1⋅σ3, as done by Shen and Karakus [14].Method 3 (lowering of equivalent shear-normal envelope).
First, based on relationships developed by previous research studies [8], a shear-normal stress envelope (line A in Figure 2) equivalent to the Hoek–Brown principal stress envelope can be determined. Next, the equivalent shear-normal stress envelope is lowered by a factor of safety (line B in Figure 2), as done by Hammah et al. [15]. More recently, a modified method to conduct the shear strength reduction has been proposed by Clausen and Damkilde [16].Figure 2
Lowering of equivalent shear-normal envelope.Method 4 (reduction of Hoek–Brown parameters).
This method selects one or several Hoek–Brown parameters (e.g., GSI,σci, and mi) to perform strength reduction. The key problem of this method is that the selection of Hoek–Brown parameters influences the factor of safety (FOS). Song et al. [17] found that the FOS value obtained by reduction of GSI and σci converges with that obtained by the limit equilibrium method.
In methods 1 and 3, the equivalent envelope is lowered by a FOS. Then, to determine the reduced Hoek–Brown parametersGSI′ and mi′, the lowered envelope should be fitted to a Hoek–Brown curve using a nonlinear regression algorithm. Therefore, the reduced Hoek–Brown shear curve is an approximate fit to the lowered envelope. It may result in some errors in determining the reduced Hoek–Brown parameters in the fitting process [18]. In method 2, instantaneous c and φ values have to be deduced in each element. Therefore, the safety factor result determined by method 2 is believed to be more accurate than that determined by method 1 [18]. In method 4, neither equivalent nor fitting is needed. Therefore, method 4 provides a more convenient way in practical use than methods 1–3.


## 3.2. Discussion on the Parameter Reduction Schemes

Two different parameter reduction schemes are compared. Both schemes belong to aforementioned method 4.Scheme 1 (reduction of σci and GSI).
The parametersσci and GSI are reduced by FOS, and the reduced parameters σci′ and GSI′ are obtained as follows:(5)σci′=σciFOS,GSI′=GSIFOS.
SubstituteGSI′ into equation (2) to obtain the reduced parameters mb′, s′, and a′:(6)mb′=expGSI−100FOS28−14DFOSmi,s′=expGSI−100FOS9−3DFOS,a′=12+16exp−GSI15FOS−exp−203.Scheme 2 (reduction of σci and expGSI/K).
The parametersσci and expGSI/K are reduced by FOS, and the reduced parameters σci′ and GSI′ are obtained as follows:(7)σci′=σciFOS,GSI′=GSI−KlnFOS,where K=28−14D+9−3D+15/3=52−17D/3.
SubstituteGSI′ into equation (2) to obtain the reduced parameters mb′, s′, and a′:(8)mb′=1FOSK/28−14DexpGSI−10028−14Dmi,s′=1FOSK/9−3DexpGSI−1009−3D,a′=12+161FOS−K/15exp−GSI15−exp−203.
The relationships of instantaneousGSI′, mb′, s′, and a′ with FOS can be used to highlight the reduction process, as shown in Figure 3. The following parameters were used for the calculation: GSI = 60 and D = 0.1; the values of FOS range from 0.8 to 5.0.
Using scheme 1,GSI′, mb′, and s′ descend sharply when FOS ranges from 0.8 to approximately 1.5. Using scheme 2, the sharp curves are significantly eased. In Figure 3(d), the nonlinear growth trend of the parameter a′ becomes approximately linear using scheme 2.Figure 3
The relationship between (a)GSI′ and FOS; (b) mb′ and FOS; (c) s′ and FOS; (d) a′ and FOS.
(a)
(b)
(c)
(d)


## 3.3. Comparison of the Strength Reduction Methods

To validate the strength reduction method used in this work, safety factor analysis was carried out for some examples described in literatures.Example 1 (rock slope).
The factors of safety for 10 m high shale rock slopes [18] with 45° slope angle were calculated with the aforementioned strength reduction methods, as shown in Figure 4(a). The results are illustrated in Figure 4(b). The parameters are as follows: Young’s modulus = 5,000 MPa, Poisson’s ratio = 0.3, unit weight = 0.025 MN/m3, σci = 40 MPa, mi = 2, and D = 0. GSI values range from 5 to 65. In the numerical model, the bottom boundary is fixed in the vertical and horizontal directions and the left and right boundaries are fixed in the horizontal direction.Figure 4
The safety factors of the rock slope. (a) The illustration of the rock slope and (b) the relationship between FOS and GSI.
(a)
(b)Example 2 (rectangular cave).
The factors of safety for 3 m wide and 4 m high rectangular caves were calculated with the aforementioned strength reduction methods, as shown in Figure5(a). The results are illustrated in Figure 5(b). The parameters are as follows: Young’s modulus = 12,000 MPa, Poisson’s ratio = 0.3, unit weight = 0.023 MN/m3, σci = 18 MPa, mi = 8, and D = 0.5. GSI values range from 40 to 100. In the numerical model, the bottom boundary was fixed in the vertical and horizontal directions and the left and right boundaries were fixed in the horizontal direction.
The results show that the FOS values derived from scheme 2 of method 4 are very close to those from the methods 1–3. However, the FOS values derived from scheme 1 of method 4 are significantly smaller than those from other methods. If method 4 is used, the author recommends that the FOS values should be determined by scheme 2 therein.Figure 5
The safety factors of the rectangular cave. (a) The illustration of the rectangular cave and (b) the relationship between FOS and GSI.
(a)
(b)


## 4. Analytical Solutions Based on the Beam Hypothesis

The antibend estimation method based on the assumption of a simply supported beam is one of the recommended methods of the Chinese technical guidance for highway foundation design and construction in a karstified area. This approach assumes that the cave roof behaves like a simply supported beam. The tensile strength is an important parameter to determine the stability of the cave roof. This paper derives the relationships between the tensile strength of rock mass and the unconfined compressive strengthσci, the rock type mi, and the geological strength index (GSI). Next, based on the assumption of simply supported beam, the analytical solution of the ultimate bearing height of the embankment is deduced. Finally, another analytical solution based on the assumption of fix supported beam is proposed. The influence of the in situ stress and roof inclination is also taken into account.


### 4.1. Tensile Strength of Rock Mass

Because of the existence of joints and fissures, the tensile strength of rock mass in engineering practice is usually lower than the measured tensile strength of the intact rock in the laboratory.A trend [19–22] is to estimate the tensile strength of rock mass by using the Hoek–Brown criterion. In this paper, the deduce process is as follows.Substitutingσ1 = 0 and σti = −σ3 into the Hoek–Brown criterion, it can be transformed as follows:(9)σtiσci1/a+mbσtiσci−s=0.For a disintegrated rock mass (the value of GSI is generally smaller than 30), its tensile strength is negligible. For a blocky rock mass (the value of GSI is generally greater than 30), it can be observed from equation (2) that the value of a ranges between 0.5 and 0.52. Therefore, it would be reasonable to assume that a ≈ 0.5 for a blocky rock mass. Substituting a ≈ 0.5 into equation (9), the following equation can be obtained:(10)σti≈mb21+4smb2−1σci.Because4s/mb2 is generally close to 0, the use of an equivalent infinitesimal substitute can simplify the equation and solve the problem easily:(11)σti≈mb22smb2σci=smbσci.Based on the Hoek–Brown 1997 criterion, Cai [20] deduced the range of the ratio of the tensile strength to the compressive strength, i.e., mi+1−1≤σti⋅σci−1<mi−1. Based on the Hoek–Brown 2002 criterion, it can be deduced that s⋅mb+1−1≤σti⋅σci−1<s⋅mb−1. Comparing equation (11) with Cai’s research results, it can be observed that the error caused by the estimation based on equation (11) is quite small. The maximum error does not exceed mb−1mb+1−1sσci, and its accuracy is adequate to meet the engineering requirements.The GSI value of the intact rock in the laboratory is known as 100. Substituting GSI = 100 into equation (2), it can be observed that s = 1.0 and mb = mi. In this case, equation (11) can be approximated as follows:(12)σti≈1miσci.Themi value ranges from 4 to 33 for some commonly encountered rocks in engineering practice [20]. According to equation (12), the tensile strength of the intact rock is 1/33 to 1/4 of its compressive strength. This result is roughly consistent with the published test results by Sheorey [23], i.e., 1/39 to 1/7.For a rock mass in engineering practice, the GSI value is always less than 100. Substituting equation (2) into equation (11), we obtain the following equation:(13)σti≈exp0.01GSI−1ησcimi,where η=2.383−D−12−D−119−11D.As seen from Figure6, the trend of η is quite nonlinear in the entire positive number domain D > 0. However, when D ranges from 0 to 1, a nearly linear trend is observed, as shown in Figure 7. In fact, according to Hoek’s definition [8], D only takes on values from 0 to 1. Using the linear approximation η ≈ 2D + 7.54 can greatly simplify the original equation η=2.383−D−12−D−119−11D. Using this linear approximation, the following equation can be used to estimate the tensile strength of rock mass:(14)σti≈1ψmiσci,where ψ=exp2D+7.541−0.01GSI.Figure 6
Trend ofη.Figure 7
Approximation ofη.By comparing equations (12) and (14), the tensile strength of rock mass is found to be ψ times lower than that of intact rock. Figure 8 shows that the value of ψ increases when lowering GSI and where construction disturbances are imposed. For typical rock masses that suffer significant disturbance due to heavy production blasting (D = 1), the value of ψ rises from 1 (GSI = 100) to 6 (GSI = 80). If the GSI value is less than 30, then the tensile strength of rock mass is at least 200 times lower than that of the intact rock, in which case, its tensile strength can be ignored.Figure 8
Trend ofψ, assuming Hoek–Brown criterion.It was encouraging to find that the value ofψ can be determined by D and GSI with a simple function. But the joy was very short because soon it was noticed that the value of ψ may have been exaggerated when compared with previous research studies. After direct tension tests on natural siltstone and medium grained sandstone samples containing joints, bedding, and mineral veins, Shang et al. [24, 25] confirmed that the tensile strength of rock mass is up to 3 times lower than that of the intact rock. The tensile strength of Precambrian granitic rock masses in Harrat Lunayyir is estimated as 1–3 MPa, averagely 1/7 times that of the intact rock samples [26], almost an order of magnitude higher than that estimated by the Hoek–Brown criterion.Following the research conducted by Aydan and Kawamoto [27], Tokashiki and Aydan [28] proposed an empirical function to estimate the value of ψ:(15)ψ=120−RMR0.2RMR,where RMR represents the rock mass rating.Set the groundwater level to 15 (the rock is completely dry) and the level adjustment to 0 (very favorable joint orientations), Hoek [29] developed the following function to link RMR and GSI:(16)GSI=RMR−5.Therefore, the value ofψ can be estimated as follows:(17)ψ=115−GSI1+0.2GSI.Values ofψ estimated by the Hoek–Brown criterion and Tokashiki are summarized in Figure 9. For conventionally encountered rock mass and under no construction disturbance (GSI ≥ 50; D = 0), the value of ψ estimated by the Hoek–Brown criterion ranges from 1 to 43 while that estimated by Tokashiki ranges from 1 to 6. In terms of numerical values, the estimation proposed by Tokashiki is more consistent with previous research studies (3 to 7).Figure 9
Trend ofψ.


### 4.2. Analytical Solution Based on Simply Supported Beam Hypothesis

Assuming that cave span isl and uniform load is q, the maximum moment MM of the cave roof is(18)MM=ql28.To prevent failure, the tensile strengthσti should be greater than the tensile stress:(19)σti≥6MMbrhr2,where hr represents the thickness of the cave roof and br represents the width of the cave roof.The load applied to the cave roof can be expressed as follows:(20)q=γsbrhs+hf+γrbrhr,where hs represents the thickness of the soil stratum over the cave roof; hf represents the construction height of the embankment; γs represents the soil weight; γr represents the rock weight.Substituting equations (18) and (20) into equation (19), the following equation can be obtained:(21)hf≤4σti3γshr2l2−γrγshr−hs.According to equation (21), the ultimate bearing height of the embankment is proportional to the tensile strength of the rock mass and the square of the thickness-to-span ratio of the cave roof.The simply supported beam hypothesis provides a worst case with least bearing capacity because it ignores the fact that the rock at the margins can provide powerful constraints on the bending moment [5]. Because of the impact of in situ stress, the tensile stress in engineering practice is significantly less than the hypothetical case in which only gravity is applied. The analytical solution based on the simply supported beam hypothesis is therefore conservative.


### 4.3. Modified Solution

Karstic bedrock is generally extremely irregular in geometry because of the subsoil dissolution at the soil/rock interface [5]. The tilted cave roof is therefore quite commonly observed. To address this situation, assume that the cave roof is a fixed-supported tilt beam and then correct the bending moment by the force method, as shown in Figure 10. The influence of in situ stress is also considered.Figure 10
Assumption of the modified solution.According to the superposition principle, the load applied to the cave roof is divided into the force applied on the cave roof and the moment applied by the fixed supports. Summing the bending moments generated by the two types of load yields:(22)Mx=MP+ML−ML+MRlcosθx,where MP is the bending moment generated by the force applied on the cave roof and is calculated by ∫Fsxdx. The parameter Fsx represents the vertical shear force. The parameter x represents the position coordinate. The parameter θ represents the tilt angle of the cave roof. ML and MR are the bending moments at the left and right supports of the cave roof, respectively. Using the force method, the following equations are obtained:(23)δ11ML+δ12MR+Δ1P=0,δ21ML+δ22MR+Δ2P=0,where δ11 and δ12 are the rotation angles of the left and right supports, respectively, of the cave roof generated by the action of the unit moment ML = 1; δ21 and δ22 are the rotation angles of the left and right supports, respectively, of the cave roof generated by the action of the unit moment MR = 1; ∆1P and ∆2P are the rotation angles of the left and right supports, respectively, of the cave roof generated by the action of the unit uniform force qx=1.The bending moment of the cave roof generated by the action of the unit momentML = 1 can be determined as follows:(24)M1¯=1−1lcosθx.The bending moment of the cave roof generated by the action of the unit momentMR = 1 can be determined as follows:(25)M2¯=−1lcosθx.The bending moment generated by the action of forceq can be determined as follows:(26)MP=∫0x12ql−qxdx=12qxl−x.Thus, the parameters in equation (23) can be determined:(27)δ11=1EI∫0lcosθM1¯2dx=lcosθ3EI,(28)δ22=1EI∫0lcosθM2¯2dx=lcosθ3EI,(29)Δ1P=1EI∫M1¯MPdx=124EIql3cos2θ2−cosθ,(30)Δ2P=1EI∫M2¯MPdx=124EIql3cos2θ3cosθ−4,where E represents the elastic modulus and I represents the moment of inertia.Substituting equations (27)–(30) into equation (23) yields the following equation:(31)ML=−112ql2cos2θ,MR=112ql2cosθ6−5cosθ.By equation (19), the total bending moment Mx can be determined as follows:(32)Mx=−12qx2−lcosθx+16l2cos2θ.The maximum momentMM of the cave roof is(33)MM=124ql21−4sin2θ.Assume that the in situ stress in the span direction isσs. Because of the presence of the cave, stress concentration occurs around it. Assume that the stress concentration factor at the centre of the cave roof is λ. Thus, equation (19) can be modified as follows:(34)σti≥6MMbrhr2−λσs.According to the elastic mechanic [30], the factor λ can be determined:(35)λ=1+2hcl,where hc represents the cave height.Substituting equation (33) into equation (34) yields the following equation, which can be used to estimate the ultimate bearing height of the embankment above the underground cave:(36)hf≤41−4sin2θσti+λσsγshr2l2−γrγshr−hs.


## 4.1. Tensile Strength of Rock Mass

Because of the existence of joints and fissures, the tensile strength of rock mass in engineering practice is usually lower than the measured tensile strength of the intact rock in the laboratory.A trend [19–22] is to estimate the tensile strength of rock mass by using the Hoek–Brown criterion. In this paper, the deduce process is as follows.Substitutingσ1 = 0 and σti = −σ3 into the Hoek–Brown criterion, it can be transformed as follows:(9)σtiσci1/a+mbσtiσci−s=0.For a disintegrated rock mass (the value of GSI is generally smaller than 30), its tensile strength is negligible. For a blocky rock mass (the value of GSI is generally greater than 30), it can be observed from equation (2) that the value of a ranges between 0.5 and 0.52. Therefore, it would be reasonable to assume that a ≈ 0.5 for a blocky rock mass. Substituting a ≈ 0.5 into equation (9), the following equation can be obtained:(10)σti≈mb21+4smb2−1σci.Because4s/mb2 is generally close to 0, the use of an equivalent infinitesimal substitute can simplify the equation and solve the problem easily:(11)σti≈mb22smb2σci=smbσci.Based on the Hoek–Brown 1997 criterion, Cai [20] deduced the range of the ratio of the tensile strength to the compressive strength, i.e., mi+1−1≤σti⋅σci−1<mi−1. Based on the Hoek–Brown 2002 criterion, it can be deduced that s⋅mb+1−1≤σti⋅σci−1<s⋅mb−1. Comparing equation (11) with Cai’s research results, it can be observed that the error caused by the estimation based on equation (11) is quite small. The maximum error does not exceed mb−1mb+1−1sσci, and its accuracy is adequate to meet the engineering requirements.The GSI value of the intact rock in the laboratory is known as 100. Substituting GSI = 100 into equation (2), it can be observed that s = 1.0 and mb = mi. In this case, equation (11) can be approximated as follows:(12)σti≈1miσci.Themi value ranges from 4 to 33 for some commonly encountered rocks in engineering practice [20]. According to equation (12), the tensile strength of the intact rock is 1/33 to 1/4 of its compressive strength. This result is roughly consistent with the published test results by Sheorey [23], i.e., 1/39 to 1/7.For a rock mass in engineering practice, the GSI value is always less than 100. Substituting equation (2) into equation (11), we obtain the following equation:(13)σti≈exp0.01GSI−1ησcimi,where η=2.383−D−12−D−119−11D.As seen from Figure6, the trend of η is quite nonlinear in the entire positive number domain D > 0. However, when D ranges from 0 to 1, a nearly linear trend is observed, as shown in Figure 7. In fact, according to Hoek’s definition [8], D only takes on values from 0 to 1. Using the linear approximation η ≈ 2D + 7.54 can greatly simplify the original equation η=2.383−D−12−D−119−11D. Using this linear approximation, the following equation can be used to estimate the tensile strength of rock mass:(14)σti≈1ψmiσci,where ψ=exp2D+7.541−0.01GSI.Figure 6
Trend ofη.Figure 7
Approximation ofη.By comparing equations (12) and (14), the tensile strength of rock mass is found to be ψ times lower than that of intact rock. Figure 8 shows that the value of ψ increases when lowering GSI and where construction disturbances are imposed. For typical rock masses that suffer significant disturbance due to heavy production blasting (D = 1), the value of ψ rises from 1 (GSI = 100) to 6 (GSI = 80). If the GSI value is less than 30, then the tensile strength of rock mass is at least 200 times lower than that of the intact rock, in which case, its tensile strength can be ignored.Figure 8
Trend ofψ, assuming Hoek–Brown criterion.It was encouraging to find that the value ofψ can be determined by D and GSI with a simple function. But the joy was very short because soon it was noticed that the value of ψ may have been exaggerated when compared with previous research studies. After direct tension tests on natural siltstone and medium grained sandstone samples containing joints, bedding, and mineral veins, Shang et al. [24, 25] confirmed that the tensile strength of rock mass is up to 3 times lower than that of the intact rock. The tensile strength of Precambrian granitic rock masses in Harrat Lunayyir is estimated as 1–3 MPa, averagely 1/7 times that of the intact rock samples [26], almost an order of magnitude higher than that estimated by the Hoek–Brown criterion.Following the research conducted by Aydan and Kawamoto [27], Tokashiki and Aydan [28] proposed an empirical function to estimate the value of ψ:(15)ψ=120−RMR0.2RMR,where RMR represents the rock mass rating.Set the groundwater level to 15 (the rock is completely dry) and the level adjustment to 0 (very favorable joint orientations), Hoek [29] developed the following function to link RMR and GSI:(16)GSI=RMR−5.Therefore, the value ofψ can be estimated as follows:(17)ψ=115−GSI1+0.2GSI.Values ofψ estimated by the Hoek–Brown criterion and Tokashiki are summarized in Figure 9. For conventionally encountered rock mass and under no construction disturbance (GSI ≥ 50; D = 0), the value of ψ estimated by the Hoek–Brown criterion ranges from 1 to 43 while that estimated by Tokashiki ranges from 1 to 6. In terms of numerical values, the estimation proposed by Tokashiki is more consistent with previous research studies (3 to 7).Figure 9
Trend ofψ.


## 4.2. Analytical Solution Based on Simply Supported Beam Hypothesis

Assuming that cave span isl and uniform load is q, the maximum moment MM of the cave roof is(18)MM=ql28.To prevent failure, the tensile strengthσti should be greater than the tensile stress:(19)σti≥6MMbrhr2,where hr represents the thickness of the cave roof and br represents the width of the cave roof.The load applied to the cave roof can be expressed as follows:(20)q=γsbrhs+hf+γrbrhr,where hs represents the thickness of the soil stratum over the cave roof; hf represents the construction height of the embankment; γs represents the soil weight; γr represents the rock weight.Substituting equations (18) and (20) into equation (19), the following equation can be obtained:(21)hf≤4σti3γshr2l2−γrγshr−hs.According to equation (21), the ultimate bearing height of the embankment is proportional to the tensile strength of the rock mass and the square of the thickness-to-span ratio of the cave roof.The simply supported beam hypothesis provides a worst case with least bearing capacity because it ignores the fact that the rock at the margins can provide powerful constraints on the bending moment [5]. Because of the impact of in situ stress, the tensile stress in engineering practice is significantly less than the hypothetical case in which only gravity is applied. The analytical solution based on the simply supported beam hypothesis is therefore conservative.


## 4.3. Modified Solution

Karstic bedrock is generally extremely irregular in geometry because of the subsoil dissolution at the soil/rock interface [5]. The tilted cave roof is therefore quite commonly observed. To address this situation, assume that the cave roof is a fixed-supported tilt beam and then correct the bending moment by the force method, as shown in Figure 10. The influence of in situ stress is also considered.Figure 10
Assumption of the modified solution.According to the superposition principle, the load applied to the cave roof is divided into the force applied on the cave roof and the moment applied by the fixed supports. Summing the bending moments generated by the two types of load yields:(22)Mx=MP+ML−ML+MRlcosθx,where MP is the bending moment generated by the force applied on the cave roof and is calculated by ∫Fsxdx. The parameter Fsx represents the vertical shear force. The parameter x represents the position coordinate. The parameter θ represents the tilt angle of the cave roof. ML and MR are the bending moments at the left and right supports of the cave roof, respectively. Using the force method, the following equations are obtained:(23)δ11ML+δ12MR+Δ1P=0,δ21ML+δ22MR+Δ2P=0,where δ11 and δ12 are the rotation angles of the left and right supports, respectively, of the cave roof generated by the action of the unit moment ML = 1; δ21 and δ22 are the rotation angles of the left and right supports, respectively, of the cave roof generated by the action of the unit moment MR = 1; ∆1P and ∆2P are the rotation angles of the left and right supports, respectively, of the cave roof generated by the action of the unit uniform force qx=1.The bending moment of the cave roof generated by the action of the unit momentML = 1 can be determined as follows:(24)M1¯=1−1lcosθx.The bending moment of the cave roof generated by the action of the unit momentMR = 1 can be determined as follows:(25)M2¯=−1lcosθx.The bending moment generated by the action of forceq can be determined as follows:(26)MP=∫0x12ql−qxdx=12qxl−x.Thus, the parameters in equation (23) can be determined:(27)δ11=1EI∫0lcosθM1¯2dx=lcosθ3EI,(28)δ22=1EI∫0lcosθM2¯2dx=lcosθ3EI,(29)Δ1P=1EI∫M1¯MPdx=124EIql3cos2θ2−cosθ,(30)Δ2P=1EI∫M2¯MPdx=124EIql3cos2θ3cosθ−4,where E represents the elastic modulus and I represents the moment of inertia.Substituting equations (27)–(30) into equation (23) yields the following equation:(31)ML=−112ql2cos2θ,MR=112ql2cosθ6−5cosθ.By equation (19), the total bending moment Mx can be determined as follows:(32)Mx=−12qx2−lcosθx+16l2cos2θ.The maximum momentMM of the cave roof is(33)MM=124ql21−4sin2θ.Assume that the in situ stress in the span direction isσs. Because of the presence of the cave, stress concentration occurs around it. Assume that the stress concentration factor at the centre of the cave roof is λ. Thus, equation (19) can be modified as follows:(34)σti≥6MMbrhr2−λσs.According to the elastic mechanic [30], the factor λ can be determined:(35)λ=1+2hcl,where hc represents the cave height.Substituting equation (33) into equation (34) yields the following equation, which can be used to estimate the ultimate bearing height of the embankment above the underground cave:(36)hf≤41−4sin2θσti+λσsγshr2l2−γrγshr−hs.


## 5. Case Study 1: Laboratory Cave Models

### 5.1. Overview

To assess the stability of underground caves under engineering imposed load, laboratory loading tests of cave models were conducted by Waltham and Swift [1]. The cave models were created by excavating in intact plasters with unconfined compressive strength (σci) of approximately 18 MPa. The cave models were centrally loaded by foundation pads of 1 m2, with no confining pressure, as shown in Figure 11. A summary of model parameters and corresponding values is provided in Table 1. Because excavation causes inevitable disturbances to the plasters [8], the disturbance factor D was assumed to be 0.5.Figure 11
Laboratory loading test (a) and numerical simulation (b).
(a)
(b)Table 1
Basic material properties.
E (GPa)
μ
γ (kN·m−3)
σ
ci (MPa)
m
i
GSI
D
1.2
0.3
23
18
8
100
0.5


### 5.2. Numerical Solutions

A total of 25 numerical models were created to simulate the laboratory loading tests of the plaster cave models. The caves were numerically simulated using FLAC3D, as shown in Figure 11. The variables covered within the numerical simulation included cave roof thickness (hr) and cave width (l). The Hoek–Brown model was employed in the numerical simulation. The material properties required for the modelling are listed in Table 1.Two strength reduction schemes were adopted to determine the relationship between load and FOS, as shown in Figure12. Scheme 1 reduces σci and GSI by FOS, until the maximum unbalance force exceeds 1 × 10−5. Scheme 2 reduces σci and expGSI/K by FOS, until the maximum unbalance force exceeds 1 × 10−5. Both schemes belong to method 4 mentioned in the section “strength reduction.”Figure 12
The relationship between load and FOS. (a) Scheme 1 and (b) scheme 2.
(a)
(b)In free field conditions (load = 0 MPa), the plaster cave models remain stable with a large FOS. As the load increases, the value of FOS gradually decreases. The load that causes FOS decreases to 1 is determined as the ultimate bearing pressure (UBP). The ultimate bearing pressures determined by these two schemes are quite close. However, scheme 2 provides a smoother load FOS curve and more sensitive FOS data.A contour map is used to describe the relationships among the ultimate bearing pressure, the cave roof thickness, and the cave width, as shown in Figure13. In the figure, the ultimate bearing pressure determined by the laboratory loading test is represented by a red dotted line and the ultimate bearing pressure determined by the numerical simulation is represented by a solid black line. The results of the laboratory loading test and the numerical simulation are found to be generally consistent. However, because of the variability of the material used in the test, the two results are completely overlapped over only small sections in the plot.Figure 13
Contour map of the ultimate bearing pressure.


### 5.3. Analytical Solutions

Two schemes were adopted to determine the analytical solution of ultimate bearing pressure.Scheme 1.
The analytical solution is based on the simply supported beam assumption.
Equation (21) can be transformed as follows:(37)UBP≤43hr2l2σti−γrhr,where UBP=γrhr+γshf+hs.Scheme 2.
The modified analytical solution accounts for the fixed support, in situ stress, and roof inclination.
The laboratory loading tests of the cave models were conducted with no confining pressure. Therefore, it is safe to assume thatσs = 0 kPa. None of the cave models had any inclined roofs. Therefore, it is reasonable to assume that θ = 0°. Substituting σs = 0 kPa and θ = 0° into Equation (36) results in the following:(38)UBP≤4hr2l2σti−γrhr.
In both schemes, the tensile strength of the intact rocks is estimated asσti=mi−1⋅σci.
In Figure14, the ultimate bearing pressure determined by the laboratory loading test is represented by a red dotted line and the ultimate bearing pressure determined by the analytical solution is represented by a solid black line. Significant gaps are observed between scheme 1 and the laboratory loading test. Scheme 2 has significantly narrowed this gap and is therefore recommended in this paper.
Data from all sources are summarized in Figure15. It can be observed again that the laboratory loading tests and the numerical simulations are generally consistent with each other. Significant gaps exist between the analytical solution (scheme 1) based on the simply supported beam assumption and the laboratory results. The modified analytical solution (scheme 2) performs fairly well in narrowing this gap.
The applied load in the analytical solutions is slightly different. In the laboratory loading tests and the numerical simulations, all models were centrally loaded by foundation pads of 1 m2 above the cave roof. In both the analytical solutions, the load was applied to the entire cave roof; thus, the stability was reduced compared to the numerical solutions. Therefore, the analytical solution could not exceed the laboratory and numerical results.
As shown in Figure15, the main factor in determining the stability is the thickness-to-span ratio (i.e., the ratio of roof thickness and cave span). A small trick was used here to exploit the relationship between ultimate bearing pressure and the thickness-to-span ratio (hr/l), as shown in Figure 16. In the new plot, the horizontal axis was converted to the square of the thickness-to-span ratio, i.e., (hr/l)2. The cave stability is found to generally increase linearly with the square of the thickness-to-span ratio.Figure 14
Contour map of the ultimate bearing pressure. (a) Scheme 1 and (b) scheme 2.
(a)
(b)Figure 15
The relationship between the ultimate bearing pressure and the thickness-to-span ratio.Figure 16
The relationship between the ultimate bearing pressure and the square of the thickness-to-span ratio.


## 5.1. Overview

To assess the stability of underground caves under engineering imposed load, laboratory loading tests of cave models were conducted by Waltham and Swift [1]. The cave models were created by excavating in intact plasters with unconfined compressive strength (σci) of approximately 18 MPa. The cave models were centrally loaded by foundation pads of 1 m2, with no confining pressure, as shown in Figure 11. A summary of model parameters and corresponding values is provided in Table 1. Because excavation causes inevitable disturbances to the plasters [8], the disturbance factor D was assumed to be 0.5.Figure 11
Laboratory loading test (a) and numerical simulation (b).
(a)
(b)Table 1
Basic material properties.
E (GPa)
μ
γ (kN·m−3)
σ
ci (MPa)
m
i
GSI
D
1.2
0.3
23
18
8
100
0.5


## 5.2. Numerical Solutions

A total of 25 numerical models were created to simulate the laboratory loading tests of the plaster cave models. The caves were numerically simulated using FLAC3D, as shown in Figure 11. The variables covered within the numerical simulation included cave roof thickness (hr) and cave width (l). The Hoek–Brown model was employed in the numerical simulation. The material properties required for the modelling are listed in Table 1.Two strength reduction schemes were adopted to determine the relationship between load and FOS, as shown in Figure12. Scheme 1 reduces σci and GSI by FOS, until the maximum unbalance force exceeds 1 × 10−5. Scheme 2 reduces σci and expGSI/K by FOS, until the maximum unbalance force exceeds 1 × 10−5. Both schemes belong to method 4 mentioned in the section “strength reduction.”Figure 12
The relationship between load and FOS. (a) Scheme 1 and (b) scheme 2.
(a)
(b)In free field conditions (load = 0 MPa), the plaster cave models remain stable with a large FOS. As the load increases, the value of FOS gradually decreases. The load that causes FOS decreases to 1 is determined as the ultimate bearing pressure (UBP). The ultimate bearing pressures determined by these two schemes are quite close. However, scheme 2 provides a smoother load FOS curve and more sensitive FOS data.A contour map is used to describe the relationships among the ultimate bearing pressure, the cave roof thickness, and the cave width, as shown in Figure13. In the figure, the ultimate bearing pressure determined by the laboratory loading test is represented by a red dotted line and the ultimate bearing pressure determined by the numerical simulation is represented by a solid black line. The results of the laboratory loading test and the numerical simulation are found to be generally consistent. However, because of the variability of the material used in the test, the two results are completely overlapped over only small sections in the plot.Figure 13
Contour map of the ultimate bearing pressure.


## 5.3. Analytical Solutions

Two schemes were adopted to determine the analytical solution of ultimate bearing pressure.Scheme 1.
The analytical solution is based on the simply supported beam assumption.
Equation (21) can be transformed as follows:(37)UBP≤43hr2l2σti−γrhr,where UBP=γrhr+γshf+hs.Scheme 2.
The modified analytical solution accounts for the fixed support, in situ stress, and roof inclination.
The laboratory loading tests of the cave models were conducted with no confining pressure. Therefore, it is safe to assume thatσs = 0 kPa. None of the cave models had any inclined roofs. Therefore, it is reasonable to assume that θ = 0°. Substituting σs = 0 kPa and θ = 0° into Equation (36) results in the following:(38)UBP≤4hr2l2σti−γrhr.
In both schemes, the tensile strength of the intact rocks is estimated asσti=mi−1⋅σci.
In Figure14, the ultimate bearing pressure determined by the laboratory loading test is represented by a red dotted line and the ultimate bearing pressure determined by the analytical solution is represented by a solid black line. Significant gaps are observed between scheme 1 and the laboratory loading test. Scheme 2 has significantly narrowed this gap and is therefore recommended in this paper.
Data from all sources are summarized in Figure15. It can be observed again that the laboratory loading tests and the numerical simulations are generally consistent with each other. Significant gaps exist between the analytical solution (scheme 1) based on the simply supported beam assumption and the laboratory results. The modified analytical solution (scheme 2) performs fairly well in narrowing this gap.
The applied load in the analytical solutions is slightly different. In the laboratory loading tests and the numerical simulations, all models were centrally loaded by foundation pads of 1 m2 above the cave roof. In both the analytical solutions, the load was applied to the entire cave roof; thus, the stability was reduced compared to the numerical solutions. Therefore, the analytical solution could not exceed the laboratory and numerical results.
As shown in Figure15, the main factor in determining the stability is the thickness-to-span ratio (i.e., the ratio of roof thickness and cave span). A small trick was used here to exploit the relationship between ultimate bearing pressure and the thickness-to-span ratio (hr/l), as shown in Figure 16. In the new plot, the horizontal axis was converted to the square of the thickness-to-span ratio, i.e., (hr/l)2. The cave stability is found to generally increase linearly with the square of the thickness-to-span ratio.Figure 14
Contour map of the ultimate bearing pressure. (a) Scheme 1 and (b) scheme 2.
(a)
(b)Figure 15
The relationship between the ultimate bearing pressure and the thickness-to-span ratio.Figure 16
The relationship between the ultimate bearing pressure and the square of the thickness-to-span ratio.


## 6. Case Study 2: The Underground Caves beneath the Changli Highway

### 6.1. Overview

Changli highway, a main road linking Nanchang and Shangli, was undergoing construction in a karstified area in the year 2015. Dozens of underground karst caves were surveyed before constructing the embankment. The span of the caves is generally 10–40 m long; the rock over the caves is generally 1–4 m thick. The thickness-to-span ratio of nearly half of the caves is less than 0.05, and the greatest ratio is only 0.35, as shown in Figure17. Local specifications for the design of highway embankments (JTG D30-2015) suggest that engineering actions must be taken in an underground cave when its thickness-to-span ratio is less than 0.8. To prevent embankment collapse caused by potential cave instability, 0.5 m thick continuous reinforced concrete slabs were built above the embankment. As a passive engineering action, although the stability of the underground caves cannot be improved, this measure could prevent the sudden collapse of the highway.Figure 17
Frequency distribution of the thickness-to-span ratio.Much concern has been paid to the cave locating at K178 + 800, which became unstable during embankment construction. Its span is approximately 27 m long; its roof is approximately 2 m thick. The unstable cave has caused significant settlement of the embankment (Figure18).Figure 18
Settlement of the embankment induced by underground cave instability.


### 6.2. Numerical Modelling

The construction process of the embankment at K178 + 800 was modelled numerically, as shown in Figure19. In the model, the height of the cave is 5 m, the thickness of the roof is 2 m, the inclination of the roof is 15°, the average thickness of the soil cover is 5 m, and the average filling height of the embankment is 6 m. The model bottom is fixed in the vertical and horizontal directions, and the left and right boundaries are fixed in the horizontal direction.Figure 19
Numerical simulation of the embankment construction process.Based on drilling data and geophysical data, there is no controlled structural plane in the rock strata. Very few cracks are found in the drilled rock. According to electrical resistivity tomography data, the resistivity contour of the rock strata is smooth, with no breaks or drastic changes. From the regional engineering geology survey in the design data, the tectonic movement in this region is not obvious, and neither faults nor folds have been found. Therefore, rock and soil stratum was simulated as the continuous material.According to drilling surveys and laboratory measurements, the underground cave is formed in the Carboniferous limestone with Young’s modulus (E) of 43,000 MPa, Poisson’s ratio (v) of 0.26, density (ρ) of 2700 kg·m−3, intact unconfined compressive strength (σci) of approximately 130 MPa, geological strength index (GSI) of 70, and intact rock material constant (mi) of 8. Because the embankment construction causes fewer disturbances to rock than blasting and tunnelling [8], the disturbance factor D was assumed as 0.3. A summary of the model parameters and corresponding values of embankement and cover is provided in Table 2. Both the rock and soil were assumed to be dry in this research study.Table 2
Basic material properties of embankment and cover.
Strata
E (MPa)
v
ρ (kg·m−3)
φ (°)
c (kPa)
Embankment
50
0.35
1900
24
10
Cover
30
0.35
1900
20
50Drilling data indicate that the underground cave contains a small amount of sludge-like filling matters. Filling matters may be beneficial to changing the stress state and improving the cave stability; however, the influence is negligible. Therefore, filling matters were not considered in the numerical simulation.The in situ stress at the site [31] could be approximately expressed as follows:(39)σH=0.03030H+2.2242,σh=0.01718H+0.8866,σv=0.02646H+0.3924,where σH is the maximum horizontal component of the in situ stress, σh is the minimum horizontal component, σv is the minimum horizontal component, and H is the depth. The direction of the maximum principal stress lies between SE-NW and ESE-WNW (approximately perpendicular to the axial direction of the karst cave).First, the natural equilibrium state before construction was simulated. A free field model, including rock strata, soil strata, and the cave, was created. Under the load of self-weight, the units in the model gradually reached equilibrium. This result shows that the underground cave would be stable under the load of self-weight alone. Next, the stratified construction process of embankment was simulated. Each time, a 1 m thick embankment layer was constructed.


### 6.3. Numerical Reproduction of the Failure Process

Before construction, a multiple point extensometer was embedded above the cave to monitor the vertical displacement of the cave roof, as shown in Figure20. The embankment construction process is shown in Table 3. After a large deformation occurred on the 34th day, observations were stopped and grouting was conducted. After grouting, the embankment and the underlying cave were far more stable.Figure 20
Numerical reproduction of the settlement process.Table 3
Construction process.
Time (d)
6
11
16
22
27
32
Embankment height (m)
1
2
3
4
5
6Two modelling schemes were adopted to simulate the settlement process, as shown in Figure20. Scheme 1 adopted the geometric model shown in Figure 19. Scheme 2 transformed the cave into an equal area rectangular cave. In the transformation, both the cave span and the cave centre were kept unchanged as well. Figure 20 shows the simulation results obtained by using the two schemes, with a comparison to the monitored data. It is observed that the simulation results by the two schemes are close and are generally consistent with the monitored data. Therefore, in the case that the cave span and roof thickness remains the same, the influence of the cave shape is not significant.The evolution of the displacement vector field in embankment construction process is shown in Figure21. The simulation results revealed that the displacement is not great (less than 5 mm) when constructing the first three embankment layers. Therefore, this section mainly discusses the displacement vector evolutionary process when constructing the 4th, 5th, and 6th embankment layers, as shown in Figure 21. In the settlement process, an inclined shear zone gradually forms in the cave roof. The cave roof slides downward along the shear zone, raising potential risks for the traction and sliding of the overlying soil. After gradual loss of the bedrock support, the overlaying soil strata begin to settle downward. Next, differential settlement of the soil strata triggers the formation of vertical shear bands. Following the settlement of the soil, the overlaying embankment begins to settle downward. Subsequently, a horizontal stratiform ripped band emerges in the embankment.Figure 21
The evolution of the displacement vector field during the embankment construction process.
(a)
(b)
(c)Both soil and embankment begin to settle downward after losing support. However, different rupture types are observed. In the soil strata, vertical shear bands are observed. In the embankment, horizontal stratiform ripped bands are observed. This observation is primarily the result of the following: first, the soil strata still bear the weight of the embankment when it settles downward; thus, no tensile stress is generated, and therefore, no ripped bands form. Then, the embankment keeps exerting pressure on the soil strata when it settles downward, while significant tensile stress emerges inside it; finally, the horizontal stratiform ripped band begins to generate and subsequently develops to the surface. This failure pattern is consistent with the surface cracks observed at the site in Figure22.Figure 22
Failure pattern observed at the site.The model reveals the failure mechanism of an underground cave under the load of embankment construction. It is therefore necessary to constrain the embankment height to prevent the failure of the underground caves. Once settlement or other recognizable precursors occurred, grouting or other engineering action should be adopted to prevent further failures.


### 6.4. Numerical Solutions

Two strength reduction schemes were adopted to determine the relationship between embankment height (Hf) and FOS, as shown in Figure 23. Scheme 1 reduces σci and GSI by FOS until the maximum unbalance force exceeds 1 × 10−5. Scheme 2 reduces σci and expGSI/K by FOS until the maximum unbalance force exceeds 1 × 10−5.Figure 23
Numerical solutions.In the free field condition (Hf = 0 m), the value of FOS is only slightly greater than 1. The value of FOS gradually decreases with the increasing embankment height (Hf). The embankment height that causes FOS decreases to 1 is determined as the ultimate bearing height (UBH). The estimation UBH = 6 m appears to be consistent with the case observed at the site.


### 6.5. Analytical Solutions

Three schemes were adopted to determine the analytical solution of the ultimate bearing height. Scheme 1 was calculated using equation (21); scheme 2 was calculated using equation (36), with θ = 15°; scheme 3 was calculated using equation (36), with θ = 0°. According to the in situ stress given by equation (39), it is reasonable to assume that σs = 2500 kPa. Other model parameters have been provided in Table 2 and the paragraphs above it.Solutions of UBH are summarized in Table4. The tensile strength of rock mass is determined by the estimation based on the Hoek–Brown criterion and the estimation proposed by Tokashiki, respectively. The UBHs obtained by the estimation based on the Hoek–Brown criterion and scheme 1 are found to be negative; i.e., the underground cave was unstable before the embankment construction. However, in practice, the settlement precursor was only observed after the embankment construction was completed. The UBH obtained by scheme 2 is very close to but slightly smaller than the numerical solution. Scheme 3 ignores the dip angle of the cave roof; thus, the results are somewhat conservative.Table 4
Solutions of UBH.
Tensile strength of rock mass
Analytical solution
Numerical solution
Scheme 1
Scheme 2
Scheme 3
Estimation based on the Hoek–Brown criterion
NEG
NEG
NEG
6.0
Estimation proposed by Tokashiki
NEG
5.5
2.0The bending moments determined by all solutions are represented in Figure24. Significant gaps are observed between the analytical solution based on the simply supported beam assumption (scheme 1) and the numerical solution. The analytical solution based on the fixed beam assumption (scheme 2 and scheme 3) performs fairly well in estimating the bending moments. Therefore, the cave stability is significantly influenced by the fixed supports. In addition, the bending moments at both ends of the top plate determined by scheme 2 are closer to the numerical solution than that determined by scheme 3. Therefore, cave stability is slightly influenced by the dip angle of the cave roof.Figure 24
The bending moment of the cave roof.


## 6.1. Overview

Changli highway, a main road linking Nanchang and Shangli, was undergoing construction in a karstified area in the year 2015. Dozens of underground karst caves were surveyed before constructing the embankment. The span of the caves is generally 10–40 m long; the rock over the caves is generally 1–4 m thick. The thickness-to-span ratio of nearly half of the caves is less than 0.05, and the greatest ratio is only 0.35, as shown in Figure17. Local specifications for the design of highway embankments (JTG D30-2015) suggest that engineering actions must be taken in an underground cave when its thickness-to-span ratio is less than 0.8. To prevent embankment collapse caused by potential cave instability, 0.5 m thick continuous reinforced concrete slabs were built above the embankment. As a passive engineering action, although the stability of the underground caves cannot be improved, this measure could prevent the sudden collapse of the highway.Figure 17
Frequency distribution of the thickness-to-span ratio.Much concern has been paid to the cave locating at K178 + 800, which became unstable during embankment construction. Its span is approximately 27 m long; its roof is approximately 2 m thick. The unstable cave has caused significant settlement of the embankment (Figure18).Figure 18
Settlement of the embankment induced by underground cave instability.


## 6.2. Numerical Modelling

The construction process of the embankment at K178 + 800 was modelled numerically, as shown in Figure19. In the model, the height of the cave is 5 m, the thickness of the roof is 2 m, the inclination of the roof is 15°, the average thickness of the soil cover is 5 m, and the average filling height of the embankment is 6 m. The model bottom is fixed in the vertical and horizontal directions, and the left and right boundaries are fixed in the horizontal direction.Figure 19
Numerical simulation of the embankment construction process.Based on drilling data and geophysical data, there is no controlled structural plane in the rock strata. Very few cracks are found in the drilled rock. According to electrical resistivity tomography data, the resistivity contour of the rock strata is smooth, with no breaks or drastic changes. From the regional engineering geology survey in the design data, the tectonic movement in this region is not obvious, and neither faults nor folds have been found. Therefore, rock and soil stratum was simulated as the continuous material.According to drilling surveys and laboratory measurements, the underground cave is formed in the Carboniferous limestone with Young’s modulus (E) of 43,000 MPa, Poisson’s ratio (v) of 0.26, density (ρ) of 2700 kg·m−3, intact unconfined compressive strength (σci) of approximately 130 MPa, geological strength index (GSI) of 70, and intact rock material constant (mi) of 8. Because the embankment construction causes fewer disturbances to rock than blasting and tunnelling [8], the disturbance factor D was assumed as 0.3. A summary of the model parameters and corresponding values of embankement and cover is provided in Table 2. Both the rock and soil were assumed to be dry in this research study.Table 2
Basic material properties of embankment and cover.
Strata
E (MPa)
v
ρ (kg·m−3)
φ (°)
c (kPa)
Embankment
50
0.35
1900
24
10
Cover
30
0.35
1900
20
50Drilling data indicate that the underground cave contains a small amount of sludge-like filling matters. Filling matters may be beneficial to changing the stress state and improving the cave stability; however, the influence is negligible. Therefore, filling matters were not considered in the numerical simulation.The in situ stress at the site [31] could be approximately expressed as follows:(39)σH=0.03030H+2.2242,σh=0.01718H+0.8866,σv=0.02646H+0.3924,where σH is the maximum horizontal component of the in situ stress, σh is the minimum horizontal component, σv is the minimum horizontal component, and H is the depth. The direction of the maximum principal stress lies between SE-NW and ESE-WNW (approximately perpendicular to the axial direction of the karst cave).First, the natural equilibrium state before construction was simulated. A free field model, including rock strata, soil strata, and the cave, was created. Under the load of self-weight, the units in the model gradually reached equilibrium. This result shows that the underground cave would be stable under the load of self-weight alone. Next, the stratified construction process of embankment was simulated. Each time, a 1 m thick embankment layer was constructed.


## 6.3. Numerical Reproduction of the Failure Process

Before construction, a multiple point extensometer was embedded above the cave to monitor the vertical displacement of the cave roof, as shown in Figure20. The embankment construction process is shown in Table 3. After a large deformation occurred on the 34th day, observations were stopped and grouting was conducted. After grouting, the embankment and the underlying cave were far more stable.Figure 20
Numerical reproduction of the settlement process.Table 3
Construction process.
Time (d)
6
11
16
22
27
32
Embankment height (m)
1
2
3
4
5
6Two modelling schemes were adopted to simulate the settlement process, as shown in Figure20. Scheme 1 adopted the geometric model shown in Figure 19. Scheme 2 transformed the cave into an equal area rectangular cave. In the transformation, both the cave span and the cave centre were kept unchanged as well. Figure 20 shows the simulation results obtained by using the two schemes, with a comparison to the monitored data. It is observed that the simulation results by the two schemes are close and are generally consistent with the monitored data. Therefore, in the case that the cave span and roof thickness remains the same, the influence of the cave shape is not significant.The evolution of the displacement vector field in embankment construction process is shown in Figure21. The simulation results revealed that the displacement is not great (less than 5 mm) when constructing the first three embankment layers. Therefore, this section mainly discusses the displacement vector evolutionary process when constructing the 4th, 5th, and 6th embankment layers, as shown in Figure 21. In the settlement process, an inclined shear zone gradually forms in the cave roof. The cave roof slides downward along the shear zone, raising potential risks for the traction and sliding of the overlying soil. After gradual loss of the bedrock support, the overlaying soil strata begin to settle downward. Next, differential settlement of the soil strata triggers the formation of vertical shear bands. Following the settlement of the soil, the overlaying embankment begins to settle downward. Subsequently, a horizontal stratiform ripped band emerges in the embankment.Figure 21
The evolution of the displacement vector field during the embankment construction process.
(a)
(b)
(c)Both soil and embankment begin to settle downward after losing support. However, different rupture types are observed. In the soil strata, vertical shear bands are observed. In the embankment, horizontal stratiform ripped bands are observed. This observation is primarily the result of the following: first, the soil strata still bear the weight of the embankment when it settles downward; thus, no tensile stress is generated, and therefore, no ripped bands form. Then, the embankment keeps exerting pressure on the soil strata when it settles downward, while significant tensile stress emerges inside it; finally, the horizontal stratiform ripped band begins to generate and subsequently develops to the surface. This failure pattern is consistent with the surface cracks observed at the site in Figure22.Figure 22
Failure pattern observed at the site.The model reveals the failure mechanism of an underground cave under the load of embankment construction. It is therefore necessary to constrain the embankment height to prevent the failure of the underground caves. Once settlement or other recognizable precursors occurred, grouting or other engineering action should be adopted to prevent further failures.


## 6.4. Numerical Solutions

Two strength reduction schemes were adopted to determine the relationship between embankment height (Hf) and FOS, as shown in Figure 23. Scheme 1 reduces σci and GSI by FOS until the maximum unbalance force exceeds 1 × 10−5. Scheme 2 reduces σci and expGSI/K by FOS until the maximum unbalance force exceeds 1 × 10−5.Figure 23
Numerical solutions.In the free field condition (Hf = 0 m), the value of FOS is only slightly greater than 1. The value of FOS gradually decreases with the increasing embankment height (Hf). The embankment height that causes FOS decreases to 1 is determined as the ultimate bearing height (UBH). The estimation UBH = 6 m appears to be consistent with the case observed at the site.


## 6.5. Analytical Solutions

Three schemes were adopted to determine the analytical solution of the ultimate bearing height. Scheme 1 was calculated using equation (21); scheme 2 was calculated using equation (36), with θ = 15°; scheme 3 was calculated using equation (36), with θ = 0°. According to the in situ stress given by equation (39), it is reasonable to assume that σs = 2500 kPa. Other model parameters have been provided in Table 2 and the paragraphs above it.Solutions of UBH are summarized in Table4. The tensile strength of rock mass is determined by the estimation based on the Hoek–Brown criterion and the estimation proposed by Tokashiki, respectively. The UBHs obtained by the estimation based on the Hoek–Brown criterion and scheme 1 are found to be negative; i.e., the underground cave was unstable before the embankment construction. However, in practice, the settlement precursor was only observed after the embankment construction was completed. The UBH obtained by scheme 2 is very close to but slightly smaller than the numerical solution. Scheme 3 ignores the dip angle of the cave roof; thus, the results are somewhat conservative.Table 4
Solutions of UBH.
Tensile strength of rock mass
Analytical solution
Numerical solution
Scheme 1
Scheme 2
Scheme 3
Estimation based on the Hoek–Brown criterion
NEG
NEG
NEG
6.0
Estimation proposed by Tokashiki
NEG
5.5
2.0The bending moments determined by all solutions are represented in Figure24. Significant gaps are observed between the analytical solution based on the simply supported beam assumption (scheme 1) and the numerical solution. The analytical solution based on the fixed beam assumption (scheme 2 and scheme 3) performs fairly well in estimating the bending moments. Therefore, the cave stability is significantly influenced by the fixed supports. In addition, the bending moments at both ends of the top plate determined by scheme 2 are closer to the numerical solution than that determined by scheme 3. Therefore, cave stability is slightly influenced by the dip angle of the cave roof.Figure 24
The bending moment of the cave roof.


## 7. Thickness-to-Span Ratio

It is necessary to determine the thickness-to-span ratio before constructing an embankment above an underground cave. The current specifications for the design of highway embankments (JTG D30-2015) in China suggest that engineering actions must be taken to treat the underground cave when its thickness-to-span ratio (hr/l) is less than 0.8.This specification might be conservative. Many embankments along the Changli highway still stand on an underground cave with a thickness-to-span ratio less than 0.8. In typical limestone karst, assuming that the rock mass surrounding the underground cave is of general quality (GSI = 60), with parameters summarized in Table5, two cases are considered. Case 1 assumes that the in situ stress is σs = 1000 kPa. For comparison, case 2 assumes that σs = 0 kPa. The validated ultimate bearing height analytical solution (i.e., equation (36)) is used to calculate the ultimate fill height. In case 1, an underground cave with a thickness-to-span ratio = 0.8 can bear an embankment of over 200 m in height. Even in case 2, it can bear an embankment of over 50 m in height. However, the local embankment height typically does not exceed 20 m.Table 5
Typical cases and the parameters.
Case
GSI
D
m
i
σ
ci (kPa)
γ
s (kN·m−3)
γ
r (kN·m−3)
h
s (m)
h
r (m)
θ (°)
λ
σ
s (kPa)
1
60
0.3
10
100000
20
27
5
4
0
1.5
1000
2
60
0.3
10
100000
20
27
5
4
0
1.5
0Waltham and Swift [1] suggested that the underground cave with a thickness-to-span ratio over 0.7 could bear a conventional foundation load of less than 2 MPa. The only drawback to Waltham’s study is that it was performed only under gravity stress. However, in shallow stratum, in situ stress might be several times the stress of gravity.The recommended method to determine the required thickness-to-span ratio in this paper is to transform the ultimate bearing height analytical solution (i.e., equation (36)) into(40)hrl≥14−sin2θγshf+hs+γrhrσti+λσs,where the unit of hr, hf, hs, and l is m; the unit of γs and γr is kN·m−3; the unit of σti and σs is kPa; the unit of θ is °. λ is dimensionless.In the case thatθ = 0°, a more simplified form could be obtained as follows:(41)hrl≥12γshf+hs+γrhrσti+λσs.Assuming that a 20 m high embankment is placed on the highway, the required thickness-span ratio is calculated using the parameters in Table5. In case 1, the analytical solution requires that the thickness-to-span ratio should be greater than 0.28. The thickness-to-span ratio = 0.28 is significantly lower than the thickness-to-span ratio = 0.8 required by specifications for design of highway embankments (JTG D30-2015). In case 2, the analytical solution requires that the thickness-to-span ratio should be greater than 0.63. The thickness-to-span ratio = 0.63 is quite close to the thickness-to-span ratio = 0.7 given by Waltham and Swift [1].


## 8. Conclusions

Significant gaps are observed between the analytical solution based on the simply supported beam assumption and the laboratory results. The modified analytical solution performs fairly well in narrowing this gap and is therefore recommended in this paper. The capacity of an underground cave to bear embankment construction is observed to be proportional to the tensile strength of the rock mass and the square of the thickness-to-span ratio of the cave roof. The tensile strength of the rock mass isψ times lower than that of the intact rock. The value of ψ is mainly determined by the geological strength index.To prevent instability of underground caves, the embankment height should be reasonably controlled. However, local construction specifications requiring that the roof thickness be greater than 0.8 times the cave span appears to be conservative. Many civil engineering works do stand on caves with much thinner roofs than recommended by the specifications. The thickness-to-span ratio derived in this paper could provide a more appropriate way in determining the stability of underground caves.



---

*Source: 3607574-2019-02-11.xml*