Numerical Verification of Empirically Designed Support for a Headrace Tunnel

In this paper, we used two empirical rock classification systems of rock mass rating (RMR) and rock quality tunnelling index (Q-system) for the support design of a tunnel in District Battagram, Khyber Pakhtunkhwa, Pakistan. Along the tunnel route, the rocks of Precambrian namely Gandaf Formation, Karora Formation and Besham Complex were exposed. During the field investigations, two shear zones were marked in the schist of Karora Formation. The discontinuities parameters collected during the field investigations, results of laboratory testing and material constants determined from RocData version 5.0 software were used during the empirical classification and numerical modelling. The support was designed for the rock mass units from RMR and Q. The quantification of the thickness of plastic zone and total displacement around the tunnel were achieved by the numerical modelling of RS2 9.0 software in both unsupported and supported conditions. The empirically designed support was installed in the model prepared in the RS2 software. According to the results, the empirically designed support when installed in models prepared in RS2 significantly reduced the plastic zone around the tunnel. The reduction in the plastic zone and displacement around the tunnel verified the support design by empirical methods. The present research concludes that empirical designed support can be used for the complex geology of Pakistan.


Introduction
Empirical methods are designed for the determination of qualities of the rock mass and support design for underground excavations. These classification systems are based on a number of case histories and are being used at various stages of the construction of the project for decades. The predefined parameters and easy tabular rating system for rock classification make these approaches popular among geoscientist and engineers [1]. Empirical rock classification systems such as Terzaghi rock quality designation and classification system [2][3][4], Rock Structure Rating [5], Rock Mass Rating [6,7] and Rock Quality Tunneling Index [8,9] were developed for tunnels and underground excavations. Among these systems, RMR and Q systems have been updated in recent years and are still in use in the tunneling industry while others are out dated and are not in use. In recent times,Palmstrom and Stille [10], Genis, Basarir [11], and Palmström [1] discussed that these classification systems have some shortcomings and emphasized upon the verification of the empiricaly design support by using any of the observational, analytical, kinematical and numerical modeling techniques. Recently, different researchers like Ozsan et al. [12]; Genis et al. [11]; Rasouli [13]; Krishna and Panthi [14]; Kaya, Bulut [15]; Panda et al. [16]; Akgün et al. [17]; Elarabi and Mustafa [18] and others in different countries and geological domains have concluded that the empirical systems should be used side by side with any of the pre-described methods at various phases of engineering projects. Larbi [19]; Kanik et al. [20]; Kanik et al. [21] highlighted that the empirical approaches do not determine the plastic zone thickness and values of total displacement. They concluded that numerical modeling should be used for the determination of these factors.
The behavior of rocks mass depends on the rock strength, modulus and in-situ stress conditions. The heterogeneity is found in the rocks behavior in different geological domains which lead to the need of verification of empirically designed support for Pakistan's geological domain. The present study focuses on the verification of proposed support system by empirical methods (RMR and Q) using numerical modeling by Finite Element Method (FEM Codes) RS 2commonly known as Phase 2 .
A case study of a small hydropower project still under construction on the Nandihar Khawar River is taken for the study of numerical verification of empirically proposed support. The project lies in the northwest of Pakistan easily accessible through Karakoram Highway (KKH), 4 km away from Thakot and 25 km from Battagram. Figure 1 shows the important localities of the project. The verification of empirically design support system through numerical modeling requires the detail information about the geology and geotechnical parameters of rock masses. The details of these parameters are given in the section 3. The methodology of the present research work is discussed in next section.

Methodology
A comprehensive methodology is defined for the current study which completes in following steps:  Geological mapping along tunnel route and surrounding areas.  Field investigation for the collection of discontinuity parameters and collection of samples for onward laboratory testing.  Assessment of ground conditions and characterization of rocks using empirical methods.  Support design for rock units using empirical approaches.  Evaluation of in situ stresses using world stress map and empirical equations.  Development of representative models for numerical analysis in RS 2 .  Running the models in unsupported and supported conditions for determination of total displacement and thickness of plastic zones.

Geology, Field and Laboratory Studies
The project site is located in the Higher Himalayas, south of Main Mantle Thrust (MMT), west of Chakesar Fault Zone and north east of Thakot shear zone. The rocks of Precambrian age, Besham Complex, Karora Formation and Gandaf Formation, are exposed in the area. Besham Complex are the basement rocks subjected to various degree of metamorphism. The ortho and para-metasediments now have resulted into the formation of Granitic Gneisses, quartz mica schists and the graphitic schists with marble at places. The Karora Formation is a sequence of marine metasediments which were laying un-conformably on top of Besham Complex rocks. The unconformity is marked by the metamorphosed pebble conglomerates that grade into a thick unit of graphitic phyllite and marble. The Gandaf Formation comprising graphitic, garnet schist, graphite slate, phyllite and schist, fine grained meta-psammite, argillite, calcite marble, tremolite marble and quartzite, underlain by Karora Formation [31]. Geological map of the project area is shown in Figure 2.

Figure 2. Geological map of the project area
Five rock mass units are encountered through the route of tunnel. Variegated Schist and Gneiss (Besham Complex), Graphitic Schist and Quartzite (Karora Formation) and Granite (Gandaf Formation) are youngest to oldest rock units exposed along tunnel route as shown in cross-section of tunnel ( Figure 3). During the field investigation and geological mapping, two shear zones were marked in the schist, so overall seven rock mass units were identified. These shear zones make the schist critical, so the tunnel route in schist is divided into three zones. The detail of rock mass units along tunnel route are given in Table 1.  Nine discontinuities surveys were conducted, along the tunnel route, to gather the discontinuities characteristics of all rock units. The Brunton Compass and google mobile application Clino (Farny, 2017) were used to record the orientations of joints along the tunnel. The dip/dip direction of each joint set was recorded on the field performa. The trend of tunnel is NW-SE with horizontal plunge. Other discontinuity parameters such as spacing, persistence, infilling and aperture, roughness and hardness were also recorded for each joint following ISRM's suggested methods [32]. The roughness of joints was measured with profilometer, while ground water conditions and weathering conditions of rock wall were observed visually. Generally, during the field investigations, it was found that the joints were rough, irregular and moderately undulating with moderately weathered walls. Through Stereonet plotting in DIPS® (RocScience Inc) four different joint sets in granite, while in all other rock mass units, three orthogonal joint sets, were identified. Schmidt rebound hammer was used for the estimation of strength of rock mass units. The rebound number value (Rn) was noted down for wall strength of each joint and different samples were also collected for the onward laboratory testing.
The rock samples collected during the field studies were tested in the laboratory for the determination of unit weights ( ) of rocks following the methods suggested by ISRM. These unit weights with Schmidt rebound numbers (Rn) were used for the estimation of uniaxial compressive strength of rocks using Miller correlation chart [33]. The spacing of joints recorded during field studies were used for the assessment of joint volumetric count (Jv) and block volume (Vb). While joint volumetric count was used for the calculation of Rock Quality Designation (RQD) of rock mass units using Palmstrom (2005) correlation. The parameters recorded during field studies, results of laboratory testing and engineering geological and geotechnical parameters ( Table 2) were used for the rock mass characterization and support estimation.

Empirical Analysis
Initially, the ground conditions were assessed using empirical systems. According to Lowson and Bieniawski [6], the underground excavation process is not always smooth as a number of problems may encounter due to unstable rock mass surrounding the opening, such as swelling, squeezing or rock burst, which threatens the safety of crew and equipment. The problems which may occur in the tunnel is determined by using empirical method of Singh et al. [34] and Goel et al. [35].
Based on case histories, Singh et al. [34] proposed an empirical method using rock mass quality (Q) and overburden (H) on rocks. For computing Q the SRF value is equal to 2.5 in this approach. Singh et al., (1992) equation is given as: Goel et al. [35] have recommended an empirical approach based on rock mass number (N).
Above Equation 2 recommends the use of N value which is obtained when Q value is determined for SRF=1 and B is width of tunnel. The graphical results obtained from these approaches are shown in Figure 4 (a and b) as suggested by Singh et al. [34] and Goel et al. [35]. The results obtained from empirical approaches shows the non-squeezing potential and minor to non-squeezing potential in each rock mass unit using Singh et al. [34] and Goel et al. [35] approaches respectively. During the empirical rock mass classification of rock masses, the ground conditions were used accordingly.

Rock Mass Rating (RMR) Classification System
The classification system was first published by Bieniawski [7], later Lowson and Bieniawski [6] updated the system in 2013. This empirical approach was developed by the experience gained in rock mechanics industry and in the light of results of case histories [6]. The basic parameters of classification system includes the jointing parameters, strength of rock (UCS), rock quality designation, presence of water, spacing of joints and major orientation of joint set. The rating values varies from 0 to 100 and it is the sum of the individual rating of each parameter. The guidelines for support design are predefined for each rock class [6].
During the current study, the tunnel excavation drive with respect to major joint set were found to be fair to very favorable. The ratings of individual parameter were given according to the guidelines. The schists are the weakest rock mass unit with the poor rock mass quality while all other rock masses are fair to good in quality. The details of ratings are given in Table 3. The support details obtained from RMR of different rock mass qualities are given in Table 5.

Rock Mass Quality Index (Q-System)
Q-system was developed for tunnels and caverns at Norwegian Institute by Barton et al., in 1974. They introduced three quotients namely block size, active stress factor and inter block shear strength. The block size is a ratio of rock quality index (RQD) and number of joints (Jn), while the ratio to joint roughness (Jr) to the joint alteration (Ja) defines the shear strength of joints. The stress factor was introduced by dividing the joint water condition (Jw) to the stress reduction factor (SRF). The individual ratings were given using the tabular guidelines and the final ratings were obtained from the product of these three quotients. The stress reduction factor (SRF) ratings were given according to the field investigation and the results of ground condition were assessed by the empirical approaches. The details of rating for each rock mass units are given in Table 4. The support estimation in Q system requires the determination of a factor of excavation support ratio (ESR) and equivalent dimension (De) according to the guidelines of Barton et al. [8]. The ESR value of 1.6 was used during the current study. The supports estimated for the rock mass units are given graphically in Figure 5 and described in Table  5.

Recommended Support System
In this research, supports of rock bolts, shotcrete and wire mesh were recommended by RMR. The support suggestion by Q system ends up to be unsupported ( Figure 5) so this system is inadequate for the current verification study. The recommended support by RMR was verified by numerical method. Recommend support with average stand up time has shown in Table 5.

Numerical Analysis
Finite element method (FEM) codes were used to assess the stresses, deformation and plastic zones in supported and unsupported conditions using a computer aided program of RocScience [36]. Many researchers have used numerical analyses that become dynamic part of geotechnical projects in the last three decades. [37][38][39][40][41][42][43] have successfully utilized the numerical techniques to sort out the rock engineering problems. Phase 2 is a computer aided two dimensional program of RocScience, (2014) which has been used for numerical verification in the current study. The modeling in phase 2 starts by adding the excavation and external boundaries. These boundaries are divided into small triangular and rectangular mesh segments. The granular behavior of rock masses follows the Hoek Brown strength and failure criteria while soils follow Mohr Coulomb strength and failure criteria in general (González de Vallejo and Ferrer, 2011). Two dimensional study of rock units was carried out by FEM using aforementioned criteria. Some input parameters such as uniaxial compressive strength (MPa) and unit weight of rocks are determined in laboratory while Hoek-Brown constants (mb, s, a) were calculated by a program RocData v. 5.0. Values of these parameters are given in Table 6. The in situ stress conditions are also basic input parameter of numerical analysis. The details of in situ stresses are given below

In Situ Stresses
In situ stresses varies with tectonic setting and geological history of the area. The inclination of sigma 1 is 30 o towards 153 o N as determined from the world stress map. Horizontal stress mainly arises due to tectonically active region while vertical stresses are related to overburden or depth of tunnel. Vertical stress is calculated by multiplying tunnel depth with density using Equation 3 of Fenner [44]. In 1998, Sengupta [45] proposed relations to calculate the field stresses for overburden less than 400 m (Equation 4 and 5) and Stephansson [46] suggested the relations for overburden less than 1000 m (Equation 6 and 7). The same correlations were used during the present study for the determination of in situ stresses. By measuring in situ stresses, it was observed that magnitude of horizontal stress is greater than vertical stress. The determined stresses are given in Table 7. These in situ stresses and basic input parameters of Phase 2 were used in elasto-plastic conditions. Automatic three noded triangular mesh was generated in the program to compute the deformation and stresses. The thickness of plastic zones, yielded elements and total displacement were determined in both supported and un-supported conditions and are given in Table 9. Fully bonded bolts of 19 mm in diameter and tensile capacity 0.1 MN were used during the analysis, the details of support properties of bolts and shotcrete are given in Table 8. Initially, the models were run in unsupported conditions. The highest displacement was observed in the schist in the range of 289-828mm and it confirms the field evidences of weak properties and because of the shear zones marked in this rock mass unit. Similarly, the thickness of plastic zone is highest in schist as well. The details of results of numerical analyses are attached in Table 9. During the analysis in unsupported condition, the software interpreted the zone of deformation formed by the shearing marked by red color cross (x) and tension marked by red color circle (o) in Figures 6 (A to N). The support recommended by empirical methods were installed in the model for the verification purpose. According to the results, the total displacement decreases drastically for all rock mass units. The installed support in RS 2 models decreases the total displacement from 828 to 273 mm for the weakest rock mass unit of schist verifies the workability of empirical support for the weakest rock mass units as well. The models of numerical analyses in unsupported and supported are given in Figure 6 (A to N).

Conclusion
During the present study, the numerical modeling technique of RS 2 was used for the verification of support proposed by empirical methods. The data collected from field and results of laboratory testing were used for the rock mass classification and determination of qualities of rock mass units along the tunnel alignment. The rocks were characterized as good to poor from RMR and fair to very poor from Q-system. The shear zones marked in the schist makes these units the most critical among others and the poor rock quality were determined from RMR while poor to very poor were determined from Q-system. The support is designed for the all the rock mass units from the both classification systems. The Q system ends up with the unsupported and concluded to be inadequate for the present study. The RMR recommended support with bolts and Shotcrete were installed in the models prepared in RS 2 software decreases the total displacement and thickness of plastic zone around the tunnel. The present study confirms the previous results of Kaya et al. (2011), Kanik et al. (2015), Kanik et al. (2018) and others that the numerical modeling is necessary for the support design and verify the support recommended by the empirical methods. The study also concludes that the empirical systems cannot predict the thickness of plastic zone around the tunnel and the displacement of tunnel inward that can be predicted well with numerical modeling such as RS 2 software.