Jinping II hydropower station was constructed on the lower half of the Yalong River bend in Sichuan Province, China. It is one of the most important cascade hydropower stations on the mainstream of Yalong River (Figure 1A). To utilize the highly natural drop for power generation, the river bend with a length of 150 km is cut by several underground tunnels. Among them, the length of the diversion tunnel is about 17 km, and the buried depth of the overlying rock mass reaches 1,500–2000 m. The circular section designed with a diameter of 12.4 m was adopted for the construction of No. 1 and No. 3 diversion tunnels using a tunnel boring machine (TBM), and a horseshoe section designed with a diameter of 13.0 m was used for the construction of No. 2 and No. 4 diversion tunnels by employing blasting. To avoid construction interference, a 60 m parallel interval between each diversion tunnel was designed as shown in Figure 1A.

1) Stratigraphic Lithology

Along the diversion tunnels, the main lithology of the strata is marble, slate, sandstone, and limestone. Along the reverse water flow direction of the diversion tunnels, they are Yantang Formation (T_2y), Baishan Formation (T_2b), Triassic Upper Series (T₃), Zagunao Formation (T_2z), and Triassic Lower Series (T₁), respectively (Figure 1B). Among them, the Yantang Formation (T_2y) mainly distributes in the Dashuigou area, and the core of the Laozhuangzi anticline is composed of marble and argillaceous limestone. The marble of the Baishan Formation (T_2b) mainly distributes in the middle of the project area, which forms the main part of the Jinping Mountains. The Baishan marble (T_2b) is stable and compact, and the thickness of the whole layer is 750–2,270 m. The Upper Triassic (T₃) mainly distributes in the main watershed, and its lithology is sandstone and slate. The Lower Triassic (T₁) stratum mainly distributes in the western part of the tunnel area, and its lithology is complex, which is composed of biotite chlorite schist, metamorphic medium-fine sandstone with thin-bedded marble, gravel or banded marble, etc.

2) Geological structure

In the project area, the folds are extremely developed and complicated, and most of them are dense folds extending near the SN direction (NNE). Generally, these folds are composed of three belts: east, middle, and west. Through the geological survey, it is revealed that compression bedding and NE-trending thrust faults are the main structural planes in the project area. Based on different structural features and distribution orientations, these planes are classified into four tectonic groups: NNW, NNE, NW∼NWW, and NE∼NEE. The main faults that the diversion tunnel passed through are F₅, F₆, F₂₅, and F₂₇, which are shown in Figure 1B. The attitude of F₅ fault is N10°∼30°E, NW∠70°, with an affecting width of 5–10 m. The occurrence of the fault F₆ is N20°∼50°E, NW or SE∠60°∼87°with the bandwidth of 1–4.2 m. The attitude of the F₂₇ fault tends to the N30∼40°W direction, which is located in the middle of the Ganhaizi Formation and distributes in the Baishan Formation (T_2b). The attitude of the fault F₂₅ is N20°E, SE∠70°with a crushing bandwidth of 1–2 m.

3) In-situ stress fields

Because of the complex geological conditions and large burial depth, the in-situ stress at the construction area of four diversion tunnels is relatively high. Using the hydraulic fracturing method, the in-situ stress field at the construction area was measured (see Figure 1C). From Figure 1C, the maximum principal stress (σ _max) reaches 72 MPa with a dip angle of 6.45–75.4°. The middle principal stress (σ _mid) reaches 34 MPa with a dip angle of 25°and is approximately perpendicular to the tunnel axis. For the minimum principal stress (σ _min), it reaches about 29 MPa in the vertical direction with a dip angle of 65°.

In the tunnel area, the lithology of the rock mass is mostly pure and brittle marble with high strength. Due to the large burial depth, the rock mass is endowed with high in-situ stress. These two factors constitute the basic conditions for the occurrence of rockbursts. Subsequently, influenced by the blasting excavation induced disturbance, the rockburst occurred frequently (Figure 2).

To analyze the evolution mechanism of rockbursts induced by blasting excavation, the rockbursts that occurred at No. 2 diversion tunnel from stake number K10+600 to K17+000 were observed and recorded, as shown in Table 1. Otherwise, the rockburst intensity at the No.2 diversion tunnel and the rockburst frequency around the tunnel cross-section is given in Figure 3. According to Figure 3, there were 75 rockbursts during the blasting excavation of No. 2 diversion tunnel from stake number K10+600 to K17+000, and the majority of the rockbursts forms were exfoliation, lump flick, and peeling. From the field observation, the surface of the rockburst crater is fresh, showing the characteristics of "shallow pit", "deep pit" and "V-shape". The crater length of the rockburst at the No.2 diversion tunnel ranges from 1 m to 72 m, and the crater depth of the strongest rockburst reached 1.5∼3.0 m. According to the length and depth of the rockburst crater, the rockburst intensity is mainly divided into four grades: slight rockburst (I), medium rockburst (II), intense rockburst (III), and drastic rockburst (IV). Among these four grades, the different rockburst grades correspond to the different sizes of rock fragments. In general, the higher the rockburst grade means the deeper crater, the larger size of rock fragments, the farther the ejection distance of rock fragments caused by rockbursts, and the louder the sound produced by rockbursts.

Combined with Table 1, Figures 1B,C and Figure 3, the frequency and intensity of the rockbursts that took place near the stake number K11+000 were more drastic than those that took place near stake number K16+000, which indicates that rockbursts tend to occur around the tunnel cross-section with high in-situ stress and buried depth. Since No. 2 diversion tunnel was excavated by the upper and lower step method, almost all of the rockbursts happened at the upper section of the tunnel. Among the total frequency of the rockbursts that happened at the No.2 diversion tunnel, 35 rockbursts took place at the right spandrel, and 17 rockbursts occurred at the left spandrel. In addition, 12, 8, and 3 rockbursts happened at the right side, top and left side of the tunnel cross-sections, respectively.

From Figures 1B,C and Figure 3, we can also find that from stake number K16+000 to K11+000, the intensity, frequency, and rockburst grade enlarge with the increase of burial depth. Furthermore, it can also be seen from the figures that the burial depth is not the only influencing factor. There is a close correspondence between the rockburst frequency and the geological structure (syncline, anticline) along the tunnel axis, which indicates that the impact of geological structures on the in-situ stress field is another important factor that induces rockbursts. To facilitate the classification of rockbursts, we regarded the sporadic rockburst and continuous rockburst as the low-grade rockburst and high-grade rockburst, respectively. Generally, the intensity and grade of rockbursts are closely related to the level of in-situ stress.

When excavating No. 2 diversion tunnel by blasting, rockbursts happened frequently, which caused work stoppage more than once and increased the lining and support costs. As a man-made engineering disaster, rockbursts are not only harmful to the stability of the rock mass near the excavation area, they also threaten the security of workers and mechanical equipment near the occurrence area. Thus, it is crucial to evaluate the rockburst risk during excavation of a deep tunnel by blasting. High in-situ stress inevitably endows high strain energy to the rock mass, which may be released in some forms by the rock mass after blasting excavation. Many studies have reported that the rockburst is a mechanical failure phenomenon caused by the energy release of the rock mass (Wang and Park, 2001; He et al., 2010; Su, et al., 2017; Su, et al., 2018; Chen, et al., 2019). Thus, the key to revealing the evolution mechanism of rockbursts is to explore the energy release regularities that result from blasting excavation of a deep-buried tunnel.

Because rockbursts are the brittle failure induced behaviour of rock masses, the mechanical property, especially the post-peak mechanical characteristic of Jinping marble should be first studied to investigate the evolution mechanism of the rockburst that took place during the blasting excavation of No. 2 diversion tunnel. Previous research has reported that the stress-strain response and the yield form of the marble do not present as a single pattern under the action of compressive stress, but are closely associated with the value of confining pressure. For instance, in the case of low confining pressure, the stress-strain curve of the marble drops rapidly after reaching peak strength. This shows that the marble presents brittle characteristics like granite under low confining pressure. Accompanied with the continuous rise of confining pressure, it shows that the stress-strain curve slowly declines and the marble after the peak has a definite bearing capacity. While the confining pressure increases to a high level, the stress-strain curve of the marble after the peak does not decline and the residual strength of the marble still stays constant, presenting perfectly plastic characteristics. Therefore, the marble after the peak presents the brittle-ductile-plastic (BDP) transition characteristic during the rising process of confining pressure (Wawersik and Fairhurst, 1970). This particular mechanical property leads to the marble after the peak, which still has a high bearing capacity under the high confining pressure and stores a large amount of strain energy, which may be abruptly released and cause rockbursts under the disturbance of the tunnel excavation. Hence, studying and describing the BDP transition characteristics of Jinping marble is not only important for revealing the energy release regularities resulting from the transient stress unloading on the excavation boundary, but also are helpful for further understanding the evolution mechanism of the rockburst that happened during the construction of No. 2 diversion tunnel.

Aiming at investigating the special mechanical properties of Jinping marble, marble samples with a burial depth of 2000 m were drilled and obtained at the east end of the auxiliary tunnel for testing (Chu, 2009). Using the MTS pressure test equipment, a triaxial compression test was performed (see Figure 4). The test result in Figure 4 indicates that the Jinping marble sample shows significantly brittle characteristics after reaching peak strength under low confining pressure (2 MPa). While the confining pressure rises to the level of 2∼8 MPa, the residual strength of the marble sample drops a little and retains a high value after the peak. It shows that the marble sample starts presenting ductility features. As the confining pressure rises to a level of 40 MPa, the residual strength of the marble sample after the peak does not reduce and shows perfectly plastic characteristics.

To reveal the energy release regularity derived from the transient stress unloading during the tunnel excavation of the Jinping II hydropower project, first and foremost, the post-peak mechanical property of Jinping marble needs to be described. In the Hoek-Brown strength criterion, the mechanical parameters m _b , s, a, etc. can be changed with the increase of plastic strain ε 3 p after the material yielding. Therefore, the hardening and softening behaviour of Jinping marble after the peak can be described by the Hoek-Brown strength criterion. In this research, the Hoek-Brown strength criterion in FLAC3D was employed to simulate the brittle-ductility-plastic (BDP) transition characteristic of Jinping marble. The yield equation is (Cundall et al., 2003): F 0 = σ 1 − σ 3 − σ c i ( m b σ 3 σ 1 + s ) a = 0 (1) here, m _b , s, and a are strength parameters relating to the quality evaluation of the rock mass GSI, the material parameters of the rock mass m _i ; σ ₁, σ ₃, and σ _ci are the first principal stress, third principal stress, and uniaxial compressive strength of the material, respectively.

In the above strength criterion, the rule is assumed that the maximum plastic strain increment Δ ε 1 p and the minimum plastic strain increment Δ ε 3 p satisfy the following relationship: Δ ε 1 p = γ Δ ε 3 p (2) here, γ is the factor associated with the stress level, and its value is updated in each calculation step with the increment of plastic strain. According to the yield stress level, four kinds of flow rules can be achieved by the Hoek-Brown strength criterion.

1) Associated flow rule

The rule for associated flow: Δ ε i p = − γ a f ∂ F 0 ∂ σ i (3)

Substituting Eq. 1 and Eq. 2 into Eq. 3, the associated flow factor γ _af can be obtained by the following formula: γ a f = − 1 1 + a m b ( m b σ 3 σ 1 + s ) a − 1 (4)

The associated flow rule is employed for the description of yield characteristics under the low confining pressure condition.

2) Isovolumetric flow rule

In the case of the higher confining pressure ( σ 3 c v ), the loading process will not terminate, and the volume of rock mass remains invariant in the subsequent loading process. Therefore, for the case of high confining pressure, this rule can describe the yield characteristic. The mathematical expression of this rule is: γ i f = − 1 (5) here, γ _if is the isovolumetric flow factor.

3) Radial flow rule

The rule is written by: γ r f = σ 1 σ 3 (6) here, γ _rf is the radial flow factor.

The radial flow rule can describe the tensile failure of rock masses under the action of the tensile stress.

4) Combination flow rule

When the confining pressure ranges from 0 to σ 3 c v , the flow rule of the rock mass after yielding should be between the associated flow rule and the isovolumetric flow rule. Hence, a combination flow rule can be used to describe the yield characteristic: γ c f = 1 1 γ a f + ( 1 γ i f − 1 γ a f ) σ 3 σ 3 c v (7) here, γ _cf is the combination flow factor.

Only the above four kinds of flow rules are still insufficient to describe the BDP transition characteristic of Jinping marble. To solve this issue, a scaling factor μ associated with the minimum principal stress σ ₃ was proposed by Cundall et al. (2003) to simulate the softening and hardening characteristic of the rock mass after yielding. By adjusting the value of the factor μ, it can be used to describe the variation characteristic of the strength parameters m _b , s, a, etc. changing with the plastic strain under the different confining pressure conditions. Using the scaling factor μ, the brittle failure of granite was successfully described by Diederichs, on this basis, the "V-shape" failure was reproduced at the Canadian underground research laboratory (Diederichs, 2007). To describe the BDP transition characteristic of Jinping marble, the scaling factor μ and 8 parameters including 4 peak intensity parameters and 4 residual strength parameters should be first determined (Zhang et al., 2010). By triaxial compression test, the peak and residual strength parameters were acquired, which are given in Table 2 and Table 3. As for the determination of scaling factor μ, in this paper, the comparison method of the numerical simulation and the actual measured result of damage zones was adopted. For example, various conditions corresponding to different scaling factors μ were assumed, and then the in-situ stress unloading caused stress change path and damage zones of the surrounding rock mass were simulated using the assumed scaling factor μ. If the low stress distribution region and the extent of damage zones obtained by numerical simulation are consistent with the range of the low velocity distribution belt in the surrounding rock mass measured by the field acoustic wave test, it means that the right scaling factor μ has been found. Due to the troublesome process, here only two assumed conditions are introduced (Table 4).

Condition one: the rock mass presents the brittle characteristic at the low confining pressure (0–2 MPa). The mechanical property transition of ductile-plastic occurs with the confining pressure reaching the level of 16 MPa.

Condition two: the ductile-plastic transition happens when confining pressure reaches 25 MPa.

Adopting the parameters given in Table 2, Table 3 and Table 4, taking the blasting excavation of No.2 diversion tunnel as an instance, the stress change path of the rock mass at a certain position of the surrounding rock mass under the disturbance of transient stress unloading on the excavation boundary was simulated, as shown in Figure 5. During the numerical modelling, the non-reflecting boundary condition was adopted, the stress dynamic adjustment processes of the rock mass elements which have a 30°angle with the x-axis (hole 4#) and have 1.0, 1.5, 2.0, 2.5, 3.0 m distance from the excavation boundary were recorded (Figure 5). In Figure 5, the fitting curve of peak strength and residual strength is also drawn. Figure 5A also presents that when the confining pressure reaches 16 MPa, the fitting curve of peak strength intersects the fitting curve of residual strength for condition one. It indicates that the confining pressure of the ductility-plastic transition is 16 MPa. While for condition two, the confining pressure of the ductility-plastic transition is 25 MPa (Figure 5B).

In Figure 5, when the stress path curve intersects the fitting curve of peak strength, it means that the rock mass element is damaged, and the larger falling of the stress path curve means there is worse damage to the surrounding rock mass. For condition one, the stress paths of the surrounding rock mass elements with the distances of 1.0 and 1.5 m from the excavation boundary both have a significant drop after reaching the fitting curve of peak strength, but the stress path of the surrounding rock mass element with 2.0 m distance from the excavation boundary does not fall. Therefore, the damage depth of the surrounding rock mass has a 30°angle with the x-axis (hole 4#) can be determined as 2.0 m. While for condition two, the stress path curve of the surrounding rock mass element with a distance of 2.0 m from the excavation boundary only has a slight drop after reaching the fitting curve of peak strength. The stress path curve of the surrounding rock mass element with a 2.5 m distance from the excavation boundary does not reach the fitting curve of peak strength, it means that the extent of damage to the surrounding rock mass is within 2.0∼2.5 m. By further drawing the stress path curve of the rock mass element with a 2.0∼2.5 m distance from the excavation boundary, it was found that the damage depth was 2.3 m at this location.

More simply, a square model (1 × 1 × 1 m) was set up to simulate the rock sample, and the numerical stress-strain curve in several confining pressure conditions was performed (Figure 6). As can be seen from Figure 6, the confining pressure plays a critical role in the mechanical characteristic of the marble after the peak. When the confining pressure is low (0∼2 MPa), the stress-strain curve falls quickly after the peak and presents the obviously brittle characteristic like granite, and there is still a significant difference between the peak strength and the residual strength. However, this difference continually declines with the rise of the confining pressure. While the confining pressure exceeds 10 MPa, the stress-strain curve does not fall quickly after yielding and shows a significantly ductile characteristic. When the numerical rock sample is under high confining pressure, the stress-strain curve does not decline after the peak and shows an ideally plastic characteristic. For condition one, the confining pressure 16 MPa is the ductility-plastic transition threshold. While for condition two, the confining pressure 25 MPa is the ductility-plastic transition threshold. Both conditions one and two can describe the BDP transition property of the marble. To deduce which condition is right, a verification using the field acoustic wave test needs to be performed.

The field acoustic wave test result of the cross-section K11+105 is given in Figure 7. From Figure 7, the damage depth of the surrounding rock mass which has a 30°angle with the x-axis (hole 4#) is 2.0 m. It can be found that the test result is consistent with the numerical simulation result of condition one in Figure 5A. On this basis, using the parameters of condition one, the damage zone of the surrounding rock mass resulting from tunnel excavation was simulated (see Figure 7). In Figure 7, we can find that except for hole 2#, there is little difference between the damage extent of simulation results and the damage extent of field test results, which proves that the simulation method adopted in this paper is correct. Because the indoor test presented in Figure 4 is hard to reproduce the real loading and unloading stress path of the surrounding rock mass, it may lead to the numerical axial stress-strain curves in Figure 6 being inaccurate. This is why the axial stress-strain curves of simulation results in Figure 6 are different from that of indoor test results in Figure 4, but the simulation result of the damage zone is the same as the field test result.

A simple assumption was made that an infinite-long tunnel is excavated with a radius of R in the initial hydrostatic stress field P ₀. Then, the stress adjustment process of the surrounding rock mass due to the quasi-static unloading of in-situ stress can be calculated by Eq. 8: σ 1 = σ θ = P 0 ( 1 + R 2 r 2 ) σ 2 = P 0 σ 3 = σ r = P 0 ( 1 − R 2 r 2 ) } (8) here, r is the distance between the rock mass and the excavation centre, σ _θ and σ _r are the circumferential stress and radial stress, respectively, and σ ₁, σ ₂, σ ₃ are the principal stresses.

After the calculation of principal stresses, the density of elastic strain energy can be counted by Eq. 9 (Solecki and Conant, 2003): U = [ σ 1 2 + σ 2 2 + σ 3 2 − 2 υ ( σ 1 σ 2 + σ 2 σ 3 + σ 1 σ 3 ) ] 2 E 0 (9) here, υ is the Poisson’s ratio, E ₀ is the elastic modulus, and U is the elastic strain energy density of the rock mass at a certain location.

Before tunnel excavation (σ ₁=σ ₂=σ ₃=P ₀), the initial rock mass is intact and its strain energy density is counted by: U 0 = ( 3 − 6 υ ) P 0 2 2 E 0 (10) here, U ₀ is the initial strain energy density of rock masses before the tunnel excavation.

When the unloading process is completed, the strain energy density of the surrounding rock mass can be obtained by substituting Eq. 8 into Eq. 9: U A = P 0 2 [ 3 − 6 υ + 2 ( 1 + υ ) R 4 r 4 ] 2 E 0 (11) here, U _A is the strain energy density of the surrounding rock mass at a certain position after stress unloading.

Comparing Eq. 10 and Eq. 11, it is not difficult to find that after the in-situ stress unloading induced by tunnel excavation, the strain energy density of the surrounding rock mass increases.

It is generally acknowledged that as a brittle material, there is a limit for the strain energy storage capacity of the rock mass. Based on the hydrostatic compressive pressure field considered in this research, the energy storage limit of the rock mass U _c for the case of three-dimensional pressure can be calculated by Eq. 12 (Xie et al., 2005): U c = σ c 3 2 E 0 ( σ 1 − σ 3 ) (12)

Under the condition of the quasi-static unloading, U _c is acquired by substituting Eq. 8 into Eq. 12: U c = r 2 σ c 3 4 R 2 E 0 P 0 (13)

From Eq. 13 and Eq. 11, it can be found that the shorter distance to the excavation centre, the weaker the energy storage capacity of the rock mass, and the more the strain energy will accumulate in the rock mass. Consequently, when it is close to the excavation boundary, the strain energy aggregated in the surrounding rock mass can easily exceed the storage limit, which will inevitably lead to the energy release and cause damage, even failure in the surrounding rock mass.

When the blasting excavation method is employed, the stress unloading process on the excavation face is rapid and transient but not a quasi-static process (Lu et al., 2012; Li et al., 2014; Zhu et al., 2014; Fan et al., 2015). The transient stress unloading induced stress field in the surrounding rock mass can be calculated by employing inversion integral, contour integration, and Laplace transform (Fan et al., 2016): σ r ( r , t ) P 0 = { 1                                                                     ⇒ ( t < ( r − R ) C p ) 1 − R 2 2 t 0 r 2 ( t − ( r − R ) C p ) − 1 π t 0 ∫ 0 ∞ Q 1 ( η , t ) d η   ⇒ ( ( r − R ) C p ≤ t ≤ ( r − R ) C p + t 0 ) 1 − R 2 2 r 2 − 1 π t 0 ∫ 0 ∞ Q 2 ( η , t ) d η                       ⇒ ( t > ( r − R ) C p + t 0 ) (14) σ θ ( r , t ) P 0 = { 1                                                                     ⇒ ( t < ( r − R ) C p ) 1 + R 2 2 t 0 r 2 ( t − ( r − R ) C p ) − 1 t 0 π ∫ 0 ∞ Q 3 ( η , t ) d η   ⇒ ( ( r − R ) C p ≤ t ≤ ( r − R ) C p + t 0 ) 1 + R 2 2 r 2 − 1 t 0 π ∫ 0 ∞ Q 4 ( η , t ) d η                            ⇒ ( t > ( r − R ) C p + t 0 ) (15) here: Q 1 = ( A r + B r ) [ cos ( ( r − R ) C p η ) − cos ( t η ) ] ( E + F ) η 2 + ( C r − D r ) [ sin ( t η ) − sin ( ( r − R ) C p η ) ] ( E + F ) η 2 Q 2 = ( A r + B r ) [ cos ( ( t − t 0 ) η ) − cos ( t η ) ] ( E + F ) η 2 + ( C r − D r ) [ sin ( t η ) − sin ( ( t − t 0 ) η ) ] ( E + F ) η 2 Q 3 = ( A θ + B θ ) [ cos ( t η ) − cos ( ( r − R ) C p η ) ] ( E + F ) η 2 + ( C θ − D θ ) [ sin ( t η ) − sin ( ( r − R ) C p η ) ] ( E + F ) η 2 Q 4 = ( A θ + B θ ) [ cos ( t η ) − cos ( ( t − t 0 ) η ) ] ( E + F ) η 2 + ( C θ − D θ ) [ sin ( t η ) − sin ( ( t − t 0 ) η ) ] ( E + F ) η 2 A r = α 1 α 2 , B r = α 3 α 4 , C r = α 1 α 4 , D r = α 2 α 3 , E = α 2 2 , F = α 4 2 , A θ = α 5 α 2 , B θ = α 6 α 4 , C θ = α 5 α 4 , D θ = α 2 α 6 α 1 = ( 2 G 2 r J 1 − η C p J 0 ) ( η C p r ) , α 2 = ( 2 G 2 R J 1 − η C p J 0 ) ( η C p R ) α 3 = ( 2 G 2 r Y 1 − η C p Y 0 ) ( η C p r ) , α 4 = ( 2 G 2 R Y 1 − η C p Y 0 ) ( η C p R ) α 5 = ( − 2 G 2 r J 1 − ( 1 − 2 G 2 ) η C p J 0 ) ( η C p r ) , α 6 = ( − 2 G 2 r Y 1 − ( 1 − 2 G 2 ) η C p Y 0 ) ( η C p r ) here, C _p is the P-waves velocity, t and t ₀ is the time and unloading time, respectively, σ _θ (r, t) and σ _r (r, t) is the dynamically circumferential stress and radial stress of the surrounding rock mass induced by transient unloading, respectively, η is the integral path, G is the Lame constant, J ₀, J ₁, Y ₀, and Y ₁ are the Bessel functions.

After calculation, the dynamically circumferential stress and radial stress of the surrounding rock mass are drawn in Figure 8A. From Figure 8A, compared with quasi-static unloading, transient unloading will enlarge the effect of radial unloading and circumferential loading in the surrounding rock mass, and result in a larger difference between the first and third principal stress. Combining this dynamic adjustment process with Eq. 12, the transient stress unloading will reduce the energy storage limit and tends to result in energy release.

Substituting the dynamic stress field of Figure 8A into Eq. 9, the strain energy density of the surrounding rock mass U _w (r, t) under the transient stress unloading can be calculated. To show the variation process of strain energy intuitively, the ratios of the strain energy density U _w (r, t) and energy storage limit U _c (r, t) to the initial strain energy density U ₀ of the surrounding rock mass at 1.0R, 1.2R, 1.4R, and 2.0R are drawn in Figure 8B. As can be seen from Figure 8B, the strain energy experiences a dynamic variation process like declining at first, raising second and stabilizing at last. In fact, before the strain energy reaches the peak value of the elastic change curve as shown in Figure 8B, it will exceed the energy storage limit and result in energy release. Obviously, under elastic conditions, the subsequent curve after the strain energy density curve intersecting with the energy storage limit curve is not the actual energy adjustment process.

On the purpose of evaluating the rockburst intensity, Cook proposed a new index, energy release rate (ERR), which can be counted by Eq. 16: E R R = ( E s + E e ) V e (16) here, E _s is the total value of released energy; E _e and V _e are the total strain energy and volume of the excavated rock mass, respectively.

Since the energy release rate (ERR) cannot reflect the intensity and location of the energy release in the surrounding rock mass, it fails to identify the occurrence position of the rockburst. To achieve this goal, an index of local energy release rate (LERR) was proposed by Jiang et al. (2010) for describing the energy release resulting from the tunnel excavation. This index is acquired by counting the difference value between the strain energy density before the material failure and the strain energy density after the material failure: L E R R i = U i p e a k − U i t r o u g h (17) here, LERR _i is LERR of the rock mass element; U _ipeak is the peak value of strain energy density of the rock mass element before the brittle failure; U _itrough is the trough value of strain energy density of the rock mass element after the brittle failure; i is the number of the rock mass element; U _ipeak and U _itrough are counted by Eq. 9.

From Figure 8B, it can also be seen that the residual value of strain energy density may be higher than the trough value of strain energy density under the elastic condition after the brittle failure occurring. Therefore, Eq. 17 should be corrected: L E R R i = U i max − U i s (18) here, U _imax is the maximum value of strain energy density for the i-th rock mass element before the brittle failure happens; U _is is the final stable value or the residual value of strain energy density after the brittle failure.

The same model (Figure 9A) and parameters in section 3.2 were adopted to analyze the energy release process due to excavation. Before excavation, adopting the ideal and elastic model, the initial in-situ stress field and strain energy of rock masses was calculated. Through a certain calculation time step, the balance was obtained and the simulation result is shown in Figure 9B. During the numerical simulation of the energy release process, the non-reflecting boundary condition, i.e. the energy absorption boundary condition mentioned in Section 3.2 was adopted. Then, based on the parameters in Table 2, Table 3, and Table 4 and the BDP model, the energy accumulation of the surrounding rock mass due to tunnel excavation was simulated as shown in Figure 9C. From the simulation result presented in Figure 9C, it was found that due to the low storage limit of the surrounding rock mass around the excavation boundary, the accumulated strain energy easily exceeds the energy storage limit. Subsequently, energy release happens and triggers off the failure of the surrounding rock mass. After the energy release process is complete, the strain energy stabilizes at the residual state, as in Figure 9D.

To reveal the dynamic process of energy release resulting from transient stress unloading, the strain energy density of the surrounding rock mass elements with the distances of 0.5, 1.0, 1.5, and 2.0 m from the excavation boundary (Figure 9A) were recorded, which is drawn in Figure 10.

Figure 10 shows that the strain energy adjustment curve has two significant amplitude values: the peak value and the trough value. For the rock mass element which is 0.5 m away from the excavation boundary (Figure 10A), the peak value of the strain energy density is smaller than its initial value. While for the rock mass element with distances of 1.0 and 1.5 m from the excavation boundary, the peak value of the strain energy density is greater than the initial value (see Figures 10B,C). This illustrates that due to transient stress unloading, the strain energy of the distant rock mass will flow to the excavation boundary and lead to the accumulated energy more easily exceeding the energy storage limit for the rock mass, which is near the excavation boundary. Thus, when the strain energy accumulates in the surrounding rock mass, it cannot reach the peak value like Figure 8B under the ideally elastic conditions, it can even be smaller than the initial value of the strain energy before excavation unloading. Therefore, it is inaccurate to use the difference value of the peak value and trough value (Eq. 17) when calculating the local energy release rate.

Figure 10 indicates that for the surrounding rock mass elements which are 0.5, 1.0, and 1.5 m away from the excavation boundary, the strain energy experiences a variation process of reducing firstly, raising secondly, then decreasing and finally stabilizing. Similar to the elastic conditions, the first reduction of strain energy will not trigger failures in the surrounding rock mass, because it is induced by the elastic unloading wave. However, the second decline of the strain energy results from the factor that the strain energy absorbed in the surrounding rock mass, which is greater than the storage limit, which will inevitably lead to the brittle failure of the surrounding rock mass. A further observation from Figure 10, it is not difficult to find that with the decrease of the distance between the rock mass and the excavation boundary, the energy release duration is shorter, and the local energy release rate is greater. These characteristics mean that the shorter distance from the rock mass to the final excavation boundary, the more violent energy release occurs.

For the deep tunnel excavation, the energy release is the inherent factor of the dynamic failure of the surrounding rock mass such as rockbursts (Hodgson and Joughin, 1966; Toksöz and Kehrer, 1972; Singh, 1988; Mikhalyuk and Zakharov, 1997; Wang and Park, 2001; He et al., 2010). Therefore, establishing a certain association of the energy release and the rockburst is crucial to predicting and controling the rockburst induced by the deep tunnel excavation. Based on the concept of LERR, the average value of the local energy release rate LERR i ¯ is proposed, which can be written by: LERR i ¯ = ∑ i n L E R R i V f (19) here, V _f is the total volume of the failure elements in surrounding rock masses.

Based on the recorded rockbursts in Table 1, the average value of LERR i ¯ was calculated. Afterwards, the relationship between the average value of LERR i ¯ and rockburst intensity was set up as shown in Figure 11. Then, using the mechanical parameters in Table 2, Table 3, and Table 4 and the in-situ stress fields in Figure 1C, the failure extent and final energy distribution of the surrounding rock mass for different rockburst grades at several stake numbers were simulated (Figure 11), and the comparison of between the numerical simulation and the measured result for the stake number K11+006 (rockburst grade IV), 10+978 (rockburst grade III), 12+644 (rockburst grade II), and 15+288 (rockburst grade I) is given in Figure 11. From Figure 11, the simulation extent of the brittle failure is consistent with the crater depth of the rockburst that occurred at the same stake number in Table 1. The simulation location of the brittle failure is consistent with the rockburst occurrence location at the tunnel cross-section in Table 1. For instance, the simulation extent of the brittle failure at stake number K11+006 is 1.5∼3.0 m, mainly at the right side and spandrel. Similarly, the field observed rockburst at K11+006 has the same characteristic. This proves the accuracy of the numerical method employed in this paper.

Using the numerical simulation method proposed above, the average value of LERR i ¯ was acquired. Then, the rockburst intensity can be predicted by the relationship curve in Figure 11. In addition, the residual state of strain energy was also simulated and drawn in Figure 11 in the form of a cloud chart to determine the location of the rockburst. After the energy release, the lowest residual value of the strain energy means the largest energy release rate and the most drastic rockburst occurred in the surrounding rock mass. Therefore, both the intensity and location of the rockburst can be predicted by employing the relationship curve in Figure 11 and the cloud chart of the residual strain energy.

The rockburst prediction is challenging work for researchers and engineers. Since the energy release plays a significant role in the formation of the rockburst, it has been the subject of much attention. The association between the rockburst intensity and the energy release rate (ERR) was developed by South African researchers, which can be expressed as: E R R = { 16 ∼ 34        ⇒ Slight rockburst  ( I ) 34 ∼ 53       ⇒ Medium rockburst  ( II ) 53 ∼ 100      ⇒ Intense rockburst  ( III ) 100 ∼ 140    ⇒ Drastic rockburst  ( IV )   ( M J / m 3 ) (20)

The magnitude in Eq. 20 is significantly different from that in Figure 11. Considering a simple case, under the burial depth of 1,000∼2000 m, an underground tunnel is excavated with the in-situ stress of 27∼54 MPa (assumption: υ=0.25, E ₀=20∼50 GPa). On this basis, the initial strain energy density of the rock mass can be estimated by Eq. 10. After calculation, the initial strain energy density of the surrounding rock mass only reaches 10.9–109.4 kJ/m³, which is similar to the magnitude in Figure 11. From this point of view, the prediction method of rockbursts proposed in Figure 11 is better than that in Eq. 20.

The fault, structural plane, and excavation method were not considered in this paper but they also play a key role in the energy release and rockburst characteristics. Thus, the prediction method proposed for the rockburst (Figure 11) still requires verification. Despite these limitations, this research provides new insights to reveal energy variation regularity and assess the rockburst risk during the excavation of the deep tunnel by blasting.