Mathematical Problems in Engineering
Volume 2020 (2020), Article ID 9418905, 18 pages
https://doi.org/10.1155/2020/9418905

Research Article

Sensitivity Analysis in the Estimation of Mechanical Parameters of Engineering Rock Mass Based on the Hoek–Brown Criterion

Zhiqiang Li,¹ Guofeng Liu,¹ Shuqian Duan,² Shufeng Pei,³ and Changgen Yan¹

¹School of Highway, Chang’an University, Xi’an 710064, China
²School of Civil Engineering, Zhengzhou University, Zhengzhou 450001, China
³College of Geosciences and Engineering, North China University of Water Resources and Electric Power, Zhengzhou 450001, China

Correspondence should be addressed to Guofeng Liu; gfliu@chd.edu.cn and Shuqian Duan; shuqian.duan@enpc.fr

Received 21 March 2020; Revised 3 July 2020; Accepted 22 July 2020; Published 30 August 2020

Academic Editor: Vittorio Zampoli

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.

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 (si). Results show that the MPs of the rock mass are more sensitive to and and are high-sensitivity 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 large-scale 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 well-developed 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 by

where 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, mi is the H–B frictional constant m for intact rock, and mb is the H–B frictional constant m for rock mass. The parameter mb is a reduced value of mi 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 near-surface 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 mi. 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 well-established parameters Joint Condition (Jcond₈₉) and Rock Quality Designation (RQD) to reduce the difference of rock mass quality estimated result between different geologists working on a single project. For mi, 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 mi, 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 mi) 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 (si) [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. 3-1, No. 4, and RS411. Therein, the f20 steeply inclined fault zone is well-developed, with a width of about 30 cm. The width of T813 is about 3∼5 cm. The thickness of the WIZs of No. 3-1 and No. 4 is 20 cm and 40∼60 cm, respectively. These are highly fractured coarse-grained 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].

Figure 1: Geological profiles along the axis of the underground powerhouse.

The surrounding rock mass of the cavern is mainly composed of cryptocrystalline basalt, porphyritic basalt, almond-shaped 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. 3-1, 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 large-scale 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.

Figure 2: The distribution of fractures and local joints along the layer I of the powerhouse.

(a)

(b)

Figure 3: In situ layout schemes of test boreholes in the cavern from chainage K0+300 to K0+330. (a) The layout of test boreholes on the K0+320 section (M-1 to M-3 are multipoint extensometers, and R-K0+320-0/1-1/2 are the test boreholes of the digital camera). (b) The planar distribution of the test boreholes.

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 pre-excavated 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 R-K0+320-0-1 and R-K0+320-0-2 (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 R-K0+320-1-1 and R-K0+320-1-2 in Figure 3(a). Herein, the digital borehole camera system (see Figure 4(a)) that can capture a 360-degree 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+320-0-1 and K0+320-0-2 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 almond-shaped basalt and relatively broken due to the cut of well-developed joints.

(a)

(b)

(c)

(d)

(e)

Figure 4: Observation results of borehole camera. (a) Digital panoramic borehole camera system. (b) An example of the panoramic image of two meters of the wall of borehole K0+320-0-1. (c) One box of the drilling core samples from the borehole K0+320-0-1. (d) The generalized structure of rock mass surrounding borehole K0+320-0-1. (e) The generalized structure of rock mass surrounding borehole K0+320-0-2.

The value of GSI can be calculated by the equation proposed by Hoek et al. [15] as expressed in (3). The Jcond₈₉ 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:

Table 1: The detailed data involved in the calculation of GSI.

Monitoring hole
Discontinuity length
Separation
Roughness
Infilling
Weathering
JCond₈₉
RQD
GSI

K0+320-0-1
1
1
6
6
6
20
50
53

K0+320-0-2
1
1
5
6
5
18
46
50

Note: the guidelines for classification of discontinuity conditions proposed by Bieniawski [27] are used to rate the factors in Table 1.

2.2.2. Measurement of σci, mi, 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 RMT150-C 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 almond-shaped basalt is within 137∼157 MPa.

Figure 5: The distribution of of the almond-shaped basalt.

The constant mi 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 MTS815-04-type 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 stress-strain curves of almond-shaped 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 almond-shaped basalt is 26, consistent with the recommended value of 25 ± 5 for basalt by Marinos and Hoek [28].

Figure 6: Typical stress-strain curves of almond-shaped basalt under confining pressure with 20 MPa.

Table 2: Results of triaxial compression tests of almond-shaped basalt.

σ₃ (MPa)
σ₁ (MPa)
(σ₁–σ₃) (MPa)
E (GPa)
ν

0
150
150
26
0.17

5
131
126
17
0.19

10
178
168
23
0.25

20
183
163
22
0.21

40
458
418
35
0.24

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 single-hole 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.

Figure 7: Rock mass wave velocity tester.

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., , mi, 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 (0-1) is given. The constants mi 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 almond-shaped basalt can cover, were chosen as the fluctuation range for mi and , respectively.

Table 3: Reference values and fluctuation ranges of IPs at the monitoring section.

IPs
GSI
(MPa)
D
mi

Reference value
53
147
0.21
26

Fluctuation range
15–65
137–157
0-1
20–30

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 (si) [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 si, (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 si represents the multiple relationships between δy and δx:

where S is the sensitivity function, si is the sensitivity factor, and xi indicate different factors, δx is the error of factor, and yi 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, xi indicates different IP (GSI, D, mi, and ), and xmax and xmin 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, mi, 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 si 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.

Table 4: E corresponding to GSI and the si of E to GSI.

GSI
15
27
34
40
46
53
65

E (GPa)
1.19
2.38
3.56
5.03
7.11
10.64
21.22

k
0
0.24
0.38
0.50
0.62
0.76
1

si
0.86
1.56
1.96
2.31
2.65
3.05
3.74

(a)

(b)

(c)

Figure 8: Results of sensitivity analysis of MPs. (a) The sensitivity curve of E to IPs. (b) The sensitivity curve of c to IPs. (c) The sensitivity curve of φ to IPs.

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 mi, the value of si has almost no change, while the sensitivity of c to is higher than mi 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 , mi, and D is ranked as  > mi > D. In the range of “0.5<k < 0.6,” the sensitivity of c changes to be ranked as  > D > mi. In the range of “0.6<k < 1,” the sequence of the sensitivity changes as D >  > mi.

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 si is D. The sensitivities of parameter φ to and mi 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 decision-makers 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 strain-softening constitutive model is adopted to describe the mechanical behavior of rock mass and WIZs.

Table 5: The MPs used for numerical simulations.

Surrounding rock
E (GPa)
c (MPa)
φ (degree)

Rock mass
10.63
4.58
46.99

WIZ
1.55
2.2
26

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+320-1-1 and K0+320-1-2 by using the wave velocity test. It shows that the simulated depths (FAI > 1) of EDZ in the direction of K0+320-1-1 and K0+320-1-2 monitoring boreholes are 1.57 m and 2.84 m, respectively. The actual EDZ depths of K0+320-1-1 and K0+320-1-2 observation boreholes are 1.9 m and 2.5 m, respectively.

(a)

(b)

(c)

Figure 9: The simulation and observation result of EDZ on the section K0+320 of the cavern. (a) The FLAC model. (b) The simulation contour of EDZ. (c) Results of in situ acoustic velocity tests.

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 stress-structure controlled collapse in some positions (see Figures 10(c) and 10(d)).

(a)

(b)

(c)

(d)

Figure 10: The response behavior of rock mass around section K0+320 of the cavern after the excavation of layer I. (a) Contour of maximum redistributed stress. (b) Contour of EDZ. (c) Photo of rock cracking and failures at the spandrel on upstream side. (d) Photo of rock cracking and failures at the arch foot on downstream side.

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 stress-induced or stress-structure controlled rock failures always occur at the excavation profile parallel to or at a small angle with the orientation of maximum far-field 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.

Figure 11: Contour of the displacement of surrounding rock after the excavation of layer I of the cavern.

Moreover, the depths of EDZ along the boreholes K0+320-0-1 and K0+320-0-2 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.

Figure 12: Results of in situ wave velocity tests from the boreholes after the excavation of layer I.

Table 6: The comparison between the simulated and observed EDZ and displacement.

EDZ (m)
Displacement (mm)

RK-320-0-1
RK-320-0-2
M1-1
M1-2
M2-2
M2-3

Actual
3.4
1.8
29.5
28.9
23.1
20.5

Simulation
3.1
1.5
24.5
20.6
20.5
17.4

Note: the actual displacement of point M2-1 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 M1-1, M2-1, and M3-1 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)

Figure 13: Results of sensitivity analysis of displacement and EDZ. (a) The sensitivity curve of displacement (M3-1) to IPs. (b) The sensitivity curve of displacement (M1-1) to IPs. (c) The sensitivity curve of displacement (M2-1) to IPs. (d) The sensitivity curve of EDZ to IPs.

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 mi 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 si is D.

4. Discussion

4.1. Comparison between the Strain-Softening and S-Shaped 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 strain-softening (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 s-shaped. Thus, the S-shaped failure criterion, considered tensile fracturing or spalling at low confinement and extensile-driven shear at high confinement, is proposed [45–47]. The two segments of the criterion are connected by a transition zone where σ₃ = /10. Therefore, it is necessary for demonstrating the applicability of S-shaped criterion in the engineering presented in this study and illustrating the differences of the prediction results based on the S-shaped criterion and SS criterion.

According the in situ stress condition and of almond-shaped basalt, the rock mass around the K0+320 section of the cavern is at high confinement. The S-shaped criterion can be applied to FLAC3D on the basis of modified Hoek-Brown constitutive model and a user-defined FISH-function [48]. The parameters involved in the simulation are shown in Table 7.

Table 7: The parameters used in the simulation.

E (GPa)
(MPa)
mb
s
a
ρ (kg/m³)

10.64
147
13.75
0.1427
0.5005
2730

The displacements of the points M1-1 and M2-1 (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+320-0-1 and K0+320-0-2 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)

Figure 14: Comparison between the actual situations of excavation response and the prediction results simulated based on S-shaped criterion and SS criterion. (a) Displacement of surrounding rock at the location of extensometer M1-1 and M2-1. (b) The depth of EDZ along the direction of boreholes K0+320-0-1 and K0+320-0-2.

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 (si). The sensitivity factors of the estimated results to the IPs under the reference values are presented in Table 8.

Table 8: The sensitivity factors of the MPs and predicted results to the IPs.

Input parameters

GSI

D
mi

Elastic modulus (E)
3.05

0.12

Cohesion (c)
0.91
0.36
0.07
0.27

Internal friction angle (φ)
0.35
0.16
0.04
0.18

EDZ
0.20
0.23
0.01
0.01

Displacement
3.66
0.10
0.01
0.05

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 high-sensitivity 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 high-sensitivity IPs (such as ±5%, 10%, ±15%, ±20%, and ±25%) are shown in Figures 15 and 16, respectively.

(a)

(b)

(c)

Figure 15: The estimation results and variations of MPs corresponding to different given fluctuations of the IPs. (a) The variations in the estimation of parameter E. (b) The variations in the estimation of parameter c. (c) The variations in the estimation of parameter φ (10.65, 4.58, and 46.99 in Figures 15(a)–15(c) are the estimation results of the MPs at the reference value of the IPs).

(a)

(b)

Figure 16: The prediction results and variations of EDZ and displacement of the rock mass corresponding to different given fluctuations of the IPs. (a) The variations in the prediction of EDZ. (b) The variations in the prediction of displacement (5.6 and 15.99 in Figures 16(a) and 16(b) are the prediction results of EDZ and displacement at the reference value of the IPs).

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 criterion-based 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 high-sensitivity parameters that affect the displacement and EDZ of the surrounding rock. The sensitivity of displacement and EDZ to the IPs is ranked as GSI >  > mi > 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 (si) 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
mi: | 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
mb: | The H–B frictional constant (rock mass)
a: | The H–B exponent
Jcond₈₉: | Joint Condition
RQD: | Rock Quality Designation
Kv: | The integrality coefficient
S: | The sensitivity function
si: | 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
Vpm: | The longitudinal wave velocity in rock mass
Vpr: | 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
xi: | The different input parameters
xmax: | The maximum value of every input parameter in its fluctuation range
xmin: | 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: | Strain-softening.

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. 2019JQ-171), 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.

E. Hoek, “Strength of jointed rock masses,” Géotechnique, vol. 33, no. 3, pp. 187–223, 1983.

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.

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.

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.

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.

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.

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.

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.

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. Carranza-torres, and B. Corkum, “Hoek-Brown failure criterion—2002 edition,” in Proceedings of the NARMS-Tac, University of Toronto Press, Toronto, Canada, 2002.

S. Sinha and G. Walton, “A progressive S-shaped yield criterion and its application to rock pillar behavior,” International Journal of Rock Mechanics and Mining Sciences, vol. 105, pp. 98–109, 2018.

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.

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.

E. Hoek, Practical Rock Engineering, A. A. Balkema, Rotterdam, Netherlands, 2000.

M. Cai, “Practical estimates of tensile strength and Hoek-Brown strength parameter mi of brittle rocks,” Rock Mechanics and Rock Engineering, vol. 43, no. 2, pp. 167–184, 2010.

C. B. Yan and G. Y. Xu, “Modification of Hoek-Brown expressions and its application to engineering,” Chinese Journal of Rock Mechanics and Engineering, vol. 24, no. 22, pp. 4030–4035, 2005.

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.

T. G. Carter and R. I. Miller, “Crown-pillar risk assessment-planning aid for cost-effective mine closure remediation,” Transactions of the Institution of Mining and Metallurgy, vol. 104, pp. A41–A57, 1995.

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.

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.

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.

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.

G. B. Baecher and J. T. Christian, “Reliability and statistics in geotechnical engineering,” Engineering Geology, vol. 30, pp. 237-238, 2005.

E. T. Brown, “Risk assessment and management in underground rock engineering-an overview,” Journal of Rock Mechanics and Geotechnical Engineering, vol. 4, no. 3, pp. 193–204, 2012.

S. K. Sharan, “An iterative method for the linearization of nonlinear failure criteria,” Advances in Engineering Software, vol. 117, pp. 89–94, 2018.

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.

Itasca, FLAC3D “Fast Lagrangian Analysis of Continua in 3-dimensions,” 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 mass-FAI and its engineering application,” Rock Mechanics and Rock Engineering, vol. 44, no. 4, pp. 401–414, 2014.

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.

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.

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.

Q. Jiang, X.-t. Feng, J. Chen, K. Huang, and Y.-l. Jiang, “Estimating in-situ rock stress from spalling veins: a case study,” Engineering Geology, vol. 152, no. 1, pp. 38–47, 2013.

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.

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.

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.

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.

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 Can-US Rock Symposiium, Vancouver, Canada, June 2007.

T. G. Carter, M. S. Diederichs, and J. L. Carvalho, “Application of modified Hoek-Brown 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.

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.

M. H. A. Korenromp, Comparative Study of the S-Shaped and Hoek-Brown Failure Criterion by Finite-Difference Modelling, Department of Geoscience and Engineering, Delft University of Technology, Delft, Netherlands, 2003.