Landslide Hazard Assessment and Runout Simulation by Utilizing Remote Sensing, Geophysical, and Geological Techniques


 Landslide events in Karakorum ranges are frequent and have already damaged local infrastructures and roads. In the hilly regions, landslide characterization and predicting its deposition pattern are essential for accurate engineering hazard assessment. To this end, numerical simulation models are commonly used tools. However, appropriate model parameters are often not available to predict and generate real landslide scenarios. This work describes the use of multidisciplinary techniques to estimate the model parameters for a slope prone to landslide and simulate the hazard level. The first important parameter, landslide boundary, and dynamics were estimated from temporal satellite images by identifying the areas with prominent deformations using the Interferometric Synthetic Aperture Radar (InSAR) technique. The susceptible subsurface strata volume and the possible landslide initiation depth were determined with the electrical resistivity method. In addition, voxel 3D electrical resistivity models were created to present the depth of the existing rupture and the nature of subsurface strata. The soil mechanical parameters were calculated during field visits and laboratory tests. The parameters adopted from different techniques helped simulate the susceptible landslide volume and initiation depth. These parameters are a critical factor in developing an accurate high-speed landslide model through numerical simulation. The applied methodology is vital to understand the dynamics of a particular slope and perform accurate engineering hazard assessment with numerical simulation. The results are essential to predict the potential deposition areas of the landslide event accurately, minimize the risk level, and take proactive mitigation measures.


Introduction
The Karakoram Himalaya is part of the south-central Asian mountain systems. It covers their highest and most rugged terrain and, with some 16,000 km 2 of perennial snow and ice, is the most heavily glacierized.
Most, if not all, of the Karakoram, was glaciated several times during the Quaternary in some of the most signi cant glacier expansions of the subtropics. They are recorded in glacially moulded landforms, old lateral-margin deposits, and erratics that lie 1000 m or more above valley oors and hundreds of kilometres beyond modern glaciers (Hewitt 1999). A large number of catastrophic landslides have occurred in these lateral-margin deposits. Occasionally, they have damaged residential buildings, roads, and power grids. Hence, precise identi cation of future landslides and their extent, volume, and activity is vital for disaster management (Kamiński et al. 2021). Currently, rainfall-induced landslides are dominant in the Karakorum ranges. With the progressive development, landslides and ash oods have caused many casualties and blocked the Karakorum highway several times (Ali et al. 2019). In most events, the landslide movement appears to be linked primarily with the topography and the surge in the bedrock saturation either by groundwater or rainfall. The other important factor is Earthquakes to be considered.
As established by new earthquakes, the Karakoram Himalaya ranges are tectonically active (Yousuf et al. 2020). In Karakoram Mountains, there exits mainly deep structural landslides with huge volumes, controlled by the relief and geological strata. One of them is the "Maicher Hill" active landslide.
The rst signi cant toppling event was reported in 1996, and the current appearance of a large crack occurred in 2018. Also, many minor tension fractures are visible, indicating the landslide creep.
Signi cant, semicircular gaps are present at the contact between granitic bedrock and Quaternary deposits. Because of the previous disaster of the Attabad landslide (Gardezi et al. 2021), it was necessary to perform a full-scale investigation of Maicher Hill to prevent future catastrophes and plan mitigation.
The hazard is de ned as the likelihood of losses instigated in the event of possible disasters. The primary consideration for landslide hazard assessment is to estimate the landslide probability, collapse pressure, and run-out distance (Fell et al. 2008). Landslides always bring high-density material with them, and its movement path is a hazard for buildings or infrastructure that lies in the run-out path (Huang et al. 2017 Initial examination of engineering properties of site topsoil was determined by geological tests and geotechnical reports (Xiansong et al. 2018). Finally, all the calculated parameters were utilized in RAMMS to simulate the high-speed landslide hazard. The study site is an elevated moraines hill body. Several historical events show that moraines have owed like characteristics and can be categorized as highspeed landslides. Therefore, the RAMMS Debris ow module can simulate high-speed landslides (Bartelt et al. 2017;Huang et al. 2017). The nal model shows the movement path, velocity, accumulation thickness, the spatial distribution of the pressure in a 3D terrain. The results are essential for local administration to plan the mitigation to avoid hazards effectively quickly.

Geology And Geomorphology
Maicher hill is a direct product of glacier transport and accumulation. It lies on the contact of Eurasian block and Kohistan Island Arc, in the MKT (Main Karakorum Thrust) Zone (Fig. 1). The Eurasian union includes granite rocks with mica, schist, and quartzite. The middle part comprises the Quaternary Geology, consisting of the thin sub-sand soil and colluvial cover from adjacent rock and boulders. The schists and quartzite in colluvial material are highly sheared and fractured. The dominant rock debris in some places is partially lled with clay and silt. Due to seasonal variations and irrigation activities, the soil cover is eroded in patches, and sinkholes/cavities are exposed, representing the fragile and cracked subsurface (Fig. 2a). The borehole record and geological investigations suggest that the depositional of the moraines is uneven (Fig. 2b).
Generally, the stability of the moraines deposits is good until the weak zones and the free faces inside the soils are not combined favourably. However, these free faces are vulnerable to external factors like rainfall and earthquakes and make the overburden slide along the weak zones (Xiansong et al. 2018).
The upper part of the hill is mainly detritus, representing rapid transport sediments by glaciers and water melting; it ts mega granular soil. The bottom of the mountain is primarily breccia soil shaped by slow transportation and accumulation. The overall structure can be termed coarse-grained soil. Five springs ow across the slope, all aligned in the nearly north-south direction, indicating cracks and increasing the role of hydrogeology for slope failure mechanism. The large blocks toppling has occurred in the past; however, the unnoticeable movement of large blocks across the slope is causing subsidence and cracks. The signi cant toppling events date back to 1903 and 1996.
After heavy rainfalls in 2018, a gap emerged in the ridge, and the village residents informed the authorities. Also, many minor tension fractures are visible, indicating the landslide creep (Fig. 3a). Large semicircular cracks are present at the contact between granitic bedrock and Quaternary deposits (Fig.   3d). The bedrock exposed in nearby localities is highly disturbed and has an inde nite trend. The thickness of the Quaternary deposits varies from place to place along the river with no exposed bedrock.

Surface movement
We prepared a surface velocity map from the measurements of extensometers. Four extensometers with Time-domain re ectometry technology (TDR) were installed previously for early warning. The mean annual recorded displacements were interpolated and shown in Figure 4a. The sliding direction is predominant eastward, connected to an increase in movement, from 0-8 cm/year at the top of the landslide to more than 20 cm/year in its lower part. Relatively high activities have occurred in lower sections as compared to the upper part. After comparison with rainfall data, some observable movements (e.g., 2.127 mm) are directly related to rainfall (Fig. 4b). The topographical position of the monitoring stations and the recoded data along these stations are plotted (Fig. 4c). The geological and geophysical data acquired on the slope, predominantly in the southern part of the hill, provided insight into the nature of the slope and landslide characteristics.

Time Series Informatory
We acquired Sentinel-1 images between January 2019 and July 2020 and processed them with the SBAS-InSAR technique. In total, forty-seven ascending orbits were examined for the targeted duration. The maximum normal baseline was set to 45% from the critical baseline, and a full 60 days temporal baseline was picked by generating 220 differential interferograms. Goldstein's lter method was adopted to improve the signal-to-noise ratio of the interferograms. The phase unwrapping is done with minimum cost ow (MCF) and the 3D-unwrapping technique by keeping the coherence threshold to 0.35 and removing the undesired interferograms. The ground collection points (GCPs) were selected from the better-unwrapped phase area with the coherence value > 0.7. Re nement, re-attening, and inversion are applied to estimate the extra height and displacement-related information. The cubic model is applied to estimate the initial displacements, a reliable and commonly used model (Ur Rehman et al. 2020). High phase temporal (365 days) and low phase spatial (1000 m) ltering are applied to remove the residual atmospheric from the displacement. Subsequently, geocoding in the line of sight (LOS) direction with a resolution of 10 m is employed.

Geophysical eldwork
The dominant geology and data acquisition plan are presented in Figure 5. Four 380m electrical resistivity pro les parallel to the observed movement were run as part of eldwork (Fig. 5b). 15 m intervals were kept between the 2D pro les, and the achieved penetration depth is approximately 80 m. Also, ten additional pro les with 180 m length were made perpendicular to the southern slope section.
The penetrating depth of these pro les is almost 40 m. The exact locations of the start and endpoints of the pro les were recorded using the Real-Time Kinematic (RTK) GPS. The dipole-dipole con guration was adopted with 6 m electrode separation.

Electrical resistivity tomography
In theory, the resistivity of the subsurface water-lled void is lower than the surrounding materials. In contrast, the resistivity of an air-lled opening is higher than that of the surrounding materials ( projections of parallel lines are arranged in 3D space with constant separation to perform a 3D survey, achieving incredible depth (Kamiński et al. 2021). The measurements accomplished in this way can be combined into a 3D three-dimensional inversion software (e.g., Earthimager3D) with the spatial dataset for processing.

Geomechanical properties
The whole landslide development process and its characteristics are the interference of uncertain factors. Physical and mechanical tests show that moraines are characterized by dense structure, high-density, to their high composition degree of compactness. According to the statistics, moderately dense morainic debris density is usually 2.0 -2.3×10 3 kg/m 3 . The density of the dense gravelly soil is usually 2.3 -2.6×10 3 kg/m 3 . The density of moraine soil with more ne particles is generally 1.9 -2.1×10 3 kg/m 3 . Moraine congeries often bear good shearing strength. The cohesive force (∁) in some cases can increase 150 kPa, and the angle of internal friction (ϕ) can be about 35° -40°; when subjected to saturation, the cohesion is reduced to 60 -100kPa, and the internal friction angle reduces by 1° -3°. Furthermore, the value ∁ is much more sensitive to water and may be reduced by half than that of ϕ

Landslide dynamics
The preliminary examination of the general environment has indicated that signi cant landslide movements have taken place. The soil creep and ground surface stretching are evident on the Miacher slope. The stretching phenomenon is most commonly observed in non-cohesive materials that do not readily form or retain minor cracks. For recognition of a potential landslide, deformation rates were calculated using the remote sensing technique. The deformation is pragmatic to be intense around the central ridge and plotted on the deformation map of the hillslope that shows average velocities within the time series for the selected CTs (Fig. 6a). The range of VLOS is between -40 and -30 mm/year, with some noticeable values observed near the moraines and bedrock contact. The visible deformation of the Maicher Hill with the help of time series corresponds to its gradients. These different deformation rates may be related to the materials involved in the displacement, primarily moraine deposits, thus having intrinsic heterogeneity. The slope orientation and ridge gradient suggest that nal landslide debris accumulation ends in the Hunza River (Fig. 6a). The deformation length is around 450 m wide and shapes a prominent scarp. The dimensions of these deformations are the basis for estimating the landslide boundary, which is used in runout analysis. The bulky active deformation on the hillslope is evident from the large fractures, tensional cracks, and surface subsidence. The detailed deformational analysis from InSAR and site-speci c sensors data are explained in Figure. 6b.

2D ERT pro les
A series of geoelectrical pro les were run across the landslide, covering the ridge's prominent scarp and lower section. Two sets of pro les arrangements were made over the sliding body. The rst arrangement contains four pro les, each 380 m long and placed perpendicular to the prominent scarp. The separation between each pro le was kept at 15 m. The second arrangement contains ten pro les with a horizontal length of 180 m and covers the terrace of the ridge. Robust inversion was applied for better visualization of the geological contacts and the slip surface. Due to the difference in elevation and depositional environment, the material composition of moraines is different. The Miacher ridge depositional environment is coarse-grained granites, and gneisses and topsoil is sub-sand soil and colluvium. As a whole, the deposits are composed of massive particles and coarse grains. The size distribution of particles in the investigated moraine deposits is quite inclusive. The giant stone (> 2000 mm) and clay particles together could be found in the moraine deposits. These particles accounted for more than 85% of the total weight of moraine deposits. The content of ne particles is less than 15%. Therefore, the moraine deposits have a signi cant inhomogeneity coe cient, and resistivity data interpretation is challenging. All the pro les consist of high-resistivity values from 2500 Ωm to 4500 Ωm, and lowresistivity deals from 5 Ωm to 2000 Ωm ( Table 2). The subsurface interpretation is based on the geological pieces of evidence and visible changes on the surface. The general trend of resistivity is high.
The high-resistivity values are associated with the occurrence of strata breaks, while the low-resistivity represents glacio-uvial moraine deposits containing glacial till, sand, clay, and water content. On the other hand, the landslide initiation zone is marked on the electrical resistivity images as an extension of the prominent scarp in the subsurface. The distinctive resistivity anomaly is evident within the same geological material. The crack is extended almost 1.4 m deep on the surface; however, its effect in the subsurface is identi ed to 35 m in ERT pro les. The resistivity pattern resembles a step faulting down the slope. To conclude, in the rst pro les arrangement (Fig. 7), the strata with values from 4 Ωm to 2000 Ωm (navy-blue) is interpreted as galcio-uvial moraine deposits with variable clay and water content. The high-resistivity values represent two types of subsurface features, 2500 Ωm to 3300 Ωm (yellow-scarlet) as main scarp dimensions and above > 3500 (dark red) as granite bedrock. The second pro les arrangement was on a ridge terrace covering the irrigation elds (Fig. 8). The low resistivity values in these sections are mainly the grassed coved sub-soil with variable clay content, uncompacted nonuniform ll, and compact ll moraine.

3D ERT models
3D models for the transverse and longitudinal pro les were created. The two models represent the prominent scarp's subsurface limits and the volume of two high-resistivity sections of the moraine, separated from granite bedrock at the contact boundary ( Fig. 9c & 9d). The prominent scarp is visible at 210 m down the slope from the moraine and granite bedrock contact boundary. The high resistivity values in this section are associated with compact moraine. Additional isosurface modelling was performed to estimate the volume of strata that lies below the prominent scarp. The isosurfaces (> 800 Ωm & 800-1500 Ωm) presents the three-dimensional extent of the mixed less compact and rugged moraine strata (Fig. 9d). To ultimately comment on the compactness and geometrical distribution of the moraine in the lower section, four horizontal resistivity maps were constructed from the ten perpendicular pro les (P-5 to P-14) depicting the cross-section of the slope (Fig. 10a). The volume of this landslide suspectable strata is calculated to be 7.6 × 10 3 m 3 (Fig. 10b). At the hillslope cross-section, the shape of the unconsolidated sliding moraine is even more pronounced at depths of -4 to -40 m (Fig. 10c). On the -4 m surface map, the less compact moraine separates the relatively high compact moraine, marked as a sliding boundary. The effect of these less sliding moraines can be observed to the depth of -40m (Fig. 10c). Geologically, they are the result of the sub-soil region dominant with variable saturated and unsaturated clay content. The remaining part of the moraine shows mostly high resistivity and becomes dense after the depth of -40 m. Their spatial distribution is vital in estimating the landslide initiation depth at the lower slope section.

Numerical Simulation results
In 2010, an avalanche of moraine deposits fell into the Hunza river by travelling a vertical distance of 1060 m and created a dammed lake. The catastrophic avalanche killed 20 people and disconnected a considerable population from the rest of the country for four months. In the event of a Maicher landslide, the potential of a new dammed at the Hunza River is possible (Fig. 11). A terrain image and digital elevation model with 1.2 m resolution was developed for Miacher hillslope using UAV. The surface boundary is marked with dominant deformations calculated through extensometers and InSAR time series analysis (see section 3.2). The mechanical properties of moraine deposits were adopted from geologic and geotechnical reports (see section 3.5). The landslide initiation depth and volume are estimated with combined 2D and 3D ERT results. The nal numerical simulation of the Miacher Hill is the RAMMS model (Fig. 11a). The simulation time is 320 s. The maximum estimated accumulation thickness is 53.4 m, and the volume is approximately 174.4×10 3 m 3 . The length of the run-out path is calculated to be 821 m in numerical simulation results (Fig. 11b). The RAMMS simulation results also contain maximum velocity, pressure, and height distribution of possible landslide events. A brief review recommends that the debris velocity reaches up to 35.42 m/s, which is directly proportional to the slope steepness (Fig. 11c). The maximum pressure in some locations can reach 2508 kPa (Fig. 11d). The accumulation body will be distributed mainly over Hunza Rive. These simulation results and the actual condition might be nearly the same if the landslide is initiated at the contact boundary of moraines and bedrock.

Conclusions
The main objective of the research was to understand the nature and the movement of the Maicher slope. Also, this research will help the administration to quickly plan the mitigation to minimize the hazards effectively and reduce the risk to local infrastructure. Remote sensing, electrical resistivity, aerial photography, geological and geomorphological investigations were used to map the surface and subsurface features of the hill slope prone to landslide. The SBAS time series InSAR technique helped demarcate the landslide boundary in rugged terrain using Sentinel-1 temporal images. The rates of motion, in some places, raise concerns about the possible development of accelerations. The deformation trend from the sensors and InSAR technique shows the spatial dynamics of changes on the surface. This information is new since, so far, this landslide has not been investigated. Therefore, developing a simulation model for future landslide activity and estimating the maximum debris accumulation is essential to understanding the hazard level. This study shows the nal RAMMS landslide model with an extended landslide surface boundary estimated from InSAR deformations points. These deformations may not represent the actual landslide boundary. For example, the movement distribution is denser at the lower slope section than the upper slope section (Fig. 6a). However, the deformations are subsentaial enough to mark them as possible landslide boundaries. Hence, the boundary parameter in RAMMS generates the maximum potential high-speed landslide hazard.
Also, we use historical boreholes data drilled in the 1960s for mining purposes. Most of them do not lie on the electrical resistivity pro les. Indeed, this factor can affect the geological interpretation of ERT pro les. However, the primary scarp model is con dently interpreted by deploying the pro les over the exposed scarp on the surface. At the lower section, the terrace has developed several surface-exposed cavities. A widespread zone with resistivity values of < 50 Ωm was recognized, indicating the transport of ne material through in ltration. Low resistivities in the compact moraines suggest that the subsurface saturation is high and strata contain clay. Therefore, the lower section has shown a considerable amount of mass movement. The lower area is highly vulnerable to move in the event of heavy long-term rainfall and or earthquakes.
The study reveals that the movement's cause is the existence of favourable geological strata at the high altitude of "Miacher Hill" (weathered and saturated moraines). Environmental changes in the form of heavy rainfall or tectonic activity can induce and enhance the movement. The InSAR practice and installed sensors provide enough evidence that the hillslope is active and can collapse in the event of intense rainfall and earthquake. Further analysis is recommended to mark the landslide boundary con dently. The run-out distance and the signi cant extent initially identi ed are crucial for disaster prevention and mitigation, such as determining the safe avoidance range and plan mitigation.    The cross-section of Miacher hill (b) the historic borehole record of Maicher hill.

Figure 3
Structural changes on Maicher Hill. (a) the horizontal stretching at the contact point between moraines granite bedrock. (b) the measured extension of the trench at the lower section on 12-12-2019, and (c) 08-07-2021 (d) the visible, prominent scarp (e) tension cracks in a walnut tree due to shear displacement at the lower ridge section. Figure 4 (a) the surface velocity map interpolated from the installed sensor data (b) maximum movement recorded against rainfall in 2020 (c) the displacement recorded at the particular location of sensors.
Page 16/19   Results of ERT measurements of four parallel longitudinal pro les over the Maicher Hill.. Geological interpretations of resistivity pro les are labelled.

Figure 8
Results of ERT measurements ten transverse pro les across the Maicher Hill. Geological interpretation of resistivity pro les is labelled at each pro le.