Analytical and numerical assessment of a preliminary support design – a case study

Abstract Underground tunneling results in redistribution of the in-situ stress state. This process can create unstable underground conditions. Tunnel instability can result in fatalities, property damages, project delays, high rehabilitation cost, and can compromise the overall functionality of the excavation. This study assesses the efficiency of the empirically recommended support design of the Pahang-Selangor Water Transfer Tunnel in Peninsular, Malaysia using convergence-confinement (CV-CF) and 2-Dimensional two-dimensional Finite Element Method (FEM). To reduce the effect of subjective uncertainties associated with the rock mass rating system, the lump rating system was replaced with the continuous rating system to classify the rock masses and provide empirical support design for the tunnel. The maximum capacity and stiffness of the support systems were analytically determined using the CV-CF method. The results of the maximum elastic deformation of the support systems were compared to the actual critical internal pressures which show that the empirically recommended support systems can satisfactorily stabilize the underground excavation. The Reliability of the support design, the plastic zones, and total deformations were determined using Phase2. The radius of plastic zones and the maximum total deformations were significantly reduced after the supports were installed. The plastic zones that remained were restricted within the reinforced zones, signifying the effectiveness of the recommended support systems.


PUBLIC INTEREST STATEMENT
A rock support system has a primary role of supporting the rock mass around a tunnel after it has been disrupted by the excavation processes. Different tools; like empirical rock mass classification systems, analytical methods, and numerical modeling are used to assess the stability and select optimal support for underground tunnels. One of the serious concerns in the use of rock mass classification systems is that they are subjective. The continuous rating system used for grading rock parameters reduces the subjective uncertainties associated with the rock mass rating (RMR) system. This method produces a controlled RMR estimate between experienced and inexperienced rock engineers. The overall results obtained in this study show that the continuous rating system, convergence-confinement method, and the 2-D Finite Element Method are useful tools for preliminary assessments of underground support designs.

Introduction
Tunnel instability can result in fatalities, property damages, project delays, high rehabilitation cost, and can compromise the overall functionality of the excavation. Tunnels are constructed in civil and mining industries to serve several purposes including; mobility of people, water transfer, conveyance of ore, irrigation, underground storage, military fortification, and hydro-electric power. Stability around an underground excavation is affected by external and internal factors. These factors include; the type of intact rock, type and nature of discontinuities, in-situ stresses, hydrological conditions, and the excavation geometries (Kanik & Gurocak, 2018;Rehman et al., 2020;Xing et al., 2019). Discontinuities such as bedding planes, joints, dykes, shear zones, and faults significantly reduce the strength of the rock mass, complicate the mechanical behavior of the rock mass, and increases the rock deformability. Instabilities may occur when the stress exceeds the rock strength or when a high differential stress is encountered. During underground excavation, redistribution of stresses causes deformation and failure around the openings which pose stability challenges in deep and shallow excavations. Underground excavation challenges can be assessed by gathering site-specific rock engineering data. Critical evaluation of data in rock engineering projects enables the selection of an appropriate excavation method and the required support system for enhanced excavation efficiency and safety. A rock support system in underground excavations has a prominent role of helping the rock mass to support itself after it has been disrupted by excavation. Different tools; empirical, analytical, and numerical modeling are used to assess the stability and select optimal support for underground tunnels. Based on experience and extensive fieldwork, several researchers have developed rock mass classification systems which are considered powerful tools for preliminary support system design in civil and mining engineering tunnels. The commonly used systems are the Q-system (Barton et al., 1974), Rock Mass Rating (RMR) (Bieniawski, 1989), Geological Strength Index (GSI) (Hoek & Brown, 1997), and Rock Mass Index (RMi) (Palmstrøm, 1996).
The advent of sophisticated computers has proved numerical methods, which have the capability of simulating ground complexities, to be an efficient tool in rock and civil engineering. Numerical modeling has been used by several researchers (Ghadimi Chermahini & Tahghighi, 2019;Ur Rehman et al., 2019;Xing et al., 2018Xing et al., , 2019Yalcin et al., 2016) to obtain the optimum support design because of its' ability to provide reliable information such as maximum displacement and thickness of plastic zone which characterizes the rock mass behavior. Based on the thickness of plastic zones and maximum deformation, Kanik et al. (2015) evaluated the effectiveness of preliminary support obtained from RMR 89 and RMR 14 . Similarly, Yalcin et al. (2016) determined the effectiveness of supports recommended by RMR and RMi classification systems. Numerical methods have the potential not only to solve complex tunneling problems, but to improve the understanding and assessment of failure mechanisms, geotechnical risks, and more efficient construction of rock reinforcement systems.
Although the empirical classification systems used by engineers make it easy to establish the preliminary support systems, they do not provide any information on the plastic yield and total maximum displacement needed to establish the optimal support system. One of the serious concerns in the use of rock mass classification schemes is that they are subjective. Field engineers with different experience levels classifying a given rock mass can result in significantly different rock mass classifications. These inherent limitations of the empirical rock mass classification system in underground support design are resolved by adopting both empirical rock classification systems and numerical modeling techniques.
Convergence-confinement is an analytical method described by Sadeghiyan et al. (2016) as a powerful tool for estimating the capacity of supports in underground excavations, the expected plastic zone extent, radial displacement, and final convergence of the opening. The method is a standard industry technique used for preliminary analysis of ground behavior after tunnel excavation. The convergence-confinement technique as indicated by González-Cao et al. (2018) is a simple and cost-effective solution for the preliminary design of excavation support. The method gives rock engineers the opportunity to estimate the load imposed on a support installed behind the face of a tunnel.
In this study, to reduce the effect of subjective uncertainties usually tagged with inexperienced rock engineers, the quality of rock masses along the Pahang-Selangor Raw Water Transfer (PSRWT) tunnel was classified using the widely accepted RMR system via the continuous rating system. Convergence-confinement method was then utilized to evaluate the capacity of the recommended rockbolt and shotcrete in stabilizing the tunnel. Finally, the effectiveness of the support design was numerically verified based on the dimensions of the plastic region and maximum deformations.

Project description and geology of the study area
The Pahang Selangor Raw Water Tunnel (PSRWT) is the sole property of the Malaysian Ministry of Energy, Green Technology, and Water and is located in Central Peninsula, Malaysia. The tunnel is earmarked for transferring raw water (approximately 1890 million litres/day) from Pahang to Selangor state to address the future water demand problems in Selangor and Kuala Lumpur states. The tunnel has a maximum raw water discharge flow rate of 27.6 m 3 /sec, a total length of 44.6 km, a diameter of 5.23 m and a longitudinal gradient of 1/1,900.
Geologically, Peninsular Malaysia is made up of four major tectonic zones, namely; the Western Stable Shelf, the Main Range Belt, the Central Graben, and the Eastern Belt. Igneous rocks with their associated metamorphic rocks are the most dominant rock types found along the tunnel. Aside the Main Range granite of igneous origin, meta-sedimentary rocks of Palaeozoic and Mesozoic age are also present.
The tunnel cuts through two major formations. The Karak Formation of Silurian-Devonian age extends from the inlet portal to chainage 3.82 km consist of metasediments such as hornfels, shale, metaschist, and phyllite. The Main Range granite is composed of coarse to very coarsegrained, porphyritic biotite granite cut by minor porphyritic differentiates and fine-grained variety of granite which represent a later phase of granite intrusion. Micro-granite, granodiorite, diorite, monzonite, granite porphyry, quartz porphyry, felsite, vein quartz, megacrysts biotite granite, megacrystic muscovite-biotite granite and equigranular tourmaline-muscovite granite are the other foundation rocks within the Main Range Granite.
The Main Range Granite extends from chainage 3.82 km to the outlet in Langat, Selangor at chainage 44.4 km. The main range granite is subdivided into three (3) sub-types, namely Bukit Tinggi Granite, Genting Sempah micro-granite and Kuala Lumpur Granite. The Kuala Lumpur Granite and the Genting Sempah Micro-granite isare separated by the Kongkoi Fault and the Bukit Tinngi Fault also separates the Genting Sempah Microgranite from the Bukit Tinggi Granite. While the Kuala Lumpur granite is megacrystic, the Genting Sempah Microgranite consistconsists of microgranodiorite. The Bukit Tinggi Granite is constituted by very coarse grainedcoarse-grained megacrystic biotite granite. The Main Range Granite formation is strongly folded due to the intrusion of granitic rocks. The contact between the Main Range Granite and the Karak formation lies approximately at chainage 3.82 km. The dips of the various rocks of the Karak formation variesvary within a short distance due to folding within the formation. The Karak fault which trends North-South cuts the tunnel at chainage 2.3 km. The rocks in the Karak formation have undergone low to medium grade metamorphism and consist mainly of mudstone, siltstone, sandstone, conglomerate, intraformational conglomerate, quartzite, schistose quartzite, quartz-mica schist, and quartzo-hornfels. Figure 1 shows the general geology along the tunnel.

Geotechnical investigations
For the entire project, scanline surveys and drilled cores according to ISRM (2015) were used to describe the discontinuities and define the rock quality designation (RQD). The discontinuity studies comprised width and infill type, degree of weathering, orientation, persistence, opening, roughness, and spacing of measured discontinuities. Dip and dip direction measurements were  processed with Dips 7.0 software based on equal-angle stereographic projection ( Figure 3) to identify the joint sets for each section. The joint sets in NATM-1 have close to very close spacing with occasional wide spacing, medium persistence, very tight to tight opening, rough to smooth planar roughness, slickenside, and moderate to high weathering characteristics. The attitudes of the identified joint sets in the various sections are listed in Table 1. Laboratory experiments including unit weight (γ), uniaxial compressive strength (σ ci ), Young's modulus (E i, ), Poisson's ratio, and triaxial test were conducted on the core samples drilled from rock blocks of all the geological units based on techniques suggested by ISRM (2015).

Classification of rock masses along the tunnel section
Rock Mass Rating (RMR) also called geomechanics classification was used to classify the rock masses along the tunnel. RMR is a supervised rock classification system, detailed and first published by Bieniawski in 1973. The RMR system relies primarily on six rock mass parameters namely: uniaxial compressive strength of rock material, RQD, spacing of discontinuities, condition of discontinuities, groundwater conditions, and orientation of discontinuities. The ratings of individual parameters are summed to give a final RMR value. The final RMR value ranges from 0 to 100. The total rating is then used to categorize the rock mass into classes. The system has five rock mass classes; I, II, III, IV, and V. In terms of rock mass competence, it decreases from class I through class V, with class V been the least competent rock mass. In tunnel excavation, these rock mass groups are used to determine the ground support system required to stabilize the tunnel walls.
The RMR system has been refined by different researchers after the introduction of the RMR 89 version. Geocontrol (2012) modified the basic RMR (RMR b ) by introducing fracture frequency defined as the number of joints per meter in the excavation face to replace RQD and discontinuity spacing. These changes eliminated the difficulty in determining the RQD from excavation faces and to obtain a good assessment of the condition of the discontinuities in the ground. Celada et al. (2014) updated the RMR 89 version after its existence for 25 years and named it RMR 14 to distinguish it from the earlier versions. The revision includes the addition of new parameters, a revised rating, and a final structure. The three new parameters include; intact rock alterability due to water (swelling), an adjustment factor for the excavation method (F e ), and an adjustment factor related to the stress-strain behavior of tunnel face. Şen and Bahaaeldin (2003) developed the continuous grading system to reduce personal judgment in the existing classical lump-rating system developed by Bieniawski (1973). Best-fit lines and equations for rating of rock quality designation (RQD), average discontinuity spacing ( � X), intact rock strength (σ), and groundwater conditions (G) were developed as given in equations 1, 2, 3, and 4 respectively. With this modification, only the ratings for discontinuity conditions and orientation of discontinuities are affected by subjective uncertainties. This method generates a controlled RMR estimate with a maximum gap of less than 10 % between experienced and unexperienced rock engineers (Şen & Bahaaeldin, 2003).
In this study, the continuous grading system was used to grade the RMR parameters. The continuous grading system only rates the RMR parameters, so the RMR 89 was used to assign the corresponding support design for the obtained rock classes. The recommended support systems are summarized in Table 2.
where r RQD represents the rating for RQD, r� X is the rating for average discontinuity spacing, r σ is the rating for the intact rock strength, r G is the rating for the groundwater conditions.

In-situ stresses
In underground excavations, in-situ stresses play an important role in the design and construction stages. Results acquired from the calculation of in-situ stresses at different places indicate that the magnitude of horizontal stress is usually greater than vertical stress at shallow depths (less than 500 m), while it occurs as hydrostatic at depths of approximately 1000 m below the surface (Brady and Brown, 2004). In-situ stresses in rocks can be obtained either directly by the use of expensive and difficult techniques such as flat-jack, over-coring, under-coring, and hydraulic fracturing or indirectly using empirical approaches and through experiences gained from other nearby underground projects (Rehman et al., 2020). In this study, the vertical stress for each section was estimated from the relation (equation 5); whereσ v = vertical stress (MPa), γ = average unit weight of rock mass (MN/m 3 ), H = overburden thickness in m.
The horizontal stress σ h ð Þ at depth H is much more difficult to determine than the vertical stresses. The equation for estimating the mean horizontal in-situ stress was formulated by Sheorey et al. (2001) and is based on the elastic constants of the rock mass, the geothermal gradient due to the cooling of the crust, and the thermal expansion coefficient. In this study, the horizontal stress for each geological unit was obtained based on the elastic constant of the rock mass using equation 6 as established by Sheorey et al. (2001).
where,where σ h = horizontal stress, H = depth of cover; E = elastic modulus; v = Poisson's ratio; γ = unit weight of the rock mass; β = coefficient of linear thermal expansion of rock (8.0 × 10 −6 /°C appears to be a reasonable representative value for different rock types but not for coal); and G = geothermal gradient (for crustal rocks = 0.0240°C/m, for coal measure rocks = 0.030°C/m).  Hoek and Brown (1980) developed the Hoek-Brown failure criterion as one of the non-linear criteria for defining the strength of rock mass. The failure criterion is generally recommended for intact rocks and moderately to heavily jointed rock masses (Hoek and Brown, 2019). The Generalised

Hoek-brown constants and rock mass deformability
Hoek-Brown constants (m b , s, and a) were determined using the following equations; where m i is the intact rock parameter and can be obtained from the triaxial tests or a published chart. GSI is the Geological Strength Index, D which ranges from 0 for undisturbed and 1 for very disturbed rock mass, is the disturbance factor that depends on the degree of disturbance from blast damage and stress relaxation (Marinos et al. 2007). D was set to 0.1, representing the situation of low disturbance from blasting. In this study, RocData 5.0 software (Rocscience Inc. 2016) was used to acquire the m i of each dominant rock material along the tunnel section.
Tunneling is known to disrupt the initial in-situ conditions, thereby altering the pre-excavation parameters (Satici & Ünver, 2015). Most rock masses at lower confinement levels reach their residual strength when strained significantly by exhibiting some post-peak strength losses. The rock mass parameters would change after excavation, so for a reliable numerical model, the residual load-bearing capacity of the rock mass is an important input parameter since it influences the tunnel stability, and attained only after noticeable plastic deformation. For this study, no postpeak behavior test was conducted. Therefore, the Residual Geological Strength Index (GSI r Þ and the Residual Hoek-Brown parameters (m br ; s r ; anda r Þwere calculated using equations 10 to 13 (Cai et al., 2007).
The deformability of rock masses can be determined directly or indirectly. The indirect methods were used since it depends on empirical methods which make use of the rock mass classification indices and the equivalent continuum approach. The equivalent continuum approach treats the rock as a continuous material and its deformability reflects the deformation properties of both the intact rock and the discontinuities (Zhang, 2017). The deformation modulus E m , and the compressive strength of the rock mass σ rm ð Þ for all the sections were obtained empirically using equations 14 and 15 developed by Hoek and Diederichs (2006) and Hoek et al. (2002) respectively.
where E m is the deformation modulus, E i is Young's modulus of the intact rock, σ ci is the uniaxial compressive strength of the intact rock.
All the parameters used for both analytical and numerical assessment of the supports are listed in Table 3.

Analytical methods (convergence-confinement)
The Convergence-confinement (CV-CF) method is a useful analytical tool for preliminary estimation of support capacities, the extent of the plastic zones, and deformations around an excavation. The method was implemented in two stages. The ground reactions around the tunnel was first analyzed using the closed-form solution developed by Carranza-Torres and Fairhurst (2000). This method is based on the Hoek-Brown failure criterion and assumes a cylindrical tunnel of radius R, exposed to a uniform far-field stress (σ 0 ) and internal pressure σ R . The uniform internal pressure and far-field stress are "scaled" to give the scaled internal pressure P i and scaled far-field stress S o , respectively, in MPa (equations 16 and 17).
where the uniform internal pressure is represented by σ R (MPa), σ o obtained from the relation (σ o = (σ v + σ h )/2) m b and s are the Hoek-Brown's constants. The applied internal support pressure on the tunnel walls is given as σ R ¼ 1 À λ ð Þσ o where λ is the de-confining rate which varies from 0 at the initial excavation stage to 1 for a fully excavated tunnel (De La Fuente et al., 2019).
The critical internal pressure, p cr i marks the elastic limit and is considered as the point where the elasticity of the rock mass ceases and transitions into the plastic behavior. The elastic behavior is maintained when σ R ≥p cr i . A plastic zone of radius R p develops around the opening if σ R <p cr i . The scaled ðP cr i ) and actual (non-scaled) critical internal pressure p cr i are obtained using the equations below: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi The relationship between the radial displacements (U e Þ and internal pressure σ R for σ R ≥ p cr i in the elastic region is given by equation 20.
The radius of the plastic region R p and the radial displacement U p that develops around the tunnel when σ R <p cr i are obtained from the equations below: The equations and methods described above were developed specifically for circular tunnels. Upon realization that most underground excavations are non-circular, Carranza-Torres and Fairhurst (2000) suggested the use of equivalent radius, R e , for non-circular tunnels such as horse-shoe cross-sections. The Equivalent Radius (R e ) which is identical to a circular tunnel of radius R was computed using equations 23 and 24 (Kabwe et al., 2020).
where A is the cross-sectional area of the horseshoe tunnel, L is the height of the tunnel and π is constant.
The input parameters used for the convergence-confinement analysis include the equivalent tunnel radius R e , σ o , σ R, m b , s, σ ci ; E m and v. The equivalent tunnel radius, R e of 2.89 m was constant for all the tunnel sections. The far-field stresses, σ o was calculated for all the sections and the uniform internal support pressure, p i was zero because the de-confining rate is taken as one (1) for fully excavated unsupported tunnel cases. The Hoek-Browns constants, intact rock strength, and the dynamic rock mass properties were changed for each section as described earlier.
In the second stage of the CV-CF analyses, the maximum support pressures and elastic stiffness based on the recommended support systems were estimated for all the outlined sections.
The maximum support pressure P max sh À � and the elastic stiffness (K sh Þ provided by shotcrete are given as: where σ cc the unconfined compressive strength of shotcrete in MPa, E c is the Young's Modulus of the shotcrete, v c is the Poisson's ratio of the shotcrete assumed to be 0.25, t c is the thickness of shotcrete which depends on the surface roughness after blasting and scaling and usually varies between 50 and 100 mm for temporary support. Singh and Bortz (1975) provided a table for the selection of typical values of σ cc and E c for both dry and wet shotcrete mixture.
For equally spaced mechanically anchored rockbolts, the maximum support pressure P max rb , and elastic stiffness K rb are given as; s c is the circumferential bolt spacing given as s c ¼ 2πR n ; where n is the number of equally spaced bolts.
For a section where two or more support systems are installed, Carranza-Torres and Fairhurst (2000) emphasized that the elastic stiffness of each support must be summed to determine the combined effects. For a case where both rockbolts and shotcrete with their maximum pressures defined as P max rb and P max sh and their elastic stiffness as K rb and K sh respectively, the combined stiffness K s is computed as K s ¼ K sh þ K rb . The computed combined effect value remains valid until a maximum possible elastic deformation is achieved by one of the two installed support systems where the combined support system is assumed to fail. The maximum possible elastic deformation is given as: The support with the lowest value of U max r defines the maximum support pressure required for the two support systems operating together. The parameters of the recommended support in Table 2 were substituted into the above equations to obtain the maximum support pressure, elastic stiffness, and the maximum possible deformation for the combined support systems.

Development of numerical model
Numerical modeling techniques are valuable tools in underground excavation designs since opening geometry, strength parameters, and in-situ stresses are considered. In civil and mining engineering, numerical modeling is used to analyze the rock mass behavior and its impact on infrastructure and support mechanisms. To determine the displacement and plastic region around an opening and to verify the empirically recommended support designs, numerical methods such as Discrete Element Method (DEM), Finite Element Method (FEM), Boundary Element Method (BEM), and Analytic Element Method (AEM) have been implemented by several researchers (Ghadimi Chermahini & Tahghighi, 2019;Kanik & Gurocak, 2018;Morris & Johnson, 2009;Rehman et al., 2020;Xing et al., 2019). For this study, Phase 2 V. 7.0, a FEM software, was used to model the deformations and plastic zones around the tunnel. Phase 2 was selected because of its plastic modeling capabilities, ability to handle multiple materials, and automatic mesh generation function. This code only allows 2-D analysis of non-linear deformations.
In rock engineering, the two commonly used numerical modeling approaches are the continuum and discontinuum methods. The choice of a continuum or discontinuum approach depends on a variety of problem-specific conditions; on the scale of the problem and the geometry of the fracture system. Continuum approach may be used if there are only a few fractures and if the opening of the discontinuity and the complete block detachment are not vital components. Discontinuum approach is suitable for moderately to highly jointed rocks where the properties of discontinuities have a bigger impact on deformation and failure modes or where there is a potential for large-scale displacement of individual blocks (Backers, 2010;Jing, 2003;Ma & Fu, 2014). Continuum numerical approach considers the rock as continuous. Though discontinuities could explicitly be added to the model, the continuum approach is only useful in instances where the length of the material's microstructure is significantly smaller than the object under consideration (Xing et al., 2017). The added joint elements only produce limited deformation without any detachment (Kulatilake & Shu, 2015). Continuum modeling assumes that the rock unit cannot be opened or broken and that the joint spacing is close or extremely close (Satici & Ünver, 2015). The majority of the discontinuities spacing in the study area varied between close to very close spacing. This qualifies the surrounding rock mass behavior to be modeled as a continuum model which assumes a continuous material throughout the body (Wyllie, 2018;Wyllie & Mah, 2017).
Due to low overburden thickness, the upper boundary was set free but the bottom, left and right sides were fixed to prevent movement in X and Y directions. Parametric studies were carried out to determine the optimal external boundary which is not influenced by the edge and corner effects, by setting different boundary limits. No edge and corner effects were observed when the excavation boundaries were set at eight times the tunnel width. To obtain reliable results, the external boundaries were set to 10 times the span of the tunnel for all sections. A six node triangular finite element mesh was automatically generated for all the sections and finer zoning was adopted for the tunnel periphery.
Since the interest of this study is to observe the failure and deformations that would occur around the excavation, the plastic material type was chosen. Elastic material type enables the estimation and plotting of strength factors but does not model the failure especially the yielded zone around the opening. The Hoek-Brown failure criterion is recommended when the rock mass is heavily jointed or contains more than two joint sets (Hoek & Brown, 2019). In this research, all the sections contain more than two joint sets (see Table 1), so the Hoek-Brown failure criterion was used to calculate the yielded elements and plastic regions around the tunnel. Tunneling is a 3-D problem, but in this study, a 2-D analysis was conducted, so the load splitting option of the software code was used to simulate the 3-D deformations before the implementation of the support system.
In the final stage of the numerical modeling, the performance of the recommended support systems by the empirical method was investigated by considering the various changes in the maximum deformations and width of the developed plastic zones.

Convergence-confinement method
From the convergence-confinement analysis, the actual critical internal support pressures p cr i calculated for all the sections are greater than the internal support pressure for the fully excavated tunnel providing a valid argument of the expected plastic zone of varying extent around the opening. The analytically determined radius of the plastic zone is 6.33 m, 15.39 m, 5.81 m, 9.84 m, 5.99 m, and 6.18 m for sections 1, 2, 3, 4, 5, and 6, respectively. The results of actual critical support pressure, maximum deformation, and the radius of the plastic zone obtained using the convergence-confinement method for each section are given in Table 4. The strain values for all the sections were empirically calculated using equation 30 (Evert Hoek & Marinos, 2000). All the calculated percentage strain values are less than 1%. According to Evert Hoek and Marinos (2000), a strain value of less than 1% implies only few instability problems are to be expected. Therefore simple tunnel support designs as recommended by the rock mass classification systems could stabilize the tunnel.
where ε is the strain, p i is the internal support pressure, p o is the vertical in-situ stress and σ cm is the strength of the rock mass.
The results of the maximum possible elastic deformation (see Table 4) were compared to the actual critical internal pressures as was similarly done by Kaya and Bulut (2019) and it was realized that the U max r of the combined supports for each section is greater than the actual critical internal pressure p cr i which signifies that the empirically recommended support systems can satisfactorily stabilize the tunnel.

Results of numerical model
The developed deformation and plastic regions were assessed based on an elasto-plastic study. The properties of the recommended support systems used are set out in Table 5. For each section, the thickness of the developed plastic regions, yielded elements, and maximum deformation was obtained for both supported and unsupported cases and presented in Table 6. The maximum total displacement values for the unsupported case for each section are very low. Nevertheless, the plastic zone thickness around the excavation for both crown and walls in all the sections indicates that the tunnel will encounter instability problems.   deformation of 0.625 cm for section 5. In numerical modeling, yielding is one of the most important criteria when the elasto-plasticity analysis is employed. Yielding condition is said to occur when the rock is loaded beyond its elastic limit.
The extent of the plastic zone around the unsupported sections pre-informs of a major instability if the tunnel is left unsupported, especially in section 2 which recorded a maximum plastic zone of 14.61 m. Similar value (i.e 16 m) was obtained for the CV-CF method. Section 2 is dominated by moderate to highly weathered shale and an estimated rock mass strength of 1.64 MPa. The analytical approach used in this research tends to produce similar results with the 2-D FEM and therefore could be deduced that both the 2-D FEM and convergence-confinement method could be reliably used in preliminary underground support designs. The differences in yielded elements, the extent of the plastic region, and the variation in deformations after the implementation of the recommended support at each section were studied and compared with the non-supported scenarios. The plastic zone defines the extent of damage when the material yields. The number of yielded elements were reduced in all the sections.
The total maximum deformation reduced averagely by 11.62 % after the recommended support systems were installed. Except for section 4, the radius of the plastic zones decreased significantly in the other sections. The behavior of the opening is considered problematic if the yielded zones extend beyond the length of the installed rockbolts (Li, 2017). Hoek and Brown (1997) indicated that to avert serious tunnel instability problems, the ratio of maximum tunnel deformation to tunnel radius must be maintained below 0.02. Section 2 produced a tunnel deformation to radius ratio of 0.08 which reflects serious instability problems. Phase 2 is known to produce small strains in numerical modeling, therefore it is limited in accommodating very large strains. Therefore, relying on the obtained magnitudes of deformation rather than the extent of plastic zones and yielded elements for stability assessment is misleading. Although the suggested supports reduced the thickness of plastic regions and yielded elements, a substantial percentage of plastic zones still remained around the openings. As observed from the supported cases in Figures 4 and Figures 5, the failure zones around the walls and crowns for all the selected sections fall inside the supported zone which confirms the stability of the excavation hence affirming the efficiency of the recommended support by the empirical rock mass rating system.

Conclusions
In this study, the stability and effectiveness of the empirically recommended support for the Pahang-Selangor Raw Water Transfer tunnel were investigated using the convergence-confinement method and numerical modeling. To minimize the subjective uncertainties associated with Rock Mass Rating system, the continuous rating system developed by Şen and Bahaaeldin (2003) was used for rating the parameters involved in the Rock Mass Rating system. The convergence-confinement method was analytically utilized as a preliminary tool to estimate the load imposed on the supports and the ground reactions around the excavation. A two-dimensional Finite Element Method was used to analyze the effect of induced stresses and the effectiveness of the empirically recommended support systems.
The main conclusions of this study can be summarized as follows: (1) The maximum elastic deformation of the combined supports (U max r ) for each section was greater than the actual critical internal pressure p cr i À � which signifies the effectiveness of the empirically recommended support design in stabilizing the tunnel.
(2) The radius of the plastic region and maximum deformations estimated by the convergenceconfinement method and 2-Dimensional Finite Element Method were very similar for all the sections along the tunnel route. Thus the two methods can be reliably used for preliminary underground support designs, especially for horse-shoe shaped tunnels.
(3) The radius of plastic regions and the maximum total deformations were significantly reduced after the installation of the recommended support systems. The failure zones that remained were restricted within the reinforced zone, hence confirming the effectiveness of the recommended support systems.
(4) From a practical standpoint, the overall results obtained in this study show that the continuous rating system, convergence-confinement method, and the 2-D Finite Element Method are useful tools for preliminary assessments of underground support designs. The convergence-confinement method used in this research is only applicable for initial or preliminary tunnel design stages.
(5) The rock mass contains joints, thus anisotropic behavior should be expected. As such, the Hoek-Brown criterion may produce deceptive results. For accurate assessments of the support system, non-isotropic models and discontinuum numerical modeling approach such as the Distinct Element Method (DEM) should be explored for future research.