Impact of Shear Zone on Rockburst in the Deep Neelum-Jehlum Hydropower Tunnel: A Numerical Modeling Approach

Rockburst is a hazardous phenomenon in deep tunnels influenced by geological structural planes like faults, joints, and shear planes. Small-scale shear-plane-like structures have damaging impact on the boundaries of the tunnel, which act as barrier and accumulate high stresses. A shear plane combined with high stress conditions is very dangerous in deep excavations. Such a shear plane exposed in the side wall of the right headrace tunnel in the Neelum-Jehlum Hydropower Project. This project is constructed in the tectonically active Himalayas under high stress conditions. The influence of a shear zone on rockburst occurrence near the tunnel is studied. The FLAC3D explicit code simulated the shear zone in the right tunnel, revealing that the stresses are concentrated near the shear zone, while no such stress concentration is present in the left tunnel. The Rock mass got displaced near this shear zone. Modeling results confirm that the presence of shear zone in side wall of the right tunnel has a major influence on rockburst occurrence. A shear slip along this plane released huge amounts of accumulated energy, causing human fatalities and property damage. A comparison of numerical simulation with empirical rockburst criteria validates the actual conditions and help us to understand the phenomenon of stress concentration near the shear zone and its impact during deep tunneling.


Introduction
Rockburst is one of the most hazardous phenomena which intensify during construction of deeper excavations. It can be caused by either brittle failure or slippage along structural planes in deep hard rock [1]. Rockburst results in a sudden or violent damage to an excavation, usually associated with seismicity due to huge energy accumulation [2]. It not only causes damage to equipment and machinery, but also to nearby structures. A lot of research work has been done to find the mechanism of rockburst in deep mining [3][4][5][6][7]. However, this phenomenon is not much studied in civil tunneling for hydropower and road projects. A few hard rock tunnels have experienced intensive rockbursts worldwide. These include the Jinping hydropower tunnels in China, the Neelum-Jhelum Hydropower Project (NJHPP) twin headrace tunnels (right and left) in Pakistan, and the Gotthard base tunnel in Europe [8][9][10][11]. Kaiser and Cai [12] identified four different factors that cause rockburst failure: geotechnical, geological, mining, and seismic. Different research works have been done to find the effect of these factors on rockburst damage [13][14][15]. This study focuses only the geological factor, which includes geological studies of the area and the influence of geological structures i.e., shear zone on the occurrence of rockburst.
The magnitude of stress on an opening usually increases with depth and geological structures have a strong impact on the orientation of in situ stress [16]. Shepherd et al. [17] have effectively reviewed mechanism of rockburst in presence of different geological structures. Usually, during deep civil tunneling, geological weaknesses like faults or shear zones are found nearby, which causes unfavorable and complex stress scenario. These stresses are usually concentrated near these structures, which results in large deformations and failure processes because rock masses with different discontinuities fail when the shear strength of the discontinuities is less than shear stress [18]. Hedley [19] worked on Ontario (Canada) hard rock mines and found that the underground mines are prone to bursting when geological discontinuities are nearby. Ortlepp and Stacey [20] explained five different mechanisms behind rockburst in deep mines: strain burst, virgin shear, face crushing, buckling and reactivated shear in prevailing faults. Rockburst occurs due to shear failure along weak planes (pre-existing faults or shear zones) within the rock mass. Durrheim et al. [21] found regional structures as a source mechanism during the investigation of 21 rockburst cases in deep South African gold mines. Ortlepp [22] has explained rockburst induced by fault slip or shear rupture in massive rock mass. The presence of faults and joints were the main reason behind rockburst events in Jinping II pilot tunnels in China [23]. Snelling et al. [24] worked at the Creighton Mine (Canada) and found that geological structures resembling shear zones experiences frequent micro-seismic and often macro-seismic events. More than 20 rockburst events associated with structural planes have been listed during the construction of the Jinping-II headrace tunnels [25]. These structure planes become barriers to redistribution of stresses in surrounding rock mass in deep tunnels. This leads to concentration of tangential stress between the excavation boundary and structural plane along with huge amount of accumulated energy [14]. Laboratory test were conducted to explain about structure type rockburst for three possible rockburst mechanisms. It was found that small-scale structural planes were main controlling factors for rockburst [10]. Manouchehrian and Cai [26] used the Abaqus2d explicit code to explain the role of discontinuities around tunnels during rockburst and confirmed that more violent rockburst occurs in the presence of geological structures.
Civil engineering tunnels are influenced by small structure planes, whereas deep mines have faults that usually are tens or hundred meters long. These small-scale structures can easily be reactivated when they are present in the vicinity of a tunnel which are usually ignored during design. Mechanism of intense rockburst induced by a nearby small-scale structures/shear zones in deep civil tunnels is still unclear and which needs further studies. A rockburst surrounding a structural plane is more destructive than one that is not surrounded by any structural plane. This plane played a very important role for the intense rockburst during the NJHPP. In this paper, the most intense rockburst of 31 May 2015, which occurred in the right headrace tunnel of NJHPP is studied. It is believed that a shear plane exposed in the wall of the headrace tunnel was a favorable condition for a rockburst to trigger. The FLAC 3D numerical simulation is used to investigate the mechanism of shear-zone-type rockburst and also compared the results with and without shear zone in twin headrace tunnels and evaluated empirically.

Project Description
NJHPP is being constructed in the north-east area of Pakistan. The construction started in 2008. The 28.5 km long headrace tunnel is constructed to divert the water from the Neelum River to an underground power house to produce 969 Megawatts (MW) of electric power. The project is in the mighty mountains of the Himalayas which are considered as hub of high stresses which causes rockburst in hard rocks and squeezing in weak rock. In addition, the presence of mixed geology and structures make worse possible conditions for tunneling. Owing to tectonics, the area is highly folded and faulted. A major fault in the region, the so-called Main Boundary Thrust (MBT) passes nearby the project area. The project site is present in an earthquake active zone and the 2005 Muzaffarabad earthquake occurred in this area and caused more than 75,000 fatalities.
The excavation of headrace tunnels in project are done with the drill and blast method and with Tunnel Boring Machine (TBM) as a single tunnel and twin tunnels, respectively. The overburden was up to 2000 m, which has the potential of rock bursting in hard rock. The TBM-excavated twin tunnels had an initial center to center distance of 33 m before the major rockburst event of 31 May 2015. This distance was increase to 66 m to further avoid any inter-tunnel pillar bursting.

Geological Settings of the Area and In-Situ Stresses
The northward movement of the Indian plate in the Himalayan orogenic system has caused several thrust faults in the region, such as the "Muzaffarabad fault" due to high horizontal stress built up [27]. The headrace tunnels pass through the areas having adverse folding and faulting with a series of anticlines and synclines structures, along with local faults and shear zones as shown in Figure 1. The bedding planes are normally perpendicular to the tunnel direction, which is a favorable condition for construction [28]. However, abnormal stress concentration caused the bedding planes to be oriented along the tunnel direction, posing a great threat to tunnel stability [29]. The entire project is excavated in the Murree formation of sedimentary origin with a stratigraphic sequence of alternative beds. These beds comprise sandstone, siltstone, and mudstone. Sandstone is the strongest rock unit, having an average UCS (uniaxial compressive strength) of 86 MPa with a stiffness of 32 GPa. This unit has thick beds which are massive and blocky, that's why most of the rockburst events occurred in this rock unit. Siltstone has medium strength and is sometimes mixed with mudstone and shale. The average UCS of Siltstone is 66 MPa with stiffness of 23.1 GPa. Mudstone is the weakest rock unit of the Murree formation and has an average UCS of 42 MPa with stiffness of 12.6 GPa. This unit has much less potential for rockburst because of its very low strength. Table 1 summarizes the engineering properties of these rock units. had an initial center to center distance of 33 m before the major rockburst event of 31 May 2015. This distance was increase to 66 m to further avoid any inter-tunnel pillar bursting.

Geological Settings of the Area and In-Situ Stresses
The northward movement of the Indian plate in the Himalayan orogenic system has caused several thrust faults in the region, such as the "Muzaffarabad fault" due to high horizontal stress built up [27]. The headrace tunnels pass through the areas having adverse folding and faulting with a series of anticlines and synclines structures, along with local faults and shear zones as shown in Figure 1. The bedding planes are normally perpendicular to the tunnel direction, which is a favorable condition for construction [28]. However, abnormal stress concentration caused the bedding planes to be oriented along the tunnel direction, posing a great threat to tunnel stability [29]. The entire project is excavated in the Murree formation of sedimentary origin with a stratigraphic sequence of alternative beds. These beds comprise sandstone, siltstone, and mudstone. Sandstone is the strongest rock unit, having an average UCS (uniaxial compressive strength) of 86 MPa with a stiffness of 32 GPa. This unit has thick beds which are massive and blocky, that's why most of the rockburst events occurred in this rock unit. Siltstone has medium strength and is sometimes mixed with mudstone and shale. The average UCS of Siltstone is 66 MPa with stiffness of 23.1 GPa. Mudstone is the weakest rock unit of the Murree formation and has an average UCS of 42 MPa with stiffness of 12.6 GPa. This unit has much less potential for rockburst because of its very low strength. Table 1 summarizes the engineering properties of these rock units.  Previous studies near the project area have witnessed high horizontal stresses due the active Himalayas [30]. In-situ stress measurements performed by over-coring in sandstone beds have also reflected high horizontal stresses in TBM tunnels. The ratio of the horizontal to vertical stress (K0) was up to 2.9 as shown in Figure 2, where the major principal stress is oriented sub-horizontally and Previous studies near the project area have witnessed high horizontal stresses due the active Himalayas [30]. In-situ stress measurements performed by over-coring in sandstone beds have also reflected high horizontal stresses in TBM tunnels. The ratio of the horizontal to vertical stress (K 0 ) was up to 2.9 as shown in Figure 2, where the major principal stress is oriented sub-horizontally and nearly perpendicular to the tunnel azimuth. The dashed curves represent hyperbolic relation between the range of extreme limits of K 0 and z (depth) [31]. The yellow dots show high overburden stresses within the limits and red dots show high abnormal stress. This high value of K 0 is due high horizontal tectonic stress which caused frequent slabbing and rock bursting. nearly perpendicular to the tunnel azimuth. The dashed curves represent hyperbolic relation between the range of extreme limits of K0 and z (depth) [31]. The yellow dots show high overburden stresses within the limits and red dots show high abnormal stress. This high value of K0 is due high horizontal tectonic stress which caused frequent slabbing and rock bursting. Therefore, hard sedimentary rock under high stresses due to tectonics along structural planes, make the condition favorable for rockburst.

Impact of Structural Planes on Rockburst in Hard Rock Tunnels
During excavation of deep tunnels, it is inevitable to encounter structural planes and most rockburst events are influenced by these planes. They can be small scale shear zones, which are usually, encountered during the excavation of hard rock tunnels in civil projects as compared to deep underground mining faces, which are mostly intercepted by large faults. Rockburst effected by these structures generally have higher intensities and produce more damage as compared to the ones that are not affected. These planes are sometimes visible on the boundary of the tunnel and sometimes present in the rock mass away from the boundary. These planes also control the boundaries of the damaged pits during rockburst [32]. The orientation of structural planes with tunnel axis and maximum principal stress is one of the important factors for rockburst intensity which can be categorized in the following different scenarios [25]. First, when the structural plane is at a large angle with the principal tangential stress and axis of tunnel, a rockburst with high intensity and deep damage pit occurs. Second, a structural plane at a tunnel wall with a large angle with the axis of tunnel along with small angle with the maximum tangential stress causes slight to moderate rockburst with moderate depth damage pit. Third, when it is parallel to the tunnel axis it causes high intensity rockburst with steep ridges and deep damage pit. Fourth, when the structural plane had a lesser angle with the axis of tunnel and passing near the tunnel boundary, causes intense rockburst with deep damage pit.
During the construction of the Jinping-II hydropower tunnel, presence of structural plane near the tunnel boundary, caused an extremely intense rockburst with huge amounts of released strain energy. Such rock bursts produce deep damage pits in the hanging wall area of the structural plane. In this type of rockburst, as the tunnel approaches the plane, large amount of concentrated energy is released, which finally resulted in the deep damage pit. Presence of rigid structural planes, subparallel to the axis of a deep hydropower tunnel, was studied by [10,33]. The huge amount of energy is accumulated and released, which resulted in very high intensity burst on 28 November 2009, hereafter referred to as the "11.28" event of rockburst in the drainage tunnel of the Jinping-II hydropower project and schematic diagram of this structural plane is shown in Figure 3a. Therefore, hard sedimentary rock under high stresses due to tectonics along structural planes, make the condition favorable for rockburst.

Impact of Structural Planes on Rockburst in Hard Rock Tunnels
During excavation of deep tunnels, it is inevitable to encounter structural planes and most rockburst events are influenced by these planes. They can be small scale shear zones, which are usually, encountered during the excavation of hard rock tunnels in civil projects as compared to deep underground mining faces, which are mostly intercepted by large faults. Rockburst effected by these structures generally have higher intensities and produce more damage as compared to the ones that are not affected. These planes are sometimes visible on the boundary of the tunnel and sometimes present in the rock mass away from the boundary. These planes also control the boundaries of the damaged pits during rockburst [32]. The orientation of structural planes with tunnel axis and maximum principal stress is one of the important factors for rockburst intensity which can be categorized in the following different scenarios [25]. First, when the structural plane is at a large angle with the principal tangential stress and axis of tunnel, a rockburst with high intensity and deep damage pit occurs. Second, a structural plane at a tunnel wall with a large angle with the axis of tunnel along with small angle with the maximum tangential stress causes slight to moderate rockburst with moderate depth damage pit. Third, when it is parallel to the tunnel axis it causes high intensity rockburst with steep ridges and deep damage pit. Fourth, when the structural plane had a lesser angle with the axis of tunnel and passing near the tunnel boundary, causes intense rockburst with deep damage pit.
During the construction of the Jinping-II hydropower tunnel, presence of structural plane near the tunnel boundary, caused an extremely intense rockburst with huge amounts of released strain energy. Such rock bursts produce deep damage pits in the hanging wall area of the structural plane. In this type of rockburst, as the tunnel approaches the plane, large amount of concentrated energy is released, which finally resulted in the deep damage pit. Presence of rigid structural planes, sub-parallel to the axis of a deep hydropower tunnel, was studied by [10,33]. The huge amount of energy is accumulated and released, which resulted in very high intensity burst on 28 November 2009, hereafter referred to as the "11.28" event of rockburst in the drainage tunnel of the Jinping-II hydropower project and schematic diagram of this structural plane is shown in Figure 3a.

Impact of the Shear Zone on the Rockburst at Neelum-Jhelum Hydropower Project (NJHPP)
In deep tunnels, rockburst is closely related to the encountered geology and rock mechanics. During excavation, geological and geo-mechanical failures are very critical. Information, either in the form of in-situ stress measurements or geological mapping is very important for revealing failure uncertainties. The shear zone can be present at different locations like on the working face, near the vault, or on the side walls of the tunnel which leads to a concentration of abnormal stress near these planes. In NJHPP, over-coring program was carried out for in situ stress measurement. Some values of K0 were out of range compared to the universally accepted depth vs K0 plot, as shown in Figure 2. According to the geological mapping of the headrace tunnel, it passes many minor faults and shear planes. Two local faults, F2 and F3, and one shear zone is expected to pass the invert level at chainage 10 + 360 m and 9 + 350 m, respectively as shown in Figure 1.
The alternate beds of sandstone, siltstone, and mudstone experience high folding and faulting under high tectonic stresses. A prominent geological anomaly encountered during the tunnel excavation. The behavior of the hard sedimentary rock of the Murree formation is further confirmed from exposure of a local shear plane (shown in Figure 3b) in the drag-folded strata shown in Figure  4 near the right headrace tunnel boundary. The over-coring, in-situ stress measurement program, also confirmed a high horizontal stress regime in this chainage range. Thus, this complex geological region resulted in unpredictable behavior of the excavation caused a series of strain burst failures due to stored strain energy before the main shear zone slip rockburst event. On 31 May 2015 at 11:35 p.m., an intense rockburst occurred with a great sound from chainage 09 + 706 to 09 + 793 at a depth of 1300 m. The maximum impact of this rockburst is received near chainage 9 + 751 received as shown in Figure 4.

Impact of the Shear Zone on the Rockburst at Neelum-Jhelum Hydropower Project (NJHPP)
In deep tunnels, rockburst is closely related to the encountered geology and rock mechanics. During excavation, geological and geo-mechanical failures are very critical. Information, either in the form of in-situ stress measurements or geological mapping is very important for revealing failure uncertainties. The shear zone can be present at different locations like on the working face, near the vault, or on the side walls of the tunnel which leads to a concentration of abnormal stress near these planes. In NJHPP, over-coring program was carried out for in situ stress measurement. Some values of K 0 were out of range compared to the universally accepted depth vs K 0 plot, as shown in Figure 2. According to the geological mapping of the headrace tunnel, it passes many minor faults and shear planes. Two local faults, F2 and F3, and one shear zone is expected to pass the invert level at chainage 10 + 360 m and 9 + 350 m, respectively as shown in Figure 1.
The alternate beds of sandstone, siltstone, and mudstone experience high folding and faulting under high tectonic stresses. A prominent geological anomaly encountered during the tunnel excavation. The behavior of the hard sedimentary rock of the Murree formation is further confirmed from exposure of a local shear plane (shown in Figure 3b) in the drag-folded strata shown in Figure 4 near the right headrace tunnel boundary. The over-coring, in-situ stress measurement program, also confirmed a high horizontal stress regime in this chainage range. Thus, this complex geological region resulted in unpredictable behavior of the excavation caused a series of strain burst failures due to stored strain energy before the main shear zone slip rockburst event. On 31 May 2015 at 11:35 p.m., an intense rockburst occurred with a great sound from chainage 09 + 706 to 09 + 793 at a depth of 1300 m. The maximum impact of this rockburst is received near chainage 9 + 751 received as shown in Figure 4.

Impact of the Shear Zone on the Rockburst at Neelum-Jhelum Hydropower Project (NJHPP)
In deep tunnels, rockburst is closely related to the encountered geology and rock mechanics. During excavation, geological and geo-mechanical failures are very critical. Information, either in the form of in-situ stress measurements or geological mapping is very important for revealing failure uncertainties. The shear zone can be present at different locations like on the working face, near the vault, or on the side walls of the tunnel which leads to a concentration of abnormal stress near these planes. In NJHPP, over-coring program was carried out for in situ stress measurement. Some values of K0 were out of range compared to the universally accepted depth vs K0 plot, as shown in Figure 2. According to the geological mapping of the headrace tunnel, it passes many minor faults and shear planes. Two local faults, F2 and F3, and one shear zone is expected to pass the invert level at chainage 10 + 360 m and 9 + 350 m, respectively as shown in Figure 1.
The alternate beds of sandstone, siltstone, and mudstone experience high folding and faulting under high tectonic stresses. A prominent geological anomaly encountered during the tunnel excavation. The behavior of the hard sedimentary rock of the Murree formation is further confirmed from exposure of a local shear plane (shown in Figure 3b) in the drag-folded strata shown in Figure  4 near the right headrace tunnel boundary. The over-coring, in-situ stress measurement program, also confirmed a high horizontal stress regime in this chainage range. Thus, this complex geological region resulted in unpredictable behavior of the excavation caused a series of strain burst failures due to stored strain energy before the main shear zone slip rockburst event. On 31 May 2015 at 11:35 p.m., an intense rockburst occurred with a great sound from chainage 09 + 706 to 09 + 793 at a depth of 1300 m. The maximum impact of this rockburst is received near chainage 9 + 751 received as shown in Figure 4.  Three workers died, and a few others injured in the 31 May 2015 event. Additionally, a large area of crown and walls fell as shown in Figure 5. Owing to the instantaneous impact of the shock wave, support system was destroyed, wire-mesh and ring beams were deformed, and TBM was destroyed. The severe damage indicates the concentration of stresses in this area due to the presence of geological structure. After this major event, the rockburst shockwave caused damage to the left wall and crown of the nearby headrace tunnel (left tunnel) already excavated by TBM 697. Three workers died, and a few others injured in the 31 May 2015 event. Additionally, a large area of crown and walls fell as shown in Figure 5. Owing to the instantaneous impact of the shock wave, support system was destroyed, wire-mesh and ring beams were deformed, and TBM was destroyed. The severe damage indicates the concentration of stresses in this area due to the presence of geological structure. After this major event, the rockburst shockwave caused damage to the left wall and crown of the nearby headrace tunnel (left tunnel) already excavated by TBM 697. After this event, a shear plane was exposed on the boundary of the tunnel, shown in Figure 3b. It is concluded here that the series of small rockburst events before the main event paved the way for the final displacement along the exposed shear zone. Those small events were due to brittle failure of hard rock. Asperities along shear zone were removed due to high stress environment, which finally caused slippage along this plane and resulted in the most hazardous rockburst event on 31 May 2015, releasing a huge amount of energy. The Fault-slip rockburst of similar occurrence was also witnessed in the headrace tunnel of the Tianshengqiao II hydropower station [34] and Jinping II drainage tunnels [33].

Numerical Simulation of Shear Zone
Rockburst involves accumulation of huge amounts of energy in the presence of shear zone. To assess the impact of the shear zone on NJHPP tunnel stability and its influence on rockburst failure mechanism around the tunnels, FLAC 3D numerical simulation is performed. From the known massive hard rock geology, it was obvious to choose the FLAC 3D continuum model to study the influence of the shear zone on tunnel stability. The field observations revealed a shear zone exposed after intense rockburst. This shear zone was nearly perpendicular to the right tunnel axis [29]. To find the reason behind the intense rockburst event, a shear zone was simulated perpendicular to the right tunnel axis. The model size was 127 m × 100 m × 100 m, as shown in Figure 6. The two blocks of 50 m each were created and named Front and Back block. The shear zone was modelled as interface zone between the two blocks. Elasto-Plastic material model with Mohr-Coulomb criteria was used for the analysis. A model with hard sandstone rock unit parameters was simulated to check the effect of shear zone on rock mass and tunnel.
The twin tunnels, each with diameter 8.53 m, were simulated according to the actual construction sequence. The center to center distance between these twin tunnels was 33 m same as discussed in project description. Three different types of mesh primitive shapes (radtunnel, radcylinder, brick) were connected together. The model has 113,600 zones. The Fine mesh was simulated around the twin tunnels to get better results near the boundary of excavation. After simulating the geometric model, reasonable boundary conditions were applied, all sides of the model were fixed, except the top side to apply vertical stresses (σzz). After this event, a shear plane was exposed on the boundary of the tunnel, shown in Figure 3b. It is concluded here that the series of small rockburst events before the main event paved the way for the final displacement along the exposed shear zone. Those small events were due to brittle failure of hard rock. Asperities along shear zone were removed due to high stress environment, which finally caused slippage along this plane and resulted in the most hazardous rockburst event on 31 May 2015, releasing a huge amount of energy. The Fault-slip rockburst of similar occurrence was also witnessed in the headrace tunnel of the Tianshengqiao II hydropower station [34] and Jinping II drainage tunnels [33].

Numerical Simulation of Shear Zone
Rockburst involves accumulation of huge amounts of energy in the presence of shear zone. To assess the impact of the shear zone on NJHPP tunnel stability and its influence on rockburst failure mechanism around the tunnels, FLAC 3D numerical simulation is performed. From the known massive hard rock geology, it was obvious to choose the FLAC 3D continuum model to study the influence of the shear zone on tunnel stability. The field observations revealed a shear zone exposed after intense rockburst. This shear zone was nearly perpendicular to the right tunnel axis [29]. To find the reason behind the intense rockburst event, a shear zone was simulated perpendicular to the right tunnel axis. The model size was 127 m × 100 m × 100 m, as shown in Figure 6. The two blocks of 50 m each were created and named Front and Back block. The shear zone was modelled as interface zone between the two blocks. Elasto-Plastic material model with Mohr-Coulomb criteria was used for the analysis. A model with hard sandstone rock unit parameters was simulated to check the effect of shear zone on rock mass and tunnel.
The twin tunnels, each with diameter 8.53 m, were simulated according to the actual construction sequence. The center to center distance between these twin tunnels was 33 m same as discussed in project description. Three different types of mesh primitive shapes (radtunnel, radcylinder, brick) were connected together. The model has 113,600 zones. The Fine mesh was simulated around the twin tunnels to get better results near the boundary of excavation. After simulating the geometric model, reasonable boundary conditions were applied, all sides of the model were fixed, except the top side to apply vertical stresses (σ zz ). Energies 2018, 11, x FOR PEER REVIEW 7 of 16 Figure 6. Rock mass model.

Interface Element
FLAC 3D provides an interface element to represent the shear zone type structures in a geologic medium. This element is a collection of triangular elements. It can be created at any location in the model. In the model, the interface element is attached to a front block at 50 m length in the right headrace tunnel to reflect the actual field conditions. After creating the interface element, the space between the front and back block was removed and back block was moved 1 m back to get attach with the front block. Both blocks were attached without attaching the shear zone area. Interface element was characterized by Coulomb sliding. Different parameters like normal stiffness (kn) and shear stiffness (ks), friction angle (°), cohesion (MPa), were assigned to interface element. FLAC [35] explained the stiffness characteristics in an explicit way. So stiffness properties of the joint are evaluated. Barton and Choubey [36] proposed following relationships by for normal stiffness (kn) and shear stiffness (ks) have been used for analysis which are discussed in FLAC manual: where: E = rock mass Young's modulus; Er = intact rock Young's modulus; kn = joint normal stiffness; s = joint spacing.
( ) where: G = rock mass shear modulus; Gr = intact rock shear modulus; ks = joint shear stiffness. For joint spacing, ref. [37] suggested different spacing values for different rock mass, i.e., for massive rock it should be greater than 10 m. Therefore, considering the value of joint spacing suggested by Palmstrӧm, the stiffness values of interface element have been calculated with the help of afore mentioned formulas. In this analysis the values of kn = 9.43 × 10 2 MPa/m and ks = 3.77 × 10 2 MPa/m for normal and shear stiffness of interface element respectively were used. Sainoki and Mitri [38] discussed that the mechanical properties of the fault which are affected by many factors and are difficult to determine due to varying material accumulation on fault surface and secondly due to scale effect during laboratory tests. As a result, different fictitious values of friction angle of 5°, 15°, 25°, 35° and c = 0 are adopted for the simulation and the sensitivity analysis has confirmed the influence of friction angle on maximum shear displacement as shown in Figure 7.

Interface Element
FLAC 3D provides an interface element to represent the shear zone type structures in a geologic medium. This element is a collection of triangular elements. It can be created at any location in the model. In the model, the interface element is attached to a front block at 50 m length in the right headrace tunnel to reflect the actual field conditions. After creating the interface element, the space between the front and back block was removed and back block was moved 1 m back to get attach with the front block. Both blocks were attached without attaching the shear zone area. Interface element was characterized by Coulomb sliding. Different parameters like normal stiffness (k n ) and shear stiffness (k s ), friction angle ( • ), cohesion (MPa), were assigned to interface element. FLAC [35] explained the stiffness characteristics in an explicit way. So stiffness properties of the joint are evaluated. Barton and Choubey [36] proposed following relationships by for normal stiffness (k n ) and shear stiffness (k s ) have been used for analysis which are discussed in FLAC manual: where: E = rock mass Young's modulus; E r = intact rock Young's modulus; k n = joint normal stiffness; s = joint spacing.
where: G = rock mass shear modulus; G r = intact rock shear modulus; k s = joint shear stiffness. For joint spacing, ref. [37] suggested different spacing values for different rock mass, i.e., for massive rock it should be greater than 10 m. Therefore, considering the value of joint spacing suggested by Palmström, the stiffness values of interface element have been calculated with the help of afore mentioned formulas. In this analysis the values of k n = 9.43 × 10 2 MPa/m and k s = 3.77 × 10 2 MPa/m for normal and shear stiffness of interface element respectively were used. Sainoki and Mitri [38] discussed that the mechanical properties of the fault which are affected by many factors and are difficult to determine due to varying material accumulation on fault surface and secondly due to scale effect during laboratory tests. As a result, different fictitious values of friction angle of 5 • , 15 • , 25 • , 35 • and c = 0 are adopted for the simulation and the sensitivity analysis has confirmed the influence of friction angle on maximum shear displacement as shown in Figure 7.

Sensitivity of Interface Element
The relationship between friction angle of shear zone and maximum shear displacement is shown in Figure 7, which reflect dependency of maximum shear displacement on friction angle in shear zone. It is important to estimate friction angle for evaluating the magnitude of shear slip in shear zone. Additionally, a parametric study of interface element was carried out to investigate the influence of its parameters on maximum displacement and maximum principal stress. The results shown in Figure 8, as the stiffness increases, the maximum displacement and maximum principal stress decrease near the shear zone. The results were similar for all the simulations carried out. The simulations for sensitivity analysis were carried out while fictitious cohesion and friction along the interface remained constant. The normal and shear stiffness were supposed to be equal and different value were used during numerical simulation. Figure 8 showed the variation in maximum displacement and maximum principal stress near the interface element which has shown inverse relation.

Initial Conditions
The FLAC 3D code is used to study the stress and deformation characteristics of rock materials during loading conditions. The shear zone was successfully simulated with this software to study the danger of rockburst. During initial conditions, the average values of the stresses were applied in three directions apart from extreme over-coring results [39]. The maximum principal stress (σxx) with a

Sensitivity of Interface Element
The relationship between friction angle of shear zone and maximum shear displacement is shown in Figure 7, which reflect dependency of maximum shear displacement on friction angle in shear zone. It is important to estimate friction angle for evaluating the magnitude of shear slip in shear zone. Additionally, a parametric study of interface element was carried out to investigate the influence of its parameters on maximum displacement and maximum principal stress. The results shown in Figure 8, as the stiffness increases, the maximum displacement and maximum principal stress decrease near the shear zone. The results were similar for all the simulations carried out. The simulations for sensitivity analysis were carried out while fictitious cohesion and friction along the interface remained constant. The normal and shear stiffness were supposed to be equal and different value were used during numerical simulation. Figure 8 showed the variation in maximum displacement and maximum principal stress near the interface element which has shown inverse relation.

Sensitivity of Interface Element
The relationship between friction angle of shear zone and maximum shear displacement is shown in Figure 7, which reflect dependency of maximum shear displacement on friction angle in shear zone. It is important to estimate friction angle for evaluating the magnitude of shear slip in shear zone. Additionally, a parametric study of interface element was carried out to investigate the influence of its parameters on maximum displacement and maximum principal stress. The results shown in Figure 8, as the stiffness increases, the maximum displacement and maximum principal stress decrease near the shear zone. The results were similar for all the simulations carried out. The simulations for sensitivity analysis were carried out while fictitious cohesion and friction along the interface remained constant. The normal and shear stiffness were supposed to be equal and different value were used during numerical simulation. Figure 8 showed the variation in maximum displacement and maximum principal stress near the interface element which has shown inverse relation.

Initial Conditions
The FLAC 3D code is used to study the stress and deformation characteristics of rock materials during loading conditions. The shear zone was successfully simulated with this software to study the danger of rockburst. During initial conditions, the average values of the stresses were applied in three directions apart from extreme over-coring results [39]. The maximum principal stress (σxx) with a

Initial Conditions
The FLAC 3D code is used to study the stress and deformation characteristics of rock materials during loading conditions. The shear zone was successfully simulated with this software to study the danger of rockburst. During initial conditions, the average values of the stresses were applied in three directions apart from extreme over-coring results [39]. The maximum principal stress (σ xx ) with a value of 60 MPa, intermediate principal stress (σ yy ) with a value of 37 MPa, and minimum principal stress (σ zz ) with a value of 35 MPa were applied. The Table 2 shows the mechanical properties (input parameters) of the rock mass applied during the numerical modeling. The simulated model was run after applying the initial conditions. As discussed in Section 1, the shear zone acts as a barrier and the stresses are usually concentrated there. The model showed that the shear zone was in a shear state and high stresses were concentrated nearby shear zone. Numerical simulation has also given promising results regarding the principal stress concentration near the shear zone before excavation and the initial displacement is set to zero in rock mass as shown in the Figure 9. The maximum and minimum principal stresses contours showed stress concentration near shear zone area on the right side of the right headrace tunnel and there is no such stress concentration in any other region of the model as shown in Figure 9a,b, respectively. The maximum principal stress (σ 1 ) and minimum principal stress (σ 3 ) have maximum value of 62.857 MPa and 38.855 MPa near right tunnel where shear zone is present nearby while such concentration is not present in left tunnel due to the absence of shear zone there. This reflects that the shear zone would be a disturbed area and would be displaced after excavation having very high chance of shear slip movement when the tunnel would be excavated through this zone.  Table 2 shows the mechanical properties (input parameters) of the rock mass applied during the numerical modeling. The simulated model was run after applying the initial conditions. As discussed in Section 1, the shear zone acts as a barrier and the stresses are usually concentrated there. The model showed that the shear zone was in a shear state and high stresses were concentrated nearby shear zone. Numerical simulation has also given promising results regarding the principal stress concentration near the shear zone before excavation and the initial displacement is set to zero in rock mass as shown in the Figure 9. The maximum and minimum principal stresses contours showed stress concentration near shear zone area on the right side of the right headrace tunnel and there is no such stress concentration in any other region of the model as shown in Figure 9a,b, respectively. The maximum principal stress (σ1) and minimum principal stress (σ3) have maximum value of 62.857 MPa and 38.855 MPa near right tunnel where shear zone is present nearby while such concentration is not present in left tunnel due to the absence of shear zone there. This reflects that the shear zone would be a disturbed area and would be displaced after excavation having very high chance of shear slip movement when the tunnel would be excavated through this zone.

Excavation Stage
The headrace tunnels in the NJHPP were both excavated by the drill-and-blast method and the TBM method. The tunnel portion affected by the intense rockburst was excavated by TBM, which has very little flexibility during construction. The behavior of the rock mass and tunnel periphery after excavation is shown in Figure 10a. The average zone state along the tunnel profile showed a shear state (shear-p) on the crown and invert, while the wall of the tunnel was in a mixed shear and tension state (shear-p and tension-p). Figure 10b showed the excavated area of both tunnels. As the excavation process approached the shear zone, the tunnel went in a state of both shear and tension (shear-p tension-p), which further intruded into the shear zone, present perpendicular to the axis of the right tunnel which is shown in Figure 10b. This indicates that the zone state extended its influence beyond the boundary of the tunnel. This unstable area, under the influence of high stresses, caused a shear slip in the shear zone area, which finally resulted in a release of huge amounts of energy in the form of the rockburst. On the other hand, such extended region is absent in left tunnel.

Excavation Stage
The headrace tunnels in the NJHPP were both excavated by the drill-and-blast method and the TBM method. The tunnel portion affected by the intense rockburst was excavated by TBM, which has very little flexibility during construction. The behavior of the rock mass and tunnel periphery after excavation is shown in Figure 10a. The average zone state along the tunnel profile showed a shear state (shear-p) on the crown and invert, while the wall of the tunnel was in a mixed shear and tension state (shear-p and tension-p). Figure 10b showed the excavated area of both tunnels. As the excavation process approached the shear zone, the tunnel went in a state of both shear and tension (shear-p tension-p), which further intruded into the shear zone, present perpendicular to the axis of the right tunnel which is shown in Figure 10b. This indicates that the zone state extended its influence beyond the boundary of the tunnel. This unstable area, under the influence of high stresses, caused a shear slip in the shear zone area, which finally resulted in a release of huge amounts of energy in the form of the rockburst. On the other hand, such extended region is absent in left tunnel.
During normal conditions, maximum principal stress is in horizontal direction, therefore, its concentration would be high at the crown and invert level. The contours of the maximum principal stress showed a normal situation i.e., the stress contour migrates from the left tunnel periphery (Figure 11a). However, due to the influence of the shear zone in the right tunnel, the contours of maximum principal stress was present near the periphery of the tunnel shown in Figure 11a. Furthermore, the concentration of maximum principal stress on the face of right hand tunnel can also be seen in Figure 11b. During normal conditions, maximum principal stress is in horizontal direction, therefore, its concentration would be high at the crown and invert level. The contours of the maximum principal stress showed a normal situation i.e., the stress contour migrates from the left tunnel periphery (Figure 11a). However, due to the influence of the shear zone in the right tunnel, the contours of maximum principal stress was present near the periphery of the tunnel shown in Figure 11a. Furthermore, the concentration of maximum principal stress on the face of right hand tunnel can also be seen in Figure 11b. The intermediate principal stress is usually ignored during 2D stress analysis, but it cannot be ignored when the maximum principal stress is in the horizontal direction [40]. At the face of the tunnel, the FLAC 3D numerical simulation showed that the intermediate principal stress was concentrated at the crown and invert level on the right tunnel as shown in Figure 12a. On the face of right tunnel, impact of shear zone on intermediate principal stress can also be seen in Figure 12b which makes conditions favorable for rock burst while there is no such visible concentration in left tunnel.
The shear zone has a significant impact on shear stress. Therefore, numerical simulation results have shown shear stress concentration near the shear zone in tunnel face area (only right tunnel is shown in Figure 13) and displacement due to induced stresses is shown in Figure 13a,b, respectively.
It is evident from the previous results (Figures 11 and 12) that there is a concentration of principal stresses near the shear plane that caused displacement in the shear zone when the excavation approaches this area as shown in Figure 13b. The shear stress has maximum concentration in the area where the structural plane is present. During this excavation stage, shear stress and shear displacement is produced in the interface that are shown in Figure 14a,b. The shear stresses are concentrated where the shear zone approaches the wall of right tunnel. The results of normal stress and normal displacement in the interface are shown in Figure 15a,b. During normal conditions, maximum principal stress is in horizontal direction, therefore, its concentration would be high at the crown and invert level. The contours of the maximum principal stress showed a normal situation i.e., the stress contour migrates from the left tunnel periphery (Figure 11a). However, due to the influence of the shear zone in the right tunnel, the contours of maximum principal stress was present near the periphery of the tunnel shown in Figure 11a. Furthermore, the concentration of maximum principal stress on the face of right hand tunnel can also be seen in Figure 11b. The intermediate principal stress is usually ignored during 2D stress analysis, but it cannot be ignored when the maximum principal stress is in the horizontal direction [40]. At the face of the tunnel, the FLAC 3D numerical simulation showed that the intermediate principal stress was concentrated at the crown and invert level on the right tunnel as shown in Figure 12a. On the face of right tunnel, impact of shear zone on intermediate principal stress can also be seen in Figure 12b which makes conditions favorable for rock burst while there is no such visible concentration in left tunnel.
The shear zone has a significant impact on shear stress. Therefore, numerical simulation results have shown shear stress concentration near the shear zone in tunnel face area (only right tunnel is shown in Figure 13) and displacement due to induced stresses is shown in Figure 13a,b, respectively.
It is evident from the previous results (Figures 11 and 12) that there is a concentration of principal stresses near the shear plane that caused displacement in the shear zone when the excavation approaches this area as shown in Figure 13b. The shear stress has maximum concentration in the area where the structural plane is present. During this excavation stage, shear stress and shear displacement is produced in the interface that are shown in Figure 14a  The intermediate principal stress is usually ignored during 2D stress analysis, but it cannot be ignored when the maximum principal stress is in the horizontal direction [40]. At the face of the tunnel, the FLAC 3D numerical simulation showed that the intermediate principal stress was concentrated at the crown and invert level on the right tunnel as shown in Figure 12a. On the face of right tunnel, impact of shear zone on intermediate principal stress can also be seen in Figure 12b which makes conditions favorable for rock burst while there is no such visible concentration in left tunnel.
The shear zone has a significant impact on shear stress. Therefore, numerical simulation results have shown shear stress concentration near the shear zone in tunnel face area (only right tunnel is shown in Figure 13) and displacement due to induced stresses is shown in Figure 13a,b, respectively.
It is evident from the previous results (Figures 11 and 12) that there is a concentration of principal stresses near the shear plane that caused displacement in the shear zone when the excavation approaches this area as shown in Figure 13b. The shear stress has maximum concentration in the area where the structural plane is present. During this excavation stage, shear stress and shear displacement is produced in the interface that are shown in Figure 14a,b. The shear stresses are concentrated where the shear zone approaches the wall of right tunnel. The results of normal stress and normal displacement in the interface are shown in Figure 15a (Figure 10b) it is clear that there is no extended shear zone in left tunnel as compared to right tunnel (which has shear zone in the vicinity). In Figure  11a maximum principal stresses are high near the boundary in right tunnel while in left tunnel it has been moved away from tunnel periphery. In Figure 12 it is also compared that intermediate principal stresses are very high in right tunnel as compared to left tunnel.

Comparison of Numerical Analysis with Empirical Rockburst Criteria
After numerical simulation, results are compared with empirical rockburst approaches of Barton and Grimstad [41,42]. The different rock and rock mass parameters were used to predict the rockburst in deep hard rock tunnel. The maximum principal and tangential stresses were extracted from FLAC 3D numerical analysis results with and without shear zone (interface). Table 3 explain the heavy rockburst criteria of Barton and Grimstad along with empirical rockburst assessment, with limit values of different stress components and rock strength parameters.

Discussion
The role of small scale shear zone as an in-situ stress-concentrating barrier is studied in deep NJHPP which has a high stress regime due to tectonics and the maximum principal stress is in the horizontal direction. Because of this high stress concentration, an intense rockburst event occurred near the shear zone during the construction of the project. A nearly vertical shear zone was exposed after the rockburst events in the right tunnel. In this section, the effect of the shear zone on the behavior of tunnels and rock mass is discussed both before the excavation as well as during construction stage. The mechanism behind the rockburst event of 31 May 2015 is specifically focused and numerically evaluated via FLAC 3D and resulted are compared with empirical rockburst criteria.  (Figure 10b) it is clear that there is no extended shear zone in left tunnel as compared to right tunnel (which has shear zone in the vicinity). In Figure 11a maximum principal stresses are high near the boundary in right tunnel while in left tunnel it has been moved away from tunnel periphery. In Figure 12 it is also compared that intermediate principal stresses are very high in right tunnel as compared to left tunnel.

Comparison of Numerical Analysis with Empirical Rockburst Criteria
After numerical simulation, results are compared with empirical rockburst approaches of Barton and Grimstad [41,42]. The different rock and rock mass parameters were used to predict the rockburst in deep hard rock tunnel. The maximum principal and tangential stresses were extracted from FLAC 3D numerical analysis results with and without shear zone (interface). Table 3 explain the heavy rockburst criteria of Barton and Grimstad along with empirical rockburst assessment, with limit values of different stress components and rock strength parameters.

Discussion
The role of small scale shear zone as an in-situ stress-concentrating barrier is studied in deep NJHPP which has a high stress regime due to tectonics and the maximum principal stress is in the horizontal direction. Because of this high stress concentration, an intense rockburst event occurred near the shear zone during the construction of the project. A nearly vertical shear zone was exposed after the rockburst events in the right tunnel. In this section, the effect of the shear zone on the behavior of tunnels and rock mass is discussed both before the excavation as well as during construction stage. The mechanism behind the rockburst event of 31 May 2015 is specifically focused and numerically evaluated via FLAC 3D and resulted are compared with empirical rockburst criteria.
During initial conditions, the right tunnel has more favorable condition for failure as compared to left tunnel where no shear zone is present in the vicinity. Stresses are concentrated near the shear zone and Figure 9 shows maximum and minimum principal stress concentration in rock mass. Near the shear zone, the maximum and the minimum principal stresses are concentrated while rest of rock mass has applied stress conditions. The contour map of maximum shear stress in rock mass after applying initial conditions have maximum value of 21.5 MPa in the center of shear zone. Therefore, after applying initial conditions, it is evident due to stress concentration that the area near stressed shear zone has very high chances of movement in shear along shear plane while such conditions were invisible in left tunnel.
During excavation stage, the zone states, principal stresses, and displacements in the shear zone are studied. The zone state is used to check whether the stresses of the rock mass zones around the tunnel satisfy yield criterion. The twin tunnels are in a shear and tension states at the walls and are totally in shear state at the crown and invert level, which is due to the high horizontal stresses. The zone states extend into the shear zone area while there is no such intrusion in the left tunnel as shown in longitudinal view of twin tunnels. The shear zone is completely in the shear state and it is supposed, huge amounts of stored strain energy is released from rockburst in this region. Because of the abnormal stress conditions, the boundary of the tunnel is in both a shear and tension state. There is no such stress concentration in left tunnel. During excavation stage, maximum principal stress is increased having the maximum value of 99.7 MPa near the right tunnel's crown and invert as shown in Figure 11a while its value is significantly low in left tunnel. In the right tunnel, the maximum principal stress has its concentration along the periphery of the tunnel while it is concentrated away in left tunnel. As a result of this high stress concentration along the periphery, increases the risk of rockburst in the right tunnel. Furthermore, the stress concentration in the right tunnel face has also confirmed the risk of rockburst. The intermediate principal stress has a high value near the crown and invert level of right tunnel and have maximum value of 65 MPa as shown in Figure 12. In case of rock burst, the intermediate principal stress also contributed to the abnormal stress concentration in the area which finally resulted in rock burst. This stress concentration has caused the maximum displacement of 97.5 mm when the right tunnel approaches the shear zone which has maximum shear stress value of 39.5 MPa while such displacement is not visible in left tunnel. The interface element reflects the shear zone in FLAC 3D . The shear and the normal stresses are concentrated in the middle of the shear zone where maximum shear and normal displacement have occurred. When the excavation is reached near the shear zone, the shear displacement reaches to the maximum value of 53.7 MPa there. Comparative study of twin tunnels has shown more stress concentration near right tunnel as compared to left tunnel due to the influence of shear zone which caused in the shear slip resulted in heavy rockburst failure. The numerical analysis result of σ 1 and maximum tangential stress are evaluated empirically with and without interface. As a result, different ratios of empirical rockburst criteria have evaluated heavy rockburst in with the interface.
Therefore, there is a close relationship between high in-situ stress, shear zone, and rockburst. The presence of a shear zone type geological structure in high stress environment makes the conditions ideal for sudden shear failure, which releases huge amounts of energy in the form of rockburst.

Conclusions
This study investigated the impact of small-scale shear zone on rockburst in deep hydropower tunnels. The effect of the structural plane on rockburst is studied theoretically first. Then, the FLAC 3D simulation is used to study this problem numerically. A FLAC 3D model with the representative rock mass and twin tunnels, is built with a Mohr-Coulomb constitutive model. The simulation procedure involves initial conditions analysis and post-excavation stage analysis. Numerical analysis results are evaluated empirically, rockburst evaluation is done with and without shear zone in twin tunnels. The relationship between high in-situ stress, shear zone, and rockburst is explored. The conclusions obtained from this study are as follows: 1.
The shear zone acts as barrier for stress accumulation, which suddenly releases in the form of a rockburst. The amount of energy stored depends on the in-situ stress conditions, the presence of a shear zone, and rock type. The hard rock, which store more energy due to its high modulus of elasticity, results in an intense rockburst in the presence of a nearby shear zone.

2.
The location of the 31 May 2015 rockburst event was near abnormal geological settings, where the bedding plane becomes horizontal due to the concentration of high stresses. These high stresses are the reason behind the aforementioned rockburst.

3.
A FLAC 3D numerical simulation validated the results that were established theoretically. In the high stress regime, the principal stresses are concentrated near the shear zone. The rock mass shows significant displacement after excavation, which reflects that the shear zone area is the most susceptible to rockburst.

4.
After excavation, impact of the shear zone is more evident on right tunnel as compared to the left tunnel. Empirical evaluation of rockburst has also witnessed heavy rockburst in the presence of interface near right tunnel which concludes, shear zone is the main reason behind the intense rockburst of 31 May 2015 in this tunnel. Additionally, the excavation in the nearby left tunnel was smooth, without any rockburst event.

5.
It can be concluded that the presence the shear zone, high stress concentration, and complex geological settings worked together and produced the dangerous condition for the 31 May 2015 event, which resulted in a great loss to human life and property.