Empirical and finite element based stability analysis of highway cut slopes in Uttarakhand Himalayan terrain, India

The present paper deals with different empirical methods and finite element method of slope stability analysis along National Highway (NH)-7, in Uttarakhand, India. The highway is only path in the hilly terrain of Lesser Himalayan for the public transport and have strategic importance due to militaries possession routes. This route is also significant due to having many holistic places, connecting to this. There was numerous landslides happened along the Highway in past due to various natural and anthropogenic activities. Hence, keeping an eye to the socio-economic development of the distant area, slope stability analysis is very crucial along the road cut sections. To identify the vulnerable locations and to collect the geotechnical data, the field investigation was carried out between Shivpuri to Byasi along NH-7 in Garhwal, Uttarakhand. Then geotechnical data was intended followed by rock mass characteristic, kinematic analysis and Qslope stability. Additionally, to review the stability results, numerical simulation (finite element method) was employed and slope mass behavior and failure mechanism of cut slopes were also evaluated. The rock mass characteristic and kinematic analysis illustrate normal and good variety of rock mass mainly wedge mode with flexural toppling of failure. The slope mass rating, continuous slope mass rating and also Qslope stability analysis showed, road cut slopes are critically stable and unstable. The results of different empirical methods shows a decent correlation between them. Further the numerical simulation analysis also evaluates that two cut slopes are unstable and other one is critically stable. This substantial empirical and numerical analysis of cut slopes provides a collective approach to stable and develop the holistic road corridor in Himalayan terrain.

areas, of which 0.04 million km 2 falls in Uttarakhand state. In Uttarakhand, Rishikesh-Badrinath (Mana) NH-7 is significant path for the transportation, public tourism and socio-economic activities. This mountainous road corridor is also a prime medium for religious activity, where pilgrims use this route to reach their respective shrine and temples located in northern proximity of Uttarakhand.
The Himalayan terrain is tectonically active, deformed and dissected by various unfavourable and oriented discontinuities. The variation in rain and temperature throughout the year, make this region more susceptible to slope failure. Also, increase in population imposes immense pressure on existing infrastructure of the Himalayan region and anthropogenic activity such as road development, tunnel and hydroelectric projects enhanced the frequency of occurrences of landslides in valley [24,42,43]. Majority of landslides occurred along NH-7 initiated from cut slopes which are poorly designed, excavated and left unsupported, affecting the road alignment and become prone to failure during rainy season [1,28,39,42]. Percolation of water through weak fracture and joint plane exposed the overburden mass resting on the cut slopes and slip it down slope [53].
Varnes [51,52] and Cruden [11] classified the earth movements downslope based on multiple processes involved such as rock and debris slide, block topple, rock fall, earth flow, avalanches. Such processes largely controlled by failure mechanism, rock type and prevailing geo-hydrological and geotechnical conditions. In recent past, several occurrences of massive landslide have caused large-scale human calamity, material and infrastructural damage and associated environmental, topography and social hazards in lesser Himalayan region [39]. The occurrence of any such hazard is devastating and the resulting effect become worse when one event triggers another, such as flash flood induced landslides of Kedarnath in June, 2013 along the entire stretch of Ganga-Alaknanda valley. Considering numerous challenges in slope stability, optimum and stable slope design is a major concern for authorities along NH-7. Therefore, it is required to predict the failure mechanism, factor of safety and vulnerable portions of cut slope to avoid any casualty due to rock mass failure [28].
In past few years several authors have performed stability analysis using conventional and numerical method of different road sections of mighty Himalayan terrain and suggested instability conditions/mechanism and general method to enhance the stability of cut slopes. Sati et al. [39] carried out stability appraisal of landslide in Uttarakhand and suggested that intensity of fresh slides is higher in the newly cut road sections, however due to reactivation and enlargement of existing landslides voluminous dimension of slides is higher as compared to new road cuts. Apart from providing cause and occurrences of many slope failures, the author has not carried out stability assessment of particular slope and recommended other agencies/department/researcher to devise methodology for the prediction of stability. Further with time many authors have done slope stability and landslide analysis in jointed rock slopes of Himalayan/mountainous region and several other part of the world [1, 30-32, 37, 38, 41, 44].
The present study was undertaken along NH-7 between Shivpuri to Byasi, in Uttarakhand, where three cut slopes were selected to understand the rock mass characteristics and slope behavior exclusively in steep jointed phyllitic quartzite rock mass. This paper presents the rock mass characteristics (GSI and RMR basic ) and slope stability evaluation through empirical methods (SMR, CoSMR and Q slope stability) with kinematic analysis. Further, numerical simulation has been done with continuum 2D elasto-plastic finite element method to acquire about failure mechanism, factor of safety (or critical SRF) and to demarcate the vulnerable portions of cut slope.

Study area
Uttarakhand state encompasses vast exposure of hilly terrain which is broadly categorized into Kumaun and Garhwal regions. The investigated area lies in Survey of India toposheet no. 53J/8 in Garhwal region. The studied cut slopes are located along NH-7 between Shivpuri to Byasi that runs parallel to the holy Ganga river valley (Fig. 1). National Highway (NH-7) passes through worst landslide prone area of Lesser Himalayan, where three locations between Byasi and Shivpuri were chosen for detailed slope stability assessment by an empirical and an advance numerical technique. The rocks of the area mainly belong to the Garhwal syncline of the Lesser Himalayan. The different meta-sedimentary litho-groups of the syncline belong to Blaini Formation (shale, siltstone, phyllite and conglomerates), Infra-Krol Formation (limestone), Krol Formation (calcareous rocks), lower and upper Tal formation (arenaceous, argillaceous, quartzite and phyllitic quartzite) ranging from Proterozoic to Cambrian in age [50]. The geotechnical parameters such as slope height with their angle, joint orientation, spacing and persistence of discontinuities, roughness, number of joints, separation (aperture), the condition of discontinuities and groundwater condition for each cut slopes were recorded carefully in the field and shown in Table 1. Further, any loosening (or dislodged block) due to blasting and excavation or tropical rainfall erosion effect due to water condition or ice wedging was also measured. The generalized views of studied cut slopes with their existing conditions are shown in Fig. 2.

Methodology
The paper focuses on the stability analysis of the three road cut sections which were selected based on their susceptibility, available in literature, landslide hazard zonation map and geological field survey. The field investigation was carried out to determine the rock masses exposed on slope surfaces, failure type, structural parameters and slope conditions. The representative rock samples were also collected from all locations, in order to determine various geotechnical properties of intact rock.

Kinematics and rock mass characteristics
The kinematic analysis is widely used to assess the potential modes of failure in jointed rock mass. The analysis is performed based on geometry of slope, material properties and angular relationships between discontinuities and slope surfaces [23,40]. The kinematic analysis of all three road cut slopes were determined using DIPS 6.0~ Rocscience Inc. [33] accordingly as Markland [27]. The rock mass rating (RMR) or Geomechanical classification was proposed by Bieniawski in (1979). Which provides quantitative, details and criterion for engineering design that refine the actual geological characteristics of rock mass. RMR system can be used to predict susceptibility of landslide and to select a method of excavation. It allows engineering geologist to demarcate the critical portion of rock mass on cut slope that could be prone to failure and accordingly suggest the required support system. RMR basic was intended from five parameters such as uniaxial compressive strength (UCS) (in MPa) of intact rock material; rock quality designation (RQD) (in %); the spacing of discontinuities (in meters); conditions of discontinuities (in meters) and groundwater condition (L/min) as Bieniawski [9]. The RQD formula was given by Palmström [29] as shown in Eq. (1).
Where volumetric joint (Jv) was estimated by the number of joints in cubic meter volume of a rock mass. In 1997, Hoek and Brown endorsed that for poor quality of rock masses (RMR < 25), RMR system was difficult to apply as RQD in most of the weak and in highly jointed rock mass is zero. The authenticity of relationship between RMR and the constants of the Hoek-Brown failure criterion begins to fail for highly fractured and weak rock masses. Therefore, Hoek et al. [19], Marinos and Hoek [25,26] developed the classification system of Geological strength index (GSI) to get over the limitations of RMR system and to provide qualitative assessment of rock masses. Further, Sonmez and Ulusay [45] modified the GSI chart of Marinos and Hoek [25,26] to obtain more quantitative numerical basis, for the evaluation of GSI. Geologically it is very sound index and their application is easy in the field. The revised GSI chart of Sonmez and Ulusay [45] include the parameters of surface condition rating (SCR) and structural rating (SR), which can be calculated on the basis of Eqs. (2) and (3) respectively.
Where Rr = Roughness rating, Rw = Weathering rating, and Rf = Infilling rating, and Jv = Volumetric joint count.

Empirical method of slope stability
The worldwide classification system Slope Mass Rating (SMR) was established by Romana, in (1985) and very useful for the preliminary evaluation of slope stability. Irigaray et al. [21] suggested that the detailed quantitative definition of the correction factors of SMR classification is one of the most important advantages over RMR system. SMR systematically describes rock mass conditions and provides some basic mandate about potentially instable modes and the required support measures. The Slope Mass Rating (SMR) is calculated by determining four correction/adjustment factors to the basic RMR [9] as in Eq. (4).
Where F 1 , F 2 and F 3 are the adjustment factors and their value rely on the interrelationship between discontinuities, affecting slope's rock mass and one additional factor F 4 is related to method of excavation.
Continuous slope mass rating (CoSMR) proposed by Tomás et al. [47] is a modified version of SMR, which rectify the discrete nature of correction factors by the continuous functions. The CoSMR allocate distinct value of each adjustment factor (F 1 , F 2 , and F 3 ) unlike a range as in SMR. CoSMR is calculated from same equation of SMR given by Romana [36]. Where, the adjustment factors F 1 , F 2 and F 3 are calculated by following continuous functions or Eqs. (5), (6) and (7) respectively.
Where A, B and C are the angles in degree and their values can be estimated from standard table given by Romana [36]. The quantitative methods RMR, SMR and continuous SMR were selected for the rock mass characterization and to suggest stable slope design. However, considerable drawback of the SMR method is to not include in the surficial discontinuity conditions and weathering effect in their final score calculation [21]. Moreover, Bar and Barton [5] stated that analytical approaches such as kinematics, limit equilibrium or finite element methods are cumbersome and it is practically impossible to assess the stability of rock cuttings and benches in real time. As in civil and mining projects, the excavation is usually too fast to perform such analysis.
Therefore, Bar and Barton [3,5] proposed the somewhat modified Q-system [6,7] known as Q slope stability, and successfully applied [4] for the stability evaluation of slopes worldwide. The Q slope is formulated by same six parameters RQD, J n , J r , J a , J w and SRF as in Q-System, but the term J w , which is now replaced by J wice , to incorporate wider range of environmental conditions (tropical intense rainfall effect, ice wedging etc.) and SRF with SRF slope which is strength reduction factor for the slope. It has the following expression as in Eq. (8).
Where RQD is the rock quality designation, Jn is the joint set number, Jr is the joint roughness number, Ja is the joint alteration number, Jr/Ja included discontinuity orientation and wedge adjustment factor (Jr/Ja) 0 , J wice is the environmental and geological condition number, and SRF slope is the strength reduction factor. Further, SRF slope factor is divided into three parts, namely, SRF a , physical condition number; SRF b , stress-strength number, and SRF c major discontinuity number.
In compare to other methods, Q slope stability classify the rock slopes in stable, unstable, and uncertain zones. It enables the engineering geologist to perform quick stability analysis and also provide the potential adjustments to slope angles for different probability of failure (PoF) to make stable slope design. This is distinguishable in field during excavation. The following Eq. (9) (For Probability of failure, PoF = 1%) was used for adjustable slope angle calculation.

Numerical slope stability analysis
Generally, the rock mass has inhomogeneous, intrinsically discontinuous, and anisotropic characteristics in nature [22]. In such media, analysis of deformation characteristics, mechanical behavior and failure mechanism is quite limited with deterministic traditional limit equilibrium method. To accomplish the limitations, numerical modelling has been applied in field of rock engineering to provide approximate solution to problems [12,14,46]. Finite element method (FEM) is most frequently used continuum numerical techniques due to their aptness to handle complex behavior of weak jointed and or heavily fractured rock mass problems [12,48]. Unlike limit equilibrium (LEM) method, FEM does not rely on any presumed location of slip surface. FEM divides the material model into several small zone and utilizes Shear strength reduction (SSR) techniques in combination with material properties to estimate factor of safety and potential failure mechanism. Recently, many personnel have attempted numerical methods to simulate the behavior of jointed/blocky rock masses and soils in Himalayan condition [15,24,28,42] and suggested the potential remedial measure for the long-term stability.
In present study, FEM analysis was performed using elasto-plastic plane strain simulator 'Phase 2 8.0' ~ Rocscience Inc. [34]. FEM is more applicable and has accessibility to handle complex problem under different conditions. To account several geotechnical parameters in analysis, material model such as Mohr-Coulomb and Hoek and Brown failure criterion have been well incorporated with FEM for the design of diverse slopes and underground excavations. Hoek [16] and Hoek et al. [18] established a good correlation between Hoek-Brown and Mohr-Coulomb failure criterion. Further, Hoek et al. [18] formulated an equation based on Hoek-Brown failure criterion and equivalent normal and shear stress of Balmer [2] to estimate equivalent Mohr-Coulomb friction angle and cohesion. In present analysis, Mohr-Coulomb failure criterion was applied for exact determination of failure mechanism and to demarcate susceptible portion of cut slope sections for instability.

Kinematics and rock mass characteristics
Stereonet plot (equal angle projection) of each slope illustrates mainly the characteristics of double plane sliding i.e. wedge (W) failure for all locations. Along with wedge failure, also potential of flexural toppling (FT) failure at location L2 was prompted (Fig. 3). These predicted failures are mainly structurally controlled, where joint J 2 is critical joint set to initiate mass failure on cut slopes. RMR basic was intended from five different rock-mass parameters based on recorded field and laboratory estimated data for all three locations. Mean values of UCS ranges from 31.8 to 36.5 MPa for phyllitic quartzite. Rating of different parameters for the RMR basic with their result are shown in Table 2, where locations 1 and 3 are under 'fair' category, while location 2 has 'good' category of rock mass. Table 3 indicate quantified GSI value, which was estimated by means of two parameters, i.e. 'surface condition rating' (SCR) and 'structure rating' (SR) as plotted in Fig. 4. The yellow color square as shown in Fig. 4 represent the quantified GSI value. The value of GSI indicate that all locations approximately fall in categories of blocky structure with fair surface condition. However, GSI value of locations 2 and 3 are slightly higher than location 1 due to blockier nature and high spacing between discontinuities. Low value of   GSI at location 1 may be due to extensive jointing (3 to 4 set) and poor surface condition of rock mass.

Empirical slope stability analysis
To calculate SMR, the rating has been assigned for each type of failure, which include the adjustment factors F 1 , F 2 , F 3 and F 4 as shown in Table 4. The result of SMR shows that location 1 is unstable while locations 2 and 3 are in partially stable condition. Similarly, continuous slope mass rating (CoSMR) was calculated according to Eq. (4) and the result shows that stability class is almost similar to the stability class of SMR (Table 5). Generally, the maximum vulnerability is obtained by CoSMR, as its final score is found to be less than the scores of discrete SMR. Moreover, SMR and CoSMR methods are most likely to give a preliminary assessment about instability condition and required support   (Table 6). Estimated Q slope values are shown in Table 7. Figure 5 shows that plot of Q slope vs. average slope angle and β angles of the investigated slope accordingly Bar and Barton [3]. Stability condition of each cut slope is marked with violet polygon, in data chart, where location 1 is unstable and rest of locations are critically stable (Fig. 5).
Bar and Barton [3] recommended that based on Q slope value, it is possible to reduce the slope angle and change in design of cut slope during excavation process based on different probability of failures. β angles for probability of failures 1% is calculated using Eq. (9) and shown in Table 7. Further, the Q slope value with average slope angle (red star) and β angles have also been plotted (violet star) for probability of failures 1%, to improve the probability and stability of cut slopes (Fig. 5).

Numerical simulation
In-situ uniaxial compressive strength (UCS) of representative rock samples were estimated in the lab from universal testing machine (UTM) as per ISRM suggested methods [49]. The Roclab program [35] was used to determine rock mass properties, where input parameters are taken as-estimated GSI (Table 3), mean value of UCS (intact rock), unit weight (0.027 MN/m 3 ), material constant and slope height. Table 8 shows the determined data for representative rock samples to assess the numerical   simulation. To analyze the investigated slopes, two dimensional plane strain condition has been considered with Gaussian elimination solver to resolve the equations. Steepest 2D profile for a slope has been incorporated in the model so that the profile with maximum vulnerability can be analyzed. The problem domain (2D slope profile) has been discretized and meshed through six nodded triangular elements to incorporate all the geometrical complexities. The base and the left boundary has been restrained in both X and Y direction as boundary conditions. The slope surface towards the road has been kept free and gravity loading, body forces and field stress were employed as initial element loading. The extracted displacement contour and its distribution pattern across FE model represents the deformation intensity in different parts of the cut slope sections. The critical SRF as from numerical simulation result implies that slopes at location L1   Fig. 6a-c represents the variation of total displacement, differential stress and maximum shear strain along the free face (apex to toe) of cut slope surfaces in the susceptible zone of failure. Further, the maximum value of these parameters was found at the toe of slope section and it was also being corroborated by field evidences.

Analyses and discussion
Bieniawski [10] suggested that it is essential to identify most critical condition of discrete geological features that mostly govern the stability of rock slopes. Whereas, Bieniawski [9] and Hoek and Brown [17] suggested that GSI can be correlated with the RMR and modified rock mass quality index Q. The result of RMR and GSI suggest that slopes has almost fair to good quality of rock with blocky structure and fair surface conditions. The calculated GSI value indicate that all locations approximately fall in categories of blocky structure with fair surface condition. However, GSI value of locations 2 and 3 are slightly higher than location 1 due to blockier nature and high spacing between discontinuities. Low value of GSI at location 1 may be due to  (Tables 4 and 5). Generally, the maximum vulnerability is obtained by CoSMR as its final score is found to be less than the scores of discrete SMR. Moreover, SMR and CoSMR methods are most likely to give a preliminary assessment about instability condition and required support measures. The average of these two rating should be considered in the design of support systems for any forthcoming project. Further, the result of SMR and CoSMR techniques is also validated from Q slope stability classification system, where all the empirical results express a well correlation to each other. In numerical simulation result, slope 1 and 3 are in unstable condition and slope 2 is critically stable. In Fig. 7 at location L1, the displacement contour of almost equal magnitude distributed along the entire length of slope face. It indicates that maximum failure zone has wide pattern near crest and gradually taper near toe of the slope, hence anticipating slight deeper zone of failure. Based on field observation and FEM result, it can be predicted that slope may undergo wedge (estimated from kinematic analysis) or possibly curvo-planar type of failure due to extensive jointing and smaller block size. At location L2 in Fig. 7 indicate that shallower damage zone is lying from mid portion to toe of slope face, clearly exhibited by maximum value of total displacement contour. Here, the slope angle (65°) is excavated in such a way that wedges are freely hanging on the roadside and discontinuity orientation across slope face also favours toppling failure (from kinematic analysis) (Fig. 3). The field data depict that at location L2, persistence of bedding joint is quite high (1-20 m) and excavation activity for further road development continuously disrupt rock masses. Due to which persistent joints are expected to get exposed, therefore enabling kinematic feasibility [13] across the slope face and anticipating thick block failure. Displacement contour pattern at location L3 (Fig. 7) shows that the most vulnerable part is distributed near middle portion of the free face which imply possibility of shallow zone of failure. However, due to high steepness of slope (75°) and intermittent joint and outward projected displacement contour with respect to slope face show potential for step path failure [20]. Further the impact of gravitational loading will be slightly higher for location 1 and 3 (H ~ 55 m) than location 2 (H ~ 50 m). Overall from empirical and numerical simulation analysis, the slope condition in present are critically stable to unstable i.e. vulnerable condition.

Conclusion
In present study, three road cut slopes have been investigated for their rock mass characteristics and stability analysis. The cut slopes are located from Shivpuri to Byasi along National Highway (NH)-7 in Pauri Garhwal region of Uttarakhand. Kinematic analysis showed wedge mode of failure in all locations with one toppling at location 2, which has also been noticed during field survey. The rock mass characteristics assessed by RMR basic and GSI technique are almost similar result and revealed that rocks are good and fair quality, with blocky structure and fair surface condition. The discrete and continuous SMR techniques also showed that slopes are partially stable with almost similar stability class. To supplement the SMR and CoSMR result, Q slope stability of cut slopes were assessed, that indicate that slope at location 1 is unstable while location 2 and 3 are partially stable, which also confirm from the result of SMR and CoSMR techniques. Further, to stable these slopes, β angles have been derived for probability of failures 1%. The numerical simulation (FEM) resulted that location 1 &3 are unstable, while location 2 is critically stable. The analysis showed that shallow to intermediate zone of failure with maximum total displacement, differential stress and shear strain found near the toe of slope section. However, the slope at location 3 is partially stable under SMR/CoSMR and Q slope system, whereas unstable in FEM model. All inclusive, the FEM and different empirical results revealed that studied road cut slopes between Shivpuri to Byasi in Uttarakhand are critically stable to unstable. Therefore, towards the stability, the sealing of discontinuities is required to install the grouting along the discontinuities surfaces. Thereafter, nets of the desired mesh or systematic bolting can be applied at vulnerable portions of slope face (as discussed) to minimize the potential threat such as any incidental rock fall or block failures.