Landsat-Based Estimation of the Glacier Surface Temperature of Hailuogou Glacier, Southeastern Tibetan Plateau, Between 1990 and 2018

: Glacier surface temperature (GST) is influenced by both the energy flux from the atmosphere above and the thermal dynamics at the ice–water–debris interfaces. However, previous studies on GST are inadequate in time series research and mountain glacier surface temperature retrieval. We evaluate the GST variability at Hailuogou glacier, a temperate glacier located in Southeastern Tibetan Plateau, from 1990 to 2018. We utilized a modified mono-window algorithm to calculate the GST using the Landsat 8 thermal infrared sensor (TIRS) band 10 data and Landsat 5 thematic mapper (TM) band 6 data. Three essential parameters, including the emissivity of ice and snow, atmospheric transmittance, and effective mean atmospheric temperature, were employed in the GST algorithm. The remotely-sensed temperatures were compared with two other single-channel algorithms to validate GST algorithm’s accuracy. Results from different algorithms showed a good agreement, with a mean difference of about 0.6 ℃ . Our results showed that the GST of the Hailuogou glacier, both in the upper debris-free part and the lower debris-covered tongue, has experienced a slightly increasing trend at a rate of 0.054 ℃ a -1 during the past decades. Atmospheric warming, expanding debris cover in the lower part, and a darkening debris-free accumulation area are the main causes of the warming of the glacier surface.


Introduction
Land surface temperature (LST), which is called the skin temperature of the Earth's surface [1], is an indispensable parameter in the physical, chemical, and biological processes at regional and global scales [2,3]. LST has been widely used in many fields of earth science due to its essential role in the mass balance and energy exchange between the atmosphere and land [4][5][6][7]. Direct field measurements of LST are still limited in regional-and global-scale research; thus, LST retrieval from multisource remotely-sensed thermal infrared data is of great importance in climatological [8], environmental, ecological [9], and meteorological studies [10].
Glaciers are recognized as indicators of climate change [11], and great changes associated with global warming have been occurring in mountain glaciers [12]. Glacier surface temperature (GST), which is strongly correlated with the ablation rate, is extremely sensitive to thermal processes related to both climatic and ice surface conditions. GST is influenced by both the energy flux from the atmosphere above and the thermal dynamics at the ice-water-debris interfaces. GST has, therefore, been used for glacier mass balance and runoff modelling and debris thickness estimations [13,14], highlighting the long-term dynamics of climatic conditions and changes of glacier surface properties.
To date, a large number of algorithms that can be used in retrieving surface temperature have been developed, including the mono-window algorithms, split-window algorithms, and multichannel algorithms [15,16]. The first generalized split-window algorithm for retrieving LST from MODIS was developed by Wan et al. in 1996 [16]. Recently, Wang et al. [17] improved the monowindow algorithm of Qin et al. [18], for which three parameters were needed (ground emissivity, atmospheric transmittance, and effective mean atmospheric temperature), and applied the algorithm to Nanjing and its vicinity in Eastern China. Meanwhile, Cristóbal et al. [19] also improved the generalized single-channel method of Jiménez-Muñoz et al. [10] to retrieve LST from Landsat 8 thermal infrared sensor (TIRS) data. Only two input parameters (near-surface air temperature, as well as integrated atmospheric column water vapor) were needed for this method. A variety of remote sensing retrieval studies based on each sensor type have been carried out on LST [20,21], sea surface temperature (SST) [22,23], and lake surface temperature (LaST) [24], which have been applied in different regions.
As for glaciers, limited investigation has been conducted on GST due to the remote locations of glaciers and the difficulty of determining ice surface emissivity. Stroeve et al. [25] used the advanced very-high resolution radiometer (AVHRR) to retrieve the surface temperatures of the Greenland ice sheet from 1989 to 1993, with the spatio-temporal variability being discussed. An attempt was made to estimate the surface temperature of Gangotri glacier, India, from advanced space-born thermal emission and reflection radiometer (ASTER) and Landsat 5 thematic mapper (TM) data, for which the retrieved temperatures and the observed temperatures showed a good correlation [26]. Additionally, using remote-sensing-retrieved data that was calculated through the single-channel algorithm, Wu et al. [27] studied the surface temperature distribution across the Qiyi glacier, China, from July 2012 to September 2013. With the development of unmanned aerial vehicles (UAVs) [28], ground-based thermal infrared images could be utilized to obtain surface temperatures for alpine glaciers, because of the UAVs' high resolution [29]. According to the literature, previous studies on GST, however, are inadequate in terms of time series investigations and mountain glacier surface temperature retrieval, and most of the studies concerning GST have concentrated on the polar regions, employing coarse spatial resolutions. Furthermore, the numbers of images in these studies are so few that the temporal variation of the GST is hard to obtain.
To address this shortcoming, this study aims to employ the modified mono-window algorithm to retrieve the surface temperature from Landsat 5 TM band 6 data combined with the Landsat 8 TIRS band 10 data for the Hailuogou glacier, southeastern Tibetan Plateau, from 1990 to 2018. The results are compared with two additional single-channel algorithms to validate their accuracy. The spatial variation coupled with temporal variation (i.e., seasonal variation as well as interannual variation) of GST are analyzed and discussed using these results.

Study Area
The Hailuogou glacier (29°34′N, 102°00′E), located on the eastern slope of Mount Gongga, in the Southeastern Tibetan Plateau (Figure 1), is a typical monsoonal temperate glacier, characterized by its high ablation rate and rapid surface velocity [30]. The length of the glacier is 13.1 km, the area is about 25 km 2 , and the elevation ranges from 2901 to 7556 m a.s.l. [31]. The glacier has a giant icefall ( Figure 1) in the middle of the ablation zone between 4900 m and 3850 m, characterized by frequent avalanches, as well as a steep glacier surface [32]. Below the icefall, most of the glacier tongue is debris-covered. The superficial debris covers 9.6% of the entire glacier, based on the latest RGI debris cover inventory [33]. The debris thickness increases from very thin (<1 cm) in the upper ablation area to more extensive cover (>0.5 m) towards the glacier terminus [34]. The local climate is dominated by the strong summer monsoon season. The annual precipitation recorded at 3000 m a.s.l. is about 1960 mm, while the mean annual temperature is 4 ℃. The glacier equilibrium line altitude (ELA) is about 4880 m [35,36]. In the context of global warming, the Hailuogou glacier has experienced remarkable terminus retreat and ice thinning, combined with an increased ablation rate. The glacier tongue has retreated by more than 2 km since the 1820s; the retreat rate was 12.7 m a -1 during 1966-1989 and 27.4 m a -1 during 1998-2008 [31]. The flow velocities of the ice tongue increase with distance from the terminus, ranging from 41 m a -1 near the terminus to 205 m a -1 in the upper glacier tongue below the icefall [37]. It has been shown that an increasing trend of glacier mass loss has emerged since 1999 [38].

Satellite Images
Both Landsat 5 TM data and Landsat 8 TIRS data were obtained from the USGS website (https://glovis.usgs.gov/), with pixel resolutions of 120 m and 100 m for the infrared thermal bands (provided by the USGS/NASA at 30-m resolution); the advanced spaceborne thermal emission and reflection radiometer global digital elevation mode (ASTER GDEM) data were obtained in the same manner. A total of 99 clear-sky images spanning 28 complete hydrologic years (from 1990 to 2018) were used to retrieve GST (Table A1) in order to study its temporal variation during this period. Note that all images were selected seasonally (March-May for the spring, June-August for the summer, September-November for the autumn, and December-February for the winter) and based on cloudfree conditions over the glacier. One Sentinel-2 image (23 August 2018) was used to map the 2018 debris cover area, due to the Landsat 8 imagery in the summer of 2018 being relatively cloudy. At the same time, ground-based observation data (including near-surface air temperature, atmospheric water vapor, and groundwater vapor pressure data) were acquired from the Gongga Mountain

GST Algorithms
In this paper, three algorithms were used to retrieve GST from Landsat thermal images. We follow the improved mono-window algorithm developed by Qin et al. [18] and compared the results with the radiative transfer equation (RTE) [39] and Jiménez-Muñoz and Sobrino's single-channel method (JM_SC) [15]. In this section, we first present the general equation of the improved monowindow algorithm. More details follow regarding the calculation of at-sensor brightness temperature, determination of ground emissivity, estimation of effective mean atmospheric temperature, and estimation of atmospheric transmittance.

The Mono-Window Algorithm
This method is widely used in estimating land surface temperatures and is developed with data specific to only one infrared thermal band [18]. The equation is expressed as follows (i = 6 for Landsat 5 TM band 6, i = 10 for Landsat 8 TIRS band 10): where Ts is the GST in degrees Celsius; Ti is the brightness temperature of Landsat 5 TM band 6 or Landsat 8 TIRS band 10; Ta is the effective mean atmospheric temperature; ai and bi are the constants used to approximate the derivative of the Planck radiance function for TM band 6 and TIRS band 10, given in Table 1; Ci and Di are the internal parameters for the algorithm, given as follows: where εi is the ground emissivity and τi is atmospheric transmittance of Landsat 5 TM band 6 or Landsat 8 TIRS band 10. Three significant parameters (Ta, εi and τi ) will be presented in the following sections.

Calculation of the At-Sensor Brightness Temperature
The brightness temperature Ti plays an important role in the process of GST retrieval from Landsat infrared thermal band data via the mono-window algorithm (Equation (1)). Firstly, the digital numbers (DNs) will be transformed into thermal radiance, and then the DN value can be converted into the brightness temperature. Considering the differences between Landsat 5 and Landsat 8, the DN values for Landsat 5 TM band 6 and Landsat 8 TIRS band 10 can be converted into thermal radiance by Equation (4a) [40] and Equation (4b), respectively.
where L6 is the at-sensor thermal spectral radiance (m W cm -2 sr -1 μm -1 ); Lmin6 and Lmax6 are the minimum and maximum detected spectral radiance, respectively; Qdn is the grey level for the analyzed pixel of the TM image; and Qmax is the maximum DN value. The values are given in Table  2.
where the L10 is the at-sensor thermal spectral radiance (W m -2 sr -1 μm -1 ) of TIRS band 10; Qdn is the grey level for the analyzed pixel of the TIRS image; M10 is the band-specific multiplicative rescaling factor for band 10; and A10 is the band-specific additive rescaling factor for band 10. The values of M10 and A10 can be acquired from the metadata file of Landsat 8 image, the values for which are presented in Table 3.  Table 3. Parameters for computing at-sensor brightness temperature from Landsat 8 TIRS band 10 data.

M10
A10 K1 (W m -2 sr -1 μm -1 ) K2 (K) 0.0003342 0.1 774.89 1321.08 With the computation of the thermal spectral radiance Li, the at-sensor brightness temperature can be calculated by using the approximation of the Planck radiance function, which is expressed as follows: where Ti is the at-sensor brightness temperature (K); and K1 and K2 are the band-specific thermal conversion constants. For Landsat 5 TM band 6 data, K1 = 60.776 m W cm -2 sr -1 μm -1 and K2 = 1260.56 K, respectively ( Table 2). For Landsat 8 TIRS band 10 data, K1 = 774.89 W m -2 sr -1 μm -1 and K2 = 1321.08 K, respectively (Table 3), which can also be acquired from the metadata file of the Landsat 8 image.

Determination of Ground Emissivity
Emissivity is related to the temperature, wavelength, surface roughness, observation angle, and composition of the surface, and is a key parameter for the surface temperature estimation. Several methods have been explored to simulate the land surface emissivity [41,42]. Owing to the fact that the NDVI(Normalized Difference Vegetation Index)-based emissivity method (NBEM) is inapplicable to this paper, we aim to calculate the GST by applying the semi-empirical model for the angular-dependent emissivity spectra of snow and ice in the 8-14 ㎛ atmospheric window [43]. Hori et al. classified the transition from snow to ice categories: fine dendrite snow, medium granular snow, coarse-grained snow, sun crust, and bare ice. They assumed that the emissivity of snow and ice surfaces can be expressed by a linear combination of emissivities for the specular and blackbody surfaces, with weighting parameters characterizing the specularity of the bulk surface. The equations are expressed as follows and more details can be found in [27,43]: where εsnow(λ, θ) is the angular-dependent emissivity of snow; λ is the wavelength; θ is the angle of exitance measured from the normal value; εbb(λ) is the emissivity of the isotropic blackbody (εbb(λ) = 1); fsp is the weighting parameter; εsp_app(λ, θ) is the apparent emissivity of the specular fraction; εsp(λ, θ) is the emissivity of an ideal smooth ice surface simulated with the Fresnel reflectance ρsp(λ, θ).
We used a band ratio approach to eliminate the impact of topography and enhance the discrimination between different snow grains. Based on three types of band ratio images ((Normalized Difference Snow Index, NDSI), Red/(Short-wave Infrared 1, SWIR1), (Near Infrared, NIR)/(Short-wave Infrared 1, SWIR1)), we utilized the K-mean unsupervised classification method together with previous field observations to classify the images into six categories: fine dendrite snow (fsp = 0.22), medium granular snow (fsp = 0.29), coarse-grained snow (fsp = 0.41), sun crust (fsp = 0.53), bare ice (fsp = 0.95), and glacier debris. In this paper, the emissivity of glacier debris has a fixed value of 0.941 [44].

Estimation of Effective Mean Atmospheric Temperature
Ground-based observations of effective atmospheric temperature (Ta) are hard to measure at the satellite pass time. Reviewing the previous literature, there are several methods that can be used to estimate the Ta. Sobrino et al. [45] developed a method that relates to the atmospheric temperature and water vapor content; this method is applicable for the places where the atmospheric profile data can be easily obtained at the satellite pass time. However, atmospheric profile data were lacking in this study. Here, we employ the method of Qin et al. [18], for which the standard atmospheric files and local meteorological data could be applied to estimate Ta. According to Qin et al. [18], the linear relations for approximation of Ta from the near-surface air temperature (To) are given in Table 4. Table 4. The linear relations for approximation of effective atmospheric temperature (Ta) from the near-surface air temperature (To) in the mid-latitude zone [18].

Estimation of Atmospheric Transmittance
Atmospheric transmittance is another essential parameter for GST retrieval. The atmospheric transmittance is influenced by the ozone, aerosol, wavelength, atmospheric chemicals, and atmospheric water vapor content (w), indicating the magnitude of thermal radiance absorption in the process of travelling to the sensor. It is commonly agreed that atmospheric water vapor content governs the atmospheric transmittance in the thermal channel. Two approaches could be used to estimate the atmospheric transmittance. The first one is to apply the atmospheric simulation program moderate-resolution atmospheric transmission (MODTRAN4) [46], while the second is to establish a function between atmospheric water vapor content and atmospheric transmittance. Here, we use the latter approach to simulate the atmospheric transmittance. Owing to the differences between Landsat 5 TM band 6 data and Landsat 8 TIRS band 10 data, the relation between atmospheric water vapor content and atmospheric transmittance can be found in Table 5. In the next step, the atmospheric water vapor content (w), which is unavailable directly, is estimated according to Yang et al. [47], using the equation: where C0 and C1 are the empirical coefficients; W is the precipitable water. For the investigate area,

  
Lat d (15) in the equations above, e is the ground water vapor pressure (hPa); H is the elevation (km); Lat is the geographical latitude (°); a0 and a1 are the empirical coefficients calculated by Equations (11) and (12), respectively.
The near-surface air temperature and ground water vapor pressure recorded by the AWS can be obtained from the Gongga Mountain Alpine Ecosystem Observation Station, Chinese Academy of Sciences (CAS) (elevation: 3000 m a.s.l.).

Validation of GSTs
Previous studies suggested that the derived temperatures using the single-channel method exhibited a bias of 1 ℃, validated against in situ observation data [19]. By comparing three different methods for LST inversion from TIRS (the single-channel method, the split-window algorithm, and the radiative transfer equation), Yu et al. [39] found that temperatures derived from the radiative transfer equation were accurate, with a RMSE lower than 1 ℃. Due to the limited samples of the ground observation data across the Hailuogou glacier, derived GSTs using the mono-window algorithm above were compared with two other single-channel algorithms (Jiménez-Muñoz and Sobrino's single-channel method and the radiative transfer equation) [15,39] for the accuracy assessment. Here, we analyzed the accuracies across the entire study area over the images acquired in different seasons and from distinct sensor types (9 December 2017; 11 January 2013; 30 August 2006; 9 May 2000) ( Table 6). The max bias was 1.3 ℃, the average bias was 0.7 ℃, and the average root mean square error (RMSE) was 0.7 ℃ when compared with the JM_SC. The max bias of the differences between derived temperatures and RTE was 3.7 ℃, while the average bias and RMSE were 1.7 ℃ and 1.7 ℃, respectively. It is noteworthy that the RTE is not accurate enough to retrieve land surface temperatures [48]. A scatter plot of the data (Figure 2) showed that the derived GSTs were in good agreement with the two other algorithms, indicating that the method mentioned in this paper can be adopted to study the spatial-temporal variations of the GSTs.

Results
Here, we analyzed the spatial-temporal variations of Hailuogou glacier. For the spatial variation, we focused on the surface temperature gradient and the anomaly. Temporal variations (i.e., seasonal variation and interannual variation) were emphasized by certain metrics, such as the average surface temperature of the Hailuogou glacier, average temperatures in the debris-covered and debris-free areas, 0 ℃ isotherm elevation, and the extent of the glacier above 0 ℃.

Spatial Distribution Patterns of GSTs
For this purpose, we analyzed the GST spatial distribution patterns across Hailuogou glacier based on the image obtained by Landsat 8 on 10 November 2018. Looking at the entire glacier, the remote-sensing-retrieved temperature distribution shows that the surface temperature decreases with increasing elevation, while the mean temperature gradient is about −0.61 ℃/100 m, which is slightly lower than the air temperature (−0.65 ℃/100 m) (Figures 3a, 4 and Figure 5). As can be seen in Figure 4, the lowest temperature is not at the top of glacier, while the highest temperature is located at the north side of glacier tongue, owing to the influences of the local topography (Figure 3b). Accordingly, surface temperatures across the accumulation area in the Hailuogou glacier show strong heterogeneity at the same elevation level (Figures 3a and 4); part of the accumulation area would exceed 0 ℃ in summer due to the exposure of bedrock ( Figure 5). In the ablation area, the surface temperatures on the north marginal area of the glacier tongue are higher than that on the south (Figure 3b,c). Notably, most of the pixels in the glacier tongue yield temperatures larger than 0 ℃ at the satellite passes, and elevations at which the surface temperature is 0 ℃ (we call this "0 ℃ isotherm elevation" in this paper) range from 3500 m a.s.l. to 5000 m a.s.l. (Figures 4 and 5).

Variations of GSTs between 1990-2018
A total of 99 cloud-free Landsat images over the glacier were utilized to study the temporal variations of surface temperatures on Hailuogou glacier. GSTs in summer are higher than in winter; the box-chart shows that GSTs change in a seasonal way and the average temperatures show a clear seasonal fluctuation (Figures 6 and 7). A comparison between the station record air temperature and satellite-derived glacier surface temperature shows that the mean surface temperature on debriscovered areas is close to the station-recorded air temperature at 3000 m a.s.l. (Figures 7 and 8a). As shown in Figure 8b, average temperatures in debris-covered and entire glacier areas both showed an increasing trend between 1990 and 2018, however the increase in the debris-covered area is more significant (Figure 8b). The 0 ℃ isothermal line elevations and extent of the glacier surface above the 0 ℃ isothermal line show some complexes in their long-term variations. The median value of the 0 ℃ isothermal line shows some variations between 4200 and 4800 m a.s.l. (Figure 9a), while the area of the glacier surface above the 0 ℃ isothermal line is more dynamic, but with an overall increasing trend for the median level since 1995 (Figure 9b). The increased extent of the glacier above 0 ℃ indicates that the ablation area of Hailuogou glacier has increased since 1995.

Errors in GST Remote Sensing
In this section, we first discussed the sources of errors in retrieving GST from Landsat thermal bands, although the remote-sensing-retrieved temperatures were in good agreement with the two other models. We believe that three major aspects would have important impacts on obtaining accurate GSTs, as follows: (1) The retrieved results would be hampered by cloud cover, which is frequently in the areas of mountain glaciers, resulting in a large number of remote sensing images being unavailable in this study. (2) It is possible that there is a mismatch in the locations of the glacier surfaces and satellite pixels. Additionally, one pixel represents the temperature of a quadrilateral (30 × 30 m), which is different from the measurements taken in the field. (3) Because natural surfaces are heterogeneous in terms of variations in emissivity, it is hard to determine the ground emissivity at the pixel scale; therefore, the derived emissivity should be taken in an effort to validate its accuracy.

Impacts of Glacier Debris Cover Expansion and Ice Surface Darkening on GSTs
The supraglacial debris plays a significant role in ice melting by controlling the surface temperature of the glacier. A debris layer of several centimeters would accelerate the ice melting owing to the lower albedo of debris cover, causing the surface to absorb more radiation compared to clean ice [34], while the debris layer would also insulate the glacier when the layer exceeded that thickness [49]. In most cases, the termini of debris-covered glaciers remain nearly stagnant, especially in Himalayan regions [50]. In the accumulation area of the Hailuogou glacier, the surface temperature shows great heterogeneity at the same elevation level because part of the area is covered by a debris layer, which is thinner than 3 cm (Figure 10). Both the debris cover in the lower part and the summer darkened ice surface in the upper part of the Hailuogou glacier have shown remarkable expansion since 1990 (Figure 11). Increasing accumulation of black carbon and dust on the glacier surface reduces the glacier albedo, increasing the net shortwave energy and resulting in higher surface temperatures.

Effect of Topography on GSTs
The local topography impacts glacier surface temperatures by restricting illumination. There are three shading areas in the Hailuogou glacier-the first one is located near the top of the glacier, the second is at the glacier tongue, and the last is at the southern margin of the accumulation area ( Figure  12). In the ablation area of the Hailuogou glacier, the topography dominates the spatial distribution. The glacier tongue surface is affected by thermal emissions from the surrounding topography. The distances between the surrounding cliffs and the glacier itself, however, can influence the magnitude of this effect [27]. Due to the reception of thermal reflection emitted from adjacent slopes, the northern margin of the glacier tongue shows a higher temperature than in the center or at the southern margin of the tongue. The southern margin of the glacier tongue is partly shaded by the north-facing topography (see the second shading area), thus exhibiting a lower surface temperature in most seasons. In the accumulation area of the Hailuogou glacier, the lowest temperature is not at the top of glacier because of the topographic shading (see the first shading area). The last shading area is not an anomaly, owing to the expansion of the exposed bedrock, which stops the occurrence of extremely low temperatures in the daytime (see the third shading area). The subsolar point changes periodically, being located in the northern hemisphere in summer and in the southern hemisphere in winter. The impact of topography in the summer is less than that in winter ( Figure 4). Since the air temperature decreases with increasing elevation, the GSTs vary in a corresponding way with the elevation change. The steeper the slope (i.e., icefalls and ice cliffs), the higher the glacier surface temperature gradient. In general, the surface temperature on the south-facing slope is high, while that on the north-facing slope is low.

Climate as a Major Factor in the Temporal Variation of GSTs
In the context of global climate change, glaciers of the Tibetan Plateau have retreated over the past 50 years, and the retreat speed in the outer edge of this region is faster than that in the inner area [51,52]. Significant warming, especially in winter, is observed in the Tibetan Plateau based on in situ meteorological data from almost every weather station, whereby the annual mean air temperature increased by 0.16 ℃ per decade during the period 1955-1996 [53] and 0.35 ℃ per decade in the period 1961-2015 [54,55]. As shown in Figures 7 and 8, the average surface temperature of the Hailuogou glacier varied consistently with air temperature, showing obviously seasonal fluctuations and an increasing trend (Figures 6 and 8). As discussed above, dynamic climate change is a major factor in the temporal variation of glacier surface temperatures.

Conclusions
An attempt has been made in this paper to utilize the modified single-channel algorithm for GST retrieval from Landsat 8 TIRS data and Landsat 5 TM data. Three essential parameters (ground emissivity, effective mean atmospheric temperature, and atmospheric transmittance) are needed for GST retrieval with this method, and determinations of the three required parameters are presented with respective details in this study. GST has been used for the glacier mass balance and runoff modelling and debris thickness estimations, so it plays an important role in the evaluating glacier ablation rate and studying the glacier hydrology. In order to fill the gap in the study of long-term glacier surface temperature changes, the primary objective of this study, therefore, was to research the spatial-temporal variations of surface temperature across the Hailuogou glacier during 1990-2018. Although the remote-sensing-retrieved temperatures were not validated by the field measurements, this paper found that the method drives reasonable GSTs and can be applied to the Hailuogou glacier when compared to the two other algorithms.
The derived temperature distribution for the entire glacier shows that the surface temperature of the Hailuogou glacier generally decreases with increasing elevation, with a temperature gradient of -0.61 ℃/1000 m. Owing to the presence of debris cover and exposed bedrock, temperatures across the accumulation area of the Hailuogou glacier show strong heterogeneity. In the ablation area, the topography is the main factor in the surface temperature distribution. For the time series of GSTs, both the glacier-wide mean surface temperature and 0 ℃ isotherm elevation show an increasing trend from 1990 to 2018, with rates of 0.054 ℃ a -1 and 12.415 m a -1 , respectively. If these trends continue, the 0 ℃ isotherm elevation will rise to about 5100 m and a total area measuring 8 km 2 will experience melting year-round in 50 years, which will result in the disappearance of the glacier landscape, a shortage of water resources, as well as frequent regional hazards. In general, the results from this study indicate that the Hailuogou glacier has become warmer and the ablation area has increased since 1990 because of the rising temperatures around the world.

Appendix A
Details of the Landsat 5 TM data and Landsat 8 TIRS data that were used to derive GSTs in this paper are shown in Table A1.