#### Abstract

The presence of a concealed karst cave above a deep highway tunnel may cause the collapse of the rock mass between the karst cave and tunnel during excavation. Rock mass collapse threatens the safety of tunnel construction personnel. A prediction method of the collapse region induced by a concealed karst cave above a deep highway tunnel is proposed on the basis of the upper bound theorem of limit analysis. An analytical expression of the collapse surface is derived from the variational principle. Using the analytical expression of the collapse surface, the shapes of the collapse surfaces are plotted for different rock mass parameters. Moreover, the minimum safe distance between the karst cave and tunnel is defined, and the computational equation of the minimum safe distance is derived. The proposed method is applied in a highway tunnel excavated in a karst terrain as a case study. Based on geological survey report parameters, the shape of the collapse surface and the minimum safe distance between the karst cave and tunnel are obtained. Finally, the collapse surface of the rock mass provided by the proposed approach is compared with that provided by numerical simulation, and the favorable result comparison shows that the proposed method is valid.

#### 1. Introduction

Karst caves are widely present in karst terrains, so karst caves are common in the adjacent regions of planned tunnels or even in the excavated regions of planned tunnels. Most highway tunnels are excavated by drilling and blasting methods. If the karst cave is too close to the planned tunnel, the rock mass between the karst cave and the tunnel may collapse under the effect of the blasting shock. Moreover, karst caves can be filled, and the karst cave fill material would rush into the tunnel instantaneously, which would threaten the personal safety of construction personnel. Even if there is no fill material or the volume of fill material in the karst cave is small, the collapse of the rock mass between the karst cave and the tunnel may damage the construction equipment, which would cause delays in the construction period and economic losses. Figure 1 shows a roof collapse induced by a karst cave located above a highway tunnel during tunnel excavation. Therefore, if we can establish a method to predict the potential collapse location and range for the rock mass between the karst cave and the tunnel, we can adopt countermeasures to avoid the occurrence of collapse during tunnel construction in karst terrains.

Because of the lack of understanding of karst characteristics, great difficulties are encountered in tunnel construction in karst terrains. Understanding the features of karst is a key factor for tunnel construction in karst terrains, and many investigators have realized the aforementioned, and great attention has been paid to this issue. Taking the Yichang–Wanzhou railway tunnel as the engineering background, Fan et al. [1] investigated the main features of karst geological disasters. By studying the influence of several unfavorable geological disasters in tunnel construction, the authors proposed treatment methods for large karst caves, and the validity of the proposed methods was verified based on monitoring data. Because a karst cave may cause geohazards during tunnel excavation in karst terrains, Li et al. [2] used ground-penetrating radar and geological drilling to predict the geometric characteristics and locations of karst caves. By applying the proposed techniques in an actual project, Li et al. [2] found that the karst cave location predicted by their technique agreed well with field observations. Subsequently, Li et al. [3] proposed a new method to study the mechanism of water and mud inrush during tunnel construction in karst terrains. Using a mud inrush computational model, Li et al. [3] investigated the integral sliding stability of the filling media on the basis of the simplified Bishop method, and their results are consistent with the solutions provided by FLAC3D. In a subsequent study, Li et al. [4] developed a method to predict the top of a concealed karst cave based on displacement monitoring during tunnel construction. To verify the validity of the proposed method, the authors compared the results derived from their method with those provided by numerical simulation. More recently, Li et al. [5] designed a true triaxial geomechanical model test to study the stability of the surrounding rock under the effect of a water-bearing cave in front of the tunnel face. To prevent water and mud inrush during tunnel construction in karst terrains, a computational model is established by Wu et al. [6] to calculate the minimum thickness of the rock stratum between the face and the cave under earthquake action. Considering that the tunnel lining and other structures are designed by a similar theory in actual projects, Tian et al. [7] employed a model experiment to study the combined influence of the surrounding rock parameters on tunnel leakage in a karst area.

Although the stability of a tunnel excavated in a karst area has drawn the attention of many scholars, previous studies were carried out on the basis of the numerical simulation technique and model experiments. However, theoretical studies of the collapse mechanism and range of the surrounding rock induced by a concealed karst cave in tunnel construction are rare. Fraldi and Guarracino [8] developed an analytical approach to determine the range and shape of the possible collapse of a cavity roof induced by excavation. Using the variational principle, the researchers obtained the equation of the collapsing block of the cavity roof on the basis of the limit analysis theorem. Because the collapse mechanism of the surrounding rock proposed by Fraldi and Guarracino [8] was constructed by an arbitrary curve, the final result of the collapsing block derived from optimal calculations was not restricted by a predefined line type. As noted by Huang et al. [9], the advantage of this approach is that the final solution of the shape of the collapsing block is consistent with the collapse observed in actual engineering. Due to the remarkable advantage of this approach, numerous scholars have used this approach to study the collapsing blocks of the roof of the tunnel and rectangular cavity under various complex geological conditions [10–17].

Considering that the variational principle is an effective method to determine the collapse region of the surrounding rock in tunnel construction, investigators used this method to study the collapse region induced by a karst cave in a tunnel excavated in a karst area. If a karst cave is present beneath a deep highway tunnel, there is a risk of collapse of the tunnel floor. In view of this problem, Huang et al. [18] derived an analytical expression of the collapse surface of the rock mass under the tunnel floor by using the variational principle and upper bound theorem. Based on this analytical solution, the range of the collapse region was obtained, and the validation of the analytical solution was proved by comparing it with the numerical simulation solution. Consequently, Yang et al. [19] employed the functional catastrophe theory to investigate the collapse region of the rock mass beneath a tunnel floor in Hoek–Brown (H–B) rock media. Considering that the failure pattern of the surrounding rock is different for different pressures of the filler material in a karst cave, Zhang and Yang [20] constructed a new failure mechanism that reflects active and passive failure patterns of the surrounding rock beneath the tunnel floor. Based on the new failure mechanism, the analytical expressions of the rock mass for active and passive failure patterns are derived, and the influence of different rock parameters on failure is obtained by parameter sensitivity analysis.

Presently, studies of the collapse region of the surrounding rock induced by a concealed karst cave focus on the rock mass located beneath the tunnel floor. However, if a concealed karst cave exists above the tunnel roof, the collapse of the rock mass above the tunnel roof occurs much more easily. In actual engineering, the pressure of the filler material in the karst cave, the gravity of the rock mass, and the blasting shock are three adverse factors that can cause the collapse of the rock pillar between the karst cave and the tunnel roof. Under extreme conditions, the three factors may work simultaneously, thereby enhancing the rock mass collapse. Therefore, it is necessary to develop a theoretical method to study the collapse mechanism of the rock mass induced by a concealed karst cave that is present above the tunnel roof.

This paper aims to investigate the collapse mode of the rock mass between the concealed karst cave and tunnel roof within the framework of the upper bound theorem. A failure mechanism is proposed based on the collapse characteristics of the rock mass induced by the concealed karst cave. Using the failure mechanism and variational principle, the analytical equation of the collapse surface of the rock mass is derived, and the range of the collapse region is obtained. To apply this method in an actual project, the Wuzhishan highway tunnel, excavated in a karst area, is selected as a case study. Based on a geological survey report and the position of the karst cave encountered during tunnel excavation, the failure mechanism in this actual engineering project is established. An analytical equation of the collapse surface of the rock mass above the Wuzhishan tunnel is derived. Furthermore, to validate the presented method, the obtained result is compared with that derived from the numerical simulation technique.

#### 2. Collapse Mechanism of a Highway Tunnel Roof Induced by a Concealed Karst Cave

A schematic diagram of the collapse mode of the rock mass under consideration is illustrated in Figure 2. A highway tunnel is excavated in a karst area, and there is a concealed karst cave just above the tunnel. Under the actions of the filler material pressure in the karst cave and construction disturbances, collapse of the rock mass between the karst cave and the tunnel roof occurs. The rock mass in a certain range would move downwards with velocity relative to the rest of the surrounding rock. Thus, the rock mass in the collapse range is called the collapsing block, and the boundary of the collapse range can be regarded as the collapse surface. Due to the movement of the collapsing block, a relative sliding movement occurs along the collapse surface, and a plastic flow region is formed in this narrow transition layer, where *e*_{i} and *n*_{i} are the tangential and normal directions of a random point along the collapse surface direction, respectively.

Because the collapse range is unknown, the collapse surface is assumed to consist of two symmetrical arbitrary curves , thus implying that the curves are not of any predefined line type, which is consistent with the rock mass collapse mode observed in actual projects. Moreover, the collapse surface extends from the tunnel roof to the bottom of the karst cave. As the collapse surface consists of arbitrary curves, the karst cave and tunnel roof locations where the rock mass begins to collapse are random. The final collapse surface of the rock mass and the location where collapse occurs are determined by theoretical calculations.

The karst cave is assumed to be a circle for mathematical convenience and the equation of the circle is . denotes the distance between the bottom of the karst cave and the roof of the tunnel and and are the half-widths of the bottom and top of the collapsing block, respectively. Furthermore, a standard internal contour for a two-lane highway tunnel is shown in Figure 3. It can be seen from the figure that the internal contour of the two-lane highway tunnel is a multicircular cross section, which consists of three arc sections. The longest arc is located at the top of the tunnel, and the same two arcs are located at the sidewalls of the tunnel. , , and are the centers of the three arcs, while and are the radii of the three arcs. Because the rock mass collapse occurs at the roof of the tunnel, the arc at the top of the tunnel would affect the shape and range of the collapsing block. Therefore, this arc is of most concern, and the equation of this arc is assumed to be .

#### 3. Calculation of the Collapse Surface of a Highway Tunnel Roof Induced by a Concealed Karst Cave

Because of the relative sliding movement between the collapsing block and the rest of the surrounding rock, energy dissipation would occur along the collapse surface. The upper bound theorem states that the upper bound solution of the equation of the collapse surface can be derived from the relation between the rate of total energy dissipation and the rate of external work. Therefore, the energy dissipation of a random point on the collapse surface should be calculated first. According to Chen [21], the rate of energy dissipation of a random point in the plastic flow region can be determined by the plastic stress/strain rate relation. Moreover, the stress/strain rate relation can be deduced from the normality condition related to the yield function. Because most rock masses in deep strata are intact rock or jointed rock masses, it is necessary to employ a suitable failure criterion to estimate the strength and deformability of the rock masses. However, the H–B failure criterion is a criterion developed to evaluate the rock mass strength for both intact rock and jointed rock masses, and the validity of this criterion has been proven by numerous scholars. Therefore, the H–B failure criterion is used here to evaluate the rock mass strength in the influence region affected by excavation. If the rock mass obeys the H–B failure criterion and its associated flow rule, the plastic potential , which is employed to calculate the stress and plastic strain rate relation, can be written in the following form:where *τ* is the shear stress, is the normal stress, and are material constants, is the uniaxial compressive strength, and is the tensile strength of the rock mass.

According to Fraldi and Guarracino [8], the rate of energy dissipation of a random point on the collapse surface is computed bywhere is the first derivative of . It can be observed in Figure 1 that the failure mechanism is symmetrical about the *y*-axis, and the energy dissipation and external work rates can be calculated based on half of the failure mechanism for mathematical convenience. The half-rate of energy dissipation along the collapse surface can be obtained by integrating (2) over interval [,], which is expressed in the following form:

The external forces in this failure mechanism include the gravity of the rock mass, supporting force, and filler material pressure. The total rate of the work in this mechanism is the sum of the rates of work performed by these external forces. The rate of work produced by the self-weight of the half-collapsing block is given bywhere is the unit weight of the rock mass. Moreover, the expressions of and arewhere is the radius of the karst cave. Assuming that the supporting pressure produced by the lining is distributed along the tunnel internal contour in the radial direction, the rate of work for the supporting pressure iswhere is the supporting pressure of the tunnel lining. Assuming that the pressure of the filler material is distributed along the karst cave in the radial direction, the rate of work of the filler material pressure is given bywhere is the pressure of the filler material in the karst cave. As stated in the upper bound theorem, the imposed loads cannot be carried by the rock massif; for any assumed failure mechanism, the rate of work conducted by the external forces exceeds the internal rate of energy dissipation. Thus, the first step to calculate the upper bound solution is to establish a virtual work equation composed of the rates of energy dissipation and external work. The virtual work equation, which is the difference between the rate of energy dissipation and the rate of external work, can be written as

Substituting equations (3)–(7) into (8), the expression of the virtual work equation is obtained:where is a functional that can be written as follows:

Because the virtual work equation includes the first derivative of , the virtual work equation can be regarded as an objective function, which is used to derive the upper bound solution of . Once the virtual work equation is established, the next step is to determine the variable that minimizes or maximizes the objective function. It can be observed that is the only variable in the objective function . Consequently, the real upper bound solution of is the solution for which the objective function reaches the extremum. Furthermore, the objective function is determined by . Therefore, the problem is converted into determining the function for which reaches the extremum. Based on the variational principle, the functional reaches the extremum when satisfies the following necessary condition:

Substituting (10) into equation (11), a linear homogeneous second-order differential equation with constant coefficients is obtained:

By solving this differential equation, the analytical expression of the collapse surface can be written aswhere and are two integration constants, and the values of these integration constants can be determined by the boundary conditions. Substituting (13) into (9), the expression of objective function is obtained:

Furthermore, the values of and are undetermined, which means that the karst cave and tunnel roof locations where the rock mass begins to collapse are still unknown. However, the values of *l*_{1} and *l*_{2} are required to derive the final expression of the collapse surface . Thus, it is necessary to establish the stress boundary condition to determine the values of and . Based on the stress distribution rule, there is no distribution of shear stress on the karst cave surface. Therefore, the shear stress of the element at the junction of the collapse surface and the karst cave surface is zero. The stress boundary condition can be expressed aswhere is the shear stress of the element at the junction of the collapse surface and the karst cave surface. According to the equilibrium differential equation, the shear stress of the element mentioned above can be written in the following form:where is the angle between the tangential and vertical directions of a random point on the collapse surface. Based on the trigonometric function, the relationship between and is given by

Substituting equations (16)-(17) and the expression of into (15), the value of *c*_{1} is obtained as follows:

From Figure 1, it can be seen that a geometric boundary condition should be satisfied. The geometric boundary condition can be expressed as

By substituting (13) and (18) into (19), an equation set that includes three unknowns, , , and , is obtained:

Therefore, we have to establish another equation to determine these unknowns. By equating the external rate of work due to external forces to the total internal energy dissipation along the collapse surface, the last equation, which is used to calculate unknowns , , and , is given by

Combining (20) and (21), a nonlinear equation that can be solved numerically is obtained. Using an iterative algorithm, the numerical solutions of , , and are obtained. Substituting the values of , , and into (13), the upper bound solution of the collapse surface can be expressed as

#### 4. The Shape and Range of the Collapse Surface for Different Parameters

The aim of this paper is to study the collapse mode of the rock mass between the concealed karst cave and highway tunnel roof. If the shape and range of the collapse surface are obtained, the collapse mode of the rock mass in this region can be examined intuitively. Based on the analytical equation of collapse surface expressed in (22), the coordinates of any point on the collapse surface are calculated. Using these point coordinates, the shape and range of the collapse surface can be drawn with the help of computer-aided drawing software. To investigate the influence of different parameters on the collapse surface of the rock mass, the collapse surfaces for parameters corresponding to *A* = 0.2–0.4, *B* = 0.6–0.8, = 1 MPa, = –, *T* = 60 kPa, = 13–19 m, = 8.21 m, = 2.5 m, and = 19–25 kN m^{−3} are illustrated in Figures 4–7. Similar to the failure mechanism shown in Figure 2, the collapse surface extends from the tunnel roof to the bottom of the karst cave, which forms an inverted funnel shape. Furthermore, the range of the collapse surface is affected by rock mass parameters , , , and . To obtain the change law of the collapse surface range as a function of these parameters, the univariate analysis method is employed in this paper. It can be observed in Figures 4–7 that greater values of and tend to induce a larger collapse surface range. In contrast, higher values of and lead to a decrease in the collapse surface range.

#### 5. The Minimum Safe Distance between a Karst Cave and Highway Tunnel

As mentioned above, the distance between a concealed karst cave and the proposed nearby highway tunnel is a major cause that induces the collapse of the rock mass in this region. Numerous investigators have used various methods to study the minimum safe distance between a karst cave and tunnel. However, presently, no theoretical equation can calculate the minimum safe distance between a karst cave and tunnel accurately. To solve this problem, the shapes of the collapse surfaces for different distances between the karst cave and tunnel are investigated. Figure 8 shows that the width of the top of the collapsing block decreases with increasing distance when the other parameters are fixed. When the distance between the karst cave and tunnel exceeds a critical value, the width of the top of the collapsing block would decrease to zero, which implies that the collapse surface and karst cave intersect at a certain point. If the distance between the karst cave and tunnel continues to increase, the collapse surface would not extend to the karst cave, and the pillar between the karst cave and tunnel will not collapse. Consequently, the critical distance between the karst cave and tunnel for which the collapse surface and karst cave intersect at a certain point can be regarded as the minimum safe distance between the karst cave and highway tunnel. Based on the above assumption, by substituting = 0 into equations (20)-(21), an equation set that includes two unknowns and is obtained:

Because (23) is a nonlinear equation set, the value of can be derived by the numerical method, and is the minimum safe distance between the karst cave and highway tunnel defined in this paper.

Based on the derivation procedure of the minimum safe distance , it is found that the value of is determined by several rock mass parameters. Thus, it is necessary to analyze the influence of these parameters on . Using (23), for different values of , , , and is plotted in Figure 9. Figure 9 shows that increases with increasing and but decreases with increasing and .

**(a)**

**(b)**

**(c)**

#### 6. Application of the Proposed Method in an Actual Engineering Project

To apply the proposed method in actual engineering, a highway tunnel located in a karst area of Guangxi Province is selected as a case study. The engineering geological conditions of this tunnel are described as follows. The Wuzhishan tunnel is a highway tunnel for the Leye-Baise highway in China. The maximum buried depth of the tunnel is 191 m, and the total length is 1230 m. This tunnel is excavated in an area where the surrounding rock has been eroded by tectonic denudation and karst is developed. The principal lithology penetrated by the Wuzhishan tunnel is intermediary weathered limestone and intermediary weathered sandstone. By taking advantage of geophysical prospecting, it was deduced that a karst cave exists above the tunnel roof in the section of KD119 + 319. The engineering geological conditions of this cross section are presented in Figure 10. Figure 10 shows that there is an irregular shape karst cave right above the planned highway tunnel, and the distance between the karst cave and the tunnel is approximately 5.7 m.

Based on the engineering geological and geometric conditions, an upper bound calculation model of this case study is constructed. The shape of the karst cave is simplified as a circle, and the radius of this circle is 3 m. Using the analytical equation of collapse surface expressed in equation (22) and the rock mass parameters listed in Table 1, the shape and range of the collapse surface for the Wuzhishan tunnel are plotted in Figure 11. From the figure, it is found that a collapse surface extends from the tunnel roof to the bottom of the karst cave, which indicates that the rock mass in this region has a risk of collapse during tunnel excavation. Furthermore, substituting the parameters listed in Table 1 into (23), the minimum safe distance between the karst cave and Wuzhishan tunnel is 5.84 m. According to the definition of the minimum safe distance, rock mass collapse can be avoided if the distance between the karst cave and Wuzhishan tunnel is larger than 5.84 m. Thus, the location of the proposed tunnel can be optimized in the planning stage to reduce the risk of rock mass collapse during tunnel construction.

#### 7. Comparison of the Numerical Simulation and Upper Bound Calculation

To verify the proposed method, the results obtained by the upper bound calculation are compared with those provided by numerical analysis using finite difference software FLAC3D. Because the internal contour of the Wuzhishan highway tunnel is multicircular, the tunnel model cannot be constructed by the modeling tool of FLAC3D directly. The tunnel model is constructed by using finite element software ANSYS, and the final model is obtained by importing the data into FLAC3D. The geometrical dimensions of this model are determined by the engineering geological conditions of the Wuzhishan tunnel mentioned above. The whole model is illustrated in Figure 12, and the dimensions of the model are as follows. The dimensions of the numerical model are 80 m in the *X* direction and 80 m in the *Z* direction. Because the failure mechanism established in the upper bound calculation is two-dimensional, the dimension of the numerical model in the *Y* direction is set to 2 m to improve the computational efficiency. The model contains approximately 8823 zones and 13462 nodes. The boundary conditions of the model were assumed as follows: the ground surface is free, and the horizontal displacement in the front and back surfaces and the vertical displacement in the bottom surface of the model are fixed. The multicircular outer contour of the Wuzhishan tunnel is composed of three arc sections. Moreover, the radius of the largest arc is 8.21 m, and the radius of the other two arcs is 5.55 m. To make the numerical model consistent with the upper bound failure mechanism, the karst cave is also assumed to be a circle, and the radius of this circle is 3 m. Based on the geological survey report, the distance between the karst cave and tunnel is approximately 5.7 m. A uniform pressure was applied to the surface of the karst cave to simulate the pressure of the filler material in the karst cave. The excavation process was simplified by assuming that the rock mass located in the excavation region was excavated instantaneously. Simultaneously, a certain uniform pressure was applied to the surface of the excavation region to simulate the supporting pressure produced by the lining of the tunnel, which can be seen in Figure 13.

Because the H–B failure criterion is used in the upper bound calculation, the H–B failure criterion was invoked in the simulation process to compare the upper bound calculation and numerical simulation under the same conditions. Furthermore, the rate of energy dissipation along the collapse surface is associated with the normal stress strain and shear stress strain, and the H–B failure criterion represented by the normal and shear stresses was employed in the upper bound calculation. However, the H–B failure criterion embedded in FLAC3D is represented in terms of the major and minor principal stresses. As the parameters of the H–B failure criterion for the two forms are different, the comparison of the collapse surfaces obtained by the upper bound calculation and numerical simulation cannot be achieved under the same conditions. To solve this discrepancy, Hoek and Brown [22] proposed a method to convert the parameters in one form of the H–B failure criterion into the parameters in another equivalent form. Based on this method, a set of equivalent parameters for the two forms of the H–B failure criterion were obtained, which are summarized in Table 2.

Using the model illustrated in Figure 12 and the parameters listed in Table 3, the rock mass deformation above the tunnel roof during the excavation process is simulated. To obtain the collapse surface of the rock mass in the limit equilibrium state, an indicator should be established to represent the zones in a certain region that enter this state. According to the manual of FLAC3D [23], the failure surface of the rock mass can be determined by the maximum shear strain value, and the maximum shear strain value can be represented by shear strain contours. Furthermore, numerous scholars have also used the shear strain contours provided by FLAC3D to define the failure surfaces of slopes [24–27]. Consequently, the contours of the maximum shear strain increment of the whole model are plotted to determine the collapse surface in the limit equilibrium state, as shown in Figure 14. It is clear that there are two shear failure belts extending from the tunnel roof to the bottom of the karst cave. The shear failure belts are symmetrical about the *z*-axis, whose shape is similar to the collapse surface derived from the upper bound calculation. To better compare the collapse surface provided by the upper bound theorem and that provided by numerical simulation, the upper bound solution of the collapse surface is plotted using red curves and superimposed on Figure 14 with the equivalent parameters listed in Table 2. Clearly, the collapse surface provided by the upper bound calculation approximately coincides with the shear failure belts. The similarity between the collapse surfaces provided by the upper bound solution and numerical simulation shows that the proposed method in this paper is valid.

#### 8. Conclusions

For a tunnel excavated in a karst terrain, excavation may induce collapse of the surrounding rock, which poses a threat to the personal safety of construction personnel. A possible prediction method of the collapse region induced by a concealed karst cave above a deep highway tunnel is presented within the framework of the upper bound theorem of limit analysis. A new failure mechanism composed of arbitrary curves is proposed to describe the collapse mode of the rock mass between the karst cave and tunnel. Based on the failure mechanism, an objective function including the equation of the collapse surface is derived from the relation between the rates of total energy dissipation and external work. The analytical expression of the collapse surface is obtained on the basis of the variational principle.

Using the analytical expression of the collapse surface, the shapes of the collapse surfaces for different parameter values are plotted. By studying the change law of the collapse surface as a function of these parameters, it is found that the collapse surface increases with increasing and but decreases with increasing and . Furthermore, a definition of the minimum safe distance between the karst cave and tunnel is proposed, and the computational equation of the minimum safe distance is derived.

Based on the geological survey report of a highway tunnel excavated in Guangxi Province, the collapse surface of the rock mass above this tunnel is calculated to verify the practicality of the proposed method. Moreover, the computed result of the proposed approach is compared with that derived from the numerical simulation technique under the same conditions. The similarity between the collapse surface results obtained by the upper bound solution and numerical simulation indicates that the proposed method is valid.

#### Data Availability

The 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 study was supported by the National Natural Science Foundation of China (Grants 51878074 and 51678071). Their financial support is greatly appreciated.