Variability of the snowline altitude in the eastern Tibetan Plateau from 1995 to 2016 using Google Earth Engine

Abstract. The glacier snowline can be used as an indicator of a glacier’s equilibrium line, which is a pivotal parameter for studying the effect of climate change on glaciers. However, the relationship between snowline altitudes (SLAs) and climatic regime, as well as the comparison between different glacier types, has received less attention. Using Google Earth Engine, we first developed an automated algorithm that employs the Otsu thresholding method to distinguish snow-covered areas from clean ice on a near-infrared band of Landsat imagery available from 1995 to 2016 and further to delineate glacier SLAs in the three regions of the eastern Tibetan Plateau (TP). The three study regions, Sepu Kangri (maritime), Bu’Gyai Kangri (continental), and western Qiajajima (continental), in the eastern TP have different climate regimes and are on a latitudinal transect from south to north. We then investigated the impacts of climatic factors on the SLA and its variability over the period studied. The results over the eastern TP indicate that (1) the SLA increased by 94, 55, and 49 m from south to north during the 22-year period, with the SLA variations of maritime glaciers being the most pronounced; (2) the southern maritime glaciers were mainly affected by precipitation, whereas the northern continental glaciers were influenced by temperature. Owing to the difference in primary climatic factors affecting snowlines, continental glaciers were found to have higher SLAs on the south slope, whereas maritime glaciers had higher SLAs on the north slope.


Introduction
Glaciers are natural water reservoirs 1,2 and are of vital importance to hydropower generation and hydrological models. 3,4 They are important indicators of climate variation 5-7 because fluctuations in climate can cause changes in their areal extent, volume, mass balance, [8][9][10] and equilibrium line altitude (ELA). 11 Mass balance and ELA, as direct reflections of climate change, integrate the competing effects of snow accumulation and ablation processes. 12 Field-based mass balance measurements are time-consuming, and it is difficult to obtain adequate representation at the scale of a mountain chain. 13 In contrast, satellite imagery allows for spatially consistent temporal reconstruction of glacier changes over recent decades. 14,15 The World Glacier Monitoring Service provides information for glaciers globally through direct field measurements and satellite images combined with geographic information science techniques. 16 The Global Land Ice Measurements from Space project was established to acquire satellitespectral images of the world's glaciers and analyze them for changes in extent and other characteristics. 17 Some glacier studies [18][19][20] have suggested that the snowline altitude (SLA) is *Address all correspondence to Zhen Li, lizhen@radi.ac.cn equivalent to the minimum elevation of perennial snow cover on a glacier near the end of the ablation period, making it a good proxy for the ELAs of mid-latitude glaciers. Regional SLA estimation using satellite imagery, when measured at the end of the ablation period, can be approximated by the lower boundary of the continuous snow-covered area (SCA). [21][22][23] Remotely sensed snow products have been applied for snowline observations; for example, Tang et al. 24,25 monitored the spatiotemporal patterns of the SLAs over High Mountain Asia using cloudremoved Moderate Resolution Imaging Spectroradiometer (MODIS) fractional snow cover data with 500 m resolution. Various methods [26][27][28] have also been developed for determining SLAs by applying various spectral band ratios to higher-resolution imagery. For instance, Hu et al. 29 and Girona-Mata et al. 12 delineated the snowlines for each Landsat scene by identifying the border between snow-covered and snow-free areas, which was determined by applying thresholds on the normalized difference snow index (NDSI) 30 in European mountains and the Himalayan catchment, respectively. However, the same threshold might not perform well in all regions because of highly variable cloud conditions, and it is essential to find an automated algorithm with the threshold adjusted to the local circumstances for each glacier individually.
Moreover, SLA and its spatial and interannual variability are known to be linked with the climate setting with sensitivity varying by region, such as in the western Himalayas 31 and in tropical and subtropical regions. 32 Further knowledge of the SLA controlling mechanisms plays a key role in assessing the potential application of the SLA to climate reconstruction and hydrological processes, especially in mountainous areas. The Tibetan Plateau (TP), referred to as the Third Pole, is one of the most sensitive areas responding to regional and global climate changes. Variations in glacier SLA in the TP have, therefore, received increasing attention. It encompasses 36,102 glaciers with a total area of 49; 688 km 2 , 33 and the glacier types in the plateau are mainly divided into continental and maritime glaciers. 34 The location of maritime glaciers is limited to the southeastern TP, whereas continental glaciers are widely distributed in other regions of the TP. 35,36 However, little is known about the connection between SLAs and climatic regimes, especially in the eastern TP, and there is a lack of comparisons between different glacier types.
In this context, this study focuses on the following two objectives. The first is to develop an automated approach that implements the Otsu threshold method 37 on a near-infrared (NIR) band to derive the snowline at a 30-m resolution from the 22-year time series of historic Landsat scenes in the eastern TP. The second is to investigate the connection between SLAs and climate regimes in different regions and describe the climatic mechanisms governing SLA dynamics for different glacier types in the eastern TP. The remainder of this article is organized as follows. Section 2 introduces the study area and datasets. The automated algorithm for extracting SLA is demonstrated in Sec. 3. In Sec. 4, the SLA and its variability during the 22-year period are presented. Section 5 discusses the linkage between the SLA interannual variability on the north and south slopes of different glaciers and the local climate conditions. Finally, Sec. 6 summarizes this study.

Study Area
To explore the relationship between local climate change and glacier SLA variation, three study areas distributed across the eastern TP were selected. These study sites include the regions of Sepu Kangri (SK; central location: 30.90°N, 93.78°E), Bu'Gyai Kangri (BK; central location: 31.82°N, 94.70°E), and western Qiajajima (WQ; central location: 33.58°N, 95.02°E), along a latitudinal transect from south to north (Fig. 1). Glaciers in the SK region are considered to be the maritime type, 38 whereas those in the other two regions are the continental type. 39

Google Earth Engine and Landsat Imagery
The Google Earth Engine (GEE) cloud computing platform 40 not only provides a large number of Landsat images in raw digital number (DN) form and top-of-atmosphere and surface reflectance orthorectified formats but also offers interactive exploration to quickly determine the spatial coverage of the study area and the temporal extent of satellite image collections. 40 As a cloud-based platform for mapping applications and data analysis, GEE enables monitoring longterm SLA variability in mountainous areas at 30-m spatial resolution. Because Landsat images taken at the end of the ablation season in the eastern TP were greatly affected by clouds, a total of 21 Landsat scenes acquired from 1995 to 2016 were finally selected for this study (Table 1). In addition to the Landsat product, high-resolution Sentinel-2 imagery under cloudless conditions was adopted to evaluate the SLA results.

Meteorological Data
Four meteorological stations within this study area record data such as daily air temperature and precipitation. Through these station records obtained from the National Meteorological Information Centre (NMIC), 41 the effects of these climatic variables on SLA can be investigated. To ensure that the dates of satellite images and meteorological measurements were in good agreement during the study period, only the mean temperature and precipitation at the four meteorological stations at the end of the ablation period (July to August) from 1995 to 2016 were analyzed. Station locations are shown in Fig. 1, and their location details are given in Table 2. These data suggested that the mean temperature at the four stations increased during the 22-year period and that the precipitation increased at each station except for Jia Li.

DEM and Glacier Outlines
The DEM scenes were used to extract SLA ranges. The ASTER GDEM, created from ASTER imagery, was released in 2009. 42 It has good spatial consistency with Landsat, with a horizontal resolution of 30 m and a vertical accuracy of ∼20 m. 13,43 The Randolph Glacier Inventory (RGI) 17,44 is a global inventory of glacier outlines freely available from the website. 45 Here, 26 glaciers larger than 1 km 2 in the three regions were selected as the research objects. Not all glaciers in the regions have a proper name; exceptions include Zuxuehui, Zhonggeimanong, Poge, and Beijia. Therefore, for the purpose of minimizing confusion, these 26 glaciers were identified using the names Glacier-A to Glacier-Z (A-Z glaciers), with their attributes shown in Table 3.

Methods
Significant differences in spectral reflectance characteristics are the basis of image classification. Both snow and ice have high spectral reflectance in the visible bands, whereas in the NIR band, the reflectance of clean ice gradually decreases and is less than the reflectance of snow. 43,46 The Otsu thresholding method is reliable in that the two classes are automatically distinguished by the threshold that maximizes the thresholding function, 47 thus realizing the distinction between snow cover and clean ice surfaces in the NIR band. 48 Subsequently, the discontinuously distributed snow that may have been more influenced by shadowing and local effects 23 was removed, and the continuous SCA of each glacier was extracted as the primary snow region. The lowest boundary of the primary snow region was considered to be the position of the glacier snowline, and the SLA of each glacier was obtained using ASTER GDEM (Fig. 2).

GEE-Based Landsat Imagery Selection
All Landsat images from 1995 to 2016, derived from the Landsat TM, ETM+ and OLI satellite instruments, were filtered on GEE through the following three criteria: 13,23 (a) Spatiotemporal criterion: in the study area, Landsat scenes of July and August from 1995 to 2016 were selected as the seasonal SCA achieved its minimum value during this period. (b) Regional cloud cover: one of the Landsat image metadata attributes, "CLOUD COVER," represents the cloud cover of the total image rather than the study area. Thus, the cloud percentage within the study area (i.e., regional cloud amount) was calculated by the "simple cloud score" function and stored in the image metadata as a new attribute, "regional cloud cover." If CLOUD_COVER > 70% or regional cloud cover >40%, the image was discarded. (c) Minimum seasonal snow extent: based on the above criteria, the NDSI with a threshold of 0.4 30 was also utilized to calculate the seasonal snow coverage for each scene, and the percentage was stored in the image metadata as a new attribute, "seasonal snow extent." The image with the minimum seasonal snow extent was finally selected as the optimal imagery to delineate the snowline.
Filtering Landsat products through these three criteria yielded only 21 images for the study area from 1995 to 2016 (detailed information in Table 1  (a) Select Landsat images, clip them using glacier outlines, and determine NIR band as the optimal feature; (b) discriminate snow from clean ice using Otsu thresholding method; and (c) delineate primary snow region and extract snowlines after negative buffering.

Feature Selection before Implementing the Otsu Method
From previous studies, the calculation of band ratios such as NIR/shortwave infrared (SWIR), red/SWIR, or NDSI has proven successful in mapping glacier extent. 49,50 However, these features map all ice and snow-covered surfaces without separating snow from clean ice, and the threshold must be chosen manually depending on the scene characteristics in different regions. 50 De Angelis 48 took a histogram of Glacier Tyndall to separate ice and snow by applying the thresholding function to NIR because the spectral reflectance characteristics differ greatly between snow and clean ice within glacierized areas in the NIR. Rabatel et al. 13 compared different band combinations and band ratios to facilitate the identification of snowlines on satellite images and showed that TM4 (NIR) with a threshold of 120 was the most appropriate for SLA retrieved from satellite images. Considering that the Otsu thresholding method performs best on the optimal feature, we compared the histogram distributions of different band combinations to find this feature. As shown in Fig. 3, the double peaks of snow and clean ice in the NIR were indeed more obvious than those of other features. Therefore, the Otsu thresholding method was applied to the NIR band to automatically distinguish SCAs from clean ice [ Fig. 2(a)].

Otsu Thresholding Method
The Otsu thresholding method uses a gray histogram in a spectral band to divide the images into target and background with a binary image as output. The segmentation threshold is adaptively obtained by calculating the between-class variance (BCV). However, this method is sensitive to the size of the target classification and noise. If many nonglacierized surroundings are considered to be the background (e.g., vegetation, land, and water) in one image, the DN value distribution is discrete, which causes the within-class variance (WCV) of the background to be relatively large. In this case, the segmentation threshold will be closer to the DN value of the background; thus, the background is likely misclassified into the target. Hence, neighboring glaciers with overlapping boundaries in RGI were merged into one validity vector mask and uploaded to GEE for clipping Landsat images [ Fig. 2(a)] to reduce the complexity of the image and prevent misclassification caused by nonglacierized surroundings.
As displayed in Fig. 2(b), the Otsu thresholding method will automatically search using DN values corresponding to the maximum BCV as the segmentation threshold to produce snow/ice classification. The key steps to implementing this method in GEE are as follows: (a) creating a list to store all thresholds; (b) computing weight and average DN values of target and background in the NIR band, respectively; (c) according to Eq. (4), looping over thresholds to calculate BCV for each; and (d) obtaining the maximum value of BCV. The value in the threshold list corresponding to the maximum BCV is used as the segmentation threshold for the image to produce clean ice and snow classification. The entire process of the Otsu thresholding method was employed for each Landsat image individually: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 2 7 2 (1) E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 6 ; 7 0 1 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 6 7 9 where n b is the pixel count of the background, n t is the pixel count of the target, N is the overall pixel count in one image, w b is the background weight, w t is the target weight, μ b is the average background DN value, μ t is the average target DN value, μ is the overall average DN value, σ 2 B is BCV, σ 2 W is WCV, and σ 2 is the total variance.

Snowline Delineation
To delineate the snowline position after the Otsu thresholding method, the largest continuous SCA of each glacier region was calculated as the primary snow region [e.g., the area surrounded by the red line in Fig. 2(c)] with the boundary of this region representing the snowline position. Moreover, the topography is known to be extremely steep in these regions, and glaciers are often disconnected from steep upper faces, which affects the distribution of snow cover and the location of the snowline. Therefore, a negative edge buffer, 51,52 with a buffer size of −200 m, removed the steep part of the glacier and ensured the stability and reliability of the snowline. Finally, snowlines were mapped for each of the 26 glaciers with selected Landsat imagery during the study period overlain on the ASTER GDEM to determine the altitude ranges of the snowlines in the three regions. To investigate whether SLA extraction performs better with the negative buffer, the SLA distribution before and after buffering was compared (Fig. 4). Figure 2(c) shows the postprocessing in GEE after the Otsu thresholding method to detect the snowline positions by the primary snow region and negative buffer.

Results
Due to the lack of any field-based measurements, we extracted SLA through manual interpretation and the Otsu thresholding method from the 10-m Sentinel-2 image (August 12, 2016) as the appropriate reference data to assess the accuracy of the SLA derived from Landsat (August 19, 2016). Although the images are not from the same date, the Sentinel-2 image, compared with Landsat, has a higher resolution available for snowline accuracy assessment. In Sentinel-2 imagery, clouds obscured partial glaciers in the BK region, and only three cloudless glaciers were available for reference in 2016 (Beijia, West Beijia, and Zhonggeimanong as a, b, c in Fig. 5). The standard deviation (SD) for the mean SLA presented here is the elevation range of individual glacier snowlines and is used to estimate the uncertainty error for each SLA. 13 For the three glaciers in the BK region, the uncertainty range of SLAs was between 35 and 59 m. As shown in Fig. 5, the position of the snowline obtained by applying the Otsu thresholding method to the Landsat scenes was close to that drawn by manual interpretation on the higher-resolution Sentinel-2 imagery, indicating that the Otsu thresholding method performs well in delineating the snowline of each glacier. Consequently, SLAs were calculated for 26 glaciers larger than 1 km 2 in the three regions: 8 glaciers in the SK region, 7 glaciers in the BK region, and 11 glaciers in the WQ region. The results (Table 3) include the mean SLAs and SDs monitored with Landsat images in the three regions and the rate of SLA change over time, which may be indicators of climate and annual climate variation.
In the period of 1995 to 2016, the mean SLAs of all glaciers in the SK region were on the rise [ Fig. 6(a)], with increases ranging from 1 to 219 m. Among them, Glacier-C had the largest glacier area, the highest elevation, and a temporal change in SLA (ΔSLA) of 219 m, the largest change in the study period. Glacier-A had the smallest ΔSLA (1 m). The highest mean SLA of all glaciers was on Glacier-H (5747 AE 18 m) in the SK region, which had the second highest elevation among all glaciers. The lowest mean SLA was on Glacier-E (5396 AE 09 m), which had the lowest elevation.
Except for the descending trend of Glacier-I in the BK region from 1995 to 2016, the mean SLAs of other glaciers were on the rise [ Fig. 6(b)]. Glacier-M, which had a slightly smaller glacier area and higher elevation, had the smallest ΔSLA (31 m). Glacier-O, with the lowest elevation, had the largest ΔSLA (150 m). The maximum and minimum values of ΔSLA were found to appear in the north and northeast of the BK region, but ΔSLA variations were different, which may be related to the elevation of different glaciers. The highest mean SLA among all glaciers in the BK region was on Glacier-I (5733 AE 34 m), and the lowest was on Glacier-O (5389 AE 23 m).
Within the WQ region, the mean SLAs of different glaciers varied from 1995 to 2016 [ Fig. 6(c)]. Glacier-R, with the smallest area, had the smallest ΔSLA (3 m). In contrast, Glacier-W, with the largest ΔSLA, descended 119 m. The highest mean SLA was on Glacier-W (5484 AE 39 m), and the lowest mean SLA in the WQ region was on Glacier-X (5389 AE 16 m). Although three glaciers (e.g., Glacier T, W, X) had negative trends, glaciers with increasing SLA (>30 m) were also found in the WQ region.

Discussion
Considering all glaciers in each region to be one group, the mean SLAs of the three regions in the eastern TP from 1995 to 2016 were 5571 m (SK), 5580 m (BK), and 5435 m (WQ) from south to north, respectively. The interannual variability in SLA exhibited increasing trends, with increases of 94 m 22 yr −1 , 55 m 22 yr −1 , and 49 m 22 yr −1 from south to north. Clearly, maritime glaciers in the SK region were more sensitive to the local monsoon climate, with the most pronounced variation in SLA (R 2 ¼ 0.56). Glaciers in the WQ region were continental glaciers and had small variations in SLA. The BK region was located in the transition zone between continental glaciers and maritime glaciers, and the SLA changed 55 m during the period studied. The relationship between the interannual SLA variability and regional climate changes in the eastern TP is worthy of further exploration (note that significance tests were conducted at the 95% confidence level). The orientations of glaciers in the eastern TP, divided into south and north slopes (the south slope is the sunny slope; the north slope is the shady slope), may influence SLA due to temperature and precipitation differences. These climatic factors were derived from data collected at four meteorological stations. Among them, the mean temperatures of the three regions during the end of the ablation period were 10.4°C, 11.7°C, and 10.5°C, and the mean precipitation totals were 125.2, 107.8, and 95.1 mm from south to north.

Effect of Climatic Factors on SLA in the SK Region
The SK region is located on the east side of Mount Nyenchen Tanglha and is obviously affected by the southwest monsoons from the Indian Ocean in summer. 53 Glaciers in this area are typically maritime in that they receive a large amount of water vapour (e.g., precipitation 125.2 mm). For the entire SK region, the precipitation was statistically significant with respect to the SLA variation [ Fig. 7(b); R 2 ¼ 0.66, P < 0.05], and the temperature very likely also had an effect on the SLA changes because of the relatively large correlation coefficient [ Fig. 7(a); For the north and south slopes, the interannual variation in SLA was related to different climatic factors. Specifically, the SLA variation on the north slope may be statistically significant in terms of temperature [ Fig. 7(c); R 2 ¼ 0.70, P < 0.05]. Both temperature and precipitation on the north slope increased in 2006, but the increase in temperature suppressed the influence of the increase in precipitation on the SLA variation, resulting in a rapid rise in SLA. In contrast, precipitation had a greater impact on the SLA of the south slope [ Fig. 7(d); R 2 ¼ 0.56, P < 0.10]. The SLA rose when the temperature decreased in 2013, which was attributable to the decrease in precipitation.
Maritime glaciers were distributed in areas with abundant precipitation and humid climates, where SLA variability was more significantly affected by precipitation. Although the highest elevation of the region ranged from 6060 to 6893 m, the mean SLA in the SK region was 5571 m. Furthermore, the primary climatic factors affecting SLA changes on the north and south slopes were different, and the mean SLA on the south slope (5524 m) was significantly lower than that on the north slope (5628 m).

Effect of Climatic Factors on SLA in the BK Region
The BK region lies in the transitional zone from maritime to continental glaciers, which are dominated by southeast and east monsoons during the precipitation season. 39 In the analysis of the transitional region, the SLA variations during the study period were much more affected by temperature [ Fig. 8(a); R 2 ¼ 0.57, P < 0.05] than precipitation. Unlike maritime glaciers in the SK region, the glacier SLAs on the north and south slopes responded to the climatic factors consistently. When the temperature and precipitation both increased in 1998, the SLA of the north slope rose [ Fig. 8(c)], revealing that the temperature had a more significant influence on the SLA than precipitation. Likewise, a tighter linkage between SLA variations and temperature was also found on the southern slope [ Fig. 8(d); R 2 ¼ 0.58, P < 0.05]. Specifically, the mean temperature on the south slope (11.9°C) was higher than that on the north slope (11.4°C), and the mean SLA (5685 m) was also higher than that on the north slope (5496 m). As an independent plateau with a maximum elevation of 6261 m, the mean SLA of the entire BK region reached 5580 m, which was the highest altitude of the snowline over the three regions in the eastern TP.

Effect of Climatic Factors on SLA in the WQ Region
The southern and southwest monsoons prevail in the southern TP at 30°N in summer and gradually weaken at 30°N to 35°N, whereas westerlies prevail north of 35°N. 54 The northernmost glaciers in the WQ region, located to the southeast of Qinghai Province, were spatially affected by westerlies, and there was a certain correlation between the SLA variations and temperature in this area [ Fig. 9(a); R 2 ¼ 0.41, P < 0.10]. Similarly, the SLAs on the south slopes were more sensitive to temperature [ Fig. 9(d); R 2 ¼ 0.45, P < 0.10] than to precipitation. From 2011 to 2016, the temperature increased annually, whereas the precipitation decreased year by year, which made the position of the snowline rise continuously to reach the maximum during this period. Although the SLA on the north slope [ Fig. 9(c)] did not pass the significance tests for temperature and precipitation, SLAs exhibited a clear change tendency under the integrated influence of these two climatic factors. For example, the decrease in temperature and increase in precipitation in 2011 caused the SLA to reach a minimum value during the study period, whereas the opposite meteorological variations led to an upward turn of SLA in 2016.
Obstructed by the extreme peaks of the Himalayas and Mount Nyenchen Tanglha, the precipitation in the WQ region was much less than those in the other two areas, and the interannual SLA variability of the continental glaciers was relatively affected by temperature. The difference in mean SLAs was slight, with 5433 m on the north slope and 5438 m on the south slope.
In addition, the spatial and temporal means of SLAs in the eastern TP from 1995 to 2016 were distributed irregularly in the north and south. The lowest and highest values were interlaced in places that may be affected by the southern terrains. Consequently, the SLA and its variability may be attributed to various factors, such as topographic differences between the regions, apart from the climate regime.

Conclusions
Based on the GEE cloud-based platform, this study develops a simple but effective algorithm that adaptively separates snow from clean ice by implementing the Otsu thresholding method on the NIR band of cloudless Landsat images to delineate glacier SLAs at the end of the ablation period from 1995 to 2016. We investigate the impacts of climatic factors on the SLAs of different types of glaciers (both maritime and continental glaciers) and their variations over the period studied. The main findings are summarized as follows:  Maritime glaciers in the SK region are more sensitive to the local monsoon climate, with the most pronounced variation in SLA (R 2 ¼ 0.56). (2) The northern continental glaciers (BK and WQ regions) are mainly affected by temperature and the southern maritime glaciers (SK region) are influenced by precipitation. Under the impact of primary climatic factors on SLA variations, continental glaciers have higher SLA on the south slope than on the north slope, whereas maritime glaciers have higher SLA on the north slope than on the south slope. The local patterns of temperature and precipitation, which influence glaciers and their snowlines, cause not only intraregional variability but also variations between regions. (3) The GEE platform and the approach proposed in this study have great potential to be developed further into an invaluable tool for obtaining more available satellite images to monitor SLA variability at 30 m resolution over longer periods or larger regions.