Rockfall hazard in the Imja Glacial Lake, eastern Nepal

In Nepal, rockfall related studies are rarely conducted and are limited to the landslide study along with few case studies on rockfall events. Rockfall problems in Nepal are more frequent in the Higher Himalayan region than Midland and Lesser Himalayan regions. In the glaciated valley and glacial lakes, rockfall and associated tsunami like huge wave are a recently initiated research. In this context, a glacial lake of the Himalaya, named as Imja Glacial Lake situated in eastern Nepal has been selected to understand the rockfall problems and their possible consequences. The lake was formed at the end part of the Imja glacier. The northern valley slope of Imja Glacier Lake, i.e. southern slope of the Island Peak has serious problem of rockfall into the lake. Rockfall simulation was performed during this research with field data. Three different Plots were defined for simulation of rockfall. Among them, Plot III seems to be the most hazardous since the detached boulders on the higher elevation can enter into the lake with the maximum velocity of 40 m/s with pre-impact energy more than 3500 kJ. This can develop a surge in the lake water that can break moraine dam. This leads to a serious threat to the downstream communities. Moreover, this research confirms that the rockfall hazard in the higher mountain region of Nepal is critical to creating flash floods due to tsunamis like surge in the lake water and breaching of the dam.


Introduction
Rockfall is a fragment of rock (a block) detached by sliding, toppling, or falling, that falls along a vertical or subvertical cliff, proceeds down slope by bouncing and flying along ballistic trajectories or by rolling on talus or debris slopes (Varnes 1978). These rock blocks can be detached by different mechanisms such as natural mechanisms like freeze -thaw cycles (McCarroll and Pawellek 1998;Matsuoka and Sakai 1999) seismic activities (Abebe et al. 2010) and by anthropogenic activities along the slope cuts or due to earth moving mechanism (Dorren 2003;Vijayakumar et al. 2011) which are slow time-dependent deformation that occur in the slope materials. The size of rock involved in a rockfall may range from pebbles to huge blocks with masses of hundreds of tones (Spang 1987) but there is no any generally accepted rock size and rock characteristic for rockfall. Therefore, rockfall can occur in any rock types ranging from granite to mudstone without any correlation between its strength and characteristics. Rockfall hazards are particularly severe in areas with heavy precipitation, frequent freeze-thaw cycles, and seismic events (Cruden and Varnes 1996).
Rockfall movement and trajectory is mostly dependent on different types of parameters such as slope geometry, slope height, slope inclination, surface irregularities, and blocks shape. It is known that, the rounded shaped massive blocks travel the maximum distance down the slope, followed by the square shaped blocks (Basson 2012).
Generally, four types of rockfall motion are known which includes free fall, bouncing, rolling, and sliding. When dealing with boulders on the top of slopes, the rockfall hazards are obvious. However, the most uncertain types of rock failure occur when a block is suddenly released from an apparently sound face by relatively small deformations in the surrounding rock mass. This can occur when the forces acting across discontinuity planes, which isolate a block from its neighbours, change as a result of water pressures in the discontinuities or a reduction of the shear strength of these planes because of long term deterioration mainly by weathering.
To understand the rockfall and its phenomenon and to give the suitable protection design, the study on the following factors are necessary; i. Topography and geology-The location of potential rockfall requires mapping to identify source areas, and the gullies in which the falls may concentrate. So, the geological studies provide the information on the likely size and shape of falls based on rock strength and on discontinuity persistence and spacing. ii. Calibration of rockfall models-Because of the complexity of rockfall behaviour, it is useful to have data on actual rockfall events to calibrate the mathematical models. Usually, falls on slopes comprising rock, talus, colluvium, asphalt, and concrete are documented to provide the calibration data. iii. Trajectory analysis-The trajectory that the rockfall follows between impacts is a parabolic path defined by gravitational acceleration, resulting in translational energy gain during the trajectory phase of the fall since trajectories define the required height of protection structures. iv. Impact mechanics-The impact process between a rockfall and the slope can be defined by the theory of impact mechanics. Application of this theory to rockfall enables calculation of changes in tangential, normal, and rotational velocities that occur during impact, and the corresponding changes in kinetic and rotational energy. v. Coefficients of restitution-The basic parameters defining the changes in tangential and normal velocities during impact are the corresponding values of coefficients of restitution. These parameters are related respectively with the friction on the contact surface and the angle at which the body impacts the slope. vi. Losses during impact-The result of velocity changes during impact are corresponding changes to the translational and rotational energies of the body. The energy changes are the result of the frictional resistance to slipping/rolling in the tangential direction, and plastic deformation of the body and slope in the normal direction. vii. Rock fall modelling-Computer programs have been developed to model rockfall behaviour and provide range values of impact-energies and trajectories for using in the design.
viii.Selection of protection structures-Methods of protection against rockfall includes construction of ditches, barriers, fences, nets, and rock sheds. Each structure has a specific range of impact energy capacity and suitability to the topography of the site. Thus, selection of the appropriate protection structure for a site involves detail study of the field data on rock fall dimensions and the analysis on frequency of occurrence.
Rockfall problems are unpredictable and incapable of being analysed with any precision that have been seen as a worldwide hazard. In Nepal, rockfall related studies are not much and are limited to the landslide study only along with some case study on rockfall events and its mitigation measures (Dahal 2016;Dahal 2018a;Dahal 2018b). Generally, in the highly elevated region i.e. in the Higher Himalayan region, maximum rockfall issues have been observed but their effects were not count since no human live there and there are no any economic activities. But on certain glacier lake area, climbers are placing base camp and rockfall problem can be a critical issue for such camping sites. The study of rockfall evaluation and consequences have been done on different approaches in different hillslopes of the world. The research on threat of rockfall issues and the occurrence of GLOFs are independently done in different aspects. Richard (2001) studied rockfalls using a mathematical model, codified for computer use, entitled GeoFall which was based on rigid body mechanics, and analyses a fall in 2D space using a new theory of impact dynamics developed by Brach (1991). These not only allowed the relative risk of rockfall between slopes to be assessed but also categorized the rock fall risk and the potential number of future rockfalls which was tested on 18 slopes at ten locations in Durham, U.K. Further, Crosta and Agliardi (2004) described the influence of different controlling factors on the dispersion of rockfall trajectories in the research which are systematically evaluated by performing 3D parametric modelling in which parametric modelling was performed at different spatial resolutions using sets of bi-planar simplified slopes characterized by different mean inclination and roughness. Also, Petje (2005) discussed the different ways of how rock masses (stones, rocks, and blocks) move along slopes and for each different way of motion (free fall, bouncing, rolling, sliding, slowing down, lubrication, fluidization) with adequate dynamic equations provided. Similarly, Leine et al. (2014) studied the dynamics of rockfall taking into account the shape of rock using the state-of-the-art methods of multibody dynamics and non-smooth contact dynamics, both of which were a fully 3D simulation technique in which the influence of rock geometry on rockfall dynamics was studied through two well-chosen numerical simulations in which the rock was modeled as an arbitrary convex polyhedron and the terrain model was based on a high resolution digital elevation model and Mikoš et al. (2006) studied the two rockfalls in the Trenta valley in NW Slovenia in which they analyzed rockfall consequences using a commercially available computer program for rockfall simulation which confirmed that the calibrated 2-D simulation models of rockfalls might be useful engineering tool when used carefully. Also, Fischer et al. (2012) discussed about the effects of changes in the highmountain cryosphere on rockfalls and rock avalanches about characteristics and distribution of recent rockslope instabilities. The lithological setting appears to influence volume rather than frequency of a slope failure. The analyses showed that not only the changes in cryosphere, but also other factors which remain constant over long periods play an important role in slope failures. Similarly, Lind (2016) performed research in a fullscale field test, to investigate how the blocks bounce where the full-scale field test was examined by filming them and then evaluated the bounce coefficient, coefficient of restitution from the individual rockfalls using photogrammetric methods. Besides, DCP test rig was used to estimate value on the coefficient of restitution and finally, during the evaluation, two rockfall simulating software (GeoRock2D and Rockfall) was used to investigate how well the results from them match the true blocks movements that was seen in the field (Lind 2016).
In Nepal Himalaya, Dahal (2012) stated that The 1255 flash flood in the Seti River of Nepal was one of the relevant example which showed that a rockfall from the Annapurna IV Range blocked the river for a while and store a huge amount of water that later converted into a massive debris flow. Meanwhile, this was occurred as a flash flood causing destructive damages on the Seti River Valley. Hence, the similar type of rockfall on the Imja valley may exist in the near future even due to small rockfall incidence which may generate waves in lake, causing collapse of the embankment, breaching of the lake etc. So, it seems an urgent for the study on the rockfall hazard assessment of the area. Nepal is facing huge loss of life and property yearly due to the frequent rockfall on the trekking routes, highways, base camps and valleys. Hence, in order to fulfill the lack of study on rockfall issues in the Glacier Lake, this research has been done.
With all these considerations, in this research a study has been performed in the Higher Himalayan region for rockfall problems and probable GLOFs hazard in the downstream. The main objective of the research is to evaluate rockfall events in the Imja Glacier Lake from surrounding peaks and their consequences.

Study area
The Imja Glacial Lake is located at 27°53′ 55" N latitude, 86°55' 20" E longitude and at an altitude of 5010 m in Everest Region of Nepal. The catchment of Imja Glacial Lake occupies the northeastern part of the Dudh Koshi sub-basin. The lake itself is located on the lower part of the glacier at the toe of its mother glaciers (snout of Imja and Lhotse Shar glaciers). The study area lies in the Sagarmatha National Park of Province No. 1 of Nepal on the south-east of Mount Everest where the Imja Valley runs east-west from the peak of Buruntse (7220 m), with Lhotse (8501 m) on the north flank and Ombigaichen (6243 m) and Ama Dablam (6856 m) on the southern part. Imja Glacier is located in the Imja River watershed in the Khumbu Region of the Nepal Himalaya ( Fig. 1).

General geology
Geologically, the Imja Glacial Lake lies in the Higher Himalayan Zone of Nepal with metamorphic crystalline basement rocks such as gneiss, schist and granite. An E-W extending regional normal fault system (South Tibetan Detachment System, STDS) passes from just north of the Imja Glacial Lake. It is also known as Qomolangma Detachment (Searle et al. 2003). Regional seismic data shows that this area lies in one of the seismically high hazard zones of the Nepal Himalaya. The Island Peak Complex is an association of mica schists, biotite gneisses, augen gneisses and quartzites, while it outcrops under Black Gneiss and associated granites on the eastern side of Island Peak. Rocks are gently dipping to the NW. The general geological map of the study area and its surroundings is given in Fig. 2.

Geohazards in Imja Glacial Lake
In the Imja Valley, different researches have been done frequently. Yamada and Sharma (1993), Sakai et al. (2003) and Somos-Valenzuela et al. (2014) conducted Bathymetric surveys of Imja Glacial Lake to examine the changes in the depth of the lake and noted increasing depth with time from past several decades. Also, Mool et al. (2001) had given notes of Imja Glacier on inventory of glaciers, glacial lakes and past Glacial Lake Outburst Flood (GLOFs) with the collaboration of International Centre for Integrated Mountain Development (ICIMOD). Similarly, Hambrey et al. (2008) carried out sedimentological and geomorphological study of Mount Everest (Sagarmatha) region, Nepal, whose main research was on the sediment, landform and dynamic context of four avalanche-fed valley Glaciers (Khumbu, Imja, Lhotse and Chukhung) in the Mount Everest (Sagarmatha) region of Nepal and Byers et al. (2015) gave notes on the geomorphological and glacial changes on Imja Glacier by using repeated photography technique and post hazard assessment of Imja, Thulagi and Tsho Rolpa glacial lakes. The research finally described about geomorphological changes of surrounding moraines with increased number of thermokrast, ice calving, piping, subsidence of water surrounding lateral moraines. As a result, Imja, Thulagi and Tsho Rolpa possess imminent threat to downstream population and mitigation steps that should be started immediately. So that, in 2016, the lake lowering project by Nepalese Army was performed in the Imja Glacial Lake by constructing an outlet and was succeed to drain over 4 million cubic meters of water from the lake.
The general engineering geological map of the Imja Glacier Lake is shown in Fig. 3 and it is clear that the source of rockfall in the lake is mainly from the southern face of Island Peak. There are many scree deposits on the southern slope of Island Peak. The base of the Island Peak has base camp of climbers to Island Peak and it lies on the bank of the Imja Glacial Lake.

Desk study
During desk study, initially different possible areas of rockfall in the southern face of Island Peak area were identified through the visual analysis of Imja Glacial Lake in Google Earth images and topographical map of the area. The possible location of the rockfall areas then named as "Plots". The expansion of lake was studied with the help of different articles and journals published. The geomorphology of the area was evaluated by the help of topographical map.
Total three Plots, namely, Plot I, Plot II and Plot III (Fig. 4) were extracted as main rockfall area based on the topography of the Island peak located on the northern side of Imja Glacial Lake.

Field study
Field visits were done to understand the geology, geomorphology and sources of rockfall in southern slope of Island Peak. The Plots identified during desk study were verified in the field. Within the Plots, the possible hanging blocks were identified the field (Fig. 5). In the analysis phase, individual slope profile on each Plot was made and trajectory paths of each hanging blocks on each Plot were simulated through GeoRock 2D software.
During the field visit, following works have been performed to gather the data.
The dimensions of each boulders deposited at each Plot were measured. The joints data was measured in each Plot area. The composition of each rock type was checked in the field. Overall walk through survey was carried. The boulders type on both the northern and southern side of lake was observed. The geomorphological estimation of the study area was done. Based on the field visit, the rock types identified and detail walk through survey, in the possible trajectory profiles the photographs were taken.

Estimation of parametric variables for modeling
The physically observed rock materials in the slope profile were defined in the certain ranges based on the control variables as frequency and roughness. The roughness and frequency values were considered as a parametric-  variables. The different parameters are defined based on the slope character and the boulders character observed in the study area. Hence, during the simulation, based on field observation and the photographs captured, four types of materials were characterized as the solid rock, rock detritus, fine debris and sand. From the photographs captured and the field verification, the roughness value and the frequency of the rock materials present in the slope profile were defined and determined for the simulation.
For the different trajectory motion on the same profile, seven rockfall simulations for each profile were estimated with a difference in the value of roughness and frequency, but the constant value of Rt and Rn. In which, R0 is the roughness value indicating the smooth surface whereas the R1, R2, R3, R4, R5, R6 are defined based on the variation on its roughness and frequency value (Table 1) determined after the rock materials observed in the slope profile in the field and photographs taken during the field.
Boulders in the depositional area were sampled in the field in four separate Plots based on the geographic separation of talus area. Statistical analysis of collected boulder dimensions was done in calculating the sample size of falling block: assuming the integrity of the falling boulder until it reached the talus. During the rockfall events, the falling blocks were assumed to disintegrate along the trajectory before coming to rest at the depositional area (talus). 2D Rockfall simulation was conducted in each assumed trajectory with parametric variation of surface roughness in the profile. The profile consisted on four lithological groups: exposed intact bedrock, rock detritus, fine debris and sand in decreasing order of altitude  respectively. A total of seven numerical simulations were done for each profile while analyzing the change in the maximum velocity at which the block falls, the maximum run-out distance it covers measured along the horizontal direction-the maximum height it gains, the maximum time it takes during each fall and the maximum preimpact energy it gains. Interpretation of collected results from numerical simulation was used in categorizing the most problematic site of rockfall, especially with consideration of the falling block reaching the Imja Glacial Lake.

Numerical modelling
For the rockfall analysis, there are different software available. Among which, during this research, GeoRock 2D, a free software was used. GeoRock 2D is the software developed by GeoStru for determining the trajectories of a boulder in free fall from a rocky ridge. The falling motion of a boulder along a rocky slope depends on many factors that are not easy to express numerically since the trajectories of the boulders depend on the geometry of the slope, on the shape of the falling boulder and on its initial velocity at the moment of detachment from the slope, and on the amount of the energy dissipated due to the impacts during the fall. The falling boulders can-slide, roll or bounce downstream depending on their shape, flattened or rounded, and on the gradient of the slope. Similarly, the energy dissipated due to impacts is generally different and varies with the characteristics of the motion and depends on the mechanical characteristics of the boulder and on the materials present along the slope (rock, soil, vegetation) that oppose in a different manner to the motion of the boulder. The simulation of falling rock is extremely challenging work because when a rock moves downslope, various modes of motion are possible: free fall, sliding, rolling, and bouncing.
Generally, there are following two analytical models in the GeoRock2D software, that are used mainly to study the phenomenon of rockfall.

i. Lumped-Mass Method and ii. Colorado Rockfall Simulation Program (CRSP)
The CRSP method is used in this research (Fig. 6,  Fig. 7, Fig. 8) in which the following parameters were initially defined for the simulation based on the field observation. The field data and the data extracted from Google Earth and Geocontext Profiler were used during the numerical modeling. The resolution of a profile is 15 m (GoogleEarth 2020) and it is enough resolution for such a long slope for rockfall simulation as in the study area (Dorren et al. 2004;Robiati et al. 2019).

Preparation of the topographic profile
For the analysis of the boulder's trajectories, first of all, it was necessary to define the topographic profile of the slope and the characteristics of the boulder. For which, the program analyzed the problem in the X, Y plane and performed any calculation which was reported to a 2D section defined in the Cartesian system (X, Y). During analysis, the profile was defined as x-axis: Distance in meters and y-axis: Elevation in meters. The individual source point of the hanging blocks and their formed topographic profile was geographically distributed so the run-out distances that come on each trajectory during simulation were also geographically distributed.

Materials characterization
The program has a database of materials that can be used to characterize each section of the slope. For each material, two restitution coefficients (normal and tangential), roughness, frequency and color were included within the simulation program. The tangential coefficient determines the reduction of the velocity component parallel to the profile section during the impact. The normal coefficient is a measure of the velocity change, normal to the profile before and after the impact. The basic material types defined during the simulation are: solid rock, rock detritus, fine debris and sand. Their respective restitution coefficient values were used as per the standard values assigned by the GeoRock2D software and the character of the materials observed from the field.

Preparation of boulder data
The falling boulder can be spherical, cylindrical, or disc, the choice of one of these options enable the corresponding geometric data to be entered during the simulation. For initial velocity, for the two components along X and Y, it can be started from equilibrium consideration of a spherical boulder of radius r placed on a plane that defines the initial velocity of the boulder (Lombardi et al. 2006) and it can be estimated by the following equation.
where, v = final velocity V 0 = initial velocity g = acceleration due to gravity α = Slope angle of detachment zone l = length of detachment zone

Computation method
During the research for the rockfall simulation, CRSP method was used and the following information were collected for each trajectory.
1. Impact points 2. Flight times between two subsequent impacts 3. Maximum height of the trajectory between two impacts 4. Flight velocity on the segment 5. Impact energy Also, in the rigid body approach, the shape of the rocks was selected from a wide variety of pre-defined smooth or polygonal shapes of the rock. The shapes were 2-dimensional and extruded for 3-dimensional size and mass. The mass of the rock determined the actual size of the rock shape, assuming that the outof-plane dimension of each rock was approximately equal to the average in-plane dimensions of the rock shape.
Since, the rigid body approach of CRSP considers the impulse reaction of the rock during the instantaneous contact period with the slope to determine the critical events of the rock (slip, sticking, reversal behavior during compression and restitution phase). The normal coefficient of restitution was used during the two phases of the contact period to calculate the terminal impulse. With the terminal impulse, the outgoing velocities at the contact point can also be calculated.
Hence, based on-site investigation for rockfall hazard mapping, altogether 35 risky hanging blocks were identified from Plot I, Plot II and Plot III which were geotagged in Google Earth in the Island Peak which had the potential to fall around and into the Imja Glacial Lake.

Analysis and results
During the rockfall simulation on different Plots, different results occurred on each Plot. Total of seven simulations (R0, R1, R2, R3, R4, R5 and R6) were performed (See Table 1). Understanding from each simulation in all plots is described below.

Plot I
According to the probabilistic distribution chart (Fig. 9) of Plot I, the average distance travel path as per its defined variable roughness and frequency is 84 m to 878.5 m. In case of smooth surface, Simulation R0, the runout distance is an average of 300 m from the source point but, only when the different roughness values were defined based on the data from the field. During simulations, three different pairs were obtained. The pair of R1 and R4, R2 and R5 and R3 and R6 were well noticed. The average distances were separated differently on each pair. These paired separated data indicated that the runout distance was more dependent on the roughness rather than on frequency. In R3 and R6, the run-out distance was less distributed from its mean value compared to R2 and R5 that was much more distributed than R3 and R6 but still, it's much more compact than R1 and R4 pair. As in R1 and R4, the run-out distance is widely distributed.
The graph R0 shows a very sharp curve having a Probability Mass Function (PMF) of 0.003 with an average travel distance of 250 m from the source. In case of surface roughness defined as R3 and R6, the average travel distance of falling rock is 500 m in the range of 84 m to 878.5 m with a standard deviation of 178.2 to 191 and a resulting Coefficient of Variance (CoV) of 0.34 to 0.36. Similarly, in R2 and R5, the average travel distance of falling boulders is 600 m in the range of 84 m to 878.5 m with a standard deviation of 207.6 to 215.8 and a resulting CoV of 0.34 to 0.35.
Likewise, in R4 and R1, the average travel distance of falling rock is 900 m from the runoff point, where the range is 84 to 878.5 m with a standard deviation of 317.8 to 297.1 and a resulting CoV of 0.36 to 0.32. Hence, it was understood that in the Plot I, the maximum rock blocks fall at an average distance of 500 m to 900 m of their total profile length.  Table 1 and it was decided as per different roughness values, i.e. 0 to 1) versus run-out distance in Plot I Simulation number versus distance explains how the run-out distance can be changed on a variation of the roughness. Also, in this work, the effect of roughness index for run-out distance was evaluated and shown in Fig. 10. Since they were dependent on each other, the run-out distance was categorized into three major categories, namely short run-out distance, medium run-out distance and long run-out distance. Based on the nature of a graph, irrespective of the source zone; the last four trajectories, 2.3, 2.6, 2.16 and 2.17 (Fig. 10) have traveled the long run-out distance that refers maximum changes have occurred in these trajectories due to change in the roughness values. However, 2.12, 2.10 and 2.8 have traveled the short run-out distance having fewer changes whereas, other trajectories have traveled the medium run-out distances. In R1 and R4, as compared to R0, the travel distance has increased at about 3 times. The graph shows three different pairs of data sets as R1 and R4, R2 and R5 and R3 and R6. In these pairs, the roughness value of R1 and R4; irrespective of frequency, have the maximum run-out distances compared to R2 and R5. Similarly, the run-out distances of R6 and R3 are comparatively less than that of R1, R2, R4, and R5.
This finally showed that in Plot I, 18 trajectories have its different source points, among which the four Plots are the most critical sites from where if the hanging blocks fall, they have the long travel distance that reach the embankment area of the Imja Glacial Lake. While others do not fall and get stuck at an elevation of 5281 m to 5176 m where the surface is covered with the rock detritus.

Plot II
Seven different trajectory data sets were prepared in Plot II. The data were used to obtain the probabilistic distribution chart and to calculate run out distances from the individual source point. Figure 11 illustrates the trajectory motion of all 7 different data sets.
According to the probabilistic distribution chart on Plot II, the average distance of travel path according to its defined variable roughness and frequency is 863 m to 1156 m. In case of smooth surface R0, the run-out distance is an average of 163 m from the source point but when we defined the different roughness values based on the data collection in the site, different three pairs were visible. Those pairs are R1 and R4, R2 and R5, R3 and R6. The average distances are separated differently on each pair. But they were not widely distributed as of Plot I. Unlike Plot I, in Plot II, the frequency plays a vital role for the run-out distances because the probability distribution chart shows that the average distances of all trajectories are same and a narrow distribution is well noticed. It indicates that the probability of distribution is wider and the probability of occurrence is less in R4 (0.3%) but in R1, the probability of occurrence is high (100%). In pair R2/R5 and R3/R6 there is no significant difference in run-out distance. Therefore, based on the simulation number (i.e. guided by roughness) versus distance Plot (Fig. 12), the maximum run-out distance was noticed in the R1 and R4. As a result, the graph with R0 shows a very wider curve having PMF of 0.0034 and an average travel distance of 164 m from its  Table 1 and it was decided as per different roughness values, i.e. 0 to 1) versus run-rut distance in Plot II Similarly, Simulation number (decided as per the roughness) versus run-out distance Fig. 12 explains how the run-out distance changed a variation of roughness values. Since they were dependent on each other, it has been divided into three major categories as short runout distance, medium run-out distance and long run-out distance. Based on the nature of a graph, irrespective of the source zone, one trajectory (Trajectory 3.5) has traveled the long run-out distance that means maximum changes have occurred in this trajectory due to change in the roughness values whereas, other trajectories have traveled the medium run-out distances. In R1 and R4, compared with R0, the travel distance has increased at about 10 times. This graph shows that trajectory 3.5 has the maximum run-out distances in R1 and R4 whereas lesser run-out distances in R3, R2, R6, and R5.
This finally showed that in Plot II, seven trajectories have their different source points, among which the only one is the most critical site from where if the hanging blocks fall, they have the long travel distance that reaches the embankment area of the Imja Glacial Lake crossing the Island Peak Base Camp area and get stuck at an elevation of 5186 m-5100 m where the surface is covered with the rock detritus.

Plot III
In Plot III, total of 10 different trajectory data sets were defined that were used to prepare a probabilistic distribution chart. Also, in this plot, each run-out distances were calculated from the individual source point. Figure 13 illustrates the trajectory motion of all 10 different data sets. According to the probabilistic distribution chart on Plot III, the average distance of travel path according to its defined variable roughness and frequency is 303 m to 625 m. In case of smooth surface R0, the run-out distance is an average of 156 m from the source point but when the different roughness values were used, in this plot also, different three pairs were obtained. The pair of R1 and R4, R2 and R5, and R3 and R6 were well noticed in Plot III (Fig. 13). The average distances were separated differently on each pair. But they were not widely distributed as of Plot I and II. Though they are in pairs but no significant variation was observed in the graph. This suggests that in Plot III, there is a little effect of roughness but there is no effect of frequency. Likewise, there is an effect of roughness but a change in roughness value is not effective in Plot III.
In Plot III, almost all the trajectories have a probability to reach into the Imja Glacial Lake (Fig. 14). The graph suggests that R0 shows a very wider curve having a PMF of 0.003 showing the average travel distance of 156 m from its source point. In the case of surface roughness defined as R3 and R6, the average travel distance of falling rock is 467 m in the range of 303 m to 625 m with a standard deviation of 101 and CoV of 0.21. In R2 and R5, the average travel distance of falling boulders is 491 m in the range of 303 m to 625 m with a standard deviation of 93 and CoV of 0.18. In R4 and R1, the average run-out distance of falling rock is 491 m from the  Table 1 and it was decided as per different roughness values, i.e. 0 to 1) versus run-rut distance in Plot III detached point. The range value of run-out distance is between 303 m and 625 m with a standard deviation of 95 and CoV of 0.20. In fact, in Plot III, almost all the hanging blocks, if fall, can enter into the Imja Glacial Lake. The evaluation of run-out distances for Plot III was avoided because simulations of all the selected hanging blocks suggest that every fallen block can enter into the lake and total run-out distance in the water was out-of-scope of this research.

Discussion
After the rockfall simulations on three different Plots of the southern slope of Island Peak, the maximum run-out distance was at an elevation of 5038 m of blocks 2.3, 2.6, 2.16, 2.17 in Plot I and block 3.5 in Plot II. In the field, all rock blocks were deposited at an elevation of 5281 m to 5176 m in Plot I and in an elevation of 5186 m to 5100 m in Plot II (Fig. 15). Other maximum blocks were accumulated around the Island Peak Base Camp area (elevation 5082 m). But in the case of Plot III, almost all blocks reached into the Imja Glacial Lake. Likewise, the simulations pointed out that the rockfalls on Plot I and Plot II may occur with the energy of about 6200 kJ -22, 500 kJ whereas on Plot III, the energy will be about 3500 kJ only with velocity of 40 m/s but rock boulder can enter into the lake.
In the present scenario, the water depth in the upstream side of glacial lake is not much high and huge ice crevasses are present at the nearby upstream area. The glacial melting process is to continue increasing the area. Therefore, in the future, it is certain that the depth of water in the lake will be higher than the present depth in Plot III area. This can cause a tsunami-like huge surge in the Imja Glacial Lake and can breach the moraine dam downstream of the lake that can create a GLOF in the downstream of Imja River.

Conclusion
This research confirmed that on northeast corner of Imja Glacier Lake, rockfall problems are prominent. Any rockfall in that area can generate a huge surge instantly except in the winter season when the lake is in a frozen state. This can bring out the critical tsunami-like surge in the lake that might lead to the breaching of the moraine dam. So, the issue of rockfall in the Imja Glacial Lake is very critical. This research also pointed out possibilities of the GLOFs in the Himalaya not only due to the natural breaching of the moraine dam but also due to possible rockfall into the lake. Many glacial lakes of Nepal are surrounded by steep rocky slopes and the consequences of rockfall in those glacial lakes also need to be evaluated.