Increase in occurrence of large glacier-related landslides in the high mountains of Asia

Globally, a large number of glaciers are retreating due to global warming. Along with climate change, glacial melting has been identified as one of the main triggers of landslide activity in high mountain areas. Evaluations of the triggered mechanism alone do not provide comprehensive insight into the overall impact of glacier accumulation and ablation on landslide-induced denudation. To investigate recent trends, we built landslide and glacier datasets for the HMA using a Landsat time-series covering the past 21 years (1998–2018). Landslides that may have been caused by major earthquakes were identified and removed, leaving an inventory that is used to explore changes that may be related to climate change. Our results show a shift in the frequency–area distribution that indicates an increasing trend of large landslides in the HMA over the last decade. A decline in glacier area is associated with the increase in landslide area.


Data and methods
Landslide inventory. Our primary dataset is a glacier-related landslide inventory, which can be very useful for risk assessment 22 . Some scholars have conducted field surveys in the HMA and generated previous landslide inventories 23 . However, these databases are mostly historical data without detailed and reliable timestamps. Remote sensing offers the possibility of addressing this data gap, by building detailed time-series of landslide inventories.
Landslides can be detected on satellite series by changes in land cover (vegetation distribution, rock outcropping exposure or soil degradation) 15,24 . Remote sensing satellite data with high resolution an extensive archive of historical acquisitions, and a generally vertical field of view that provides an almost planimetric perspective, are particularly valuable for generating time-series data. Sentinel and Landsat are sensors that have both used by scholars to study landslides 25,26 . Although the spatial resolution of Landsat is slightly coarser than that of Sentinel, Landsat provides a longer observation time. To avoid biasing our time series, with more landslides being identified from the more recent high-resolution imagery, we use only Landsat imagery.
For our study, Landsat imagery was obtained from the U.S. Geological Survey (USGS) EarthExplorer website. Images were obtained for each year for the period from 1998 to 2018. The Landsat data include images from the Thematic Mapper (TM), Enhanced Thematic Mapper (ETM+) and Operational Line Imager. The spatial resolution of these data is 30 m. The entire study area is covered by 17 Landsat path/row scenes. For each scene, cloud-free images at the end of summer (September through October) were selected. Additional images were selected during the winter season of each year so that each year's imagery includes one imagery at beginning or end of the year.
Landslides were manually identified by photo-interpretation of the scenes, and by comparing imagery between dates. Landslide scars and deposits resulting from landslides can be identified on Landsat imagery from contrasts in landform and land cover with the surrounding areas, and confirmed by lack of these distinctive features in the summer image of the previous year. The multi-temporal database was used to constrain the date of the landslide, which was labeled with the year the landslide is first visible. The digitization of the landslide polygons was performed in ENVI software. SRTM3 digital elevation data is used to calculate the aspect and slope of a landslide.
In the process of mapping landslides, we tried to map all the affected area of each landslide. However, the failed areas of some landslides were not clearly visible because of shadows or a lack of contrast between ground disturbed by the landslide, and undisturbed snow/ice or bare rock. This leads to the possibility that the landslide area may include the failure area and the runout area, or it may include only the runout area or the accumulation area.
We classified the types of landslides according to the criteria of Varnes (1998) 27 . Due to the limitation of Landsat image resolution, it is difficult to distinguish the type of movement. Therefore, we also used higher spatial resolution imagery available in Google Earth to supplement our Landsat interpretation. The landslide material types were rock, debris and moraine, as illustrated in Fig. 1.
Because we focus on the possible effects of glacier retreat on landslides, we removed from our inventory landslides that could potentially have been triggered by earthquakes and that were more than 50 km from glaciers. Keefer 28 has proposed a typical maximum distance for the occurrence of seismically-induced landslides from the epicenter; a distance that varies as a function of the magnitude of the earthquake and represented earthquakes (Fig. 2). We therefore used the USGS earthquake hazards website to identify all major earthquakes near or within the study area during the 20-year interval studied, and then removed from our inventory any landslide that could have occurred at the time of the earthquake and was located within the threshold of distance. We chose the disturbed landslide because it is the most sensitive to earthquakes, which means that the earthquake has a large impact distance on such landslides. We found four earthquakes (Table 1) and an associated four landslides that met the magnitude and distance criteria shown in Fig. 2. These landslides were removed from the inventory.
Glacier mapping. The second major dataset that we generated for this study is a 20-year time series of maps of glacial extent for the study area. As with the landslide inventory, we used Landsat remotely sensed data. The Landsat sensors TM, ETM+ , OLI capture the distinctive spectral reflectance of snow and ice in the VNIR (visible and near-infrared) band and SWIR (shortwave infrared) band that is used to separate it from the surrounding terrain. Typical Landsat spectral measures include spectral band ratios such as TM3/TM5, TM4/TM5, and the NDSI (normalized difference snow index). In this study, we used the Landsat NDSI index to obtain glacial boundaries. The NDSI index is defined as: www.nature.com/scientificreports/  www.nature.com/scientificreports/ where Green and SWIR represent the values of the green and shortwave infrared bands for Landsat series images. We set 0.4 as the extraction threshold, which is commonly used by scholars to extract ice from remote sensing images 29 . Spectral ratio methods are effective in delineating clean ice. However, seasonal snow, cloud and frozen lakes may be misclassified as glaciers and require manual editing of the classification. we therefore first classified ice using the NDSI. Then, we clipped our data based on the Randolph glacier inventory 30 , to reduce snowfall and cloud mass errors. Finally, manual interpretation was used to remove some errors associated with glacial lakes and around the edges of images.

Probability-area distribution of landslides.
A number of complex natural phenomena exhibit powerlaw frequency-area relationships, including earthquakes, which are considered a classic example of such phenomena. Landslides are thought to be another natural hazard that exhibits power-law frequency-area relationships under a wide variety of circumstances 31 . Some scholars question whether small landslides follow the power distribution, because small landslides sometimes evolve into larger landslides. There is also a problem that smaller landslides may be difficult to identify, or many individual small slides appear to be one large slide 32 . The main focus of this paper is medium to large landslides, so these concerns are not directly relevant to our work.
The probability-size distribution of landslides in the study area is given by where p(A L ) is a probability density function, δN L is the number of landslides with areas between A L (km 2 ) and A L + δA L , δA L (km 2 ) is based on a log scale and N LT is the total number of landslides in an inventory. The frequency density of landslides, f (A L ) , is given by: The frequency-size distribution of landslides in the study area was compared with those for various event magnitudes proposed by Malamud et al. (2004) 33 to assess the nature of a landslide event. The p value represents the probability density of landslides with area A L (km 2 ) as follows: where ρ is a parameter controlling the power-law decay for medium and large landslides, a (km 2 ) is a parameter controlling the location of the maximum probability distribution, s (km 2 ) is a parameter controlling the exponential rollover for small landslides, and Ŵ(ρ) is the gamma function of ρ . Malamud et al. 33 also proposed a magnitude scale for a landslide event, m L , as follows: The combination of Eqs. (2)-(5) provides the frequency density of landslides linked to the magnitude scale of a landslide event. www.nature.com/scientificreports/ For our study, the power-rate distribution was obtained using the R language.

Results and discussion
Factors influencing identification of landslides. Landslides can be identified in alpine and glacial regions by (1) high contrast compared to surrounding snow and ice, (2) the influence of landslide on river, and (3) lobate forms typical of rock-avalanche deposits 15,34 . Local variations in tone, texture or pattern, and the presence of lineaments can also be used to infer slope instabilities 35 . However, there are some factors that affect us to see the certain landslide in the identification process. The factors affecting landslide identification include landslide characteristics, snow cover, and the quality of remote sensing image.
Landslides are a complex movement process. Landslides include failure area, transportation area and accumulation area. Due to the influence of image quality (resolution, clouds, etc.) and the geological characteristics of the alpine region (glacier development, lack of vegetation), it is difficult to distinguish some slides which didn't cause optical characteristic variation in the image. They include some small area debris flow and moraine failure, and rock avalanches that have not formed significant accumulation. In our identification process, we found that the influence of landslide on river is beneficial to our identification work.
The influence of external factors mainly comes from the influence of snow cover. According to our identification process, although a landslide occurs in winter and is covered with snow (Fig. 3A,B), it can be seen after the snow melts as the area is large or the accumulation is large (Fig. 3C). Coe et al. 15 mentioned that the accumulation of landslides greater than 0.5 km 2 would be obvious.
The quality problems of Landsat series satellites mainly come from cloud cover, seasonal snow cover and so on. We try to ensure a late summer and cloudless image every year to identify landslides. But sometimes it can't meet two conditions. Hence, we selected multiple images from the same year to identify the landslide. In particular, there is a problem with the black stripe error in the Lansat7 after it broke down in 2013. We think it is inevitable that we will miss some of the landslides that occur within these bands or most of them. At the same time, we can be sure that some large areas of the landslide, even if the impact of the strip, we can identify and determine its occurrence.
Activity of landslides and climate change in the HMA. A total of 127 landslides were detected in the Landsat images of the study area, covering the period from 1999 to 2018 ( Fig. 4 and Supplementary Information). The landslides are mainly concentrated in the Karakoram Mountains, eastern part of the Pamir Mountains, western Himalayas and south of the Hindu Kush. Based on the source area the landslides were divided into three categories. A total of 72 rock landslides were identified, 45 debris landslides, and 10 moraine landslides. Table 2 and Fig. 5 summarize the landslide characteristics. The average area of the 127 landslides was approximately 58 ha, with 6.35 landslides occurring per year. The average elevation of the landslides is 2,966 m, and the average slope is 10.88° at the center point of the failure area. Most of the mapped landslides were at elevations more than 2,000 m above sea level (Fig. 5A). Landslide source areas had more slope directions (aspects) between 90° and 270° (i.e., southwest-to southeast-facing slopes, Fig. 5B).
The landslides of rock type are widely distributed in the HMA, and concentrated in the western Kunlun Mountains, although they also occur frequently in the Hindu Kush and the western Himalayas. Debris landslides are mainly distributed in the West Kunlun Mountains and the Hindu Kush, with additional occurrences in the Pamirs and western Himalayas. As expected, morainic landslides occur at higher elevations, including the Karakorum Mountains, the western Himalayas, the Kunlun Mountains, and the Pamirs, at an average elevation of 4,139 m. We also analyzed the distribution of large landslides, defined as landslides that affected areas greater www.nature.com/scientificreports/ than 2 km 2 . A total of three large landslides were found in the Pamirs, four in the western Himalayas, and two in Karakorum. A graph of annual landslide area is show in Fig. 6. The Thompson Tau method 36 is used to find the outliers with the critical probability alpha value is 0.01. The result indicates that the landslide area points in 2003, 2010 and 2016 are outliers, as shown in Fig. 6. After removing the outliers, linear-fit line (R 2 = 0.26) shows that the annual landslide area in the HMA are increasing (Fig. 6).   www.nature.com/scientificreports/ The area of landslides was greatest in 2016. The average annual temperature was unusually high during this year, and the glacial area had notably declined in the previous 3 years (Fig. 8E). To provide additional insight into the impact of climate change on landslides, we compared the annual landslide area data with annual El Niño data for the study years 37 . We found that the landslide area outliers tend to coincide with El Niño years. The one exception is 2006, an El Niño year with no associated peak in landslide area. However, 2006 was associated with generally lower average temperatures in Fig. 8E. Overall, however, the time-series of landslide area, and particularly the outliers, appears to show periodicity, and the trend in landslide area is increasing. Scientists have hypothesized that as El Niño becomes more frequent and stronger, alpine glaciers are also being affected 38 . Warmer temperatures caused by El Niño can accelerate melting of glaciers, and the increase of rainfall also increases the conditions that favor the occurrence of landslides.

Shifts in landslide frequency-area distribution.
To investigate whether landslide characteristics have changed over the past 20 years, we divided the landslide data set into two datasets: landslides that occurred before 2009, and those that occurred in 2009 or later. Table 3 shows that both the number of landslides and the average area of landslides increased in the second 10 year period. The difference in the average slope of the landslides was not significant.
The power distribution of landslide data in the study area over the entire 20 years of the landslide inventory has a decay factor ρ of approximately 1.14, which is in the normal range of previous studies. Although there is some variability in the findings of previous research, most prior landslide data sets follow noncumulative powerlaw frequency statistics, and the range of ρ values is approximately 1.5 ± 0.5 32 . However, it is important to note that previous studies have obtained different power distributions, fitting different tail attenuation coefficients. Furthermore, the data sets used in these studies come from different environments and the associated earthquakes are assumed to have been triggered by different factors, including rainfall, earthquakes and snowmelt.
The power-rate distribution for the landslides separated into the two 10-year periods is shown in Fig. 7. Both periods show a distinct power distribution. The decay factor in the second 10 year period (ρ = 0.85, 95%CI: 0.70-1.01) is smaller than that of the first 10 year period (ρ = 1.32, 95%CI: 1.02-1.61). This indicates that the probability of large and medium-sized landslides in the second 10 year period is higher than the first period. The confidence intervals of the two values of ρ do not overlap, and thus the difference is statistically significant.

Relationship between landslides and glaciation.
For the entire study area, we used the Randolph glacier outlines as a mask to minimize misclassification due to factors such as seasonal snow. However, seasonal snow can also lead to errors in the calculations. In order to evaluate the error associated with the glacier delineation, we compared automatic extraction of glacial area with a semi-automatic approach that was regarded by scholars as a relatively accurate method 39 . The semi-automatic method uses the ratio threshold method to generate approximate results, when are then modified by hand, using visual interpretation by an expert interpreter.  www.nature.com/scientificreports/ We selected data from 3 regions for data verification ( Table 4). The calculated difference is less than 10%, and the average difference is about ± 6.3%. For our study, we think the error of the glacier area results is acceptable. The glacier area in the HMA has shown a downward trend over the past 20 years. Before 2005, the glacier area decreased and increased repeatedly, which may be related to the monsoon climate. After 2005, the area of glaciers decreased significantly. The trend of glacial degradation over the past decade was significantly stronger than that in the previous decade.
For a comparison of the glacial extent between years, we chose the remote sensing image from 2002 as the basic reference because the imaging times for the mosaic image in 2002 were not very different. We selected the glacier area in four subsequent years (2006, 2009, 2013 and 2017) to compare with the reference to analyze the change in glaciers in the study area ( Fig. 8A-D). We found that except for the Karakoram Mountains, which exhibited only small areas of glacial retreat over the past 2 decades, the retreat in the other areas was more extensive. An anomaly of the Karakoram Mountains was found due to the relatively high elevation 40 . In addition, climate anomalies and the strong summer monsoon climate in the region have brought additional moisture to the Karakoram Mountains, leading to increased snow 41 . Figure 8E summarizes the total area of glacial extent of the 20-year period studied. The figure also shows the mean annual temperature, averaged over the region. MAAT was calculated based on MODIS surface temperature product MOD11A2, which has a spatial resolution of 1000 m and a temporal resolution of 8 days.
The HMA glacial area has had a downward trend over the past 20 years (Fig. 8E). Looking at the record in more detail, the glacial area decreased prior to 2005, and then increased, possibly in response to variations in the associated region's monsoon climate. After 2012, the area of glaciers decreased notably. The glacial area showed little no overall change in the first decade, and a notable decline in the second decade. The temperature record does not show evidence of warming over the period 2000 to 2018. The year 2000 was unusually warm compared to the subsequent 18 years. Even excluding 2000, there is no clear evidence of a warming trend over the remaining years.
However, glacial area can be assumed to be an environmental proxy that integrates climate over time and space. Glacial area could therefore be a useful predictor of slope stability. Glacier area does indeed show a negative association with landslide area (Fig. 9), indicating landslide area is increasing as glacial area declines. Glacial retreat is an indicator of climate change, and our results support the hypothesis that warming of the region is associated with an increase in landslide occurrence. However, the small number of samples and the scatter indicate notable uncertainty, and the association therefore requires further research.  www.nature.com/scientificreports/

Conclusion
In this work, we generated an inventory of 127 landslides in the HMA using Landsat data covering the 2 decades between 1999 and 2018. Annual glacial area maps for the same 20 years were also generated. In addition to the impact of image resolution, the rock surface and ice cover make landslide identification difficult in the HMA, which may lead to the number of landslides being underestimated. Nevertheless, the landslide and glacial area maps provide detailed spatial and temporal information. Based on our custom inventories, we studied the interaction between landslides and glaciers in the study area. The following major conclusions can be drawn from this work: www.nature.com/scientificreports/ 1. In the HMA, landslides are widespread. Over the past 2 decades, the area affected by landslide disasters has increased. Both the number of landslides and the affected area was larger during the most recent 10 years (2009-2018) compared to the previous 10 years (1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008). The area of landslides was anomalously high in three out of four El Niño years. 2. Landslides in the HMA follow power-law (fractal) frequency statistics. The attenuation coefficient of the 10 years from 2009 onwards decreased compared to the previous 10 years, indicating that the probability of occurrence of large landslides has increased. 3. Glacial area in the HMA has shown an overall downward trend over the past 20 years. Except for the Karakoram Mountains, most glaciers throughout the study area have been retreated. As the glacial area decreases, the area of annual landslides has increased.
In summary, retreat of glaciers in the HMA of 20 years appears to be associated with more frequent landslides, and larger landslides. A weak negative correlation between annual landslides and glacial areas is evident. Although local influences may have triggered landslides in the HMA, glacial retreat may be a useful proxy for aspects of climate that control slope stability. If climate warming continues, the area of glaciers will further decrease. As a result, the probability of large landslides in the HMA may continue to increase.
Received: 23 June 2020; Accepted: 5 January 2021 Figure 9. Graph of landslide area and glacier area. GA represents glacier area, and LA represents landslide area, and ln represents the natural log. The black dotted line is the linear trend line.