Rock glacier inventory of the western Nyainqêntanglha Range, Tibetan Plateau, supported by InSAR time series and automated classification

The western Nyainqêntanglha Range on the Tibetan Plateau reaches an elevation of 7,162 m and is characterized by an extensive periglacial environment under semi‐arid climatic conditions. Rock glaciers play an important part of the water budget in high mountain areas and recent studies suggest that they may even act as climate‐resistant water storages. In this study we present the first rock glacier inventory of this region containing 1,433 rock glaciers over an area of 4,622 km. To create the most reliable inventory we combine manually created rock glacier outlines with an automated classification approach. The manual outlines were generated based on surface elevation data, optical satellite imagery and a surface velocity estimation. This estimation was generated via InSAR time series analysis with Sentinel‐1 data from 2016 to 2019. Our pixel‐based automated classification was able to correctly identify 87.8% of all rock glaciers in the study area at a true positive rate of 69.5%. In total, 65.9% of rock glaciers are classified as transitional with surface velocities of 1–10 cm/yr. In total, 18.5% are classified as active with higher velocities of up to 87 cm/yr. The southern windward side of the mountain range contains more numerous and more active rock glaciers. We attribute this to higher moisture availability supplied by the Indian Monsoon.

likely to increase over the coming decades due to global glacial retreat. 12 Rock glacier inventories have been used to estimate the regional lower permafrost limit, 13,14 though subsurface ice may still be found in favorable conditions at much lower elevations. 15 Rock glaciers have displayed accelerated motion rates over the past decades in the European Alps, northern Norway, the Carpathian Mountains, and the Tien Shan. [16][17][18] This can potentially lead to hazardous mass wasting processes due to destabilizing rock glaciers. 19,20 The potential of rock glaciers to act as climateresistant water storages combined with the importance of the Tibetan Plateau (TP) as a source of fresh water to $1.65 billion people in Asia 21 and the amplified warming effect of climate change on the TP 22 emphasize the necessity to study rock glacier distribution on the TP.

| Rock glacier morphology
The precise definition of what constitutes a rock glacier and its potential origins have been discussed but remain contested. Rock glaciers are described as the visible expression of cumulative deformation by long-term creep of ice/debris mixtures under permafrost conditions. 23 Rock glaciers may develop in periglacial, paraglacial, or glacial environments. 24 They may exist in areas where permafrost is unlikely outside the rock glaciers themselves when gravity-driven creep extends the landform beyond the expected periglacial zone. 25 They are typical landforms of high mountain environments and can extend up to multiple kilometers in length and several hundred meters to a kilometer in width.
Rock glacier debris often originates from avalanche processes and rock falls, which can be triggered by continuous headwall weathering, sudden heavy precipitation, or earthquakes. 3 Rock glaciers may also develop from frozen moraines or debris-covered glaciers. 10,26,27 The latter often include glacial ice within the rock glacier body. Rock glaciers developing without a connection to a glacier form ice within the debris by freezing rain or melt water and/or by burying accumulations of snow and ice. 28

| Rock glacier activity
Rock glaciers that contain ice are categorized as active if they display motion and as inactive if they are unmoving. Unmoving rock glaciers without ice in their bodies are referred to as fossil or relict rock glaciers. 1 Active rock glaciers often display horizontal velocities of centimeters to decimeters per year but may reach velocities of several meters per year. 16,18 Velocities are influenced by changes of ground temperature as well as moisture availability, which may accelerate or decelerate a rock glacier on a decadal scale. [29][30][31][32] Rock glaciers also display strong seasonal variations in their surface motion in many cases, with higher velocities in summer and autumn compared to winter and spring. [33][34][35]

| Creating rock glacier inventories
The spatial coverage and comparability of rock glacier inventories are lagging behind inventories of glaciers on a global scale. 12 This is due to a number of reasons including: (a) the necessity of highresolution satellite data to identify rock glaciers and to distinguish them from surrounding landforms; (b) the difficulty in judging their status of activity and their ice content; and (c) the subjectivity of identifying and outlining rock glaciers and the dependency on the analyst's experience to distinguish them from similar landforms such as debris-covered glaciers and ice-cored glacial moraines. 26,36 The availability of freely accessible global map providers such as Bing Map, Google Maps, or Zoom Earth are alleviating issue (a). However, remote areas often display lower temporal and spatial data coverage compared to urban centers and the positional accuracy may vary. Microwave remote sensing data and Interferometric Synthetic Aperture Radar (InSAR) techniques have been applied successfully to determine rock glacier activity in some large-scale inventories in recent years (Table 1). This addresses issue (b). InSAR techniques allow monitoring surface motion on the order of a few millimeters to decimeters per year depending on the technique used, the wavelength and the revisit time of the satellite constellation. 37 Issue (c) has led to problems in the comparability of rock glacier inventories of different regions due to varying inventorying methodologies. Identifying a landform as a rock glacier and subsequently manually creating its outline is subjective. This leads to large variabilities in both the number of rock glaciers and their surface areas when comparing the rock glacier inventories created by different analysts. 36 The recently formed action group on "Rock glacier inventories and kinematics" of the International Permafrost Association (IPA) aims to address these issues. They explore the feasibility of developing widely accepted standard guidelines for inventorying rock glaciers on a global scale, including information on their kinematics. 38,39 Furthermore, varying techniques to detect rock glaciers automatically based on remote sensing data have achieved promising results. [40][41][42][43][44] In this study we follow the practical guidelines of the IPA action group to generate a manual inventory and compare it to automated classification results to work towards a global rock glacier inventory.

| Purpose of this study
We present an inventory of actively moving rock glaciers within the western Nyainqêntanglha Range on the southeastern TP ( Figure 1). We combine a manual approach with an automated classification algorithm to generate an exhaustive rock glacier inventory. Following the IPA guidelines ensures the comparability with future rock glacier inventories following these guidelines.
Combining the manual approach with a semi-automated classification increases the reliability of the inventory. This inventory is the first of its kind in this mountain range and will help to fill the gaps in a global rock glacier inventory. We use Sentinel-1 satellite data of 2016-2019 processed with InSAR time series techniques to estimate the interannual surface velocity and infer rock glacier activity. In the following we introduce our study area and related work of previous studies. We then explain the data and methods we used to generate the inventory. Finally we present the inventory followed by a discussion on the limitations of our methods, potential implications regarding the regional climate, and finally our conclusions.
2 | STUDY AREA: THE WESTERN NYAINQÊNTANGLHA RANGE

| Climate
The western Nyainqêntanglha Range is located in the southeastern center of the TP (Figure 1) with an elevation of $4,500-7,162 m. The median elevation of the study area is 5,369 m with an interquartile range of 5,109-5,632 m. The study area is about 230 km in length T A B L E 1 Rock glacier inventories with velocity information based on differential interferometric synthetic aperture radar (DInSAR) techniques  56 The mean annual air temperatures at those meteorological stations during the same time periods were 1.6, À0.6, and À5.9 C respectively. The contribution of increased glacial meltwater to the increase in lake level was estimated to be 10%

| InSAR processing
Interferograms are a spatial representation of the phase difference between two SAR acquisitions and can be used to determine LOS surface displacement. 37 The phase stability, so-called coherence, is often used to represent the quality of an interferogram. A variety of differ- and is included in SARScape as an option of the SBAS processing F I G U R E 3 The workflow we followed to generate the manual rock glacier inventory and the automated classification results [Colour figure can be viewed at wileyonlinelibrary.com] chain labeled "disconnected blocks". Intermittent SBAS interpolates time periods where the coherence of a pixel drops below the chosen coherence threshold in some interferograms. This may occur due to heavy rainfall, snow cover, or large displacement. Intermittent SBAS results in a greatly improved spatial coverage compared to the original SBAS algorithm for partially vegetated study areas. 72,73 In our study area it helps to interpolate time periods in spring and autumn when the transition between snow-free summers and frozen ground in winter causes decorrelation. For our analysis we chose a coherence threshold of 0.2, which is lower than the threshold of 0.3 used by others for intermittent SBAS methods. 72,73 The lower threshold is justified in our opinion, as good spatial coverage is most important for this study. Pixels need to exceed the coherence threshold in at least 60% of all interferograms during our 3-yr observation period. A pixel may therefore contain up to 40% interpolated information, though the median value of all pixels on rock glaciers is only 2%.
We use Sentinel-1 Level-1 single look complex data from the interferometric wide swath mode to perform our InSAR time series analysis. The wide swath mode has a ground resolution of 20 m in azimuth and 5 m in range. 74 We used a multilooking factor of 4 in range

| Removal of large-scale phase delay
It is necessary to carefully select stable reference points in the study area to achieve InSAR results with a high accuracy. 37 Poorly selected reference points can introduce seasonal artifacts and interannual shifts leading to misinterpretation of displacement patterns. No longterm data from the global navigation satellite system is available for this area to act as reference points. We therefore selected probable stable areas on the outskirts of the mountain range for our reference points. These 385 reference points were placed along both the northern and the southern side of the mountain range. We only selected points without visible signs of displacement and with high coherence throughout all interferograms. We avoided areas near the valley bottom during the selection of our reference points, as they probably experience heaving and thawing related to freezing and thawing of the ground. 64, 65 We applied an atmospheric high pass filter of 100 days and a low pass filter of 1,200 m to reduce atmospheric artefacts.

| Downslope projection
We projected the LOS surface velocity of both geometries in the direction of the steepest slope. 75 The following areas were masked prior to the projection: (a) areas of glaciers where free ice is visible; (b) areas with a slope <2 , as these areas are generally too flat to support downslope motion 76 ; and (c) areas with a low slope and strong vertical displacement signal (based on a decomposition of ascending and descending time series results). This displacement is probably not gravity-driven and therefore not directed along the steepest slope.
(d) We masked wetlands and the immediate surrounding of rivers and lakes. Changing soil moisture in these areas may lead to a false identification of displacement. 77 To perform the downslope projection, we calculated a sensitivity coefficient with a value between 0.2 and 1 and divided the LOS velocity by this coefficient to receive the downslope velocity. A sensitivity coefficient of 1 represents slopes where the angle between the LOS vector and the downslope direction is zero. The sensitivity of the satellite to displacement in that direction is therefore very high. A sensitivity coefficient of 0.2 represents areas where the angle between the LOS vector and the downslope direction is large. The sensitivity of the satellite to displacement in that direction is therefore very small.
Ascending and descending data sets have different LOS and their respective sensitivity coefficients are therefore also different when observing the same area. We included the spatial distribution of this sensitivity coefficient in the Supporting Information. After projecting both ascending and descending data sets into the downslope direction, we combined both data sets according to the following criteria: 1. We only use the downslope projection of the more sensitive geometry in areas where the sensitivity coefficient of the more sensitive geometry is at least 50% larger than the sensitivity coefficient of the other geometry and is larger than 0.2. September. 79 The irregularity of these dates may introduce a bias but they were necessary to produce an NDVI map with near total coverage of the study area. We masked snow-covered areas as well as clouds and their shadows prior to the calculation of the NDVI. The NDVI data sets were then merged into one raster with pixels present in multiple NDVI data sets receiving their mean value. The free ice glacier outlines were derived from a cloud-free Sentinel-2 acquisition of January 30, 2018 with a supervised maximum-likelihood classification and checked manually.
We use the 12-m-resolution TanDEM

| Manually generating the rock glacier inventory
We followed the guidelines of the IPA Action Group on "Rock glacier inventories and kinematics" to identify, outline, and classify all actively moving rock glaciers in the study area. 38

| Supervised maximum-likelihood classification
Manually outlining rock glaciers is inherently subjective and dependent on the experience of the analyst. Multiple studies in recent years have therefore evaluated different classification algorithms and varying input parameters to automatically identify and classify rock glaciers. [40][41][42][43][44] In addition to the manual approach described above, we also performed a pixel-based supervised maximum-likelihood classification. We masked glaciers, flat areas, wetlands, and the immediate surrounding of rivers and lakes (see section 4.3). We used the ENVI software to perform the automated classification based on the surface features described in We also selected 20 nonrock glacier areas, to assess the number of pixels falsely classified as rock glacier. We evaluated different feature combinations to determine the most suitable selection to automatically classify rock glaciers in our study area. All features were resampled to the 12 Â 12-m resolution of the TanDEM-X DEM. For the initial evaluation we performed the classification with the following five features: elevation, slope, NDVI, slope variability, and height error. Elevation, slope, and NDVI scored highly in a similar study by Brenning et al. 40 Slope variability was to represent surface roughness.

| Manual rock glacier inventory
The rock glacier inventory contains a total of 1,433 rock glaciers (Table 3; Figure 5). The IPA guidelines classify rock glaciers with a

| Maximum-likelihood estimation
The best result was achieved by combining the following nine features: downslope velocity, elevation, slope, NDVI, slope variability, height error, TPI with 200 m radius, and Sentinel-2 bands 2 and 3. We give an overview of these features in The entire inventory displays a TP of 69.5% and an FP of 10.4% (example area in Figure 7). In total, the automated classification predicts a combined rock glacier area of 428.4 km 2 . This is equivalent to an overestimation of 243.0% when compared to the 124.9 km 2 of manually outlined rock glaciers. After checking the manual outlines against the automatic classification, we identified 18 additional rock glaciers. To assess the relative importance of the features to each other, we repeated the classification with only eight of the nine features (leaving out a different feature every time). The resulting feature ranking is presented in Table 5. Vegetation growth, snow fall, and soil moisture variations between two SAR acquisitions may cause decorrelation or misinterpretation of InSAR data. 16,72,73,77 The high altitude of our study area means that large vegetation, such as trees, are not a problem. Actively moving rock glaciers are sparsely vegetated, and we therefore consider vegetation to have little impact on our InSAR results. 78 The winter climate in our study area is semi-arid. 53 Optical satellite imagery shows that snow fall occurs mainly between autumn and spring but our field work has shown that it may also occur in summer. We observe the lowest overall coherence during the spring and autumn months, which we attribute to the transition between frozen and Varying types of atmospheric properties may delay the microwave emitted from the satellite, leading to major uncertainties in InSAR analysis. 86 The two dominant types we observe in our study area are turbulent delay and seasonal stratified delay. Turbulent delay is associated with short-term variations in the atmosphere and can be considered random in time. We observe the effect of turbulent delay in our individual interferograms, where they lead to false identification of surface displacement of up to half a phase cycle. As turbulent delay is random in time, its effect on longer observation periods becomes negligible. We therefore consider the effect of turbulent delay on our time series results to be relatively small. Seasonal stratified delay may lead to an apparent seasonal cycle in the surface displacement if there is a large elevation difference between the reference points and the area of interest. 87 This does not affect the interannual surface velocity and we therefore expect its effect to be negligible for this study.

| Rock glacier inventory
In  southern side of the mountain range, which contains 91.3% of all active rock glaciers. The majority of the moisture available in this region is delivered by the Indian Monsoon from the south and southeast. 49 The southern (windward) side receives most of this moisture, which results in larger glaciers compared to the northern side. 54 It is likely that the higher moisture availability also facilitates the higher number and activity of rock glaciers.

| Semi-automatic classification
The classification algorithm used in this study is a supervised maximum-likelihood classification. Initial testing had shown that this algorithm performs better than other classification algorithms of the ENVI environment in our study area. We did not assess the performance of other classification algorithms outside the ENVI environment. ENVI is a popular remote sensing environment and easy to use without in-depth machine-learning knowledge. However, this also precludes the use of more sophisticated machine-learning algorithms.
The algorithm used in this study performs the classification based on individual pixel values. Other algorithms, such as convoluted neural networks, take into account more complex patterns and textures. 42,43 These algorithms would probably perform better to separate rock glaciers from other debris-covered landforms, such as debris-covered glaciers.
The following features are used for the automated classification: but it also imposes the limitations associated with InSAR onto the automated classification. We discuss these limitations in detail in section 6.1 and will therefore focus here on the most important limitations: InSAR sensitivity and spatial data coverage. InSAR techniques have poor sensitivity towards displacement with a strong north or south component. Rock glaciers moving mainly into those directions are therefore probably excluded from the classification if their overall velocity is slow or if they are inactive. High spatial coverage of InSAR surface velocity data is necessary to apply our approach to a study area. This can be difficult to achieve in high mountain areas, where heavy snow fall may lead to decorrelation of InSAR in winter. 16 Using exclusively snow-free summer acquisitions or exploiting the low coherence associated with moving landforms could be an alternative in those regions. 88 Our automated classification approach is likely to work well in other regions with dry winter climates, such as the central and northern TP. Stable winter conditions in the northern Tien Shan make this region another potential candidate. 18 The supervised maximum-likelihood classification produced a TP