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

#### 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 *m*_{i}, 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:where *σ*_{1} and *σ*_{3} are the major and minor principal stresses, respectively, and compression is taken as positive. The strength parameters *m*_{b}, *s*, and *a* are functions of the geological strength index GSI:where *m*_{i} 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, equivalent *c*-*φ* 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:where and is the upper limit of the confining stress. The value of can be determined by searching the maximum *σ*_{3} in the case.

Subsequently, the equivalent *c*-*φ* 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.

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

*Method 4 (reduction of Hoek–Brown parameters). *This method selects one or several Hoek–Brown parameters (e.g., GSI, *σ*_{ci}, and *m*_{i}) 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 parameters and , 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 and are obtained as follows:Substitute into equation (2) to obtain the reduced parameters , , and :

*Scheme 2 (reduction of σ _{ci} and ). *The parameters

*σ*

_{ci}and are reduced by FOS, and the reduced parameters and are obtained as follows:where .

Substitute into equation (2) to obtain the reduced parameters , , and :The relationships of instantaneous , , , and 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, , , and 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 becomes approximately linear using scheme 2.

**(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/m^{3}, *σ*_{ci} = 40 MPa, *m*_{i} = 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.

**(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 Figure 5(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/m^{3}, *σ*_{ci} = 18 MPa, *m*_{i} = 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.

**(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 *m*_{i}, 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:

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:

Because is generally close to 0, the use of an equivalent infinitesimal substitute can simplify the equation and solve the problem easily:

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., . Based on the Hoek–Brown 2002 criterion, it can be deduced that . 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 , 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 *m*_{b} = *m*_{i}. In this case, equation (11) can be approximated as follows:

The *m*_{i} 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:where .

As seen from Figure 6, 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 *η* ≈ 2*D* + 7.54 can greatly simplify the original equation . Using this linear approximation, the following equation can be used to estimate the tensile strength of rock mass:where .

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.

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 *ψ*: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:

Therefore, the value of *ψ* can be estimated as follows:

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

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

Assuming that cave span is *l* and uniform load is *q*, the maximum moment *M*_{M} of the cave roof is

To prevent failure, the tensile strength *σ*_{ti} should be greater than the tensile stress:where *h*_{r} represents the thickness of the cave roof and *b*_{r} represents the width of the cave roof.

The load applied to the cave roof can be expressed as follows:where *h*_{s} represents the thickness of the soil stratum over the cave roof; *h*_{f} 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:

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.

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:where *M*_{P} is the bending moment generated by the force applied on the cave roof and is calculated by . The parameter represents the vertical shear force. The parameter *x* represents the position coordinate. The parameter represents the tilt angle of the cave roof. *M*_{L} and *M*_{R} are the bending moments at the left and right supports of the cave roof, respectively. Using the force method, the following equations are obtained: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 *M*_{L} = 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 *M*_{R} = 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 .

The bending moment of the cave roof generated by the action of the unit moment *M*_{L} = 1 can be determined as follows:

The bending moment of the cave roof generated by the action of the unit moment *M*_{R} = 1 can be determined as follows:

The bending moment generated by the action of force *q* can be determined as follows:

Thus, the parameters in equation (23) can be determined:where *E* represents the elastic modulus and *I* represents the moment of inertia.

Substituting equations (27)–(30) into equation (23) yields the following equation:

By equation (19), the total bending moment can be determined as follows:

The maximum moment *M*_{M} of the cave roof is

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:

According to the elastic mechanic [30], the factor *λ* can be determined:where *h*_{c} 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:

#### 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 m^{2}, 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.

**(a)**

**(b)**

##### 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 FLAC^{3D}, as shown in Figure 11. The variables covered within the numerical simulation included cave roof thickness (*h*_{r}) 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 Figure 12. Scheme 1 reduces *σ*_{ci} and GSI by FOS, until the maximum unbalance force exceeds 1 × 10^{−5}. Scheme 2 reduces *σ*_{ci} and by FOS, until the maximum unbalance force exceeds 1 × 10^{−5}. Both schemes belong to method 4 mentioned in the section “strength reduction.”

**(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 Figure 13. 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.

##### 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:where .

*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:In both schemes, the tensile strength of the intact rocks is estimated as .

In Figure 14, 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 Figure 15. 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 m^{2} 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 Figure 15, 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 (*h*_{r}/*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., (*h*_{r}/*l*)^{2}. The cave stability is found to generally increase linearly with the square of the thickness-to-span ratio.

**(a)**

**(b)**

#### 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 Figure 17. 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.

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 (Figure 18).

##### 6.2. Numerical Modelling

The construction process of the embankment at K178 + 800 was modelled numerically, as shown in Figure 19. 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.

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 () 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 (*m*_{i}) 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.

Drilling 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: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 Figure 20. 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.

Two modelling schemes were adopted to simulate the settlement process, as shown in Figure 20. 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 Figure 21. 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.

**(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 Figure 22.

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 (*H*_{f}) 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 by FOS until the maximum unbalance force exceeds 1 × 10^{−5}.

In the free field condition (*H*_{f} = 0 m), the value of FOS is only slightly greater than 1. The value of FOS gradually decreases with the increasing embankment height (*H*_{f}). 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 Table 4. 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.

The bending moments determined by all solutions are represented in Figure 24. 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.

#### 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 (*h*_{r}/*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 Table 5, 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.

Waltham 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)) intowhere the unit of *h*_{r}, *h*_{f}, *h*_{s}, 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:

Assuming that a 20 m high embankment is placed on the highway, the required thickness-span ratio is calculated using the parameters in Table 5. 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.

#### 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 there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

This research was funded by the Technology Program of Jiangxi Transportation Hall (Grant number 2015C0022) and the Fundamental Research Foundation of the Central Universities (Grant number 310821175026).