Stability assessment of landslide-prone road cut rock slopes in Himalayan terrain: A finite element method based approach

Srd Prsd Prdhn, Triq Siddique,b,*

a Department of Earth Sciences, Indian Institute of Technology Roorkee, Roorkee, 247667, India

b Department of Geology, Aligarh Muslim University, Aligarh, 202002, India

Keywords:Landslides Numerical modeling Finite element method (FEM)Slope stability

A B S T R A C T Large-scale slope destabilization could be aggravated due to swift urbanization and ever-rising demands of geoengineering projects such as dams, tunnels, bridges and widening roads. National Highway-58 connects Delhi to Badrinath in India, which passes through complex geomorphological and geological terrain and often encounters cut slopes susceptible to slope failures. In the present investigation, a detailed geotechnical appraisal is conducted along the road cut slopes from Rishikesh to Devprayag in the Himalayas.Twenty vulnerable road cut slopes were demarcated for detailed slope stability analysis using Phase2D finite element modeling simulator. Nonlinear generalized Hoek-Brown (GHB) criterion was adopted for stability analyses. Out of 20 slopes, five slopes (S6, S7, S18, S19 and S20) are unstable with factor of safety (FoS) less than or equal to 1, and thus needs immediate attention. The FoS values of four slopes(S2,S9,S13 and S17)lie between 1 and 1.3,i.e.marginally stable,and slopes S1,S3,S4,S5,S8,S10,S11, S12, S14, S15 and S16 are stable. Mohr-Coulomb (MC) criterion was also adopted to compare the slope stability analysis with GHB criterion. The FoS calculated from GHB criterion is close to that using MC criterion for lower values of FoS whereas for higher values, the difference is marked.For the jointed rock in the Himalayan region, the nonlinear GHB criterion gives better results as compared to MC criterion and matches with the prevailing field conditions. Accordingly, some suggestions are proposed to strengthen the stability of cut slopes.

1. Introduction

Landslide involves downward and outward movements of slope-forming materials due to gravitational force by a variety of motions like falling, sliding, flowing and any combination of the above (Cruden and Varnes, 1996). Frequent slope failure along natural and engineered slopes is of significance as it threatens the lives, hampers the socio-economic growth, and deteriorates the habitat (Leroueil, 2001). Mountainous terrain like Himalayas is characterized by a variety of landslides.In such precarious regions,local residents, tourists and pilgrims are under sustained threat of various natural hazards like seismicity, landslides, floods, cloudburst, and forest fires. The anthropogenic intervention to the extensively fragile and complex behavior of Himalayan terrain is fundamental for occurrence of large- and small-scale landslides in this region.The stability of road cut slopes in the region degraded notably due to inadequate excavations during road construction and widening projects(Siddique et al.,2017).Extensive mechanical excavation and improper blasting generate secondary fractures within the rock mass. However, it has been recommended that in such problematic terrain conditions, controlled blasting should be performed (Mondal et al., 2016).

Enormous researches on the slope stability in precarious Himalayan terrain have been undertaken (e.g. Devkota et al., 2013;Kundu et al., 2016; Kumar et al., 2017; Sah et al., 2018). The consequences of landslides became much hazardous if they occur along transportation corridors or near townships. Uttarakhand state of India is well known for its ecological diversity and richness along with massive tourism and religious activities.Natural hazards like frequent seismicity, cloudburst, landslides, floods, lightning,and forest fires are of the major concern among authorities and stakeholders.The occurrence of any such event is catastrophic and the impacts are manifold when one event triggers another. Many such incidences have been reported over the past years, for example the one occurring in Kedarnath valley of Uttarakhand in 2013. The area witnessed extremely heavy rainfall with massive cloudburst due to the fusion of westerlies with the Indian monsoonal cloud system (Nair and Singh, 2014) that caused massive flash floods and consequently large-scale landslides(Dubey et al.,2013).As per Map the Neighbourhood in Uttarakhand(MANU) reported by Ahmad et al. (2015), Wadia Institute of Himalayan Geology (WIHG) investigated the damage occurrence within the Bhagirathi-Ganga-Nayar valley during the disaster. The result indicates that 1034 landslides were either initiated or reactivated during the calamity. The study being conducted here is a part of this valley.

National Highway-58 (NH-58) has exceptional importance due to tourism and pilgrimage activities and acts as an artery for socioeconomic development of the region. For restoring environment with sustainable socio-economic development, routine geotechnical assessment is prerequisite. Several highways are under continual threat of a variety of slope failures.From several decades,frequent mass movement phenomena are the major hindrance to ongoing traffic. Enormous endeavor has been made by various geotechnical researchers to assess the stability grade of cut slopes along NH-58 in Uttarakhand (Siddique et al., 2015; Vishal et al.,2017). Slope stability analysis along NH-58 has also been performed by finite element (FE) modeling technique (e.g. Kanungo et al., 2013; Pain et al., 2014; Sarkar et al., 2018). These case studies highlighted the cause of the problem and some of them also suggested remedial measures to reinforce the cut slopes but still,there are certain zones along the highway that require comprehensive assessment. Road cut slopes from Rishikesh to Devprayag along NH-58 are highly vulnerable which were evaluated by rock mass classification and kinematic tool by Siddique et al. (2017).They marked many potential zones along the transportation corridor and also suggested general guidelines to enhance the stability of slopes. Furthermore, stability analysis of debris cut slopes in the same stretch has been undertaken by deterministic and sensitivity approach (Siddique and Pradhan, 2018).

In the present investigation, highly landslide-prone area along the road section from Rishikesh to Devprayag was taken for detailed slope stability assessment by an advanced numerical tool(Fig. 1). Due to highly fractured conditions of the rock mass,curvilinear generalized Hoek-Brown (GHB) failure criterion has been adopted for the study. Critical factor of safety (FoS) was calculated for each slope and accordingly remedial measure is suggested. Moreover, the widely used linear Mohr-Coulomb (MC)criterion was also used to compare the results.

2. Study area

Uttarakhand comprises hilly terrain with rugged topography that encompasses 13 districts which are categorized into two units:Garhwal and Kumaon regions. Differential meteorological phenomena and extreme weather conditions in the state are reported due to its geographical location and subtle geomorphology.It is also characterized by extensive glacier ice caps in upper regimes.Along with some seasonal and perennial rivers, two major rivers, i.e. the Ganga and the Yamuna, originate from these extensive cover of glaciers. Beautiful landscapes, unique ecosystem and holy shrines attract tourists and pilgrims across the globe.

2.1. Geological framework of the area

The Himalayan orogenic belt is the ramification of continentcontinent collision of Indian and Eurasian plates (Gansser, 1964;Dewey and Burke,1973; Yin, 2006). The Himalayan orogenic belt is bounded by some major structural features, i.e. Indus-Tsangpo suture zone (ITSZ) in the north, Chaman fault in the west, Sagaing fault in the east and main frontal thrust(MFT)in the south(LeFort,1975; Yin, 2006). The Himalayas is subdivided into sub-Himalaya(Siwaliks), Lesser Himalaya (Himachal), Higher Himalaya (Himadri) and Tethys Himalaya (Gansser,1964; LeFort,1975).

Fig.1. Investigated road cut slopes on the geological map of study area (modified after Siddique and Pradhan, 2018).

According to Geological Survey of India toposheets, the study area lies in 53J/8 and 53J/12, i.e. in the Lesser Himalaya which is predominantly comprised of meta-sedimentary sequence (LeFort,1975) along with some volcanics, meta-volcanics and gneiss(Frank et al., 1995). The Lesser Himalayan sequence has been formed due to contraction in multiple phases and structurally bounded by main boundary thrust(MBT)in the south and by main central thrust (MCT) in the north (Valdiya, 1983). Furthermore,extensive work of certain researchers proposed that the geological formations of the Lesser Himalayan sequence can be categorized into two major subdivisions, i.e. Garhwal and Kumaon Lesser Himalaya (Valdiya, 1980). The study area lies under Garhwal Himalayas.It comprises Damtha,Tejam,Jaunsar,Mussoorie,Sirmur,Ramgarh and Almora groups of rocks.These groups of rocks can be further divided into outer and inner Lesser Himalayan sequences which encompasses different formations(Valdiya,1980;Yin,2006).The investigated slopes lie along NH-58 in Garhwal syncline of the outer Lesser Himalaya, i.e. running parallel to the holy river the Ganga (Fig. 1). The area comprises of meta-sedimentary rocks(shale, siltstone and conglomerates of Blaini Formation; limestone of Infra-Krol Formation;calcareous rocks including limestones and dolomites of Krol Formation; argillaceous, arenaceous, siliceous and calcareous rocks of Tal Formation; quartzite of Nagthat Formation; and sandstone of Chakrata Formation) of Proterozoic to Cambrian in age(Valdiya,1980; Jiang et al., 2003).

Seismic vulnerability and its ongoing tectonism are important components to understand the behavior of the rock mass. Such continuous adverse activities steadily decreased the inherent strength characteristics of the rock mass and an extreme event may trigger large-scale landslides in the region. As manifested by frequent seismic events, the entire Himalayan mountain chain is inherently fragile. Relatively, greater Himalayan sequence witnessed much frequent seismicity as compared to the Lesser Himalaya. Interestingly, the Lesser Himalayan sequence is extensively deformed terrain that remained calm for a long time in the aspect of higher magnitude earthquake. It is quite notable that cumulative and progressive nature of stress accumulation within the Lesser Himalayan sequence is due to extensive horizontal stresses owing to sub-surface Delhi-Haridwar, Faizabad and Monghr-Saharsa ridges that are underlying extensions of Delhi-Aravalli, Vindhyan and Satpura rocks, respectively (Valdiya,1992).This could be one of the major contributors accounting for extensive seismic events that are capable of triggering substantial landslides. Valdiya (2002) assessed the major events during the evolution of the Himalayas and proclaimed that spasmodic rise of marginal Lesser Himalayan sequence and suggested that extensive contraction in the Siwaliks range at 1.6 million years ago caused rampant landslides.Sati et al.(1998)studied the impact of regional tectonic setting on a variety of landslides and suggested that most slopes prone to failure are parallel to the regional trend of the Himalaya, i.e. WNW-ESE. It is due to the factor that the extensive horizontal stress is perpendicular to the regional trend of the Himalayan orogeny. Consequently, the sections parallel to this trend would be sheared and shattered for greater distance, thereby forming a significantly large unstable zone.Mithal(1988)assessed litho-tectonic landslides in Garhwal-Kumaon region of the Lesser Himalaya and suggested that dynamic forces due to developmental activities in the form of excavation are continuously endangering the geo-environmental conditions, particularly along NH-58,Rishikesh to Badrinath highway.

According to Kumar and Dhaundiyal (1980), the rocks in the area had experienced significant crustal stresses that resulted in fan folding due to which the central portion of Garhwal synform represents a doubly plunging anticline.Among the major thrust faults in the Himalayas,MBT is the most proximal to the studied section.Being the plane of underthrusting of Indian plate under the extensive Himalayas,MBT is geodynamically active(Valdiya,1983).Apart from MBT, there exist few major faults in the Garhwal synform itself. These are Pulinda, Bedasini, Fatehpur, Bonga, Singtali,Duwadhar and Maidan faults (Kumar and Dhaundiyal, 1980;Valdiya, 1980; Sati et al., 2011). Among these, Duwadhar and Singtali faults intersect road cut slopes at distinct locations along the highway.However,care must be taken as Singtali fault is being named by Garhwal and Binj thrust in the literature.Singtali thrust often cross-cuts the road section several times due to which slopes near Kaudiyala are very much fragile. Such adverse geological and structural conditions are responsible for inherent degradation in the stability of cut slopes along the highways. Such hazard-prone zones must be evaluated on a routine basis to ensure better safety and swift sustainable development of the region.

3. Methodology

The vulnerable slopes along NH-58 from Rishikesh to Devprayag can be identified in available literature, landslide inventory and reconnaissance field survey.A detailed field survey was conducted to record geometrical, geological, structural, geotechnical, and hydrological parameters with respect to slope stability analysis.Furthermore, laboratory experiments were carried out for determination of geomechanical properties of the intact rock. The data obtained from the field were integrated with the laboratory results which were used to perform numerical modeling for stability analysis.

3.1. Numerical modeling

Numerical modeling technique is widely used in resolving various issues related to geotechnical engineering projects. Swift progression in computational efficiency enabled the researchers to understand the geomechanical response of slope forming material under static and dynamic loadings. In a broader sense, numerical methods can be classified as continuum,discontinuum and hybrid methods. In continuum method, the whole material is treated as continuous, i.e. uniformly distributed throughout the slope;whereas in discontinuum modeling, the slope forming material is treated as heterogeneous mass(Jing and Hudson,2002;Jing,2003).Continuum methods include finite element method (FEM), finite difference method (FDM), finite volume method (FVM), boundary element method (BEM) and meshless methods, while discontinuum methods include discrete element method (DEM) and discrete fracture network method (DFM) (Jing and Hudson, 2002;Jing, 2003). According to Jing (2003), hybrid method is applied to understanding the flow and stress-related problems. The use of above-mentioned numerical methods depends upon site-specific conditions (Jing, 2003). In the present study, continuum modeling has been used, because it is probably best suited for weak and jointed rock mass.

3.2. Finite element modeling (FEM)

FEM is an elasto-plastic method widely accepted in geotechnical research area. It does not require any pre-assumption for the location of slip surface and the outcomes are nonlinear and iterative in nature (Griffiths and Lane, 1999; Rocscience, 2001). It is much economic and time-effective with iterative capability which makes it a much reliable method to make the best possible and efficient solutions to the problem. In FEM, failure occurs naturally along the zone where prevailing shear stresses overcomes the shear strength of the material without any pre-assumption(Griffiths and Lane, 1999). For slope stability assessment via numerical methods, FoS is probably the easiest and quickest way to predict the risk due to failure.Such quantitative assessment enables user or reader to have a quick overview at a glance. It may be termed as the ratio of the actual shear strength of the slope forming material to the minimum shear strength required to resist failure.In shear strength reduction(SSR),large numbers of simulations are performed for a series of trial of FoS (Matsui and San,1992). The shear strength of the material is reduced during successive iterations until the failure occurs. In each trial, shear strength parameters,viz.cohesion(c)and angle of internal friction(φ),are reduced as per Eqs.(1)and(2).Thus,critical strength reduction factor(SRF)is computed which is equivalent to FoS.

where c is the cohesion, φ is the friction angle, and f is the safety factor.strength (UCS) of the intact rock; and m and s are the material constants.

The parameter m is equivalent to the friction of the rock while s is related to the degree of fracturing that relates to the cohesion of the rock mass (Eberhardt, 2012). Since the development of the HB criterion, it has been updated and refined several times to reduce the limitations from the gained experience. Some of the major attempts for modification have been discussed and compiled by Hoek and Marinos(2007).The concept of GHB criterion was introduced.It involves the replacement of rock mass rating(RMR)by geological strength index (GSI) system and also the concept of disturbed and undisturbed rocks was removed that allows the users to reduce the GSI values by careful judgment of prevailing field conditions.Along with this, several terms mb, s and a were introduced for good and poor quality rock masses having GSI values ＞25 and ＜25, respectively. The GHB criterion is given as

3.3. Failure criterion of rock mass

In highly fractured rock masses, the failure surface principally runs along discontinuities and partially through the intact rock.At low normal stresses, individual fragment or block may move or rotate due to low cohesion;but at high normal stress levels,friction is diminished as a response of crushing.As compared to intact rock,the failure mechanism in the jointed rock mass is very complex.The failure in such jointed rock slopes may occur by shearing of the rock mass, sliding of blocks along discontinuities and/or rotational and translational movement of individual intact rock blocks.The major and widely used failure criteria for numerical modeling of rock mass are MC criterion (Coulomb,1776; Mohr,1900), Hoek-Brown(HB) criterion (Hoek and Brown, 1980), Ramamurthy criterion(Ramamurthy et al., 1993), Ramamurthy-Arora criterion(Ramamurthy and Arora,1994), GHB criterion (Hoek et al.,1995),and modified MC criterion(Singh and Singh,2012).Coulomb(1776)assumed that shear strength of rock is the function of cohesion and angle of internal friction and failure envelope is linear. But the failure mechanism changes continuously in jointed rock mass and the shear strength envelope is curvilinear,especially at low range of normal stresses. Therefore, it is inadequate to directly use geomechanical properties of intact rock for any rock slope design.The degree of fracturing and properties related to the joint surface must be taken into consideration.From the numerous experimental data,it was observed that the failure envelope in jointed rocks is curvilinear rather than straight line and the curve shows concavity towards the axis of normal stress(Hoek et al.,2013;Singh,2019).HB failure criterion for rock mass is based upon several empirical relationships that characterize the stress conditions related to the failure of intact rock and rock mass. It has been derived from the crack theory of Griffith (1920, 1924) in brittle rocks (Hoek, 1968)and subsequently modified as per the field and laboratory observations (Marsal, 1967, 1973; Jaeger, 1970). Unlike the MC linear criterion,HB is an empirically derived criterion which relies on the nonlinear increase in peak shear stress with confining pressure(Eberhardt,2012).In numerical modeling of the jointed rock mass,nonlinear shear strength envelopes have been ignored for a long time (Barton, 2013). The original HB failure criterion was introduced by Hoek and Brown(1980)(Eq.(3))in an attempt to develop an empirical relationship that could be scaled in relation to geological data (Hoek and Marinos, 2007).

where mb= miexp[(GSI 100)/28, s = exp[(GSI 100)/9and a = 0.5 for GSI＞25(good quality);s = 0 and a = 0.65 GSI/200 for GSI＜25(poor quality).miis the material constant that depends upon the type of rock,mbis the reduced value of miwhich accounts for the strength by reducing effects of jointed rock mass that relies upon GSI values and disturbance factor, and s and a are the curve fitting parameters which are being determined by using GSI and D values.

Several changes were incorporated in the above criterion by Hoek et al. (2002) and the concept of disturbance factor (D) was introduced to consider the damage caused due to blasting.Further,it was conceived by researchers that there exists a hiatus between good and poor quality rock masses that was proposed in GHB criterion of 1995 edition.By considering this fact,an attempt has been made for smoother transition among good and poor quality rock masses and the equations of GHB criterion (1995 edition) were updated by Hoek et al.(2002).Probably,GHB(2002 edition)is the most updated version of HB criterion derived. The GHB failure criterion (2002 edition) for slopes is written as

where σ′cmis the UCS of rock mass,γ is the unit weight of the rock mass, and H is the height of the slope.

Hoek and Diederichs(2006)proposed an empirical relationship to calculate the deformation modulus of rock mass (Erm). It is designated as generalized Hoek-Diederichs criterion:

Furthermore, it has been realized that many numerical modeling tools do not incorporate HB criterion,and they use simple MC criterion.To overcome this problem,a window based program called ‘RocLab' was also developed to calculate equivalent shear strength parameters(Hoek et al., 2002).

3.4. Failure criterion of joints

where σ′1and σ′3are the major and minor effective principal stresses at failure, respectively; σciis the uniaxial compressive

In jointed rock mass, stability is greatly influenced by discontinuities especially at shallow depth or at low stress levels. In such slopes, failures particularly shear failures tend to occur along the zone or plane of the least resistance, i.e. discontinuity. The shear strength along such planes primarily relies upon shape and roughness of asperities, degree of alteration, matching along the fractured surface, type and thickness of infilling material (Hoek,2006). The well-known model for estimating the shear strength of discontinuities is the MC criterion:

τf= cj+σntanφj(7)

where τfis the shear strength at failure,cjis the cohesion of joint,σnis the effective normal stress, and φfis the friction of joint.

However, there are certain limitations in MC criterion. The failure envelope is linear and also the tests are too expensive to be conducted on each site,and for shear testing,jointed surface is to be allowed to fail. Patton (1966) recognized the importance of roughness and performed a series of direct shear tests on sawtooth triangular joints. As a result, he suggested the bilinear failure criterion.The effective normal stress acting on a particular discontinuity surface depends upon its orientation, depth, and weight imposed due to overburden along with density and hydrological conditions (Hoek, 2006). According to Barton (1973)and Barton and Choubey (1977), dilation induced roughness significantly contributes to the nonlinear behavior among normal stress (σn) and shear stress (τ). Roughness of potential sliding surface can be termed as first order and second order asperities.The former are major undulations and measured at large scale and later are small ripples and bumps at a very small scale. However,Barton (1973) suggested that the importance of first and second order asperities largely depends upon the magnitude of the normal load. At low normal stress levels, second order asperities are significant while first order asperities play a major role only at high normal stresses. By considering all the major aspects, the shear strength of discontinuity surfaces within a jointed rock mass is the coupled effect of surface irregularities or asperities,strength, normal stress, and shear displacement along the potential sliding surface (Wyllie and Mah, 2004). Barton (1973)proposed an extension of MC and Patton's models by incorporating sliding and shearing simultaneously. He replaced the constant friction angle by a function of normal stress,roughness,and strength of sliding surface. The Barton's shear strength criterion can be expressed as

where φbis the basic friction angle on unweathered surface;JRC is the joint roughness coefficient which can be estimated by comparing measured roughness with standard profiles by Barton and Choubey (1977) or it can be determined by measuring the amplitude of roughness by straight edge method (Barton, 1982);and JCS is the joint wall compressive strength that can be estimated by Schmidt hammer rebound test by Miller (1965) or by the standard field method suggested by International Society for Rock Mechanics (ISRM,1981).

Some researchers conducted point load test suggested by Broch and Franklin(1972),but it may mislead the results as often strength is slightly lower along the joint surface particularly in weathered joints.Thus it was recommended to conduct Schmidt hammer test in the field to obtain much more realistic results.While performing Schmidt hammer test,the care must be given to the orientation of hammer during rebound measurements,and corrections suggested by Barton and Choubey (1977)must be applied.

Furthermore,the stiffness of rock joints defines the deformation under both normal and tangential loads.The normal stiffness may be defined as the normal stress per unit closure of the joint while shear stiffness of a joint is the ratio of peak shear stress to the shear displacement (Barton,1972).

The normal stiffness of joint is represented as (Barton,1972):

where Knis the normal stiffness of joints,Emis the modulus of rock mass,Eiis the modulus of intact rock,and L is the mean spacing of joint.

The shear stiffness of joint is represented as (Barton,1972):

where Ksis the shear stiffness of joints,Gmis the shear modulus of rock mass, and Giis the shear modulus of intact rock.

3.5. Geological strength index (GSI)

GSI system is a major input in numerical modeling and has been used to calculate the deformation and strength parameters of the rock mass. This classification tells about the rock mass and does not count the orientation of discontinuities.Nevertheless, it helps to overcome the restrictions of expensive and timeconsuming laboratory investigations of intact samples only. It has been widely used for estimating strength characteristics in a large number of international tunneling projects. GSI values depend upon the blockiness of the rock mass and existing conditions of discontinuity surfaces. Although GSI values can be estimated by existing relationships among various rock mass classification (RMR and Q system), for much accurate results, it has been recommended to use published GSI chart (Hoek, 2006).Although Hoek (2006) suggested that single value of GSI should not be assigned, it is good to consider a range of GSI values, but as the GSI is directly linked to the empirical equations of GHB criterion and deformation modulus of rock mass,the subjectivity of the system ought to be removed. To surmount these limitations, Sonmez and Ulusay (1999, 2002) proposed two parameters, i.e. structure rating (SR) and surface condition rating (SCR).SR depends upon the degree of fracturing or blockiness of the rock mass at a particular site. It has been calculated by the empirical relationship among SR and volumetric joint count (Eq.(11)) and SCR relies upon roughness, weathering and infilling material.This parameter is being calculated by the algebraic sum of ratings (Eq. (12)).

SR = 17 lnJv+79.8 (11)

SCR = Rr+Rw+Rf (12)

where Jvis the volumetric joint count,Rr is the roughness rating,Rw is the weathering rating, and Rf is the infilling rating.

4. Slope stability assessment along NH-58

In the present study,slope stability assessment of road cut rock slopes was conducted by coupling various parameters and integrating a variety of information from the geotechnical field and laboratory investigations. During reconnaissance stages of the investigation, an initial field survey was conducted along NH-58 from Rishikesh to Devprayag, to demarcate hazardous zones, and 20 vulnerable cut slopes have been identifeid for detailed geotechnical assessment (Fig. 2). A variety of meta-sedimentary rock formations belonging to different geological formations have been illustrated in Table 1.

Fig. 2. Field photographs depicting the condition of investigated slopes along NH-58.

The investigated road cut rock slopes are heterogeneous and dissected by 3-4 sets of discontinuities. It was evidenced during field surveys that failures in such slopes are controlled by unfavorably oriented discontinuities. Recurrent nature of structurally controlled failure in the form of wedge and planar failure was observed at Slope S13 (Fig. 3). Pre- and post-failure conditions at Slope S13 have been depicted in Fig. 3a, b and d in 2016,2017 and 2018, respectively. A massive failed block of 5-6 ft(1 ft = 30.48 cm) is shown in the inset view (Fig. 3c) and even much larger block has been witnessed in successive failure(Fig. 3d). Slope S13 lies in proximity to thrust zone near Kaudiyala. Massive structurally controlled failures in such slopes indicate the validity of results obtained by the present numerical modeling. Among 20 investigated slopes, there are many other critical slopes which are likely to experience similar mass failures in the near future. Apart from extensive or mass failures,small block failures in the form of wedges are significantly common in the investigated area. It is often noted that even though the overall stability of any slope is reasonably fair but such small wedges may form occasional rock falls due to the presence of multiple sets of discontinuities. Such issue was propounded at Slope S15, where the overall slope is stable but small-sized wedges have often destroyed the road, roadside garders and wall (Fig. 4). The site for wedge initiation (encircled)at S15 is depicted in Fig. 4a and the impact of such small block failures is illustrated in Fig. 4b,c.

Table 1 Lithology and geological formation at corresponding studied slopes.

The geological and geotechnical field data were collected for characterizing identified slopes via GSI system.The GSI values of 20 slopes are calculated (Table 2) and also illustrated in GSI chart(Fig. 5).

Fig.4. Impact of occasional wedges on road at Slope S15 near Kaudiyala:(a)Encircled portion highlighting the zone of wedge initiation; (b) Damaged roadside garders and walls; and (c) Inset view showing the damage to the road and associated structures.

Geomechanical properties required for numerical modeling were determined during field and laboratory investigations(Table 3). From each slope, undisturbed and representative rock chunks were collected and NX-size (54.7 mm in diameter) core samples were prepared to determine the UCS as per the guidelines suggested by ISRM (Bieniawski and Bernede, 1979). The laboratory investigation also incorporates determination of unit weight of intact samples as per the standards of Bureau of Indian Standards codes (IS 1122,1974). The rock constant (mivalue) for each slope was determined by identifying rock type and thin section under a petrological microscope. Disturbance factor was determined in the field as per the guidelines suggested in the GHB criterion. Poisson's ratio plays a significant role in elastic deformation of rocks and rock masses subjected to static or dynamic loading. Among many widely used mechanical properties of rocks in rock engineering practices, the importance of Poisson's ratio has not been appreciated much due to very narrow range among different rock materials (Gercek, 2007). Laboratory experiments on the Poisson's ratio of rocks are rarely conducted (Vásárhelyi and Kovács, 2017). There are several indirect methods available in the literature to estimate the Poisson's ratio of rocks by using UCS values and RMR classification system. Furthermore, on the basis of GSI values and rock material constant (mi), Vásárhelyi (2009) proposed a chart to calculate Poisson's ratio of rock mass (νrm). This chart based approach has been used in the present study to determine the Poisson's ratio of rock mass(see Table 3).Schmidt hammer hardness test is

Fig.3. Field photographs depicting recurrent failure at Slope S13 near Kaudiyala:(a)Pre-failure condition of slope with encircled probable zone observed during initial field surveys and measurement in 2016;(b)Structurally controlled mass failure observed during successive field survey in 2017;(c)Inset view showing the massive failed blocks;and(d)Much larger failure occurring in 2018.

Table 2 Quantification of GSI for studied cut slopes.

Note: SR = 17.5 lnJv+ 79.8.a non-destructive and in situ technique which is often employed in many rock mechanics and rock engineering practices. It has been widely used to determine several properties such as elastic modulus, hardness, surface smoothness and strength of rocks. In the literature, there are many empirical relationships between UCS and Schmidt hammer rebound values(Hr).Schmidt hammer rebound values for each joint set have been measured during field survey (Fig. 6b) as per the guidelines suggested by ISRM(1981). From each slope, ten Schmidt hammer rebound values for each joint set were taken and the median Hrvalue was used to determine joint wall strength. Then, the equation proposed(Eq. (13)) by Katz et al. (2000) was used to calculate UCS(Table 4) of each joint wall, because the equation was proposed on more or less similar lithology with an appreciably good regression of 0.92. Moreover, Hrand UCS of Eq. (13) are 24-73 and 11-259 MPa, respectively. The calculated UCS or JCS was used as input in simulating jointed rock mass by Barton-Bandis(BB) failure criterion.

UCS(MPa) = 2.208e0.067Hr(13)

Fig. 5. GSI of investigated slopes along NH-58 from Rishikesh to Devprayag.

Table 3 Geomechanical inputs used for numerical modeling.

Furthermore,modulus of elasticity(Table 4)was determined by following empirical relationship proposed by Yagiz (2009):

Fig. 6. Measurement of geotechnical parameters during field survey: (a) Roughness profile of joint wall by using Barton comb; and (b) Schmidt hammer hardness of the jointed surface.

Table 4 Schmidt hammer hardness,Young's modulus and joint wall compressive strength of different joint sets at each investigated slope.

Shear strength of discontinuities is largely influenced by the roughness or asperities of the joint wall. Planar and smooth joints offer least shear strength and are susceptible to sliding.The roughness of rock joints plays a more important role particularly at low-stress levels, while at higher stresses, asperities are sheared posing smooth surface. The measurement of surface roughness in the laboratory is amenable, and it can be measured in field using different sized plates (Hencher and Richards, 2015). Barton comb was used during field survey to record the roughness profile of prevailing discontinuities at each slope in the study area (Fig. 6a). Later, to determine JRC values(Table 5), the recorded profiles have been compared with standard profiles proposed by Barton and Choubey (1977).Similarly, roughness friction also plays a significant role in guiding shear strength of joints in rocks. Tilt test has been conducted to measure the residual friction of rock joints(Table 5) for numerical modeling by BB failure criterion. The degree of fracturing in the rock mass was determined by the spacing of joints. Joint spacing is the perpendicular distance between two joints in a particular set. For each joint set, joint spacing (Table 6) has been measured carefully during the field survey by meter scale. Further, modulus of rock mass (Table 6)has been determined by an empirical relationship (Eq. (6))proposed by Hoek and Diederichs (2006). Joint stiffness is the one of the fundamental properties that is widely used in numerical modeling of the jointed rock mass. It was measured by direct shear field testing in field (Barton and Choubey, 1977;Bandis et al., 1983). The measurement of in situ joint stiffness is expensive and time-consuming. Barton (1972) and Rocscience(2001) recommended a method for determination of normal and shear stiffnesses of joints based upon the deformation properties of the rock mass and intact rock (Eqs. (9) and (10)).The normal stiffness for each joint set was calculated by Eq. (9)(Table 6). According to Singh and Goel (2002), the normal stiffness may vary from 10 to 30 times its shear stiffness. In the Himalayan region, the shear stiffness has been estimated as onetenth of the normal stiffness (Pain et al., 2014). As the shear modulus is not available, a similar approach has been followed here to estimate the shear stiffness of joints.

In rock slope engineering projects, FEM has been widely applied under diverse slope conditions. SSR analysis by FEM has been significantly employed to determine the FoS. As discussed earlier, the shear strength failure envelope is reduced systematically until the deformations are unacceptably large or solution does not converge. Due to its popularity and applicability under a variety of diverse conditions, SSR method has been adopted in the present study. A two-dimensional (2D) planestrain simulator ‘Phase2D' has been used for modeling. The models have been generated for all 20 slopes. A variety of information related to slope stability has been collected and synthesized from extensive field investigations, rigorous calculations and laboratory experiments. From the investigated cut slopes, the input parameters obtained from rock mass and intact rocks are showing widespread data across a range. Thegeometry of the model of each slope was generated by coupling the data obtained from laser inclinometer and Brunton compass.The plane strain analysis is being performed with metric units and Gaussian elimination solver. Stress analysis is conducted by 500 maximum iterations with a tolerance of 0.001. The modeling is being performed under gravity loading by discretization with a 6-noded graded triangle. For each slope, the bottom of the model has been fixed and the cut slope has been made free in both x- and y-directions, while for y-axis of the model, displacements and stresses have been restricted in the xdirection. An example of the generated model has been shown in Fig. 7.

Table 5 Residual friction angle and joint roughness coefficient at different joint sets in investigated road cut slopes along NH-58.

Fig. 7. The input model of Slope S17 along NH-58.

Table 6 Intact rock and rock mass moduli along with spacing and normal and shear stiffnesses at different joint sets.

Table 7 MC shear strength parameters calculated by GHB parameters.

After successful generation of the model, selection of an appropriate failure criterion is the utmost to have reasonably fair outcomes. The nonlinear GHB criterion has been employed for rock mass, while BB criterion has been followed for joints.Furthermore, for comparative analysis, MC criterion was also applied to the same slopes by using equivalent MC parameters.Few decades ago, certain numerical techniques do involve input parameters of nonlinear GHB criterion. Later, Hoek et al. (2002)derived empirical relationship from hundreds of case studies of slopes and underground projects and also introduced a window based program called ‘RocLab' to calculate equivalent shear strength parameters (cohesion and friction) from GHB parameters (Hoek and Marinos, 2007). The equivalent shear strength parameters were determined by using RocLab program(Table 7). Shear strain contours by GHB and MC criteria of each slope were represented in Figs. 8 and 9, respectively. These strain contours are significant in identification of the most prominent mode of failure.

Fig. 8. Shear strain contours and SRFs by GHB criterion of investigated slopes along NH-58.

The critical SRFs have been determined by using GHB and MC criteria,respectively(Table 8).The outcomes obtained by GHB and MC criteria have been compared and represented in Fig.10.It can be evidenced from the comparison that for slopes having SRF close to 1,the results are nearly identical;while SRF determined by GHB is sufficiently greater than 1,and SRF determined by MC increases by approximately 1.5-2 times. Thus it can be inferred that for fair heterogeneous rocks,MC criterion gives higher stability.Therefore,linear MC criterion can only be applied for extensively fractured or jointed rock mass, which can be assumed homogeneous. As nonlinear GHB criterion is available,MC criterion should be avoided in the jointed rock mass. Furthermore, a linear relationship between SRFs determined by GHB and MC criteria has been suggested by this study (Fig.11).

The zig-zag and planar distribution of shear strain contours indicates that structurally controlled failures are prominent in most of the slopes. It is also quite notable that the locations of shear strain by GHB and MC criteria are nearly matching with each other,but the thickness of slip surfaces by MC criterion is slightly less and thinner in contrast to the slip surfaces by GHB criterion. There are certain slopes in the investigated area,having critical SRF less than 1. Such slopes require immediate treatment by coupling rock bolting, grouting, and shotcreting. Retaining walls along with proper drainage system are also required for improving the geotechnical grounds of such slopes.The slopes having SRF of 1-1.3 are marginally stable and need an implementation of proper remedial measures.In most of the investigated slopes,structurally controlled failures have been evidenced from prevailing field conditions and shear strain pattern observed from simulation work.The sealing of discontinuities by employing grouting will increase shear strength along the discontinuities.For further safety,nets can be installed at particular slope facets to reduce the threat such as any occasional block failures or rock fall.The overall stability of the slopes having SRF greater than 1.3 is acceptable in terms of mass failure.But,in the case of jointed rock mass,such slopes are under continual threat of generating occasional wedges. Such adverse conditions were evidenced from Slope S15, where free-falling wedges destroyed the roadside garders and walls on either side of the road. To overcome such rock fall, nets of the desired mesh should be installed.The slope stabilization process is being done by considering several factors such as stability grade, capital to be invested, availability of raw material, and availability of dumping sites.

Fig. 9. Shear strain contours and SRFs by MC criterion of investigated slopes along NH-58.

Table 8 Stability grade and critical SRF determined by GHB and MC criteria.

Fig.10. Graphical representation of outcomes determined by GHB and MC criteria.

Fig.11. Linear relationship between SRF determined by GHB and MC criteria.

5. Conclusions

Slope stability assessment by numerical modeling technique is a significant part in formulating or achieving the safe and sound design in various rock slope engineering practices. The vulnerable slopes in the study area have been evaluated and stability grade has been quantified by means of critical SRF. Due to the adverse orientation of discontinuities, structurally controlled failures are prominent in the region.Moreover,meteorological,geological and geotechnical factors are the major factors accounting for frequent failures along the highway. Extensive rainfall during monsoon season also reduces the stability of cut slopes to great extent. The investigated patch has been dissected by major thrust faults. Most of the investigated slopes lie proximal to such sections.Due to such reasons, slopes near Kaudiyala and Shivpuri are very prone to failure. One of the slopes near Kaudiyala (Slope S13) had experienced mass failure during the period of the study. The slope was evaluated as unstable (critical SRF = 1.22). The outcomes, viz.critical SRF and shear strain contour distribution within the slopes,are being supported by existing field conditions.On the basis of SRF,road cut rock slopes in the investigated patch is categorized and demarcated as unstable, marginally stable and stable in terms of mass failure.Furthermore, some general guidelines for prevention and stabilization of critical slopes are also suggested to reduce the continual threat of landslides along the highway. As some occasional block failures have been reported even from stable slopes,mitigation for such rock fall should be adopted. MC criterion is more suitable for homogenous rock or debris/soil slopes.However,extensively jointed rock mass can be treated as homogeneous and the criterion can be adopted for stability study. In case of moderately heterogeneous jointed rock mass, nonlinear GHB criterion is well applicable.

Declaration of Competing Interest

The authors confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.

Acknowledgments

The authors are grateful to NRDMS Division, Department of Science and Technology, Government of India for providing financial assistance for field investigations.The authors are also thankful to Rock Sciences and Engineering Laboratory, Indian Institute of Technology Bombay for simulation work.