Numerical modelling for geotechnical assessment of rock mass behaviour and performance of support system for diversion tunnels using optimized Hoek-Brown parameters

Purpose. Empirical and numerical methods play a vital role in assessing rock mass behaviour quantitatively and qualitatively to design underground structures/caverns and support systems. This research aims to assess and evaluate the rock mass be-haviour for safe, stable, efficient, and economical design of support system for underground structures especially tunnels in diverse rock mass conditions. Methods. In this research, such empirical design methods as Rock Mass Rating (RMR), Q-system and GSI were used to characterize and classify the rock mass environment along the tunnel for the preliminary design of twin tunnels and support systems. The geomechanical parameters, Hoek-Brown failure criterion, and its variants for assessing rock mass behaviour were optimized using multiple regression of Stewart, generalized and globalized variant of nonlinear regression method. The rock mass was classified for the selected section A-A. The excavation method and support system for the said section were designed based on the results obtained from empirical modelling. 2D elasto-plastic finite element method (FEM) was used for numerical analysis of rock mass behaviour and performance of the designed supports in section A-A. Findings. The major rock type encountered in the diversion scheme comprises gabbronorite (GN) and Ultramafic Association (UMA). Based on the quantification of RMR, Q-system, and GSI, section A-A’s rock mass ranges from very poor to poor. From the numerical analysis for the said rock mass environment both RMR and Q system support recommendations are equally efficient to support the rock mass surrounding the tunnel. However, keeping in view the yield zone, especially in the crown, the rock bolt’s length should not be less than 5 meters. Based on the analysis of results, both the tunnels are at a safe distance from each other. Originality. In this research, the design input parameters for numerical modeling were optimized by using different techniques to eliminate the chances of error in evaluating rock mass behaviour and designing an optimum support system in the said rock mass environment. Practical implications. The assessment of rock mass behaviour and the design of optimum support systems in heterogenous conditions is quite challenging and requires thorough investigation through different design techniques. This research provides a refined meth-od to be used for the safe, stable, and economical design of tunnels.


Introduction
The design and construction of unground structures involve certain potential risks due to the nature and characteristics of their spatial variation, rock mass behaviour, and level of knowledge. The success of an underground project can be achieved through advance and effective geotechnical investigation, adoption of the effective design method, effective ground stabilization, and monitoring techniques [1]- [3]. In the preliminary stage of execution of any under-ground civil and mining project, limited data about subsurface geology, ground hydrology, strength & stiffness of rock mass, and response or behaviour of rock mass to excavation is available [4]- [6]. Empirical design methods have success stories in the design of underground structures, both soft and hard rocks [7]- [13]. At the primary stage of tunnelling projects, empirical design methods, especially rock mass classification systems, can solve rock engineering problems [14]- [16]. Among these classification systems, Rock Mass Rating (RMR) and Q-system are internationally accepted design methods commonly used in the field of tunneling [4], [6], [17]. Although the empirical methods provide acceptable design for underground structures, these methods do not evaluate the response of excavation, rock mass behaviour, and effectiveness of support system in detail. An empirical analysis of tunnel in rocks, modeling of rock masses is challenging due to anisotropy, heterogeneity, non-elastic and nonlinear nature of rock mass, and requirement of quality input data [18]- [20]. Furthermore, the design aspects i.e. shape and size of tunnel/exaction and support sequence, make the modeling more complex [21]- [23]. To evaluate the performance of support structures, stress redistribution, and stress deformation around tunnels, the empirical design methods are aided by numerical methods to produce more viable, authentic, safe, and economical design for excavation and supports [24]- [26].
Due to the low cost, time-efficient nature, convenience, and availability of user-friendly codes, numerical methods have got more attention in civil engineering and rock engineering for the solution of complex geometries tunnels and rock conditions [27]- [29]. Moreover, the addition of numerical analysis minimizes the risk uncertainties in the design. However, selecting a method out from available numerical methods depends on many factors including the nature of the problem, the capability of a method to solve a problem, and the simplicity of the codes available. The numerical methods give an optimum mathematical solution to a problem based on engineering judgment and rock mass behavior.

Materials and methods
In this study, the geological and geotechnical data of the Diamer Basha Dam site, especially the diversion scheme, is acquired from different investigative reports and literature. The ground conditions are established by collecting geological, hydrological, and mechanical data obtained from different reports and research papers related to the site. To characterize the rock mass of the dam site and obtain quality inputs to both empirical and numerical methods, the data is scrutinized and analyzed through statistical tools. RMR and Q-systems are used as design tools to design excavation and supports for diversion tunnels, whereas GSI is used to obtain inputs for the numerical tool. 2D Finite Element tool Phase 2 is used to evaluate the performance of the design support system and the response of ground to excavation.

Diversion scheme
The project consists of two parallel diversion tunnels with a centre-to-centre distance of 50 meters on the Right Bank (RB) of the dam and a diversion canal (DC), passing on the RB. Diversion Tunnel No. 1 (DT-1) is proposed to be used later as a flushing tunnel for the RB intake. The Diversion Tunnel No. 2 (DT-2) shall be plugged after its operation is finished. In this study, only tunnels of the diversion scheme will be analyzed. The DT-1 is 782 m long, while the DT-2 is 911 m long. These are D-shaped tunnels with 15.5 m width and 15.5 m height.

Location and geology
Diamir Basha Dam is proposed on River Indus, between the Khyber Pakhtunkhwa province and Gilgit Baltistan in Pakistan. The dam site is approximately 315 km upstream from the existing Tarbela Dam, 165 km downstream from Gilgit city, and 40 km from Chilas as shown in Figure 1  For feasibility, NEAC Consultant performed a comprehensive range of geotechnical studies, including exploration drill holes, geological mapping, geo-physical survey, two exploratory adits. The Diamer Basha Dam project is situated within the Ju-rassic-Cretaceous island arc in northern Pakistan, known as the Kohistan Arc. Due to the Indo-Pakistan Plate subduction below the Eurasian Plate is mainly produced between > 130 to 55 Ma [30].
The arc of Kohistan is an ancient island arc, created along with the arc of Ladakh as an eastern continuation, divided by the Nanga Parbat's syntax. They were sutured along with MKT (Main Karakorum Thrust, Shyok Suture) to the Eurasian Plate's previously accreted Karakoram block (80 ~ 90 Ma), and along with MMT in the south and east with the Indian Plate.
The major rock type encountered an area where the proposed diversion scheme comprises Gabbronorite (GN) and Ultramafic Association (UMA). The Ultramafic magmas have been injected into the Gabbronorite while they were still in molten/partially molten form. The rock body is made up of intermixing of the two types of magma.

Physical & geotechnical properties of the rocks
At the outcrops, the GN seems to be fresh and hardly weathered. Some part of the rock is just discoloured and mechanically weathered up to a depth of few centimetres. The intensity of this weathering seems to be proportional to the content of mafic minerals. These parts' surfaces show a light brown to light orange colour from iron oxides, which is derived from the weathering of Fe-bearing minerals. Below that, the rock shows occasional slightly discoloured pyroxenes and is otherwise fresh and hardly weakened. The aver-age unit weight of the rock is 29.2 kN/m 3 .
Weathering and alteration affect the UMA rocks more than GN due to their minerals' less chemical and physical resistance. Usually, the classic pyroxenites and websterites are showing a red to rust-coloured staining on exposed surfaces. The grain boundaries are sometimes weakened to the extent that the rock can be easily crushed by hand. However, this zone of decomposition is diminishing after a few centimeters or more in some cases. The UMA is heavier than GN, with an average unit weight of 32.3 kN/m 3 .
Seven boreholes in the area along the axis of diversion tunnels grouped the rock mass along the tunnel into three distinct geotechnical units. Data from boreholes and core logging suggest that the rock mass with high RQD values and mild to no weathering is of good quality massive rock. GN rock was observed in the area from the inlet up to 482 m along the tunnel axis in geotechnical unit -1. According to the results from the other two boreholes, the rock mass from 482 to 630 m along the alignment is expected. The region of joint GN is located sub-parallel to the tunnel at and around the elevation of the diversion and attack. Two boreholes are drilled in Geotechnical unit 2. Ultramafic rock is located at 633 m from the Diversion Tunnels inlet. With moderate to strongly weathered rocks, the rock mass is more intensively joined together. The rock mass mostly consists of joints and therefore shows low RQD. Geotechnical Unit-3 starts from chainage 800 m; inlet up to the outlet of the tunnel goes in GN rock. As both the tunnels are parallel to each other and hence due to similarities in geological conditions, the discus-sion is on only Diversion Tunnel-1. Great variations have been observed in mechanical properties within the same geotechnical unit, especially strength, Poisson's ratio, and elastic modulus. To investigate the variation in the data and to obtain more realistic values of the parameters, the collected data is re-analyzed.
The geotechnical data was statistically analyzed along the tunnel axis in detail. Analysis of the results is presented in Table 1. The results show that the average value of water absorption and porosity for UMA is 0.7% and 2.22 respectively, while for GN Their value is 0.22% and 0.65%, which shows that UMA is more porous and permeable than SGN. Similarly, for UMA, the average uniaxial compressive intensity with a standard deviation of 32.92 is 85.80 MPa. The mean value of Brazilian tensile strength for UMA's is 5.60 MPa with a standard deviation (SD) of 1.02, and for GN, the mean value of compressive and tensile strength is 124.10 and 8.20 MPa with an SD of 34.10 and 1.79, respectively, which demonstrates that GN rocks are highly resistant compared to UMA. To obtain the behaviour of rock at different stress levels and estimating rock and rock mass parameters, Hoek-Brown failure criterion and its variants, the Multiple regression of Stewart, Generalized and Globalized variant was applied to laboratory tests data. The failure criterion parameters are optimized using Excel add-in "Solver" as shown in Figures 2 and 3

Figure 3. Curve fitting of Hoek & Brown failure criterion and its variants for data to intact UMA
The comparison of various variants is carried out based on residual. It is observed from the analysis that the globalized variant optimally described the behaviour of both the intact GN and UMA rock at different stress levels.
The various parameters gained from different fitting techniques are presented in Table 2 and 3. After analysing data, the most probable values of physical and mechanical properties are summarized in Table 4 and  Table 5, respectively.

Rock mass classes long tunnel axis
The geological section along the alignment of DT1 is shown in Figure 4. Based on the analysis of geological and geotechnical data, the rock mass along the axis of the diversion tunnel was classified into three zones, based on the dominating range of RMR values. The rock mass from inlet up to 482 m (Ch-00 to 742) is a massive rock of good quality, with medium to widely spaced joints GN. According to the RMR system, "Good" rock mass is the dominating class in this geotechnical unit. From 482 m (Ch-742) up to 633 m (Ch-893) the rock mass is of "Fair" to "Poor" quality with weak and closely jointed zones at places. From 633 m (Changes 893) up to outlet (Ch-1143), the tunnel passes through UMA, and the quality of the rock mass decreases as the rock shows many fractured zones and weathered joints with a 15 m thick zone of weakness, cutting the tunnel at a low angle and "Poor" rock class is predominant in this zone.
For numerical analysis, section A-A is selected along the alignment of the tunnel. The RMR value of rock mass is calculated from geological and geotechnical, and borehole data as presented in Table 6.   GSI and Q values were estimated from RMR values. Classification of rock mass for the selected section based on RMR, GSI, and Qthe system is given in Table 7.

Classification system
Section A-A System rating Rock mass quality RMR-system 37 Poor Q-system (Q = exp((RMR-44) / 9) 0.46 Very poor GSI (GSI = RMR + 5) 42 Poor The support system recommended for the said section using RMR and Q-system is presented in Table 8.

Rock mass properties
Some rock mass properties are needed to be used as input for numerical analysis of excavation design. UCS, UTS, deformation modulus and Hoek-Brown parameters for rock mass mb, s, and a are the most significant input parameters.
The strength parameters of intact rock and GSI in RocLab software version 1.1 are used to obtain rock mass parameters for the said section along the tunnel alignment (see Figure 4) summarized in Table 9.

Classification system
Rock mass quality  In design, the deformation module is the most representative input parameter, especially in numerical analysis. Two pre-existing models as proposed by [31] and [32] are used to estimate the deformation modulus of the rock mass along the tunnel's alignment. The deformation modulus values obtained from these models are given in Table 10.

In-situ stresses
In this research, the following Equation 1 was used for the estimation of vertical stresses: where: γthe unit weight, MPa; Hthe height of overburden, m. The ratio between the horizontal stress and vertical stress is called constant represent by K and mainly depends on the depth of overburden. Determining the value of this ratio by calculation is more practical. The theoretical approach to determining horizontal stress from vertical stress is, however, easy to use. Horizontal stress is shown to depend on the constant elasticity of the rock mass [7]. The following Equation 2 is used for horizontal tension stresses: where: υpoison ratio; βindicate the coefficient of thermal expansion, and mostly its value for rocks is 8·10 -6 /ºC; Ermthe rock young modulus, MPa; Gthe rock thermal gradient, ºC/m. However, the below-stated relationship is adopted in this research for estimation of the horizontal stress as: ( The overburden stresses along the tunnel axis are summarized in Table 11.

Numerical analysis of support and stability
Based on geoengineering interpretation and consideration, together with the subsurface investigation and tunnel face observation, the rock mass along the tunnel is classified into three distinct zones varying from good to poor rock conditions. In the zone where the predominant rock mass belongs to the poor class, section A-A (Fig. 4) is selected for numerical analysis to analyze the stability of the tunnel and validate the support systems designed empirically.
Software Phase 2 is used to analyze the applicability of empirical design methods, i.e. RMR and Q, to determine the induced deformation around different sections and investigate rock and recommended support interaction.
The parallel diversion tunnels DT1 and DT2 are modelled in five stages. In stage-1, the virgin field conditions are validated, and in the stage, the rock mass behaviour was studied after excavation and support installation. The model is simulated using an elasto-plastic constitutive model with Generalized Hoek-Brown failure criteria for RMR and Q support systems. Figures 5, 6, 7, 8 and 9 show the numerical simulation results using Q support systems in Phase 2 .

Figure 13. Yield zone and elements at stage-5 (RMR Support)
The Results from both RMR and Q support from simulated models are presented in Table 12 and 13 respectively.
The numerical analysis results validate that both the RMR and Q systems suggest a safe support system for the given rock mass condition. No yield elements are observed in rock bolts as well as shotcrete as recommended by RMR and Q system. The RMR support model (stage 4) shows that the yield zone in the crown is 4.6 m, while for the Q Support model, the Yield zone in the crown is 5.1 m.

Conclusions
The major rock type encountered an area where the diversion scheme is proposed, comprised of Gabbronorite (GN) and Ultramafic Association (UMA). The geological and geotechnical data indicate that in GN rock mass, the variation in the mechanical properties is less than UMA rocks, probably due to the inertness to chemical weathering of GN rock. Near the surface and at depth below the water table, the intact UMA rock pieces have low strength and stiffness properties.
Based on the quantification of RMR, Q-system, and GSI the rock mass quality at section A-A ranges from very poor to poor. For the said rock mass environment, both RMR and Q system support recommendations are equally efficient to support the rock mass surrounding the tunnel safely. No yield elements are observed in rock bolts as well as shotcrete as recommended by RMR and Q system. The RMR support model (stage 4) shows that the yield zone in the crown is 4.6 m, while for the Q support model, the yield zone in the crown is 5.1 m. However, keeping in view the yield zone, especially in the crown, the rock bolt's length should not be less than 5 meters. Based on the analysis of results, both tunnels are at safe distances from each other. This research provides a refined method to be used for the safe, stable, and economical design of tunnels.