Numerical analysis of rockfall and slope stability along the Karakorum Highway in Jijal-Pattan

, ,3,4, , , , ,

1a.Key Laboratory of Mountain Hazards and Earth Surface Process; 1b.Institute of Mountain Hazards and Environment, Chinese Academy of Sciences (CAS), Chengdu 610041, P.R. China;2. University of Chinese Academy of Sciences, Beijing 100049, P. R. China; 3. CAS Center for Excellence in Tibetan Plateau Earth Sciences, Beijing 100101, P. R. China; 4. China-Pakistan joint Research Center on Earth Sciences, Islamabad, Pakistan)

Abstract: Along the Karakorum Highway (KKH), the key route for the China-Pakistan Economic Corridor, there are many rockfalls and unstable slopes, usually caused by tectonic movement and rainfall on the fractured rocks and slopes. This paper presents a numerical investigation of the rockfall and slope stability along the Karakorum Highway in Jijal-Pattan, Northern Pakistan using DIPS, GeoRock 2D and SLIDE, focusing on rockfall and slope stability along the KKH to develop countermeasures. Along the KKH, two major sections susceptible to rockfalls were selected to investigate the mechanism of rockfall and slope instability. The stereographic projection analysis following four sets of joints indicates that both sections are prone to plane failure and wedge failure. Based on the limit equilibrium theory, under static loading, the slope for Section 1 showed a stability coefficient of 0.917, representing its instability, and the slope in Section 2 has a stability coefficient of 1.131 depicting its slight stability. However, under the seismic condition, the stability coefficients of the slopes were lower than 1 for both sections, which indicates their instability. The results by GeoRock 2D reveal that in Section 1 the fallen rock mass attained the bounce height of 33 m, and in Section 2 it attained a bounce height of 29 m. The fallen rocks in Section 1 have the total kinetic energy of 1 135.099 kJ with a velocity ranging from 0.5 m/s to 44 m/s, while in Section 2 the fallen rocks have a velocity ranging from 0.5 m/s to 40.901 m/s with a damage capacity of 973.012 kJ. This study showed the rockfalls and landslides along the KKH have great damage potential.

Keywords:Karakorum Highway; steep rock slope; stereographic projection; slope stability analysis; dynamic process

1 Introduction

Rockfalls and landslides are widely known hazards in mountainous areas. Rockfalls usually include the quick movement of rock boulders in the form of falling, bouncing, or rolling[1]and are a great threat to people, their livelihoods, environmental services and resources, infrastructure, and economic, social and cultural assets[2]. One of the major causes of rock slope failure is the construction of roads without proper geological and geotechnical engineering investigation of the natural rock slopes[3-5]. In addition, it is not possible to continuously monitor the rock slope, particularly in the rainfall season. The threat of rockfall exists whenever humans or nature disturb the natural balance of the rock slope[6-8].

The forces triggering rockfalls are usually earthquakes, temperature fluctuation, and neotectonic activity[9-12]. The study of rockfalls along highways is of interest to many researchers. Singh et al.[5]investigated rockfall activity along the Luhri hydro-electric project on the Sutlej River in Himachal Pradesh, India, and performed a kinematic analysis to assess the failure mode. RocFall v4.0 was used to study the trajectories and energy dissipation of the falling rock. It was observed that the main reasons for the rockfall were weak rock mass and rainfall. They described that the fall of rock blocks was a potential threat to human lives and infrastructure. Singh et al.[4]also analyzed the stability of the road cut cliff face along SH-121, in Maharashtra, India through rockfall analysis and finite element modelling. They reported that the study area was prone to rockfall hazards, particularly in rainfall events, due to the steep and highly jointed slopes along the roads.

Slope failure is another great threat along the KKH. Slope failure is the result of forces such as increased destabilization or seismic events, external load, undermining, increased water pressure in rock cracks, hairline cracks and frost wedging, mining and loss of capillary pressure[13]. Slope stability analysis can be carried out using the limit equilibrium method, numerical modelling techniques and kinematic analysis. Kinematic analysis is suitable for identifying slope failure types using discontinuities and joint orientations[14-16]. Slope stability issues can be minimized with in-depth monitoring and analysis[17]. Akram et al.[18]carried out the stability evaluation of a slope in Balakot, Pakistan, which is one of the seismically active regions in Pakistan. These researchers performed kinematic analysis using limit equilibrium methods to assess the failure modes of slopes and to evaluate the stability of slopes under different conditions. It was concluded that the slope failure modes were plane, wedge, and toppling, with less likely chances of circular failures. The above-mentioned studies were focused mainly on the assessment of rockfall due to slope orientation, rock joint condition, and dynamic stability. However, in the Jijal-Pattan area, the rockfall results from earthquakes along with weak rock conditions, steep slopes, and the lack of geotechnical engineering investigations. The rockfall assessment in such areas should be carried out by performing slope stability analysis in both static and dynamic conditions along with kinematic analysis to assess the failure modes of the rockfalls.

The northern part of Pakistan is comprised of high mountain ranges with a history of rockfalls due to seismic activity, particularly in the area between Jijal and Pattan[19]. The only mode of transportation in such a mountainous area is by road, but recurring rockfalls and landslides lead to damage to the infrastructure, residents, and travelers. The Jijal-Pattan road is an important part of the China-Pakistan Economic Corridor (CPEC). However, due to complex tectonic conditions and multiple seismic events, this section of the road is characterized by highly fractured and jointed rocks. Further, ill-considered rock cuts for infrastructure development in this area facilitate rockfalls and landslides. Nonetheless, there has been little study of the rockfalls and landslides in this area, and it is crucial to investigate the mechanism of the rockfalls and landslides here due to the threat to human life and infrastructure.

This research aims to reveal the stability of the slope and the extent of the threat from rockfalls and landslides along the road from Jijal to Pattan through field investigation and numerical studies. Along and across the slopes, DIPS was initially used to obtain the geological orientation and perform kinematic analysis of the major planes[20]. GeoRock 2D software was used to display the rockfall analysis based on the kinematic analysis[21]. The software SLIDE 6.0 was used for the numerical slope stability analysis[22]. The key contribution of this study is the usage of three different models to assess the slope stability and rockfall risk in the Jijal-Pattan area to fill the research gap.

2 Geo-location of the research area

Northern Pakistan is linked to Western China through the Karakorum Highway (KKH), which forms a part of the China-Pakistan Economic Corridor (CPEC). The rising of the Himalayan, Karakorum, and Hindu Kush Mountains represent the collision of the Indian and the Eurasian plates and the Kohistan Island Arc[23]. The study area is the Lower Kohistan District (Jijal-Pattan) along the Karakorum Highway in the Khyber Pakhtunkhwa Province, Pakistan. The Lower Kohistan District extends from latitude 34°54′ to latitude 35°52′ north and from longitude 72°43′ to longitude 73°57′ east. It borders the Ghizer and Diamer districts on the north and northeast, the Manshera District on the south-east, the Battragram District on the south, and the Shangla and Swat districts on the west. The geology of the study area mainly contains sedimentary rocks, igneous rocks, and metamorphic rocks. Along the KKH, highly active rockfall areas have been identified. The lithology of the study area consists of the Besham group, the Jijal Complex, the Kamila amphibolite, and the Chilas Complex.

From Pattan to Kamila and north along the Indus River, the Kamila amphibolites are well exposed. South and north of Kamila lie large gneissic and huge granitic bodies consisting of sheet-like intrusions. During the field visit it was observed that the section is dominated by amphibolies of the Besham group, which is of Cretaceous age intruded by younger granodioritic gneiss with little schists at the top of the slope.

The northern part of Jijal along the KKH lies in a highly vulnerable zone. It consists of the highly fragmented Jijal Complex ultra-mafic rocks. It is extremely jointed and locally sheared because the study area is located in the hanging wall of the Main Mantle Thrust (MMT)[19].

This seismic area is only 3 km away from the site of the earthquake in Pattan (Magnitude=6.2, Depth=22 km) on 28 December 1974[24]. In this part, the topography is steep, mostly with slope angles of more than 50°, even up to 90°. The area is located in the monsoon region, where the annual average precipitation is over 400 mm[25]. A large number of rockfalls have occurred due to heavy rainfall and biological weathering, blocking the roads for weeks. The surface of the slope is moderately weathered, which has produced clay with medium vegetation cover. Due to the dominant weathering, the slope surface is covered by rock fragments ranging in size from pebbles to boulders. Fault closeness, biological weathering, strong seismicity, fractured rock mass, heavy rainfall, and steep topography, are all responsible for the rockfalls in this region.

Geological Cross-Section of the slopes: Section 1 is located along the Karakorum highway between the Pattan Tehsil and the Mali Dhera Kohistan District. During the field survey, it was observed that the section is dominated by amphibolites of the Besham group. Cretaceous age intruded by younger granodioritic gneiss with little schists at the top of the slope were also observed. The section is moderately jointed with an almost 100 m-high slope facing N30E and a dip angle ranging from 67° to 80°. The Cretaceous amphibolite sheet intrusion sub-parallel to the fabric and banding is very common throughout the Besham group. Most parts of the slope were covered by fallen rocks, indicating the high risk of rockfall impacting the asphalt road. The geological cross-section of both sections is given in Fig.1. Section 2 is situated a few kilometers away from Dubair, Kohistan. Cretaceous amphibolite dominates in this highly fractured section and a well-defined joint system was observed in this section. The slope face dips in the N58E direction with an average dip angle of 70°. The surface of the slope is moderately weathered, producing clay that favors medium vegetation cover. Due to the dominant weathering, the slope surface is covered by rock fragments ranging in size from pebbles to boulders.

Fig.1 Geological cross section of (a) Section 1 and (b)Section 2

These amphibolites are coarse-grained with a fracture filling of quartz. The overall dip and strike of the granodioritic gneiss are 60° and N82E, respectively. A geological cross-section is given in Fig.1.

3 Material and methods

The design of the cliff is an iterative process, and no principles are defined throughout all areas[26], therefore, every survey is important. The rockfall prone areas were identified during a field survey. The slope height, slope angle, block size, block shape, joint spacing and biological weathering of each section were identified and measured during the field study. The open joints, blocks overhanging the KKH and biological weathering were found to be vulnerable to instability. However, only a few man-made cut slopes are located along the KKH. The Rocscience software SLIDE 6.0 and DIPS were used to analyze the profile of each section, GeoRock 2D software based on kinematic analysis was used to display the orientation of major planes along and across the slope to investigate the slope failure. In this study, two different vulnerable sites were chosen and analyzed kinematically for rockfall by numerical analysis, which are explained below.

The slope may be naturally formed or man-made.The man-made slope includes excavation/cut for construction, borders of embankments, dams, canals etc. Many factors contribute to slope failure, including 1) forces due to the seepage of water, 2) gravitational forces, 3) sudden lowering of the water table adjacent to the slope, 4) earthquake forces, 5) reduction in strength of the material and 6) a non-engineered cut. Slope failure occurs in several modes. Cruden and Varnes[27]classified slope failure into five major categories: fall, slide, topple, spread, and flow.

Different methods are available to compute the slope stability for rock and soil. Due to the advancement of computer technology, a number of slope stability tools exist for both rock and mixed rock-soil slopes[28].

Kinematic analysis shows the orientation of rock discontinuities (joints, fault, bedding, etc.) is the leading factor influencing the stability of rock slopes[26]. Different failure modes are associated with the orientation of discontinuities for plane failure, wedge failure, toppling failure and circular failure[5-6,29-30]. Kinematic analysis using stereographic projection gives the geometry of discontinuities and analyzes the result to predict the type of failure. Scanline survey was used in this study to find the parameters of the rock discontinuities for stereographic analysis. These parameters include the type of discontinuities, persistence, aperture, property of infilling, spacing, roughness, water condition and lithology.

A rockfall is the movement of a rock or boulder sliding, toppling or falling along a steep or sub-vertical slope, which proceeds down a steep hill both bouncing and flying or rolls downwards over debris slopes or talus[1]. Various geometrical, geological, geotechnical and climatic influences lead to the initiation of significant rockfall incidents in mountainous regions. In this study, the widely used GeoRock 2D software was used to display the orientation of major planes along and across the slope to analyze the slope failure.

The most commonly used factors of the environmental condition of the study area[21, 31-32]are shown in Table 1.

Table 1 Boulder properties of Section 1 and Section 2

The minimum and maximum mass of a rock boulder in Section 1 is 12 kg and 4 486 kg, respectively, while the minimum and maximum mass of a rock boulder in Section 2 is 9.8 kg and 5 972 kg, respectively, as shown in Table 1. Section 1 is a highly weathered rockmass jointed with some large spacing and size blocks in the face zone, while Section 2 is highly fractured with a thin cover of weathered overburden. The rockfall simulation technique calculated the trajectory, kinetic energy, velocity motion, bounce height and run-out distance of the falling blocks based on the theory of collision and laws of motion[33].

The numerical analysis of Section 1 and Section 2 was carried out using the commercial software Rocscience SLIDE 6.0. Since the strength of the rock mass is controlled by discontinuities, the Hoek-Brown failure criterion is used in the analysis. The safety factor was calculated based on the limit equilibrium method (LEM). The limit equilibrium method is commonly used in geotechnical engineering problems related to seepage and the stability of slopes. It uses the perfectly plastic Mohr-Coulomb criterion to model soil stress-strain behavior.

4 Results

4.1 Kinematic slope stability analysis

Based on stereographic analysis, four sets of joints were observed. The types of critical discontinuity planes that are prone to fail are listed in Table 2. The planar sliding analysis results for Section 1 show that only 3/7 poles are in the critical zone, having 42.86% probability of plane failure. When considering the individual sets, 2/3 poles are in the critical zone in set 1 with 66.67% probability of plane failure. For set 2, 1 pole has a 100% probability of failure (Fig.2(a)). For Section 2, out of 6 poles, only 2 poles are in the critical planar sliding zone, with 33.33% probability of plane failure. However, for set 2, all the poles are in the critical state having 100% probability of plane failure (Fig.2(d)). The results revealed that set 2 in both sections is more susceptible to planar failure (Table 2). Nagendran et al.[34]mentioned in their research, the overall plane failure was 8.66%, where the probability of failure was 15.40% (set 1), which is in line with our results. Section 1 wedge sliding analysis presents that 11/15 inter-sections are under the critical condition, with 73.33% probability of failure. For Section 2, the critical inter-sections are 10/15, with 66.67% probability of failure (Fig.2(b) and (d)), which is comparable to the results of Tiwari et al.[35]and Sazid et al.[36]. This indicates that the percentage of critical intersections in these analyzed sections is very high and more prone to wedge failure (Table 2).

Table 2 Dominant joint set data at two rockfall sites

4.2 Rockfall slope stability analysis

Ritchie[37]proposed that falling blocks achieve numerous types of motion, depending on the slope structure and the mechanical properties of the blocks. In free fall movement, there is hardly any interaction with the slope surface, while the falling mass continually interacts with the surface in other motions such as rolling, sliding and bouncing and each impact changes their impact and energy.

Fig. 2 Plane failure and wedge failure of

Impact positions of the falling rock bodies (X(m) andY(m)), the falling rock mass (kg), first strike point and run-out distance of falling rocks in both sections are described. The slope height is 103.87 m in Section 1 while in Section 2 it is 123 m as shown in Fig.3, which causes the blocks to achieve greater velocity and even extensive momentum. In all situations, rock movement starts with sliding. The falling rocks strike the rock slopes and bounce and roll several times before stopping or resting at the asphalt. The estimated maximum bounce heights for both sections plotted with the run out distance, are shown in Fig.4(c) and (d).

Fig.3 Slope face, trajectory, motion and run out the distance of the falling body at (a) Section 1 and (b) Section 2

The results showed that the rock mass in Section 1 attained a maximum bounce height of 33 m, while in Section 2, it attained a maximum bounce height of 23 m. The damage capacity of the rockfall based on the translational velocity and kinetic energy values obtained through rockfall analysis, shows that the damaging impact of the rockfall is as high as 1 135.099 kJ with a falling velocity ranging from a minimum of 0.5 m/s to a maximum of 44 m/s in Section 1, as given in Fig.4(a), (e) and Table 3.

The values of the damage capacity of the fallen rock mass in Section 2 are estimated to be 973.012 kJ, with a maximum velocity of 40.901 m/s and a minimum velocity of 0.5 m/s, as shown in Fig.4(b), (f) and Table 3. The parabola height, kinetic energy and velocity of falling rocks are greatly influenced by slope height, slope angle and the

Table 3 Statistic computations of Section 1 and 2

Fig.4 Energies at each strike of the fallen rock mass (a)(b),the trajectories and their parabolic heights attained by the fallen rock bodies(c)(d),and velocity pattern of the fallen rock bodies(e)(f) at Section 1 and Section 2.

weight of the falling boulders. Similar results were obtained by Choi et al.[31]. Perret et al.[38]divided kinetic energy into three intensity groups to determine the hazardous zones. The highest intensity zone had a kinetic energy of more than 300 kJ, which is achieved in this study in both sections, as shown in Table 3. The medium intensity zone has a kinetic energy of 30～300 kJ and the low- intensity zone has a kinetic energy of <30 kJ. Dorren[12], Perret et al.[38]and Basson[39]suggested that a descending block reaches a velocity of 5～30 m/s and eventually it stops underneath a slope of 30°.

These analyses showed that the majority of the fallen rocks stopped at the roadside after losing most of their energy, with few falling further down to the valley floor. So the chance of impacting the commuters on a mountain road is very high.

4.3 Numerical analysis

Limit equilibrium method (LEM) parameters including cohesion, angle of friction, unconfined compressive strength (UCS), unit weight, the Geological Strength Index (GSI) value and loading type were used for the purpose of calculation, as shown in Table 4. Considering the seismic loading of 0.24g, Seismic Zone 2B[40]parameters were followed. These analyses followed Janbu methodology, and 638 possible slope slice surfaces were considered. For Section 1, 25 sets of different critical parameter slices were calculated such as base cohesion, base friction angle, shear stress, and shear strength, with the numeric model shown in Fig.5 (a) and (c). For Section 2, 550 slices were considered for the calculation of these parameters, and 19 slice sets were defined. Graphical representation of Section 2 under static and dynamic conditions, is shown in Fig.5 (b) and (d). These analyses show that the stability of the slope is directly dependent on the safety factor, which is 0.917 for Section 1, showing that it is unstable in static conditions, while the safety factor for Section 2 is 1.131, showing that it is slightly stable under static conditions (Table 5, Fig.5 (a) and (b)), whereas the safety factor is defined from the numerical analysis under dynamic conditions for Section 1 and Section 2, which is 0.539 and 0.784, respectively (Table 5, Fig.5 (c) and (d)). Exposure to any seismic activity, whether by artificial blasting or natural, such as an earthquake, can induce movement down the slope under gravity and both the slopes can fail at any time producing socio-economic disaster in the area. A similar trend of static condition to the dynamic condition was acquired by Akram et al.[18]

Table 4 Average values of the parameters

Note:UCS is unconfined compressive strength; GSI is Geological Strength Index;sandaare constant values which depend upon the characteristics of the rock mass; mb is a reduced value (for the rock mass) of the material.

Table 5 Derived parameters obtained from the lab analysis

Fig.5 Numerical slope stability model without seismic loading of (a) Section 1 and (b) Section 2 and 13-5.tifwith seismic loading of (c) Section 1 and (d) Section 2

5 Discussion

At each site, hundreds of dip/dip directions were calculated by Brunton compass while dominant joint set data were measured by pole density. In a major part, four sets of slightly weathered joints were identified at each site at diverse persistence and frequency, as listed in Table 2. The results of the kinematic studies are used in this study to evaluate the failure mode in both sections. Rocscience Software DIPS was used to display the major planes along and across the slope and to analyze the data for failure types, as shown in Fig.2.

Rockfall slope stability analysis shows that the blocks were triggered to fall due to the steep rock slope face, reaching the road and causing undesirable consequences. The trajectories of their fall and their endpoints for Section 1 and 2 are shown in Fig.3. Most of the rock boulders reached the bottom of the slope due to the absence of benches. Their trajectories are decided by the collision of the rock boulder on the face of the slope.

Additionally, the orientation of the blocks is determined by the properties of the slope. Upon the first impact, a large amount of energy is lost and the blocks are separated into smaller sections. Most of them may stop at their first impact, and some may move hundreds of meters downhill to the valley floor.

Numerical slope stability analysis analyzes the equilibrium between driving forces and resisting forces and takes into account material parameters like cohesion, angle of friction, USC, unit weight and GSI value, as shown in Table 4. Limit equilibrium analysis is used for numerical slope stability analysis to define the critical surface on which the movement of rock or soil occurs or is expected to occur. The critical surface is based on the minimum factor of safety. The limit equilibrium method (LEM) is based on a common approach "resisting forces/driving forces"[28].

Based on the analysis of rockfall events,the failure characteristics of the slopes, and the energy of the falling blocks, appropriate structural countermeasures are proposed to improve the stability of the selected sections with the aim of avoiding geological hazard. The suggested structural countermeasure is an anchored rock mesh system with rock bolts as additional support. The anchored rock mesh system consists of a steel mesh anchored by bolts, which covers the rock surface and restrains the movement of small rock blocks on the slope. Additionally, the selected sections are also susceptible to slope failure; therefore, rock bolting is recommended to increase the safety factor to prevent the sliding of the slope.

6 Conclusion

The Karakorum Highway in Pakistan is not only a very important route for the business trade between China and Pakistan, but also important for domestic trading. The area along the Karakorum highway is very vulnerable to slope failure and rockfalls that could put people at risk and result in significant finacial losses. The kinematic analysis showed that the two study sites are highly jointed with a dip/dip direction of 67/32, indicating 100% susceptibility to sliding under gravity. Rockfall analysis showed that the fallen rock mass in section 1 has attained a maximum bounce height of 33 m whereas the fallen rock in in section 2 has attained a maximum bounce height of 23 m. The damage capacity of the fallen rock in section 1 was probably 1135.099 kJ, with the velocity varying from 0.5 m/s to 44 m/s, while in section 2, it was 973.012 kJ, with a minimum velocity of 0.5 m/s and a maximum velocity of 40.901 m/s. These analyses also showed that the majority of the fallen rocks stopped at the road, having lost most of their energy, with few falling further down to the valley floor. Therefore, the chance of falling rocks impacting commuters is high. Based on the numerical analysis, the stability of the slope directly depends on the safety factor, which was 0.917 for section 1, showing that it is unstable in the static condition, while the safety factor of section 2 was 1.131, showing that it is slightly stable under the static condition. However, the safety factors for both sections were less than 1 under the dynamic condition, which means that the slopes are unstable and can slide any time. To avoid the hazards of rockfall and landslides, engineering countermeasures are proposed.

Acknowledgements

The authors would like to acknowledge the financial support from The Strategic Priority Research Program of the Chinese Academy of Sciences (Grant No. XDA20030301) and "Belt & Road" international cooperation team for the "Light of West" program of CAS.