Changes in the End-of-Summer Snow Line Altitude of Summer-Accumulation-Type Glaciers in the Eastern Tien Shan Mountains from 1994 to 2016

Abstract

For summer-accumulation-type glaciers, the glaciological literature is lacking studies on determining the snow line altitude (SLA) from optical images at the end of the summer as an indicator of the equilibrium line altitude (ELA). This paper presents a workflow for extracting the SLA from Landsat images based on the variation in the albedo with the altitude in the central line area of glaciers. The correlation of >0.8 at the 99% confidence level between the retrieved SLAs with ELAs derived from the interpolation of ground-based, mass balance measurements indicated that the workflow can be applied to derive the SLA from end-of-summer satellite data as an indicator of ELA. The ELA was under-estimated by the calculated SLA. The relationship between the end-of-summer SLA and the ELA depends on the intensity of glacier melting. Subsequently, the workflow was applied to the seven glaciers in the Eastern Tien Shan Mountains, and a time series of the SLA was obtained using 12 end-of-summer Landsat scenes from 1994 to 2016. Over the whole study period, a mean SLA of 4011.6 ± 20.7 m above sea level (a.s.l.) was derived for the seven investigated glaciers, and an increasing SLA was demonstrated. The increase in SLAs was consistent for the seven glaciers from 1994 to 2016. Concerning the spatial variability, the east–west difference was prominent, and these differences exhibited a decreasing trend. The average SLA of each glacier is more influenced by its morpho-topographic variables. The interannual variations in the average SLA are mainly driven by the increasing summer air temperature, and the high correlation with the cumulative summer solid precipitation reflects the characteristics of the summer-accumulation-type glaciers.

1. Introduction

Mountain glaciers are recognized as high-confidence climate indicators due to their sensitive and identifiable reaction to even minor climatic changes [1]. Glaciers have a real impact on local hazard situations, regional water cycles, and the global sea level [2]. To deepen our knowledge of glaciers’ response to climate fluctuations in mountain ranges up to a regional scale and to understand the contribution of glaciers to the water resources in high-mountain cryosphere watersheds, we need to explore climate–glacier relationships.

Glacier mass balance is a critical parameter that responds most directly to climate, which can be measured using glaciological or geodetic methods. Due to the commitment of human power and time, as well as the challenges posed by the harsh natural environment, glacier mass balance observations are limited worldwide (about 450 glaciers), which has restricted the representativeness for whole mountain systems [3]. The geodetic method usually determines variations in the glacier-wide volume for time scales of multiple years to decades [3], which cannot satisfy the need to know annual mass changes on a regional scale in geoscientific and socioeconomic applications [4].

Consequently, an alternative method is needed to indicate whether the glacier mass of a mountain system is growing or shrinking. The measurement of the annual snow line altitude (SLA) on many glaciers at the end of the summer is one of the most commonly used methods due to the annual end-of-summer SLA being an optimal indicator of the equilibrium line altitude (ELA). For a given glacier, the ELA has a stable relationship with the mass balance. These relationships have been confirmed for temperate maritime glaciers (e.g., [5,6,7,8]) and glaciers in the outer tropics [9]. Remote sensing has provided a unique opportunity to allow the spatially consistent temporal reconstruction of SLA at a regional scale since the 1970s [10]. Consequently, changes in the ELA and mass balance can be reconstructed using remote sensing images (e.g., [2,5,6,7]), and the climate–glacier relationship can be studied in remote areas without data from direct observations (e.g., [10,11,12]). For example, in the Western Alps, variations in the average SLAs are mainly driven by the summer temperature and winter precipitation [13]. However, for inland glaciers located at high and middle latitudes and with a summer precipitation maximum, using end-of-summer SLA as a proxy of the ELA is debateable. Therefore, the climate–glacier relationship is still poorly understood.

Although manually mapping the snowline is usually highly accurate, for repeat applications over large regions, this method is time-consuming. A method based on various band algorithms, combined with a selected threshold, has been widely applied in the study of the snow line over large regions [14,15,16,17]. In addition, spectral mixture analysis, the maximum likelihood classification of principal components, and unsupervised classification have also frequently been used for identification of the snow line (e.g., [18,19]). The most obvious disadvantage of these methods is that the same threshold does not apply to all glaciers because of the strong impact of terrain and glacier surface conditions (i.e., the presence of impurities or metamorphic snow) on spectral reflectance. Hence, it is vital to find a threshold that is adjusted to local conditions, i.e., for each glacier individually for every satellite. Guo et al. [20] manually determined the threshold for each snowline by analyzing the variation in albedo with altitude. Rastner et al. [4] developed a method for automatically selecting a threshold to discriminate glacier ice and snow based on a histogram of the reflectance frequency in the near-infrared waveband.

In the current study, we expanded upon the research of Rastner et al. [4] and Guo et al. [20] by analyzing the spatial variations in broadband albedo derived from satellite imagery. A threshold was automatically identified to separate snow from ice for individual glaciers in each satellite image. The workflow was validated on two glaciers in the Eastern Tien Shan Mountains by comparing ELAs derived from the interpolation of ground-based, mass balance measurements. Then, based on this workflow, we investigated the SLA time series for seven glaciers over 1994–2016 using 12 Landsat snapshots acquired in late summer, discussing the relationship between the measured ELAs determined using the glaciological method and the end-of-summer SLAs derived from optical images, and the impacts of climatic variables and morpho-topographic variables on SLAs. Our aims were (1) to examine the representativeness of the end-of-summer snowline from optical images as a proxy of the ELA for the summer-accumulation-type of glacier and (2) to deepen our knowledge of the relation between glaciers and climate for summer-accumulation-type glaciers.

2. Study Site

The Tien Shan Mountains lie in the hinterland of Eurasia. We focused on seven glaciers near the Tianger Peak, located in the Eastern Tien Shan Mountains (Figure 1). The glaciers were selected based on the following factors: (1) Available meteorological data (Daxigou station is the only weather station that can provide continuous ground-based observational records since 1990, so the glaciers selected for this study had to be located near the weather station); (2) the glacier area, as, according to the most recent Tien Shan Glacier Inventory, about 90% of the Tien Shan glacier is less than 1.5 km² in size [21], thus, the upper limit of the investigated glaciers had an area of 1.5 km²; (3) the glacier aspect: We selected glaciers with as different aspects as possible in this study; and (4) available images: During the study period, each of the study glaciers on every useable image were affected by clouds and scan line corrector (SLC) failure as little as possible. To obtain up-to-date glacier outlines, a Landsat-8 image acquired on 23 August 2016 was selected. The glacier outlines were manually delineated, and a 6.11 km² total area was calculated. The glacier outlines obtained in 2016 were used for all subsequent analyses. The seven glaciers are shown in detail in

Table 1. Overview of all seven study glaciers, including their area, elevation range, and mean aspect. UG1, Urumqi Glacier No. 1. The elevation range and mean aspect for each glacier were extracted from ASTER GDEM V2 with a 30 m spatial resolution (downloaded from http://earthexplorer.usgs.gov/, accessed on 1 March 2021). The aspect of the glacier matched the average aspect of the area of the glacier.

Glacier Label	Area (km2)	Elevation Range (m a.s.l.)	Aspect (°)
A/East branch of UG1	0.97	3770–4245	12.6
B/West branch of UG1	0.59	3860–4450	63.4
C	0.95	3700–4430	326.2
D	0.38	3790–4290	55.1
E	1.23	3650–4343	332.6
F	1.46	3695–4343	328.4
G	0.53	3785–4240	357.9

.

Urumqi Glacier No. 1 (UG1), consisting of two small valley glaciers, was included in the investigated glaciers. In this study, the east branch of UG1 was labeled A, and the west branch was labeled B. Among the seven investigated glaciers, continuous long-term mass balance measurements have been carried out at Urumqi Glacier No. 1 using the glaciological method since 1959 [22]. Mass balance data and related ELAs on A and B glaciers were used to evaluate the representativeness of the end-of-summer SLA computed from the satellite images as an indicator of the ELA computed from field measurements.

The region is characterized by a continental climate, which is cold and dry in winter and warm and humid in summer [23]. The prevailing westerlies are the primary atmospheric circulation regime, bringing moisture-laden air. This area is affected by other weather disturbances, such as the Siberian anticyclonic circulation and the cyclonic disturbances of the west wind circulation [24]. The annual air temperature averages −5 °C and the annual total precipitation is 460 mm at a nearby meteorological station (DXG station, Figure 1) located at 3549 m above sea level (a.s.l.). Over 78% of the annual precipitation falls as snow, hail, and sleet between May and August at higher elevations [25]. Therefore, unlike maritime-type glaciers and glaciers in transitional climates in polar, temperate, and tropical regions, for the glaciers in the Eastern Tien Shan, most of the accumulation and ablation primarily occurs from May to August. These glaciers are classified as summer-accumulation-type glaciers [26,27].

3. Data and Methods

3.1. Data

3.1.1. Optical Satellite Data

To obtain as high a spatial resolution as possible and enable a long-term comparison, while minimizing the uncertainties arising from different satellites and the spatial resolution, Landsat satellite missions were chosen. Specifically, Landsat-5 Thematic Mapper (TM), Landsat-7 Enhanced TM Plus (ETM+), and Landsat-8 Operational Land Imager (OLI) images were used as a basis to obtain snow lines with a spatial resolution of 30 m. Images were manually selected from late summer to show the higher SLA compared to the images acquired from other periods of summer. Unfortunately, in several years, the images were not available due to the underlying glacier being completely occluded by temporary snowfalls occurring in late summer or cloud cover. These cases mainly occurred in the 1990s, as well as in 2000, 2001, 2004, 2007, 2008, and 2013. Finally, 12 satellite images (path/row: 145/030) were acquired over the 1994–2016 period from USGS (https://earthexplorer.usgs.gov/, accessed on 1 March 2021) (

Table 2. Landsat data were used to retrieve the snow line altitudes (SLAs). The mass balance measurement data were employed to calculate the equilibrium line altitudes (ELAs).

Landsat Date YYYYMMDD	Landsat Sensor	Mass-Balance Measurements Date	Landsat Date YYYYMMDD	Landsat Sensor	Mass-Balance Measurements Date
19940824	TM5	19940831	20100804	TM5	20100911
20020814	TM5	20020901	20110807	TM5	20110906
20030817	TM5	20030830	20120902	ETM+	20120901
20050814	ETM+	20050901	20140831	OLI	20140901
20060801	ETM+	20060901	20150826	ETM+	20150902
20090724	ETM+	20090917	20160804	OLI	20160902

).

The major drawbacks of Landsat data are the scan line corrector (SLC) failure in the ETM+ data post May 2003, and the saturation problem over snow-covered areas in the TM and ETM+ data, resulting in missing data. Compared with the image in 2016, the percentage of areas with missing data was less than 5% for individual glaciers. This value is relatively small, and the areas with missing data were neglected in this study. For Landsat TM/ETM+, the saturation problem exists over snow-covered areas, but it is not present in the Landsat OLI data. The retrieval of albedo values from the Landsat images is limited to the bands in 0.52–0.6 and 0.75–0.9 μm. Generally, the band in 0.52–0.6 μm can be easily saturated. When this occurs, we applied narrow-to-broadband conversion that only depends on the albedo in the 0.75–0.9 μm band (Equation (5)).

3.1.2. Digital Elevation Model

Digital elevation model (DEM) data were necessary in this study for (1) topographic correction of broadband albedo retrieval; (2) determining the altitude of the snow line; and (3) computing the morpho-topographic variables used—aspect and elevation—which were calculated as the arithmetic mean of the elevation and aspect of each pixel of the DEM included within the glacier’s outline, respectively. In this study, we used 30 m spatial resolution ASTER GDEM V2 (downloaded from http://earthexplorer.usgs.gov/, accessed on 1 March 2021) for the full time series. The vertical and horizontal accuracies were 17 and 26.5 m, respectively, in the study area at a 95% confidence level [28].

3.1.3. Mass Balance Data

Mass balance and related ELA and accumulation area values determined using the glaciological method were available over the 1994–2016 period for glaciers A and B. The field measurements were carried out twice a year around 1 May and 31 August, and the mass balance year is defined as the period from 1 September to 31 August in the following year. More than 40 ablation stakes were evenly distributed on UG1, covering different altitude zones, with intervals in altitude of about 50 m (Figure 1). The number and location of stakes varied slightly from year to year. The point mass balance was acquired by measuring the thickness and density of the annual layer in the accumulation zone and the distance between the top of the stake and the ice surface in the ablation zone. Then, the location of measured points was marked on the topographic map of a glacier with 50 m contour intervals. The 50 m contour interval was chosen to ensure that we had at least one stake in each interval along the longitudinal profile of the glacier. For individual altitude intervals, the mass balance was calculated as the average of all point mass balances within the interval. By multiplying each interval surface area with the corresponding value of the mass balance, the specific mass balance was obtained, and it is reported as a function of altitude, showing the vertical mass balance profile (VBP). The ELA is defined as the altitude where the VBP is zero. The glacier accumulation zone is above the equilibrium line. The ratio of the area of the accumulation zone to the area of the glacier is the accumulation–area ratio (AAR). In this study, the ELAs and AARs from 1994 to 2016 were provided for glaciers A and B by the world glacier monitoring service (WGMS) [29,30].

The ELA is a theoretical line that defines the altitude at which annual accumulation equals the ablation. Therefore, the accuracy of ELA estimates depends on field measurements and mass balance calculations. The measurement and calculation of the mass balance contain various sources of systematic and random errors that are difficult to quantify. These uncertainties include the measured accuracy of stake readings and thickness of the annual layer, snow/firn density measurements, the spatial distribution of the measured sites’ network, the extrapolation for unmeasured areas, internal accumulation, and calving, as well as any changes in the glacier area that are usually not subject to annual measurement [31,32]. For Urumqi Glacier No.1, the uncertainty in the glaciological mass balance range is between ±0.11 and ±0.13 m w.e. [33]. ELA uncertainty was calculated from the standard error of the linear regression of the specific mass balance with altitude [9].

3.1.4. Meteorological Data

To investigate the impact of climatic variables on the SLA interannual variability, meteorological data from 1994 to 2016 from the Daxigou (DXG) station were used (Figure 1). This is the nearest weather station with continuous ground-based observational records for the whole study period. These data were supplied by the Xinjiang Meteorological Information Center. The climatic variables used in the analysis were the mean daily air temperature and solid precipitation. Here, solid precipitation events were recognized according to the air temperature, and the following three sets of temperature thresholds were used in this study. Kang and Atsumu [34] first reported that the temperature threshold of solid precipitation was 2.8 °C on Urumqi Glacier No. 1 based on field observations from 1986 to 1989. In addition, in energy and mass balance simulation research on Urumqi Glacier No. 1, the threshold temperature of 1 °C is usually used to determine the amount of solid precipitation [35]. In this study, the length of the summer period was from 1 May to 31 August, and the length of the winter period was from 1 September to 30 April of the next year; thus, one year usually ranged from 1 September to 31 August of the next year. Notably, the air temperature and solid precipitation were not extrapolated from the automatic meteorological station (AWS) to the (mean) glacier elevation based upon the altitudinal lapse rate. Given the microclimate differences, it was unreasonable to apply a constant altitudinal lapse rate for all study glaciers. However, the altitudinal lapse rates for individual glaciers were difficult to quantify due to limited meteorological measurements of the study glaciers.

3.2. Methods

3.2.1. Broadband Albedo Retrieval

To improve the accuracy of snow–ice differentiation, the Landsat images used for SLA measurements were calibrated to broadband albedo as a method of standardization. The calibration procedures included geolocation, radiometric calibration, atmospheric correction, topographic correction, anisotropic correction, and narrow-to-broadband conversion.

For the process of deriving surface albedo and its parameter selection, we referred to Yue et al. [36] for Urumqi Glacier No. 1. Here, only the main retrieval steps are briefly described.

We performed a correction of reflectance for atmospheric interference using FLAASH in ENVI 5.3 software. For topographic correction, an improved topographic correction approach based on the C-factor correction was applied [37]. The improved expression is

$$L_H=(L_T-L_{\min })\times \frac{[\cos \beta -\cos {\alpha }_{\min }]}{[\cos \alpha -\cos {\alpha }_{\min }]}+L_{\min },$$

(1)

where L_T is the reflectance of the inclined surface, L_H is the reflectance of the horizontal surface, L_min is the minimum reflectance of the satellite image, α is the local solar incident angle, and β is the solar zenith angle.

Anisotropic correction converts the directional reflectance into spectral albedo using the bidirectional reflectance distribution function (BRDF). The anisotropic factor (f) was calculated for ice by Greuell and de Ruyter de Wildt [38] and for snow by Reijmer et al. [39], shown as Equations (2) and (3), respectively:

For ice:

$$f_{ice}=a_0+a_1\times \cos \theta +a_2\times {\theta }^2\times \cos \phi +a_3\times {\theta }^2\times {(\cos \phi )}^2,$$

(2)

For snow:

$$f_{snow}=b_0+b_1\times {\theta }^2+b_2\times {\theta }^2\times \cos \phi +b_3\times {\theta }^2\times {(\cos \phi )}^2,$$

(3)

where f is the anisotropic reflection factor, φ is the relative view azimuth angle, and θ is the satellite zenith angle with respect to the inclined surface variables. a_i and b_i are coefficients used in two previous studies. These were determined by performing a regression analysis on all individual samples using the criterion of the minimum residual standard deviation in f, and can be found in the studies of Greuell and de Ruyter de Wildt [38] and de Ruyter de Wildt [38].

Finally, we calculated the broadband albedo from the narrow albedo in bands 3 and 5 for OLI, and bands 2 and 4 for TM or ETM+, following the parameterization in Knap et al. [40]:

$$\alpha =0.726\times {\alpha }_{2or3}-0.322\times {\alpha }_{2or3}{}^2-0.051\times {\alpha }_{4or5}+0.581\times {\alpha }_{4or5}{}^2.$$

(4)

Satellite images in band 2 easily saturate because fresh snow has a high reflectivity at visible wavelengths. When this occurs, another formula is used:

$${\alpha }^{\prime }=0.782\times {\alpha }_{4or5}+0.148\times {\alpha }_{4or5}{}^2.$$

(5)

As mentioned, the BRDF parameterizations differ for ice and snow, but whether an area is ice or snow is not known beforehand. Here, we calculated broadband albedo with both BRDF parameterizations. Then, we used the mean value with both BRDF parameterizations as the final albedo value of the pixel.

Because the procedures take into account all essential processes that substantially affect the relationship between the reflectance properties of the glacier surface and the satellite signal, good agreement between the albedo measured in the field and the retrieval albedo from Landsat images was found for Morteratschgletscher, Switzerland [41]; Urumqi Glacier No. 1, China [36]; Forni Glacier, Italy [42]; Haut Glacier, Switzerland [43]; and five glaciers in western China [44]. Mean differences between Landsat-retrieved albedo values and albedo values derived from radiation measurements for the time the satellite passes over were all within 0.05 (

Table 3. The overview of uncertainties in the Landsat-retrieved albedo values at different glaciers.

Glacier Name.	Landsat Sensor	Field Measurement Method	Uncertainty
Morteratschgletscher Glacier	TM and ETM+	measured incoming and reflected shortwave radiation for a fixed point	0.03
Urumqi Glacier No. 1	ETM+ and OLI	Analytical Spectral Device (ASD) Spectroradiometer for 17 ground points	0.02
Forni Glacier	ETM+	albedo measured by CRN1 net radiometer for a fixed point	0.01
Haut Glacier	TM	Kipp and Zonen CM 7B albedometer for a fixed point	0.02–0.05
Rongbuk Glacier, Qiyi Glacier, Laohugou No.12 Glacier, Dongkemadi Glacier, and Palongzangbu No.4 Glacier	TM and ETM+	albedo measured by CRN1 net radiometer for a fixed point	0.048

). Here, an average uncertainty of 0.05 is regarded as the uncertainty for all investigated glaciers in this study.

Furthermore, Naegeli et al. [45] reliably classified snow and ice over a melting glacier surface based on the derived shortwave broadband albedo map. These results suggested that the values of albedo derived from Landsat images are reliable.

3.2.2. Snow-Cover Mapping for Individual Glaciers

Because of the strong optical contrast between bare-ice and snow-covered areas, the delineation of bare-ice surfaces vs. snow-covered surfaces is feasible based on surface albedo values if the threshold value for snow or ice is defined [43]. As a single threshold applied to the surface albedo map of all the study glaciers did not provide reliable results, a workflow of automatically selecting a threshold for individual glaciers is proposed (Figure 2).

First, the centerline and a 200 m wide centerline area were extracted. Here, the centerlines were determined by maximizing the distance to the glacier margins and respecting downhill flow, according to Le Bris and Paul [46]. Second, contours were extracted from the registered ASTER GDEM at intervals of 50 and 5 m. Third, the layer of contours of 50 m was overlaid on the derived albedo map for the central-line area, and the albedo in the center-line area of the glacier at intervals of 50 m was then obtained (Figure 3a). The change in albedo with the elevation at intervals of 50 m is shown in Figure 3b, where the error bar represents the standard deviation at each elevation interval. The sharp transition between ice and snow occurred over a short distance (e.g., [20,41,43,47]). Therefore, the point with the largest standard deviation should be given special attention as it indicates that the changes in the surface properties are most drastic in this altitude zone. The albedo for this altitude was considered to be the glacier- and image-specific albedo threshold discriminating snow and ice and is henceforth termed α_thre. Therefore, the glacier surface was classified according to the α_thre. Areas with α ≥ α_thre were classified as snow regions. The albedo snow line (α_s) threshold in each year is shown in

Table 4. Albedo threshold (α_thre) of the glacier ice region and snow-covered region, taking the east branch glacier of UG1 (A) as an example.

	1994	2002	2003	2005	2006	2009	2010	2011	2012	2014	2015	2016
A	0.58	0.30	0.59	0.48	0.41	0.39	0.30	0.21	0.17	0.48	0.33	0.39

, using the east branch glacier of UG1 (A) as an example.

Snow lines were delineated as the lower boundary of the near-continuous snow-cover area across the entire glacier width (all the results of the snow-cover area and SLAs in the supplementary document). The fragmentary snow patches were ignored, which were more attributed to local effects, such as drifting snow, avalanches, or hollows. However, the area of these snow patches is included in the snow-cover area in Section 4.1.2 and Section 4.2. The ratio of the area of the snow-cover area to the area of the glacier is the snow-cover ratio (SCR).

3.2.3. Snow Line Altitude Determination

Each delineated snowline was overlaid on the extracted layer of 5 m contours to identify the altitude range of the snowlines. The snow line was not necessarily parallel to the contours, and usually crossed several contours; especially near the glacier banks, shadows from surrounding slopes resulted in a lower snowline altitude, such as on the eastern side of glacier A. Therefore, the mean values are provided for the annual SLA of each glacier. Values discussed throughout the text are the average SLA along its delineation for individual glaciers.

3.3. Uncertainty

The SCRs are subject to uncertainties arising from (1) the application of a static glacier outline, (2) the missing data resulting from the SLC-off in ETM+ data, and (3) the snow–ice mixed pixels. The area decreased by 0.28 km² in 1994–2016 for the whole Urumqi Glacier No. 1, representing about 15.2% of the area in 1994 [48]. When the glacier was divided by 50 m intervals, the area changes for different altitude belts showed that decreases mainly occurred below 4085 m a.s.l. for Glacier B and below 3935 m a.s.l. for Glacier A, corresponding to the middle and lower parts of the ablation area [48]. In the current study, the snow-covered area was almost all located above 4000 m a.s.l. for Glacier A and 4100 m a.s.l. for Glacier B. The error resulting from the variation in the area of the upper part of the glacier was less prominent. However, the shrinkage of the area in the middle and lower parts of the glacier caused an overestimation of SCRs. We assumed that the area of the other five glaciers steadily reduced by 15.2% during 1994–2016, resulting in an overestimation of less than 10%. For the uncertainty (2), as described in Section 3.1.1, we assumed that the area of missing data was all snow cover; SLC failure resulted in an underestimation of the snow area by a mean of 3%, ranging from 1.6% to 5%. For the uncertainty (3), in this study, Landsat images had a spatial resolution of 30 m. During the ablation season, even within a pixel, it is probable that snow and ice exist together, especially for pixels near the snow line. As this uncertainty is hard to quantify, no exact number is given here.

Uncertainty in the derived SLA mainly stemmed from (1) the vertical error of the GDEM, (2) the vertical error resulting from the horizontal uncertainty in the location of the snowline, and (3) the time gaps in the acquisition of Landsat images and GDEM. The first two uncertainties can be quantified. The error due to (1) was about 17 m (P < 0.05) in this study area according to Tachikawa et al. [28]. The uncertainty due to (2) was determined by multiplying the pixel size by the slope of the glacier in the vicinity of the SLA. However, for (3), the uncertainty was hard to quantify, despite being quantifiable, as it is more ambiguous for each study glacier and year. Hence, the impact of (3) was not calculated in the present study, but the impact of (3) is discussed in Section 5.4. The resulting total uncertainty of each SLA was added using error propagation (the root of the quadratic sum of the different independent errors), and varied between 17.7 and 24.4 m, depending on the year and glacier (

Table 5. Mean and uncertainty (SD) for SLAs on each glacier in the Eastern Tien Shan Mountains extracted from a 2000 ASTER GDEM. The mean SLA, displayed at the bottom right, is the individual glacier mean for each year. The uncertainty of the mean SLA was calculated as an average of population statistics. Labels (A, B, C, D, E, F, G) correspond to glaciers in Figure 1.

	A		B		C		D
Year	Mean	Std	Mean	Std	Mean	Std	Mean	Std
1994	3950	18.5	4067.5	18.6	3957.5	20.6	3940	19.9
2002	4097.5	21.2	4127.5	19.4	4110	21.5	4110	24.0
2003	3947.5	18.9	4072.5	18.4	3922.5	19.4	3925	19.1
2005	4037.5	20.1	4097.5	20.3	4065	21.6	3962.5	20.0
2006	4017.5	20.2	4087.5	18.7	4007.5	21.0	3965	21.0
2009	3902.5	18.6	3992.5	19.7	3917.5	19.4	3910	19.4
2010	4110	21.1	4135	19.6	4055	20.9	3997.5	23.1
2011	4095	21.6	4290	21.1	4192.5	22.1	4092.5	24.1
2012	4150	23.6	4230	23.4	4180	22.1	4077.5	23.3
2014	3960	19.5	4105	18.7	4002.5	21.1	3970	20.6
2015	4080	22.0	4122.5	19.5	4060	21.5	3990	22.1
2016	4120	22.3	4142.5	20.2	4075	22.2	3997.5	21.3
Mean	4039	20.6	4122.5	19.8	4045	21.1	3995	21.5
	E		F		G		Mean SLA
Year	Mean	Std	Mean	Std	Mean	Std	Mean	Std
1994	3827.5	19.8	3877.5	18.3	3872.5	20.2	3927.5	19.4
2002	3950	21.1	3962.5	17.8	NA	NA	4059.6	20.8
2003	3815	21.4	3855	18.1	3890	20.0	3918.2	19.3
2005	3940	21.2	3912.5	17.7	4010	22.4	4003.6	20.5
2006	3857	20.1	3907.5	17.7	4020	23.7	3980.3	20.3
2009	3800	20.7	3832.5	18.0	3837.5	18.4	3884.6	19.2
2010	4022.5	22.6	3965	18.8	NA	NA	4047.5	21.0
2011	3977.5	22.0	3982.5	21.4	4172.5	24.0	4114.6	22.3
2012	4232.5	24.4	3992.5	19.1	4202.5	23.4	4152.1	22.8
2014	3888	20.2	3965	17.7	3927.5	24.1	3974	20.2
2015	3923	21.9	3985	19.7	4022.5	20.2	4026.1	21.0
2016	4015	21.7	3985	19.9	4022.5	22.9	4051.1	21.5
Mean	3937	21.4	3935	18.7	3998	21.9	4011.6	20.7

).

4. Results

4.1. Comparison of Measured Values and Retrieved Values

4.1.1. Snow-Cover Ratio vs. Accumulation Area Ratio

The snow-cover ratios (SCRs) for glaciers A and B are shown along with the measured accumulation area ratio (AAR) in Figure 4. The root-mean-square error between the retrieved SCRs and the measured AARs was about 18% for glacier A and about 15% for glacier B, with determination coefficients (r²) of 0.71 and 0.53 for glaciers A and B at the 99% confidence level, respectively. Throughout the study period, similar multi-year trends in the two variables were demonstrated, with SCR in the majority of cases being slightly higher than the AAR for the two glaciers. The average AAR values for glaciers A and B were 33.6% and 45.2%, respectively, whereas they were 42.3% and 49.6% for the SCR, respectively.

4.1.2. Snow Line Altitude vs. Equilibrium Line Altitude

A comparison of the derived SLAs and measured ELAs for the two test glaciers is shown in Figure 5. In most cases, the SLA values slightly underestimated the ELA for both glaciers. For glacier A, the mean SLA for all years was about 48.7 m lower, and for glacier B, it was 10.1 m lower. Nevertheless, a significant correlation was observed at the 99% confidence level (r = 0.888 for glacier A and r = 0.808 for glacier B). The root-mean-square error between the retrieved SLAs and the measured ELAs for the whole study period was about 61 m for glacier A and about 43 m for glacier B.

The results of the comparison presented above suggested that the end-of-summer SCRs derived from satellite images overestimate the AARs; correspondingly, the SLA computed from snow-cover mapping underestimates the ELA. Therefore, the retrieved end-of-summer SLA (or SCRs) did not numerically approximate the measured ELA (or AARs) for the individual glaciers in specific years. Consequently, the mass balance reconstructed on the basis of the retrieved SLA may be unreliable. However, the significant correlation indicates that the snow-covered area on the glacier and the SLA derived from satellite images at the end of the summer provide a reliable indicator of the accumulation area and the ELA.

4.2. Changes in the End-of-Summer Snow-Cover Ratio

The SCR at the end of summer is a high-confidence indicator for the mass balance due to a high correlation of the glacier mass balance with the snow-cover ratio on a glacier. The data for glacier G may be unreliable due to only a few snow patches being identified in 2002 and 2010. The six (out of seven) glaciers with the most complete time series were used to analyze the temporal trend in snow-cover ratio change throughout the Tianger Peak from 1994 to 2016.

The annual SCR for individual glaciers is shown in Figure 6a. Overall, the SCR exhibited large variability, ranging from 13.2% in 2016 for glacier D to about 87.0% in 2011 for glacier F over the entire study period. All glaciers recorded their maximum SCR in 2009; however, their minima occurred at different times. Three out of the six glaciers had their minimum in 2011: Glacier A reached its minimum in 2012 and the remaining two glaciers recorded their minima in 2016. During the whole study period, the mean SCR for glacier E was the highest, followed by C, D, F, B, and A. In general, the SCR for individual glaciers decreased from 1994 to 2016, but this decreasing tendency was only statistically significant for glaciers E and F at the 95% confidence level. Five out of the six glaciers recorded the largest decrease in SCR between 2009 and 2011, with the SCR decreasing by 34% to 63%, and glacier E experienced the second largest reduction after 2015–2016 over the same period. By 2016, the SCRs for all glaciers were less than 50%.

Spatially, the east–west differences in SCR were remarkable (Figure 7). For both western and eastern slope glaciers, the mean SCR decreased over the study period. Nevertheless, western-slope glaciers were almost always larger during this period than their eastern-slope counterparts. The mean western slope SCR was 8.9% larger than that of the eastern slope from 1994 to 2016. The maximum east–west difference in SCR was 18.5% in 2012, whereas the SCR was almost equivalent for the western and eastern slopes in 2014.

In addition to individual SCRs, we calculated the regional mean SCR for all glaciers from 1994 to 2016 (Figure 6b). For the seven studied glaciers, the mean SCR over the whole period was about 50.5%. The difference between extreme years was as high as 55.7%, with the lowest value appearing in 2012 (25.4%) and the maximum in 2009 (81.2%). In general, a reduction in SCR was observed from 1994 to 2016. An analysis of the annual variation revealed four declines in the regional mean SCR. The largest reduction in the mean SCR from 2009 to 2012 was a difference of 55.7%, and the second largest reduction occurred during the period of 1994–2002, with the mean SCR decreasing by 29.5%. This value again decreased to 23% from 2014 to 2016, which is less than half of that in the period from 2009 to 2012. The smallest of the four declines occurred in 2003–2005, with a value of about 20%. However, the average SCR temporarily increased by 32% from 2002 to 2003 and by 29% from 2012 to 2014.

4.3. Multi-Year Variations in the End-of-Summer Snow Line Altitude

Individual glaciers’ SLAs are shown in Figure 8a. Individual glaciers’ SLAs ranged from 3800 ± 20.7 m a.s.l. (glacier E, 2009) to 4290 ± 21.1 m a.s.l. (glacier B, 2011). The seven studied glaciers showed increased SLAs from 1994 to 2106. Notably, the increasing trends for individual glaciers were not statistically significant at the 95% confidence level, except for glacier F (P = 0.033). The oscillations in SLA (ΔSLA) ranged from 160 m for glacier F to 365 m for glacier G. The SLA of glacier G (standard deviation: 115 m) and glacier E (standard deviation: 114 m) demonstrated the most drastic changes from 1994 to 2016, whereas glacier F (standard deviation: 53.7 m) demonstrated the least change among all glaciers from 1994 to 2016.

To investigate variations in the SLA for the entire range, the mean SLA of the seven glaciers was used to characterize the prevailing behavior of the range (Figure 8b). A mean SLA of 4011.6 ± 20.7 m a.s.l. was derived for all the glaciers and the whole study period, ranging from 3884.6 ± 19.2 m a.s.l. in 2009 to 4152.1 ± 22.8 m a.s.l. in 2012. The rise of the mean SLA was presented from 1994 to 2016. The largest increase in the average SLA was recorded from 2009 to 2012, which was one time greater than the second largest increase from 1994 to 2002. The other two smaller increases occurred from 2003 to 2005 and from 2014 to 2016, with increases of 85 and 77 m, respectively. In general, before 2010, the annual SLAs were lower than the average SLA of 4011.6 m a.s.l for the whole study period (except in 2002, when it was about 48 m above the average). Since 2010, the annual SLAs have consistently been above the mean SLA until the end of the study period (except in 2014, when it was about 37.6 m below the average). Although no significant tendency in SLA was observed before or after 2010, the average SLA increased by about 98.6 m from the first subset (1994–2009) to the second subset (2010–2016). Hence, the year 2010 can be considered a breakpoint in the time series. However, given that the time series was short for the two subsets, the results should be interpreted with caution. The long-term monitoring records show that Urumqi Glacier No. 1 has experienced steady negative mass changes since 1959, and the maximum mass loss was recorded in 2010 (–1441 mm w.e. in the east branch glacier and −1115 mm w.e. in the west branch glacier). In that year, the measured ELAs exceeded the summit of the east branch glacier (A) and the west branch glacier (B). Therefore, positive glacier mass balances may have caused the higher SLAs after 2010.

In addition, the spatial difference demonstrated that the mean SLAs from glaciers facing west were always lower than those facing east. The average SLA was 80 m lower on the western slope than the eastern slope, ranging from 117.5 m in 2003 to 17.5 m in 2012 (Figure 9). Although the mean SLA for both eastern and western slope glaciers increased from 1994 to 2016, the east–west differences in the mean SLA exhibited a decreasing trend. The east–west differences in the mean SLAs, as well as their temporal changes, are not easily explained. One possible reason is associated with morpho-topographic (e.g., the average altitude and area) and/or microclimate (e.g., the differences in temperature, precipitation, and humidity) differences of the glaciers.

4.4. Changes in Climatic Conditions

Figure 10 shows the change in summer temperature and summer precipitation over the 1994–2016 period for Daxigou station. The average summer temperature for the whole period was 3.8 °C. The interannual variability in the average summer temperature was high, with a standard deviation of 0.55 °C. The difference between extreme years was as high as 1.9 °C, with the lowest value measured in 2003 (2.8 °C) and the highest value measured in 2015 (4.7 °C). In addition, over the study period, the summer temperature time series showed an average increasing trend of 0.04 °C year⁻¹, assuming a linear trend, which was statistically significant considering a risk of error of 5%. Especially during the period from 1994 to 2002, when the SLA was lacking due to Landsat scenes being unavailable, the increasing trend in summer temperature was more significant, averaging 0.15 °C year⁻¹ considering a risk of error of 1%.

Regarding summer solid precipitation, when the threshold used to define solid precipitation was 2.8 °C, the mean annual summer solid precipitation was about 183 mm, ranging from 283 mm in 1996 to 79.6 mm in 2001. For the solid precipitation threshold of 1 °C, the mean annual summer solid precipitation decreased to 100.4 mm, ranging from 44.8 mm in 2011 to 183 mm in 2000. Regardless of the solid precipitation threshold used, no significant trend was observed over the study period.

5. Discussion

5.1. Relationship between the Equilibrium Line Altitude Calculated Using the Glaciology Method and the Snow Line Altitude Delineated from Remote Sensing Images

In the present study, the end-of-summer SLAs delineated from satellite images were, in the majority of cases, lower than the measured ELAs. A similar result was reported by Wu et al. [49] using the high-spatiotemporal-resolution HJ-1 satellite images for Urumqi Glacier No. 1 in 2009–2010. The underestimated ELA appeared in studies like those of Rabatel et al. [9] and Guo et al. [50] for the glaciers of the outer tropics and Qilian Mountains, respectively.

Our hypothesis is that the image-derived SLA at the end of the ablation season was the highest SLA: The possible errors were negligible. For cold glaciers where superimposed ice is accreted to the glacier, the relationship between the end-of-summer SLA from optical images and the ELA falls into one of four scenarios:

(1) Over the surface of the melting glacier at the end of summer, the present year’s accumulated snow can completely cover the superimposed ice and the firn region from the previous year. In other words, superimposed ice and firn are not exposed on the surface of the melting glacier at the end of summer. In this situation, the equilibrium line and image-derived snow line are essentially the same line, and the image-derived SLA theoretically equals the ELA. This situation probably occurs in years when the mass balance is greater than zero, e.g., 2008/2009;

(2) With the enhancement of ablation, a superimposed-ice zone is exposed at the surface of the melting glacier at the end of summer, while the firn region is still covered by the present year’s snow. In this situation, the lower boundary of the superimposed-ice zone is the equilibrium line, while the image-derived snow line is the boundary of the present year’s accumulated snow. Hence, in theory, the ELA is lower than the image-derived SLA. This situation can happen in both positive and negative mass balance conditions;

(3) In the case of more intense ablation, the firn accumulated in the previous year begins to melt. The image-derived snow line is actually the lower boundary of the firn region from the previous year, called the firn line, whereas the equilibrium line (which represents the snowline from the present year’s accumulated snow) is above the firn line. The ELA is higher than the image-derived snow line SLA. This situation is likely to occur in negative mass balance years;

(4) As a severe case of the third scenario, for example, on glacier A in 2010 and 2011, the measured ELA exceeded the summit of the glacier, while the image-derived SLA was 4110 and 4095 m a.s.l., respectively. In reality, snow cover from the present year has melted, and the image-derived SLA is still the firn line. Therefore, there is no theoretical relationship between SLA and ELA.

Considering this, the relationship between the end-of-summer SLA and the ELA is complex and depends on the intensity of glacier melting. Urumqi Glacier No. 1 has experienced steady mass loss since 1959, and an accelerating trend in mass loss has been observed. As observed by Li et al. [51], the superimposed-ice zone diminished, and part of the superimposed-ice zone even disappeared with increasing summer temperatures. As a consequence, for the first and second scenarios, the probability of occurrence is decreasing, while the probability of the third scenario happening is increasing, and the fourth scenario may even occur. With these considerations in mind, the image-derived SLA being, in most cases, lower than the measured ELA is understandable.

5.2. The Impact of Climate on Snow Line Altitude Changes

Snow line altitudes are sensitive to changes in climatic variables, i.e., air temperature, precipitation, and humidity. Generally, increasing (decreasing) solid precipitation results in a decrease (increase) in the snow line altitude, and an increasing (decreasing) temperature results in an increase (decrease) in the snow line altitude. For example, in 2009, the mean air temperature in summer was much lower than the multi-year average of 3.8 °C for 1994–2016, coupled with unusually high summer solid precipitation, thus resulting in the anomalously low SLA in 2009 (Figure 10). During the whole study period, the highest correlation coefficient was found between the mean summer air temperature and the annual mean SLA, reaching r = 0.828 with a risk of error lower than 1% (

Table 6. The correlation coefficient between the annual mean snow line altitude (SLA) for the study glaciers and the daily temperature (T) and the cumulative solid precipitation (S-P) from Daxigou weather station at different time periods. * represents the solid precipitation threshold of 2.8 °C, and # represents the solid precipitation threshold of 1.0 °C.

	Period	Annual			Summer			Winter
	Climatic Variables	T	S-P*	S-P#	T	S-P*	S-P#	T	S-P*	S-P#
SLA vs.	Correlation coefficient	0.373	−0.491	−0.318	0.829	−0.661	−0.587	−0.120	0.112	0.81
	Significance level (p)	0.232	0.105	0.155	0.001	0.019	0.045	0.710	0.728	0.713

). Regarding solid precipitation, regardless of the solid precipitation threshold used, the annual average SLA and the cumulative summer solid precipitation were negatively correlated, considering a risk of error of 5% (Table 6). This was different from a significantly negative correlation being found between the annual average SLAs and winter precipitation over the other mid-latitude mountain glaciers [8,13]. In the present study, the stronger correlation between the annual average SLA and the cumulative summer solid precipitation reflects the characteristics of the summer-accumulation-type glaciers.

In addition, the sensitivity of SLA to climatic variables was estimated by establishing a linear regression between SLA anomalies and mean summer temperature anomalies (Figure 11a) and cumulative summer solid precipitation anomalies (Figure 11b). Due to the correlation being higher when the threshold of solid precipitation is 2.8°C, the sensitivity of SLA to solid precipitation was only observed for the case of 2.8°C. Here, the anomalies were computed from the multi-year average values. The results suggested that the SLA increased by 104 m when the mean summer temperature increased by 1 °C, while the sensitivity of the SLA to summer solid precipitation was 97 m/100 mm. The result of the SLA sensitivity to summer temperature was higher than the value reported in the literature, which was 61.7 m °C⁻¹ for Urumqi Glacier No. 1 [52]. The difference mainly resulted from the following: (1) The time series of the data: The study period was 1959 to 2008 in the literature, while for our study, it was 1994 to 2016; (2) in the literature, researchers computed the sensitivity of ELA to summer air temperature, while in terms of the value, the ELA was not equivalent to SLA for this study, although SLA was an indicator of ELA; and (3) Dong et al. [52] only focused on Urumqi Glacier No. 1 (namely, glaciers A and B in our study), while we considered seven glaciers, including glaciers A and B. Other reported values of SLA sensitivity to summer temperature from the mid-latitude mountain glaciers are 172 [53], 100 [54], and 115 m °C⁻¹ [13]. These values are close to our estimate.

The above results indicate that the end-of-summer SLA is highly sensitive to changes in air temperature and solid precipitation, especially in the summer. This implies that the regional mean SLA derived from optical images at the end of the summer can be used as a regional climate proxy in remote regions.

5.3. The Impact of Morpho-Topographic Variables on the Snow Line Altitude

As shown in Section 4.3, although there was an adjacent location and similar climatic conditions for the seven glaciers in the Eastern Tien Shan Mountains, the differences in the average SLA of the glaciers studied and their changes were distinct in 1994–2016. This suggests that other factors may have been involved, such as morpho-topographic variables, including the altitude and surface area.

Scatter plots of the average SLA vs. the average altitude and area are presented in Figure 12a,b. In Figure 11, each point represents a glacier. Figure 12a shows that the SLA increased with an increase in the average altitude. On the contrary, an increase in the glacier area resulted in a decrease in SLA (Figure 12b). Glacier D had the smallest surface area, while its mean SLA was not the highest of the seven glaciers, which occurred for glacier B, and appears to be associated more with a higher mean altitude. As a result, the altitudinal effect was stronger than the surface area effect. Accordingly, the morpho-topographic effect may be the principal reason for the east–west differences in the mean SLA noted above. The bigger glacier surface area and lower mean altitude caused the lower mean SLA on the western side of the ridge.

5.4. Uncertainty Stemming from the Time Gaps in the Acquisition of Landsat Images and GDEM

Uncertainty is introduced by using the ASTER GDEM V2 from 2000 for the full study period. This is because the DEM only represents the glacier surface from 2000, while the surface elevation decreased continuously. Therefore, we evaluated the uncertainty due to the general surface lowering on Urumqi Glacier No. 1 (glaciers A and B in this study). The results of the field measurement suggested that the cumulative mass balance of Urumqi Glacier No. 1 was −11.6 m w.e. during 2000–2016 [29,30]. This implied that the thickness of Urumqi Glacier No. 1 reduced by an average of 11.6 m from 2000 to 2016, if the change in surface area and the effect of the glacier dynamic process are neglected. In addition, the ice thinning continuously decreased with an increasing altitude; between 3950 and 4200 m a.s.l., it was about 0.4–0.5 m year⁻¹ [55]. Most of the SLA was located between 3950 and 4200 in this work, which is an altitudinal range for which the error resulting from the variation in surface altitude was less prominent (up to +3 m before 2000 and as low as −8 m after 2000). Hence, the vertical error associated with the surface elevation changes was neglected in this study.

5.5. Advantages, Drawbacks, and Robustness of the Albedo Method Applied to Mountain Glaciers

We have presented a workflow to determine the SLA from end-of-summer satellite images based on variation in the albedo with the altitude in the central line area of the glacier (we temporarily call it the albedo method). Although Guo et al. [20] proposed and applied a similar method to the Western Himalayas, the threshold in their study was manually determined; thus, their method was more subjective and required more human resources and time, resulting in the method being unsuitable for extension to a larger scale, such as a whole mountain system. In this study, we tried to select the point with the largest standard deviation in the albedo change as the threshold to enable the automatic discrimination of glacier ice and snow. Raster et al. [4] presented a fully automated approach to determine snow cover on glaciers and derive the SLA from multi-temporal Landsat data and a DEM. Their approach works very well when snow and ice facies have a good contrast, resulting in a bimodal histogram, but the histograms were not always exactly bimodal, for example, Figure 3b, d in their study. In these cases, there was no clear minimum detectable in the histogram, resulting in a larger difference between the automatically- and manually-derived SLAs. Moreover, the drawback of this approach was that the reflectance values from Landsat images are limited by the top-of-atmosphere reflectance at near-infrared wavelengths. Less attention has been paid to the glacier surface features, including light-absorbing impurities and liquid meltwater, which substantially influence snow albedo in the spectral region, where reflectance by snow is the strongest (λ = 0.3–0.7 μm) [56]. In the current study, the threshold was automatically determined based on the variation in the derived shortwave broadband albedo with the altitude for individual glaciers, which considers the impacts of glacier surface features on the snow cover as much as possible. As a consequence, the albedo thresholds for the years 2011 and 2012 were lower than for other years, being about 0.21 and 0.17, respectively, which were within the range of 0.15–0.40 recommended by Cuffey and Paterson [57] for impurity-rich firn snow. The lower albedo thresholds implied that the glacier experienced strong ablation and firn snow was exposed over the glacier, resulting from the warmer summer temperatures and low solid summer precipitation.

Regarding the robustness of the albedo method applied to Urumqi Glacier No. 1, Zhao et al. [58] evaluated the performance and accuracy of using visual interpretation, the normalized difference snow index, a single-band threshold, and the albedo method to extract end-of-summer SLA in comparison with the annual ELA. Their results showed that the albedo method was more consistent with the annual ELA obtained from field measurements, except for visual interpretation.

The successful application of the derived SLA from end-of-summer satellite data as a proxy for the ELA for Urumqi Glacier No. 1 indicated the potential for this workflow to be extended to a larger scale within the whole mountain system. Nevertheless, the SLAs derived from the end-of-summer Landsat images underestimated the ELAs; thus, the mass balance deduced from the derived SLAs may be imprecise. Although precise values for the glaciers in terms of ELA data were not obtained using this method to derive quantitative mass balance values, the method has merit in resolving long-term trends. Especially in remote areas without data from direct observations, an attractive advantage of the albedo method is its lower subjectivity and fewer logistical requirements. As a consequence, the magnitude and variability of SLAs on glaciers distributed throughout the mountain system could provide new insights into spatial contrasts in the behavior of glaciers, as well as their relationship with local and global climates.

However, the SLA derived from Landsat images was the transient snow line, while the ELA is the theoretical mean line on a glacier where the annual mass accumulation equals the annual mass loss during one mass-balance year. In fact, optical remote-sensing derived the highest SLA as the indicator of the ELA. Hence, in previous studies, including this study, images were preferentially selected at the end of summer, mainly with the expectation of obtaining a higher SLA compared to the images acquired from other periods of summer. Temporary snowfall is likely to occur at the end of summer, and the glacier was fully snow covered. Clouds can hide glaciers at the time of image acquisition. Landsat data exhibit the long revisit time of 16 days. These factors mean that the number of useful images from the end of summer is often insufficient for capturing the highest SLA. In this regard, it can be expected that applying the albedo method to images with a shorter revisit time (e.g., Sentinel-2A/B) will improve the future data availability drastically. In addition, as reported by Wu et al. [49], optical remote-sensing images only distinguish snow (or firn) from ice, but cannot distinguish snowfall in the present year from that accumulated previously. This limitation will confuse the firn line from previous years with the snowline from the present year, while the equilibrium line should refer to the snowline from the present year’s accumulated snow. This will result in an underestimation of the ELA.

Although the snow cover and the related SLA retrieved by the albedo method agreed with the measured values in most cases, in some years (such as in 1994), larger differences arose between the remote-sensing-derived and measured results. These larger differences may be related to the threshold. Notably, the variability in albedo thresholds was high for individual glaciers, ranging from 0.17 to 0.58. In this regard, there is a need for synchronous and accurate albedo field measurements, especially near the glacial snow line, to verify the snow–ice threshold derived from remote sensing images using the albedo method. In addition, due to the missing data in the 1990s, as well as in 2000, 2001, 2004, 2007, 2008, and 2013, the trend in SLAs, as well as the correlativity with air temperature and solid precipitation, should be interpreted with caution. In order to more reliably study the relationship between glaciers and climate, it will be significant to obtain longer-term and complete SLA data.

6. Conclusions

This paper has presented a workflow for extracting the SLA from end-of-summer satellite images. The threshold between the snow-covered regions and glacier-ice regions was automatically selected based on the variation in the albedo with the altitude in the central line area of the glacier. Subsequently, a layer of 5 m interval contours was overlaid on the snow-cover map to determine the snow line altitude.

Comparing the SCRs and SLAs retrieved from the end-of-summer satellite images with the AARs and ELAs obtained from field measurements, agreement was observed at the 99% confidence level. In terms of the value, the end-of-summer SCRs derived from satellite images overestimated the AARs; correspondingly, the SLAs computed from snow-cover mapping underestimated the ELAs. The relationship between the end-of-summer SLAs and the ELAs is complex and depends on the intensity of glacier melting. In summary, the performance of the workflow is reliable for studying the SLA derived from end-of-summer satellite data as an indicator of ELA.

The SCRs and SLAs were calculated using the workflow above for seven glaciers located in the Eastern Tien Shan Mountains over 1994–2016. The temporal variations in the SLAs and SCRs were presented. In parallel, the impacts of morpho-topographic and climate factors on the SLAs were investigated using bivariate correlation analysis.

Over the 1994–2016 period, the average SCR for the seven glaciers was about 50.5%. In general, the mean SCR decreased from 1994 to 2016. During the same time period, a mean SLA of 4011.6 ± 20.7 m a.s.l. was derived for all of the studied glaciers, and a rise of the SLA was shown from 1994 to 2016. The temporal variability in the SLAs and SCRs was homogeneous among the seven glaciers. Concerning the spatial variability, an east–west difference was prominent, with more rapid variations on the western side of the ridge.

The spatial differences in the snowline position were mainly dominated by the morpho-topographic variables of glaciers, namely the average altitude and surface area. The temporal variations in the average SLA of the glaciers studied were attributed to the air temperature and solid precipitation, especially in summer. As expected, the highest correlation coefficient was found between the SLAs and the summer air temperature. However, a higher correlation between the annual average SLAs and the cumulative summer solid precipitation was present, which is different from winter-accumulation-type glaciers and reflects the characteristics of the summer-accumulation-type glaciers. For the seven glaciers in the Eastern Tien Shan Mountains for 1994–2016, the rising SLA was mainly driven by an increasing summer air temperature.