Research Article  Open Access
Zhiqiang Li, Guofeng Liu, Shuqian Duan, Shufeng Pei, Changgen Yan, "Sensitivity Analysis in the Estimation of Mechanical Parameters of Engineering Rock Mass Based on the Hoek–Brown Criterion", Mathematical Problems in Engineering, vol. 2020, Article ID 9418905, 18 pages, 2020. https://doi.org/10.1155/2020/9418905
Sensitivity Analysis in the Estimation of Mechanical Parameters of Engineering Rock Mass Based on the Hoek–Brown Criterion
Abstract
Geological strength index , disturbance factor (D), material constant , and uniaxial compressive strength of the intact rock are essential input parameters of the Hoek–Brown criterion. Mechanical parameters of the engineering rock mass, including elastic modulus , cohesion , and internal friction angle estimated by the H–B criterion, and the predicted excavation response of surrounding rock, including the displacement and excavation damage zone based on the MPs, are of high relevance with the four IPs of the H–B criterion. In this paper, the deep and huge underground cavern excavated in basalt from a hydropower station under construction in the southwest of China is used to analyse the sensitivity of the IPs on the MPs, the displacement, and EDZ of the surrounding rock mass. Firstly, the H–B criterion is applied to estimate the MPs, among which the IPs are obtained from a series of in situ and laboratory tests, including borehole camera observation, wave velocity test, uniaxial and triaxial compression tests, and so on. Secondly, the sensitivity relationships between IPs, MPs, and prediction results of displacement and EDZ are established and described quantitatively by the sensitivity factor (s_{i}). Results show that the MPs of the rock mass are more sensitive to and and are highsensitivity parameters affecting the displacement and EDZ. Finally, the variations in the estimated MPs and associated prediction results concerning excavation response, which are caused by the uncertainties in the determination of the IPs, are further quantified. This study provides a straightforward assessment for the variability of the rock mass parameters estimated by the H–B criterion. It also gives a valuable reference to similar geotechnical engineering for the determination of rock mass parameters in the preliminary design.
1. Introduction
The study on the mechanical properties and behaviors of rock mass is the basis of rock engineering practices due to its inherent uncertainty [1]. Moreover, it has aroused great attention to rock mass mechanical parameters (MPs) such as elasticity modulus , cohesion , and friction angle as essential elements to describe the properties and behaviors of rock mass in engineering [2–5]. On one hand, the reliable parameters are highly demanded in the preliminary engineering design and stability assessment phase. On the other hand, reliable acquisition of parameters is necessary to ensure construction safety and the engineering economy. Hence, the acquisition of more reliable MPs of largescale rock mass before engineering construction has always been the keystone and difficult problem in rock mass engineering practice [1].
Many methods have been proposed to solve the aforementioned difficulties. The method of converting rock parameters obtained by laboratory tests into engineering rock mass parameters, based on a rock mass classification system, has been increasingly recognized by specialists in the geotechnical engineering [6–9]. This method fully considers the integrated influence of joints and fissures in rock mass, groundwater, and size effects. In these methods, the Hoek–Brown (H–B) criterion is now considered to be one of the most welldeveloped methods, with the influence of structural characteristics on rock mass strength comprehensively reflected. The H–B criterion for rocks was firstly proposed by Hoek and Brown [10], according to the statistical analysis of laboratory triaxial tests and field tests for hard rocks, and it can be expressed bywhere and are major and minor effective principal stresses at failure, is the uniaxial compression strength of the intact rock, m is the H–B frictional constants, and s is the H–B material constants. The constant m is comparable to the frictional strength of the rock and s indicates how fractured the rock is and is related to the rock mass cohesion [11].
The H–B criterion has been frequently modified many times in the past few decades. The most commonly used expression (2002 edition) [12] is as follows:where a is the H–B exponent, m_{i} is the H–B frictional constant m for intact rock, and m_{b} is the H–B frictional constant m for rock mass. The parameter m_{b} is a reduced value of m_{i} which accounts for the strength reducing effects of the rock mass conditions defined by geological strength index (GSI) [11]. GSI is one kind of rock mass classification system. D is a new index introduced to account for nearsurface blast damage and stress relaxation in the rock mass.
According to (2), the critical point of using the H–B criterion is to determine its input parameters (IPs), including GSI, D, , and m_{i}. The difficulty of criterion application and the accuracy of rock mass parameters calculation could be determined by the accuracy of input parameters, to a great extent [13]. For GSI, a quantitative evaluation method that can consider the distribution rate of the discontinuous surface, roughness, weathering, and infilling properties was proposed in [14]. Hoek et al. [15] proposed a quantification method of the GSI chart based on two wellestablished parameters Joint Condition (Jcond_{89}) and Rock Quality Designation (RQD) to reduce the difference of rock mass quality estimated result between different geologists working on a single project. For m_{i}, a quantified calculation method was built by Hoek [16] based on the statistical analysis of the results of a set of triaxial tests. Cai [17] established the relationships between m_{i}, the crack initiation strength of the uniaxial compression test, and the direct tensile strength based on Griffith’s theory of strength. The relationship between D and the integrality coefficient (kv) calculated by the wave velocity of the rock and rock mass was established by Yan and Xu [18]. A disturbance rating for open pit mine slopes developed by Rose et al. [19] are based on the geologic structure, rock mass conditions, groundwater, in situ stresses, slope geometry, poor blasting, etc. More quantification proposed in these methods could reduce the subjectivity of IPs’ selection and also improve the calculation accuracy of rock mass MPs. However, there is uncertainty in the determination of IPs for the complexity of natural rock mass [20, 21]. The errors caused by this uncertainty are inevitable in practical engineering. There are two issues that should be taken seriously in engineering practice: (1) how the error of each IP (GSI, D, , and m_{i}) acts on the estimated results of rock mass MPs (E, c, and φ) and the predicted results of displacement and excavation damage zone (EDZ) of engineering rock mass and (2) how it influences the calculation accuracy of MPs, and what range of error could be acceptable with few effects on the preliminary engineering design.
Based on a deep and huge underground cavern from a hydropower station under construction in Southwest China, the sensitivity of the H–B criterion’s IPs on the estimated results of rock mass MPs as well as the predicted response results of the surrounding rock is presented in this paper. Firstly, the H–B criterion’s IPs are determined before the cavern excavation according to a series of necessary laboratory and field tests. Then the sensitivity function (S) and sensitivity factor (s_{i}) [22] are introduced to analyse the sensitivity between the MPs of rock mass and IPs within a specific fluctuation range. After that, the sensitivity relationships between the IPs and the corresponding displacement and EDZ of engineering rock mass are established. The research results provide an explicit theory of interpretation on the sensitivity relationships between the H–B criterion’s IPs and rock mass MPs, displacement, and EDZ, which has important guiding significance for the acquisition of rock mass MPs in the preliminary engineering design.
2. Sensitivity Analysis of the MPs Estimated by the H–B Criterion to the IPs
2.1. Engineering Background
The research conducted in this paper is based on an underground powerhouse cavern with a span of 34 m and a maximum height of 88.7 m at a hydropower station located in the southwest of China. The dimension of this cavern is 453 m × 34 m × 88.7 m (length × span × height) and the cavern has a depth of 420∼800 m in the horizontal direction and 420∼540 m in the vertical direction, with the axis direction of N10°W. It is a typical underground cavern of high side walls, large span, and great burial depth [23]. The main geological structures involved in this large cavern are shown in Figure 1. Several weak structural belts are distributed, mainly including the f20 steeply inclined fault zone, the large fracture of T813, and the weak interlayer zones (WIZ) of No. 3, No. 31, No. 4, and RS411. Therein, the f20 steeply inclined fault zone is welldeveloped, with a width of about 30 cm. The width of T813 is about 3∼5 cm. The thickness of the WIZs of No. 31 and No. 4 is 20 cm and 40∼60 cm, respectively. These are highly fractured coarsegrained materials, mainly composed of mud and debris with different sizes and shapes. The weak interlayer zone of RS411 outcrops in the upper part of the cavern arch, about 5 m thick, with densely developed cracks [24, 25].
The surrounding rock mass of the cavern is mainly composed of cryptocrystalline basalt, porphyritic basalt, almondshaped basalt, breccia lava, and tuff in the Upper Permian Emeishan Formation, with good integrity, mostly in massive and subblocky structure. The geostress is mainly dominated by tectonic stress, with horizontal stress higher than the vertical stress. The maximum principal stress involved in the cavern is within 22∼26 MPa, with the orientation of N0°∼20°E and the dip angle of 2°∼11°. The intermediate principal stress is about 14∼18 MPa, and the plunge is 0°∼8°. The minimum principal stress is approximately vertical, and its value is generally 13∼16 MPa. The initial strength/stress ratio of the rock is 3∼5. However, in the case of considering the redistributed stress after excavation, the maximum stress tends to increase to 50∼60 MPa, and the strength/stress ratio would become 2∼3. Overall, the underground cavern groups are located in a high stressed area [23, 26].
According to the engineering investigations, geological discontinuities involving the rock mass surrounding the cavern from chainage K0+300 to K0+330 are rather developed. Cut by several sets of large fractures and the weak interlayer zones of No. 3 and No. 31, rock mass in this section tend to be more fractured than that of other sections in the cavern (see Figure 2). The existence of these largescale discontinuities could sharply weaken the mechanical strength of rock mass and influence the stability of the cavern. Therefore, it is of considerable significance for accurately determining the MPs of the rock mass around this section in the design phase of the cavern engineering since it directly affects the excavation and support schemes, construction safety, and engineering cost. The cavern was excavated layer by layer from top to bottom by drilling and blasting method. The pilot tunnel of layer I was excavated firstly and then the expanding excavation towards two sides was conducted. The detailed excavation sequences for the first three layers at the K0+320 section are presented in Figure 3.
(a)
(b)
The following work will focus on the sensitivity and error analysis on the estimation of MPs of the rock mass surrounding the cavern from chainage K0+300 to K0+330. It will provide an illustration on how to evaluate the rationality of the rock mass MPs estimated by the H–B criterion and the related excavation response results predicted based on the estimated MPs.
2.2. Determination of the IPs
2.2.1. Estimation of GSI
Fifteen monitoring sections are arranged along the underground cavern, among which six sections are used to monitor the cracking process and associated EDZ of surrounding rock. The monitoring holes are symmetrically set on the upstream and downstream side of the cavern. The remaining sections are set to monitor the displacement of surrounding rock with the excavation, each of which is equipped with three sets of multipoint extensometers. As shown in Figure 3, there are two monitoring sections arranged in the zone from K0+300 section to K0+330, namely, sections K0+320 and K0+330 for monitoring both the EDZ and displacement of the surrounding rock. It needs to be mentioned that some small auxiliary caverns are designed and preexcavated around the large underground caverns in underground hydropower engineering, such as the anchorage caverns and drainage caverns, which are very helpful for setting the preexcavated observation boreholes towards the concerned cavern, such as the boreholes RK0+32001 and RK0+32002 (with 110 mm in diameter and approximately 26.5 m in depth) shown in Figure 3(a). In addition, some boreholes for EDZ observation could be flexibly set after the excavation of the pilot tunnel, such as the boreholes RK0+32011 and RK0+32012 in Figure 3(a). Herein, the digital borehole camera system (see Figure 4(a)) that can capture a 360degree panoramic image of the borehole wall is used as an efficient way to recognize the geological structure planes that cut through the boreholes, as well as the cracking process of surrounding rock. The highest circumferential accuracy of the borehole camera is 0.1∼0.2 mm. Details on the principle and components of the borehole camera system can be found in work by Liu et al. [23]. The generalized structure of rock mass surrounding the boreholes K0+32001 and K0+32002 are, respectively, interpreted based on the observations results from the borehole camera and the drilling cores as shown in Figures 4(b) and 4(c), for example, and shown in Figures 4(d) and 4(e). It indicates that the rock mass around the cavern roof at the section K0+320 is mainly composed of almondshaped basalt and relatively broken due to the cut of welldeveloped joints.
(a)
(b)
(c)
(d)
(e)
The value of GSI can be calculated by the equation proposed by Hoek et al. [15] as expressed in (3). The Jcond_{89} involved in (3) was estimated based on the above interpretation results of the generalized structure of rock mass surrounding the boreholes. The value of the RQD was calculated according to the coring results. The final calculation result shows that the average value of the GSI of rock mass surrounding the cavern roof at the section K0+320 is 53 (55 and 50 for the two boreholes, respectively). The detailed data involved in the calculation are shown in Table 1:
2.2.2. Measurement of σ_{ci}, m_{i}, and D
Core samples from the boreholes of the concerned monitoring section K0+320 were partially processed into standard rock specimens for laboratory mechanical tests. The uniaxial compression test was performed on the RMT150C multifunctional test system to obtain the uniaxial compressive strength (). In the test, the axial pressure was applied at a fixed loading rate of 0.001 mm/s, and the distribution of of almond basalt was obtained when 20 sets of rock tests are repeated (see Figure 5). obeys the normal distribution, in which the mean value is 147 MPa, and the variance is 24.33. According to the 95% confidence interval, the average of the almondshaped basalt is within 137∼157 MPa.
The constant m_{i} could be calculated by the equation proposed by Hoek [16]:where n is the number of the triaxial tests and σ_{1i} and σ_{3i} are the major and minor principal stresses.
The triaxial compression test was performed on the MTS81504type rock mechanical triaxial testing system. During the test, the axial pressure was applied at a fixed loading rate of 0.02 mm/min. The confining pressure is first applied to a predefined value, and then the axial pressure is applied until the failure of specimen. The confining pressures adopted in the conventional triaxial tests were 0 MPa, 5 MPa, 10 MPa, 20 MPa, and 40 MPa, respectively. The typical stressstrain curves of almondshaped basalt are shown in Figure 6. The obtained mechanical parameters, including the strength, elastic modulus, and Poisson’s ratio, are shown in Table 2. Afterward, the calculation result indicates that the f almondshaped basalt is 26, consistent with the recommended value of 25 ± 5 for basalt by Marinos and Hoek [28].

The parameter D could be calculated by the following equation proposed by Yan and Xu [18]:where indicates the longitudinal wave velocity (LWV) in the rock mass and indicates the LWV in rock specimen.
The in situ acoustic velocity testing of the rock mass was carried out in the monitoring section K0+320 by the wave velocity tester (see Figure 7) through the singlehole method. The test boreholes located in the cavern roof as shown in Figure 3 could also be used for acoustic velocity testing at any time. During the test, the water was used as a coupling medium between the probe and the borehole wall, and the wave velocity data were collected every 20 cm along the borehole. The measured average LWV of the rock mass around the borehole is 4206 m/s. Meanwhile, the average LWV of rock specimens obtained in the laboratory wave speed tests is 4773 m/s. Afterward, the value of D can be obtained as 0.21.
2.3. Sensitivity of the MPs Estimated by the H–B Criterion to the IPs
The values of the IPs used to calculate the MPs of the rock mass surrounding the cavern from chainage K0+300 to K0+330 were determined as stated above. However, there is always uncertainty in the determination of IPs [20, 21]. For the complexity of natural rock mass, this uncertainty is prevailing in the analysis of geotechnical engineering. Carter and Miller [20] divided this uncertainty into two types based on its causes. Besides, they also established the relationship between uncertainties due to natural variability and data inadequacies as a function of data accuracy and cost of site investigations. Baecher and Christian [29] and others have described these sources of uncertainty as being aleatory or epistemic. Aleatory uncertainty is the irreducible randomness or variability associated with phenomena that are naturally variable in time or space, even when the system is well known. Epistemic uncertainty, on the other hand, arises from limitations in our fundamental knowledge or understanding of some aspects of a problem [30]. During the estimation of the MPs, the uncertainties arise from a series of data collection processes, such as the changes in the location of the observation boreholes, the errors of field tests, the changes in the number of rock specimens for laboratory tests and the test errors, and the subjective errors in the artificial analysis. There are errors in the determination of IPs, while these errors are always limited. In practice, the values of the IPs of the H–B criterion inevitably fluctuate within a certain range, and the estimation results of MPs will inevitably fluctuate accordingly. People are more concerned that what degree of variation range in the estimation results of MPs would be caused by the fluctuation range of IPs.
For example, the average value of GSI at the K0+320 section is estimated as 53 by using the method introduced above. Perhaps the result will be 45, 55, or other values close to 53 if another method is used by other experts in the same condition. However, the result is unlikely to be 30 or 75, which is far from 53. The errors in the determination of other IPs (i.e., , m_{i}, and D) also have the same characteristics. To simulate the effect of this error on the calculation results of MPs, displacement, and EDZ, IPs determined by actual investigations and tests above were taken as the reference value, and then the fluctuation ranges of IPs should be artificially given. However, in order to understand more about the sensitivities of the estimation results to these IPs, the fluctuation ranges of IPs are considered as much as possible, as shown in Table 3. The GSI and D are two parameters related to the rock mass. For GSI, its fluctuation range is determined as 15–65 in which the basic H–B criterion can be better used [15]. For D, the complete range (01) is given. The constants m_{i} and are related to the intact rock. According to the studies by Marinos and Hoek [28] and the results of laboratory tests, the ranges of 20–30 and 137–157 MPa, which the almondshaped basalt can cover, were chosen as the fluctuation range for m_{i} and , respectively.

Based on the H–B criterion, a series of estimation results of MPs corresponding to the IPs in the fluctuation ranges are obtainable, and then the sensitivity analyses of MPs to IPs can be carried out. The sensitivity function (S) and sensitivity factor (s_{i}) [22] used to describe the sensitivity relationships between IPs, MPs, and the predicted results of displacement and EDZ are defined as in (6) and (7)). To illustrate the physical meaning of s_{i}, (7) is transformed into the expression, as shown in (8). The numerator of (8) is defined as the variation of characteristic δ_{y} and expressed as (9). The denominator of (8) is defined as the error of factor δ_{x} and expressed as (10). Hence, it can be found that s_{i} represents the multiple relationships between δ_{y} and δ_{x}:where S is the sensitivity function, s_{i} is the sensitivity factor, and x_{i} indicate different factors, δ_{x} is the error of factor, and y_{i} indicate different characteristics, δ_{y} is the variation of characteristics, and F represents the functional relationship between y and x.
In this study, the IPs are considered as different factors, and MPs are regarded as different characteristics. represents the functional relationship between MPs and IPs. Among the MPs, the parameters c and φ can be calculated according to (11) and (12) and the iterative procedures proposed by Sharan [31]. Equation (13) proposed by Hoek et al. [12] is applied to calculate the rock mass parameter E. Then, the functional relationships between MPs and every IP can be established. After that, the sensitivity of each MP of the rock mass to every IP is deeply analysed according to (6)–(10). Herein, the dimensions of various IPs are normalized through a regularization transformation by (14) in order to make a better comparative analysis on the sensitivity of each parameter:where and are linearization parameters for the material, is the minor principal stress, and represents the functional relationship between major and minor principal stress:where k indicates the normalized result of the IP, x_{i} indicates different IP (GSI, D, m_{i}, and ), and x_{max} and x_{min} indicate the maximum and minimum value of every IP in its fluctuation range, respectively.The sensitivity of E to GSI is taken as an example to illustrate the process of sensitivity analysis. The parameter E corresponding to different values of GSI can be calculated by (13) [12] as shown in Table 4 when the remaining IPs (D, m_{i}, and ) are kept unchanged at the reference value. Then the functional relationship between E and GSI is established based on data fitting. The sensitivity function (S) of E to GSI is obtained according to (6) and s_{i} of E to GSI can be calculated by S accordingly (see Table 4). Similarly, the sensitivity analysis of other MPs corresponding to each IP is conducted. Figure 8 shows the variation curves of the calculated MPs corresponding to the whole fluctuation range of each IP.

(a)
(b)
(c)
The elastic modulus of rock mass is an important parameter in any analysis of rock mass behavior [32]. It was calculated by (13), which only involved GSI and D. Figure 8(a) shows that the sensitivity of E to GSI is more significant than D and it increases linearly with GSI. And the sensitivity of E on D is also gradually increasing as D increases.
Figure 8(b) indicates that the estimated result of parameter c is most sensitive to GSI, and the sensitivity tends to be higher at an increasing rate with the increase of GSI. The sensitivity of c to D shows the same significant rising tendency. In the fluctuation ranges of and m_{i}, the value of s_{i} has almost no change, while the sensitivity of c to is higher than m_{i} according to the comparison. Besides, there is a significant change in sensitivity priority over the entire variation range for four IPs. When k changes in the range of “0 < k < 0.5,” the sensitivity of c to , m_{i}, and D is ranked as > m_{i} > D. In the range of “0.5<k < 0.6,” the sensitivity of c changes to be ranked as > D > m_{i}. In the range of “0.6<k < 1,” the sequence of the sensitivity changes as D > > m_{i}.
From Figure 8(c), when k changes in the range of “0 < k < 0.75,” the sensitivity of parameter φ to GSI is the highest. In the range of “0.75 < k < 0.1,” the IP with the highest s_{i} is D. The sensitivities of parameter φ to and m_{i} have almost no change over the whole fluctuation ranges of k.
As described above, within the given fluctuation ranges of IPs in Table 3, the estimation of rock mass deformation parameter (E) and strength parameters (c and φ) are more sensitive to the GSI and D. Therefore, it is vital to determine the reasonable value of GSI and D when estimating the rock mass MPs by using the basic H–B criterion.
3. Sensitivity Analysis of the Displacement and EDZ to the IPs
3.1. Numerical Simulation on the Excavation Response of Rock Mass
In the design stage of excavation and rock support for such a large cavern, numerical simulation analysis is essential for the decisionmakers to predict the probable response behaviour of rock mass induced by adjacent excavations. According to the MPs of rock mass estimated above, the deformation and failure behaviour of the surrounding rock around the concerned cavern section under excavation had been simulated in advance by using the FLAC3D program (Fast Lagrangian Analysis of Continua [33]). The MPs of rock mass and WIZ are shown in Table 5, among which the parameters of rock mass are calculated by (11)–(13) based on the reference value of IPs in Table 3. The WIZ is a permanent deformation zone with different thickness and spacing in hard stratified rock masses in areas where regional geotectonic movement has been historically active. The MPs of WIZ are determined by the recommended values given by Duan [24] and Duan et al. [25]. The commonly used strainsoftening constitutive model is adopted to describe the mechanical behavior of rock mass and WIZs.

Failure approaching index (FAI), proposed by Zhang et al. [34], is a practical index and widely used to describe the mechanical response status of surrounding rock. The higher the FAI is, the more cracks that occurred in the surrounding rock mass. FAI < 1 means that the rock mass is at an elastic state. FAI = 1 indicates that the stress reaches the peak strength of rock mass, FAI = 2 indicates that the rock mass is completely destroyed, and FAI ＞ 1 means the surrounding rock is at a plastic state. Therefore, FAI could be used as a straightforward index to assess the stress concentration area and the damage degree of the surrounding rock mass. More details on the index of FAI could be found in the work by Zhang et al. [34]. The numerical prediction result of EDZ concerning the cavern section K0+320 after the excavation of the central pilot is shown in Figure 9(b). Figure 9(c) shows the actual investigation result of EDZ along the boreholes K0+32011 and K0+32012 by using the wave velocity test. It shows that the simulated depths (FAI > 1) of EDZ in the direction of K0+32011 and K0+32012 monitoring boreholes are 1.57 m and 2.84 m, respectively. The actual EDZ depths of K0+32011 and K0+32012 observation boreholes are 1.9 m and 2.5 m, respectively.
(a)
(b)
(c)
The excavation response of the surrounding rock around the section K0+320 of the cavern after the excavation of layer I was simulated. The maximum principal stress redistributed in the rock mass is shown in Figure 10(a), and it indicates that the stress concentration after the excavation of layer I occurring at the spandrel on the upstream side and the arch foot on the downstream side of the cavern. Figure 10(b) shows the EDZ of the surrounding rock which was evaluated by FAI, and it can be seen that the strongly damaged zones were consistent with the stress concentration locations. During the in situ investigation work over the period of the excavation from chainage K0+300 to K0+330, it was found that the surrounding rock at the spandrel on the upstream side as well as the arch foot on the downstream side was subjected to obvious cracking from the surface of the cavern wall to the inside gradually and even developed into spalling damage and stressstructure controlled collapse in some positions (see Figures 10(c) and 10(d)).
(a)
(b)
(c)
(d)
Based on the field observation, the actually damaged zones were in accordance with the simulation result. It should be noted that the asymmetric feature of the EDZ on the two sides of the cavern was related to the initial principal stress direction on the cross section. As shown in Figure 9(a), the orientation of the maximum principal stress is at a certain angle from the horizontal plane. Meanwhile, the above analysis also concluded that stressinduced or stressstructure controlled rock failures always occur at the excavation profile parallel to or at a small angle with the orientation of maximum farfield principal stress in the cross section of the tunnel [3, 35–38]. Figure 11 shows the simulated displacement of surrounding rock after the excavation of layer I, and it is observed that the position of the larger displacement is at the spandrel on the downstream side, different from the EDZ. That is because the simulation in this study was carried out based on the continuum analysis, and the displacement along the unloading direction of maximum principal stress, therefore, would be bigger than other positions.
Moreover, the depths of EDZ along the boreholes K0+32001 and K0+32002 were determined by the wave velocity tests (see Figure 12), and the displacements of surrounding rock could be recorded by the multipoint extensometers. The comparison between the simulation results and actual investigated results of the excavation response concerning several specific observation points around the roof of the cavern is shown in Table 6, and the high agreement between the two indicates that the rock mass MPs estimated above and the reference value of IPs are basically reasonable.
 
Note: the actual displacement of point M21 is not effectively recorded. 
3.2. Sensitivity of Displacement and EDZ to IPs
The predicted displacement and EDZ are closely related to the MPs and therefore affected by the error in the estimation of IPs. Hence, it is necessary to analyse the sensitivity between the prediction results (displacement and EDZ) and IPs, similar to the sensitivity analysis procedures as illustrated above. First of all, several key points where the displacement or EDZ is of concern around the roof of the cavern should be selected to generate the data samples used for sensitivity analysis. Herein, the dataset concerning the displacement of point M11, M21, and M31 that were monitored by the multipoint extensometer as shown in Figure 3 and the largest depth of the EDZ are employed. Then, a series of prediction results of the displacement and EDZ corresponding to the IPs in the fluctuation ranges (see Table 4) are obtained. After that, the sensitivity analyses of the predicted displacement (EDZ) to the IPs can be carried out according to (6) and (7). The obtained sensitivity curves are shown in Figure 13.
(a)
(b)
(c)
(d)
The sensitivity curves of the predicted displacement at different points to the IPs have the similar sensitivity regularities (see Figures 13(a)–13(c)). It shows that the predicted displacement is most sensitive to GSI and insensitive to m_{i} and . The sensitivity of the displacement increases with D. The sensitivity of the predicted EDZ to IPs is more complicated than that of the displacement (see Figure 13(d)). The sensitivity of EDZ to GSI increases first and then decreases, finally increasing over the entire range of GSI. The sensitivity of EDZ to D continue to increase over the entire range of D. When k changes in the range of “0 < k < 0.5” and “0.9 < k < 1.0,” the sensitivity of EDZ to GSI is the highest. In the range of “0.5 < k < 0.9,” the parameter with the highest s_{i} is D.
4. Discussion
4.1. Comparison between the StrainSoftening and SShaped Failure Criterion
Numerical simulation is a technique of major importance in various technical and scientific fields [39]. For rock engineering, such simulation is essential for studying the fundamental processes occurring in rocks and for rock engineering design [40]. On the one hand, the prediction results of excavation response (displacement and EDZ) from numerical simulations could be used to illustrate the rationality of the reference IPs in comparison with actual results. On the other hand, the numerical simulation could be employed to predict the response of surrounding rock under excavation corresponding to different IPs so as to construct the dataset for sensitivity analysis of related parameters. Hence, the efficiency of the numerical simulation must be guaranteed. Choosing a proper constitutive model that can describe the mechanical behavior of rock mass correctly is of importance to the rationality of simulation results. In above analysis, the strainsoftening (SS) criterion is chosen to carry out the numerical simulation for its effectivity at characterizing the mechanical behavior of deep hard rock [41–44]. According to the literatures, it is suggested that the failure envelop for the entire confinement range of brittle rocks and rock masses is distinctly sshaped. Thus, the Sshaped failure criterion, considered tensile fracturing or spalling at low confinement and extensiledriven shear at high confinement, is proposed [45–47]. The two segments of the criterion are connected by a transition zone where σ_{3} = /10. Therefore, it is necessary for demonstrating the applicability of Sshaped criterion in the engineering presented in this study and illustrating the differences of the prediction results based on the Sshaped criterion and SS criterion.
According the in situ stress condition and of almondshaped basalt, the rock mass around the K0+320 section of the cavern is at high confinement. The Sshaped criterion can be applied to FLAC3D on the basis of modified HoekBrown constitutive model and a userdefined FISHfunction [48]. The parameters involved in the simulation are shown in Table 7.

The displacements of the points M11 and M21 (see Figure 3) located in the multipoint extensometers are recorded during the simulations. Besides, the depth of EDZ along the direction of test boreholes K0+32001 and K0+32002 are extracted. The comparison between the predicted results based on above two criteria and the actual results are shown in Figure 14. It is found that the prediction results of displacement and EDZ based on the two criteria are largely consistent. The prediction results of excavation response simulated based on the two criteria are nearly close; that is, both criteria could be adopted in the numerical analysis in this study.
(a)
(b)
4.2. Variations in the Estimated MPs, Calculated Displacement, and EDZ of the Rock Mass Caused by the Fluctuations of IPs
The relation of the IPs to the MPs estimated by H–B criterion and the predicted displacement (EDZ) of the rock mass have been established based on the sensitivity analysis and described quantitatively by the sensitivity factor (s_{i}). The sensitivity factors of the estimated results to the IPs under the reference values are presented in Table 8.

The errors caused by the uncertainty in the determination of IPs are inevitable in practical engineering, which would influence the accuracy of the estimated MPs and predicted excavation response of the rock mass. What sort of variations in the estimation and prediction results would be made by the errors of IPs and what range of variations could be acceptable with few effects on the engineering design and preliminary risk estimation are all worthy of further discussion.
Taking the sensitivity of MPs, displacement, and EDZ of the rock mass to two highsensitivity IPs under the reference values as an example, the estimation results and variations of MPs and predicted displacement and EDZ corresponding to the different fluctuations of the two highsensitivity IPs (such as ±5%, 10%, ±15%, ±20%, and ±25%) are shown in Figures 15 and 16, respectively.
(a)
(b)
(c)
(a)
(b)
Figure 15 shows the variations in the estimation of the MPs caused by the fluctuations in different degrees in the determination of each IP. For a certain fluctuation in the determination of IPs, the larger the associated variation of the MP is, the more sensitive the estimated result is to the IP. For the parameter E, 5% fluctuation in determination of D causes 0.6% variation in the estimation of E, while 5% fluctuation in determination of GSI could cause 14% variation in the estimation of E. Namely, GSI is a very important IP needing to be paid more attention in the determination process when using the H–B criterion to estimate elastic modulus. For the parameter c, GSI is still the sensitive IP that mostly affects the estimation result compared with other IPs, but the variation in the estimation result is relatively smaller than that of parameter E under the same fluctuation condition in GSI, and it means that the accuracy of parameter c is more likely to be assured than parameter E in the estimation work. By contrast, the parameter φ can be estimated with a high accuracy, and large variations in estimation result are produced only in the case of great fluctuation in the IPs as shown in Figure 15(c).
The variations in estimated results of the excavation response caused by the fluctuations in different degrees in the determination of each IP are shown in Figure 16. It can be seen that the fluctuations in the determination of GSI and could cause the variations in the estimation of the EDZ, but the influences are limited, and the high accuracy of the EDZ estimation therefore could be assured. For the displacement, GSI is still a highly relevant IP that could significantly cause large variations in the prediction result, and that is because the estimation of E is highly related to GSI as discussed above. In engineering practice, the deformation behaviour induced by excavation in hard rock is generally small and cannot be the main issue during construction. However, for soft rock mass characterized by the large deformation, especially under high geostress condition, the prediction of surrounding rock displacement is an important work, which responds to a need for the accurate MPs of the rock mass, and the GSI should be carefully determined in this scenario if the estimation method of MPs based on the H–B criterion is employed.
It should be noted that the above discussion points on the accuracy of the parameters estimation are drawn under the specific rock mass conditions presented in this study. Actually, the sensitivity of estimation results to the IPs under different reference values is different, which can be observed from Figures 8 and 12. In other words, the accuracy level of each parameter in the estimation is not fixed and mainly depends on the geological conditions, rock mass structure, the strength of intact rock, excavation disturbance, etc. Therefore, for specific engineering condition, such sensitivity analysis process needs to be conducted once again to evaluate the accuracy of rock mass parameters in the H–B criterionbased estimation.
5. Conclusions
The Hoek–Brown criterion is used to estimate the rock mass parameters of a deep underground cavern at a hydropower station located in the southwest of China. The reference values of IPs involved in the H–B criterion are obtained from the in situ tests and laboratory tests that were carried out in the concerned section of the cavern and then used to estimate the MPs of the rock mass. Based on the estimated results, the excavation response, including the displacement and EDZ of the surrounding rock, was predicted through numerical simulations. The estimation results were validated by the in situ investigations, illustrating the reasonability of the reference values of the IPs determined before. The sensitivity analysis of the MPs and excavation response behavior of the rock mass to the IPs of H–B criterion was deeply carried out on this basis. The results show the following:(1)Different parameters of the rock mass to be estimated have different sensitivity to the input parameters of the H–B criterion. Generally, the deformation parameter (E) and strength parameters (c and φ) are more sensitive to the GSI and D. In order to acquire reasonable MPs of the rock mass in the preliminary design phase of the engineering according to the H–B failure criterion, the key is that appropriate technique should be applied that take into account the determination of GSI and D.(2)The GSI and are highsensitivity parameters that affect the displacement and EDZ of the surrounding rock. The sensitivity of displacement and EDZ to the IPs is ranked as GSI > > m_{i} > D at the reference values concerning the cavern engineering present in this study. However, it should be noted that different reference values will lead to different conclusions.(3)The sensitivity of estimation results to the IPs under different reference values are different, and variations in the estimation of MPs and associated predicted excavation response of the surrounding rock, which are caused by the uncertainties in the determination of IPs, could be quantified via an error analysis. It should be noted that the accuracy of the estimation result of each parameter are usually not fixed and mainly depend on the geology conditions, rock mass characteristics, and so on.(4)The sensitivity factor (s_{i}) in the estimation of MPs based on the H–B criterion is of great value to engineering practice, which can be used to estimate what range of error could be acceptable with few effects on the preliminary engineering design.
Abbreviations
GSI:  Geological strength index 
D:  Disturbance factor 
m_{i}:  Frictional constant 
:  Uniaxial compressive strength 
IP:  Input parameter 
H–B:  Hoek–Brown 
MP:  Mechanical parameter 
E:  Elastic modulus 
c:  Cohesion 
φ:  Internal friction angle 
EDZ:  Excavation damage zone 
:  The major effective principal stresses at failure 
:  The minor effective principal stresses at failure 
m:  The H–B material constant 
s:  The H–B material constant 
m_{b}:  The H–B frictional constant (rock mass) 
a:  The H–B exponent 
Jcond_{89}:  Joint Condition 
RQD:  Rock Quality Designation 
Kv:  The integrality coefficient 
S:  The sensitivity function 
s_{i}:  The sensitivity factor 
WIZ:  Weak interlayer zone 
n:  The number of triaxial tests 
σ_{1i}:  The major principal stresses 
σ_{3i}:  The minor principal stresses 
:  Poisson’s ratio 
V_{pm}:  The longitudinal wave velocity in rock mass 
V_{pr}:  The longitudinal wave velocity in rock 
LWV:  Longitudinal wave velocity 
x:  The different factors 
y:  The different characteristics 
F:  The functional relationship between y and x 
k:  The normalized result of input parameter 
x_{i}:  The different input parameters 
x_{max}:  The maximum value of every input parameter in its fluctuation range 
x_{min}:  The minimum value of every input parameter in its fluctuation range 
FAI:  Failure approaching index 
δ_{x}:  The error of factors 
δ_{y}:  The error of characteristics 
SS:  Strainsoftening. 
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 no conflicts of interest.
Acknowledgments
This research was supported by the Basic Research Program of Natural Science from Shaanxi Science and Technology Department (Grant no. 2019JQ171), the Fundamental Research Funds for the Central Universities, CHD (Grant nos. 300102210110), the National Natural Science Foundation of China (Grant no. 51909241 and 41927806), and the Natural Science Foundation of Shaanxi Province of China (Grant no. 2018JQ5001).
References
 E. T. Brown, “Estimating the mechanical properties of rock masses,” in Proceedings of the Southern Hemisphere International Rock Mechanics Symposium, Perth, Australia, September 2008. View at: Google Scholar
 E. Hoek, “Strength of jointed rock masses,” Géotechnique, vol. 33, no. 3, pp. 187–223, 1983. View at: Publisher Site  Google Scholar
 E. Hoek, P. K. Kaiser, and W. F. BAawdn, Support of Underground Excavations in Hard Rock, A. A. Balkema, Rotterdam, Netherlands, 1995.
 P. Marinos and E. Hoek, “Estimating the geotechnical properties of heterogeneous rock masses such as flysch,” Bulletin of Engineering Geology and the Environment, vol. 60, no. 2, pp. 85–92, 2001. View at: Publisher Site  Google Scholar
 K. Hong, E. Han, and K. Kang, “Determination of geological strength index of jointed rock mass based on image processing,” Journal of Rock Mechanics and Geotechnical Engineering, vol. 9, no. 4, pp. 702–708, 2017. View at: Publisher Site  Google Scholar
 M. Cai, P. K. Kaiser, H. Uno, Y. Tasaka, and M. Minami, “Estimation of rock mass deformation modulus and strength of jointed hard rock masses using the GSI system,” International Journal of Rock Mechanics and Mining Sciences, vol. 41, no. 1, pp. 3–19, 2004. View at: Publisher Site  Google Scholar
 V. Marinos and T. G. Carter, “Maintaining geological reality in application of GSI for design of engineering structures in rock,” Engineering Geology, vol. 239, pp. 282–297, 2018. View at: Publisher Site  Google Scholar
 Q. Zhang, X. Huang, H. Zhu, and J. Li, “Quantitative assessments of the correlations between rock mass rating (RMR) and geological strength index (GSI),” Tunnelling and Underground Space Technology, vol. 83, pp. 73–81, 2019. View at: Publisher Site  Google Scholar
 M. Cai, P. K. Kaiser, Y. Tasaka, and M. Minami, “Determination of residual strength parameters of jointed rock masses using the GSI system,” International Journal of Rock Mechanics and Mining Sciences, vol. 44, no. 2, pp. 247–265, 2007. View at: Publisher Site  Google Scholar
 E. Hoek and E. T. Brown, “Empirical strength criterion for rock mass,” Journal of Geotechnical Engineering Division, vol. 106, no. 9, pp. 1013–1035, 1980. View at: Google Scholar
 Ulusay and J. A. Hudson, The Complete ISRM Suggested Methods for Rock Characterization, Testing and Monitoring: 1974–2006, ISRM Turkish National Group, Ankara, Turkey, 2007.
 E. Hoek, C. Carranzatorres, and B. Corkum, “HoekBrown failure criterion—2002 edition,” in Proceedings of the NARMSTac, University of Toronto Press, Toronto, Canada, 2002. View at: Google Scholar
 S. Sinha and G. Walton, “A progressive Sshaped yield criterion and its application to rock pillar behavior,” International Journal of Rock Mechanics and Mining Sciences, vol. 105, pp. 98–109, 2018. View at: Publisher Site  Google Scholar
 H. Sonmez and R. Ulusay, “Modifications to the geological strength index (GSI) and their applicability to stability of slopes,” International Journal of Rock Mechanics and Mining Sciences, vol. 36, no. 6, pp. 743–760, 1999. View at: Publisher Site  Google Scholar
 E. Hoek, T. G. Carter, M. S. Diederichs et al., “Quantification of the geological strength index chart,” in Proceedings of 47th US Rock Mechanics/Geomechanics Symposium, San Francisco, CA, USA, June 2013. View at: Google Scholar
 E. Hoek, Practical Rock Engineering, A. A. Balkema, Rotterdam, Netherlands, 2000.
 M. Cai, “Practical estimates of tensile strength and HoekBrown strength parameter m_{i} of brittle rocks,” Rock Mechanics and Rock Engineering, vol. 43, no. 2, pp. 167–184, 2010. View at: Publisher Site  Google Scholar
 C. B. Yan and G. Y. Xu, “Modification of HoekBrown expressions and its application to engineering,” Chinese Journal of Rock Mechanics and Engineering, vol. 24, no. 22, pp. 4030–4035, 2005. View at: Google Scholar
 N. D. Rose, M. Scholz, J. Burden et al., “Quantifying transitional rock mass disturbance in open pit slopes related to mining excavation,” in Proceedings of the XIV International Congress on Energy and Mineral Resources, Seville, Spain, April 2018. View at: Google Scholar
 T. G. Carter and R. I. Miller, “Crownpillar risk assessmentplanning aid for costeffective mine closure remediation,” Transactions of the Institution of Mining and Metallurgy, vol. 104, pp. A41–A57, 1995. View at: Google Scholar
 J. A. Hudson and X. T. Feng, Rock Engineering Risk, CRC Press, Boca Raton, FL, USA, 2015.
 W. S. Zhu and M. C. He, Stability of Rock Surrounding in Complex Condition and Dynamic Construction Mechanics, Science Press, Beijing, China, 1995.
 G. F. Liu, X. T. Feng, Q. Jiang et al., “In situ observation of spalling process of intact rock mass at large cavern excavation,” Engineering Geology, vol. 226, pp. 52–69, 2017. View at: Publisher Site  Google Scholar
 S. Q. Duan, Study on the Deformation and Failure Mechanisms and Overall Stability of Rock Mass with Interlayer Staggered Zones in Large Underground Caverns under High Geostress, Institute of Rock and Soil Mechanics, Chinese Academy of sciences, Wuhan, China, 2017.
 S. Q. Duan, Q. Jiang, D. P. Xu et al., “Experimental study of mechanical behavior of interlayer staggered zone under cyclic loading and unloading condition,” International Journal of Geomechanics, vol. 20, no. 3, pp. 52–69, 2020. View at: Publisher Site  Google Scholar
 X.T. Feng, S.F. Pei, Q. Jiang, Y.Y. Zhou, S.J. Li, and Z.B. Yao, “Deep fracturing of the hard rock surrounding a large underground cavern subjected to high geostress: in situ observation and mechanism analysis,” Rock Mechanics and Rock Engineering, vol. 50, no. 8, pp. 2155–2175, 2017. View at: Publisher Site  Google Scholar
 Z. T. Bieniawski, Engineering Rock Mass Classifications, Wiley, New York, NY, USA, 1989.
 P. Marinos and E. Hoek, “GSI: a geologically friendly tool for rock mass strength estimation,” in Proceedings of the International Conference on Geotechnical and Geological Engineering, pp. 19–24, Lancaster, UK, November 2002. View at: Google Scholar
 G. B. Baecher and J. T. Christian, “Reliability and statistics in geotechnical engineering,” Engineering Geology, vol. 30, pp. 237238, 2005. View at: Google Scholar
 E. T. Brown, “Risk assessment and management in underground rock engineeringan overview,” Journal of Rock Mechanics and Geotechnical Engineering, vol. 4, no. 3, pp. 193–204, 2012. View at: Publisher Site  Google Scholar
 S. K. Sharan, “An iterative method for the linearization of nonlinear failure criteria,” Advances in Engineering Software, vol. 117, pp. 89–94, 2018. View at: Publisher Site  Google Scholar
 E. Hoek and M. S. Diederichs, “Empirical estimation of rock mass modulus,” International Journal of Rock Mechanics and Mining Sciences, vol. 43, no. 2, pp. 203–215, 2006. View at: Publisher Site  Google Scholar
 Itasca, FLAC3D “Fast Lagrangian Analysis of Continua in 3dimensions,” Version 5.0, Itasca Consulting Group, Inc., Itasca, Minnesota, 2012, Manual.
 C. Q. Zhang, H. Zhou, and X. T. Feng, “An Index for estimating the stability of brittle surrounding rock massFAI and its engineering application,” Rock Mechanics and Rock Engineering, vol. 44, no. 4, pp. 401–414, 2014. View at: Google Scholar
 O. Stephansson, T. Savilahti, and B. Bjarnason, “Rock mechanics of the deep borehole at Gravberg, Sweden,” August Aimé Balkema, vol. 2, pp. 863–870, 1989. View at: Google Scholar
 R. S. Read, “20 years of excavation response studies at AECL’s underground research laboratory,” International Journal of Rock Mechanics and Mining Sciences, vol. 41, no. 8, pp. 1251–1275, 2004. View at: Publisher Site  Google Scholar
 C. D. Martin and R. Christiansson, “Estimating the potential for spalling around a deep nuclear waste repository in crystalline rock,” International Journal of Rock Mechanics and Mining Sciences, vol. 46, no. 2, pp. 219–228, 2009. View at: Publisher Site  Google Scholar
 Q. Jiang, X.t. Feng, J. Chen, K. Huang, and Y.l. Jiang, “Estimating insitu rock stress from spalling veins: a case study,” Engineering Geology, vol. 152, no. 1, pp. 38–47, 2013. View at: Publisher Site  Google Scholar
 J. F. Sigrist, Numerical Simulation, an Art of Prediction 1, Wiley, New York, NY, USA, 2019.
 J. Li, “A review of techniques, advances and outstanding issues in numerical modelling for rock mechanics and rock engineering,” International Journal of Rock Mechanics and Mining Sciences, vol. 40, no. 3, pp. 283–353, 2003. View at: Google Scholar
 V. Hajiabdolmajid, Mobilization of Strength in Brittle Failure of Rock, Department of Mining Engineering, Queen’s University, Kingston, Canada, 2001.
 V. Hajiabdolmajid, P. K. Kaiser, and C. D. Martin, “Modelling brittle failure of rock,” International Journal of Rock Mechanics and Mining Sciences, vol. 39, no. 6, pp. 731–741, 2002. View at: Publisher Site  Google Scholar
 V. Hajiabdolmajid, C. D. Martin, and P. K. Kaiser, “Modelling brittle failure of rock,” in Proceedings of the Fourth North American Rock Mechanics Symposium, Seattle, WA, USA, August 2000. View at: Google Scholar
 J. W. Zhou, W. Y. Xu, M. W. Xu et al., “Application of rock strain softening model to numerical analysis of deep tunnel,” Chinese Journal of Rock Mechanics and Engineering, vol. 28, no. 6, pp. 1116–1127, 2009. View at: Google Scholar
 M. S. Diederichs, J. L. Carvalho, and T. G. Carter, “A modified approach for prediction of strength and post yield behaviour for high GSI rock masses in strong brittle ground,” in Proceedings of the 1st CanUS Rock Symposiium, Vancouver, Canada, June 2007. View at: Google Scholar
 T. G. Carter, M. S. Diederichs, and J. L. Carvalho, “Application of modified HoekBrown transition relationships for assessing strength and post yield behaviour at both ends of the rock competence scale,” Journal of the Southern African Institute of Mining & Metallurgy, vol. 108, pp. 325–338, 2008. View at: Google Scholar
 P. K. Kaiser, B. Kim, R. P. Bewick, and B. Valley, “Rock mass strength at depth and implications for pillar design,” Mining Technology, vol. 120, no. 3, pp. 170–179, 2011. View at: Publisher Site  Google Scholar
 M. H. A. Korenromp, Comparative Study of the SShaped and HoekBrown Failure Criterion by FiniteDifference Modelling, Department of Geoscience and Engineering, Delft University of Technology, Delft, Netherlands, 2003.
Copyright
Copyright © 2020 Zhiqiang Li et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.