Variations in Winter Surface Temperature of the Purog Kangri Ice Field, Qinghai–Tibetan Plateau, 2001–2018, Using MODIS Data

: In the context of global warming, the land surface temperature (LST) from remote sensing data is one of the most useful indicators to directly quantify the degree of climate warming in high-altitude mountainous areas where meteorological observations are sparse. Using the daily Moderate Resolution Imaging Spectroradiometer (MODIS) LST product (MOD11A1 V6) after eliminating pixels that might be contaminated by clouds, this paper analyzes temporal and spatial variations in the mean LST on the Purog Kangri ice ﬁeld, Qinghai–Tibetan Plateau, in winter from 2001 to 2018. There was a large increasing trend in LST (0.116 ± 0.05 ◦ C · a − 1 ) on the Purog Kangri ice ﬁeld during December, while there was no evident LST rising trend in January and February. In December, both the signiﬁcantly decreased albedo ( − 0.002 a − 1 , based on the MOD10A1 V6 albedo product) on the ice ﬁeld surface and the signiﬁcantly increased number of clear days (0.322 d · a − 1 ) may be the main reason for the signiﬁcant warming trend in the ice ﬁeld. In addition, although the two highest LST of December were observed in 2017 and 2018, a longer data set is needed to determine whether this is an anomaly or a hint of a warmer phase of the Purog Kangri ice ﬁeld in December.


Introduction
Land surface temperature (LST), which is an essential variable of the terrestrial energy budget, is a key variable of the surface-atmosphere interface that drives mass and energy exchange between the atmosphere and land surface [1][2][3][4][5], thereby playing a critical role in hydrological, meteorological, environmental, and climate change studies [6][7][8][9][10]. Compared with traditional in situ LST observations that are limited to a small area with relatively homogeneous land-cover types, the thermal infrared satellite remote sensing data available since the 1970s can offer surface temperature information with large spatial coverage and a long time series (depending on the revisit cycle of the satellites). These satellite data have become one of the most important data sources for assessing climate change at regional and global scales [11].
Currently, the sensors most commonly used for LST retrieval include Advanced Very High Resolution Radiometer (AVHRR), Moderate Resolution Imaging Spectroradiometer (MODIS), Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM+), and Thermal Infrared Sensor (TIRS) [12][13][14][15]. MODIS thermal infrared data have been available since 2000 with~1 km spatial resolution, and view the entire surface of the Earth every one to two days [16]. Due to its high spatiotemporal resolution, MODIS LST has been widely used to study spatial and temporal variations in regional surface temperature, especially for polar regions and alpine mountainous areas with harsh environments and sparse observation networks. Many studies have been performed to estimate temperature variations on the Antarctic and Greenland ice sheets using MODIS LST, revealing a warming trend and extreme melting events (LST extremes) in recent years [17][18][19][20]. For the Qinghai-Tibetan Plateau, the MODIS LST product has been applied to the lower-altitude regions (<5000 m) with meteorological observation station data in terms of near-surface temperature estimation, delimitation of permafrost extent, and lake surface temperature analysis [21][22][23][24].
The accuracy of the MODIS LST product on snow/ice surfaces has been extensively studied in the Arctic regions, revealing substantial discrepancies among different regions with diverse altitude and climate conditions, as well as different versions of MODIS LST products [25][26][27][28]. However, the ground-based validation of MODIS LST products is comparatively rare in the high-altitude glaciers of the Qinghai-Tibetan Plateau, possibly due to the very sparse meteorological station network and the relatively small area (mostly <1 km 2 ) of mountain glaciers that are generally represented as mixed pixels of MODIS data. Although Zhang et al. [29] validated the errors of MOD11A1 and MYD11A1 for the Xiaodongkemadi Glacier and the Palongzangbo Glacier on the Qinghai-Tibetan Plateau, their results were based on only two mixed pixels (non-glacier parts accounting for 28% and 32%), which may result in spurious accuracy evaluation of MODIS LST products on mountain glaciers. Currently, it is generally believed that the errors in the MODIS LST product on snow or ice surfaces mainly result from failures in determining cloud contamination, especially of thin clouds. Therefore, when the MODIS LST product is applied on snow or ice surfaces, it is essential to eliminate the pixels affected by cloud via statistical methods [30].
Since the 1980s, mountain glaciers on the Qinghai-Tibetan Plateau have generally suffered from a retreating and thinning trend [31,32]. The role of the EDW effect (elevation-dependent warming) in high-altitude mountainous areas [33,34] is still in question owing to the sparse station distribution and to the errors of both remote sensing surface temperature products and climate models. With the global warming in recent decades, rising temperature has led to intensified glacier ablation on the Qinghai-Tibetan Plateau and consequently unstable runoff of downstream rivers, which might jeopardize the livelihood of hundreds of millions of people in the region [35][36][37]. Moreover, the warming of high-altitude areas is directly related to the occurrence of mountain glacial disasters, such as glacier debris flows, glacial lake outburst floods, and glacier avalanches [38][39][40][41]. Therefore, understanding the warming rate of high-altitude glacial areas has a practical significance for water resources management and safety in the downstream area of glaciers [31]. Although Wu et al. [4] retrieved high-precision mountain glacier surface temperatures from Landsat ETM+ data on the Qinghai-Tibetan Plateau, it is difficult to analyze the interannual variations and trends in the LST retrieved from Landsat thermal infrared data because of 16 day revisit time of the satellite.
As the largest glacier on the Qinghai-Tibetan Plateau of China, the Purog Kangri ice field is an ideal area to study the variations in glacier surface temperature using MODIS LST products. An ice core drilled in the Purog Kangri ice field revealed an evident warming since the 20th century [42]; recent data from optical remote sensing satellites show that the Purog Kangri ice field has shrunk by 15.29 km 2 (about 3.6%) from 1992 to 2014 [43], and multi-source Digital Elevation Model (DEM) data indicate an increasing thinning rate of the glacier from −0.049 ± 0.200 m·a −1 from 2000 to 2012 [44] to −0.317 ± 0.027m·a −1 from 2012 to 2016 [45]. All these findings reflect the fact that the high-altitude mountainous areas of the Qinghai-Tibetan Plateau are influenced by global warming, but there remains a lack of quantitative knowledge of the specific rate of warming. It is predicted that the warming rate in the Qinghai-Tibetan Plateau will increase until the end of the 21st century, with more pronounced increases in winter than in summer [46]. The rapid wintertime warming could lead to a reduced glacier cold storage and an earlier glacier ablation season [47]. In addition, the glacier LST will stop increasing after reaching 0 • C during the ablation season, which will affect the calculation of the trends of summer glacier LST. Therefore, this paper attempts to analyze the interannual variations in winter LST on the Remote Sens. 2020, 12, 1133 3 of 17 glacier from 2001 to 2018 using MODIS LST products after eliminating cloud-contaminated pixels, in order to understand the trend of winter LST on the Purog Kangri ice field in the 21st century.

Study Site
The Purog Kangri ice field (33 • 44´N-34 • 03´N, 89 • 00´E-89 • 20´E) in the north-central part of the Qinghai-Tibetan Plateau (Figure 1a) covers an area of 422.58 km 2 and is the largest ice field in China [48,49]. The ice field is relatively flat over its central part, and is surrounded by more than 50 glacier tongues extending outwards along the mountain valleys, with the longest tongue reaching the foothills (Figure 1). The elevation of the ice field is between 5345 m and 6450 m, and the mean equilibrium-line altitude (ELA is the altitude at which glacier surface accumulation and ablation equal in a year, i.e., the glacier mass balance is equal to 0) during 2001 to 2012 was 5748 ± 37 m [50]. Due to the high altitude and relatively large area, there have only been minor changes in the ice area since the Little Ice Age; in particular the area of the southeastern and western parts is relatively stable, while there is a more obvious shrinkage in the northern part [48].
Remote Sens. 2020, 12, x FOR PEER REVIEW 3 of 18 after eliminating cloud-contaminated pixels, in order to understand the trend of winter LST on the Purog Kangri ice field in the 21st century.

Study Site
The Purog Kangri ice field (33°44´N-34°03´N, 89°00´E-89°20´E) in the north-central part of the Qinghai-Tibetan Plateau (Figure 1a) covers an area of 422.58 km 2 and is the largest ice field in China [48,49]. The ice field is relatively flat over its central part, and is surrounded by more than 50 glacier tongues extending outwards along the mountain valleys, with the longest tongue reaching the foothills (Figure 1). The elevation of the ice field is between 5345 m and 6450 m, and the mean equilibrium-line altitude (ELA is the altitude at which glacier surface accumulation and ablation equal in a year, i.e., the glacier mass balance is equal to 0) during 2001 to 2012 was 5748 ± 37 m [50]. Due to the high altitude and relatively large area, there have only been minor changes in the ice area since the Little Ice Age; in particular the area of the southeastern and western parts is relatively stable, while there is a more obvious shrinkage in the northern part [48].

MODIS LST Products and Processing
The MODIS/Terra and Aqua Land Surface Temperature/Emissivity Daily L3 Global 1 km SIN Grid product is used in this study, which was downloaded from the NASA DAAC (Distributed Active Archive Center) at EROS (Earth Resources Observation and Science) Data Center. Both MODIS/Terra and Aqua pass over the Purog Kangri ice field twice a day, and their overpass times are ~10:30 (MOD11A1_day data)/22:30 (MOD11A1_night data) and ~1:30 (MYD11A1_night data)/13:30 (MYD11A1_day data), respectively. The MOD11A1 and MYD11A1 products are produced with the generalized split-window algorithm that uses MODIS bands 31 and 32 [51] (10.78-11.28 μm and 11.77-12.27 μm, respectively). In 2017, NASA provided a new MODIS LST product (V6), which included the elimination of cloud-contaminated data in MODIS Level 2 and 3, an update of the coefficients of the split-window algorithm, and the adjustment of emissivity values based on ground classification [52]. The detailed information of MODIS LST products and their interactions with other MODIS products can be found in Wan [53,54].

MODIS LST Products and Processing
The MODIS/Terra and Aqua Land Surface Temperature/Emissivity Daily L3 Global 1 km SIN Grid product is used in this study, which was downloaded from the NASA DAAC (Distributed Active Archive Center) at EROS (Earth Resources Observation and Science) Data Center. Both MODIS/Terra and Aqua pass over the Purog Kangri ice field twice a day, and their overpass times are~10:30 (MOD11A1_day data)/22:30 (MOD11A1_night data) and~1:30 (MYD11A1_night data)/13:30 (MYD11A1_day data), respectively. The MOD11A1 and MYD11A1 products are produced with the generalized split-window algorithm that uses MODIS bands 31 and 32 [51] (10.78-11.28 µm and 11.77-12.27 µm, respectively). In 2017, NASA provided a new MODIS LST product (V6), which included the elimination of cloud-contaminated data in MODIS Level 2 and 3, an update of the coefficients of the split-window algorithm, and the adjustment of emissivity values based on ground classification [52]. The detailed information of MODIS LST products and their interactions with other MODIS products can be found in Wan [53,54]. The MOD11A1 and MYD11A1 tile h25v05 provides complete coverage of the Purog Kangri ice field, and data from this tile for the three winter months (December, January, and February) from 2001 to 2018 are used here. The MODIS LST product is reprojected onto a geographic projection using the MODIS Reprojection Tool (MRT) from the original sinusoidal projection, and a total of 212 pixels within the Purog Kangri ice field (according to the Second Glacier Inventory Dataset of China [55]) were used for further study, excluding mixed pixels of glacier and bare rock (Figure 1b). Each grid of MODIS LST product has a quality control (QC) flag ranging from 0 to 3 indicating average errors of <1, 1-2, 2-3, and >3 K. Only pixels which passed quality control at the QC 0/1 level were included for analysis. Cloud contaminated pixels (2) or pixels contaminated for another reason (3) were omitted. The monthly mean available days of four MODIS LST products after quality control (QC) (Figure 2) show that the MOD11A1_day data (red bars) is more appropriate to analyze the variations in LST, while the insufficient availability of the other three data makes it difficult to retrieve a convincing LST trend, the reason for which is also discussed in Section 5.2.
Remote Sens. 2020, 12, x FOR PEER REVIEW 4 of 18 The MOD11A1 and MYD11A1 tile h25v05 provides complete coverage of the Purog Kangri ice field, and data from this tile for the three winter months (December, January, and February) from 2001 to 2018 are used here. The MODIS LST product is reprojected onto a geographic projection using the MODIS Reprojection Tool (MRT) from the original sinusoidal projection, and a total of 212 pixels within the Purog Kangri ice field (according to the Second Glacier Inventory Dataset of China [55]) were used for further study, excluding mixed pixels of glacier and bare rock (Figure 1b). Each grid of MODIS LST product has a quality control (QC) flag ranging from 0 to 3 indicating average errors of <1, 1-2, 2-3, and >3 K. Only pixels which passed quality control at the QC 0/1 level were included for analysis. Cloud contaminated pixels (2) or pixels contaminated for another reason (3) were omitted. The monthly mean available days of four MODIS LST products after quality control (QC) (Figure 2) show that the MOD11A1_day data (red bars) is more appropriate to analyze the variations in LST, while the insufficient availability of the other three data makes it difficult to retrieve a convincing LST trend, the reason for which is also discussed in Section 5.2. The cloud identification algorithm for the MOD11A1_day products still misses some pixels that are affected by cloud contamination, which is primarily related to thin clouds, cloud shadow contamination, or large view angles [56]. To further improve the data quality and reduce the impact of other sources of uncertainty on the data, following Williamson et al. [30], a 3 × 3 grid cell moving window was used to compare values between each central pixel and its 8 surrounding pixels; the central pixel was excluded if the difference between the central value and the mean (or mean plus standard deviation) of all the surrounding pixels exceeded 5 °C (or 3 °C). After screening for "abnormal" data, the three months have different numbers of valid days with sufficient remaining MODIS pixels. To determine the number of valid days required to adequately represent the monthly mean LST value, the pixel with the greatest number of valid days was selected to calculate the RMSE between the LST averages of randomly selected days and the monthly averages based on the entire dataset. The sample size for each number of days is 10,000. From Figure 3, the RMSE shows a rapid decline when the number of days increases from 1 day to 6 days, and then begins to level off. When the number of days is set to 6 as the threshold for the pixels to be used in this study, the RMSE is around 1. This provides acceptable accuracy for monthly average temperature, while permitting more pixels to be studied. Therefore, in this study, those pixels with at least 6 valid days per month The cloud identification algorithm for the MOD11A1_day products still misses some pixels that are affected by cloud contamination, which is primarily related to thin clouds, cloud shadow contamination, or large view angles [56]. To further improve the data quality and reduce the impact of other sources of uncertainty on the data, following Williamson et al. [30], a 3 × 3 grid cell moving window was used to compare values between each central pixel and its 8 surrounding pixels; the central pixel was excluded if the difference between the central value and the mean (or mean plus standard deviation) of all the surrounding pixels exceeded 5 • C (or 3 • C). After screening for "abnormal" data, the three months have different numbers of valid days with sufficient remaining MODIS pixels. To determine the number of valid days required to adequately represent the monthly mean LST value, the pixel with the greatest number of valid days was selected to calculate the RMSE between the LST averages of randomly selected days and the monthly averages based on the entire dataset. The sample size for each number of days is 10,000. From Figure 3, the RMSE shows a rapid decline when the number of days increases from 1 day to 6 days, and then begins to level off. When the number of days is set to 6 as the threshold Remote Sens. 2020, 12, 1133 5 of 17 for the pixels to be used in this study, the RMSE is around 1. This provides acceptable accuracy for monthly average temperature, while permitting more pixels to be studied. Therefore, in this study, those pixels with at least 6 valid days per month are used to calculate the monthly average LST.
Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 18 The number of pixels available (more than 6 days of data) from 2001 to 2018 is shown in Figure  4. In some years, excessive cloud cover has reduced the number of pixels, such as 2008 with only 6 pixels available, which may cause large uncertainty in the LST analysis. Therefore, the years of 2005, 2008, 2010, 2012, 2014, and 2016 with too few pixels are excluded in the analysis of LST variations. After data screening, there are 26 pure snow/ice pixels with six days or more of data in December, January, and February together. The analysis of the LST trend in this paper is carried out on those 26 pure ice pixels (shown in Figure 5).

MODIS Albedo Products
In this study, the MODIS/Terra Snow Cover Daily L3 Global 500m SIN Grid product (MOD10A1 V6) was employed to analyze the trend of albedo; this is available from the NASA DAAC at the National Snow and Ice Data Center (NSIDC). The MOD10A1(V6) contains snow extent, snow albedo, fractional snow cover, and quality assessment data at 500 m resolution and the albedo algorithm within the MOD10A1 algorithm suite was developed by Klein et al. [57]. The MOD10A1 product is generated when the viewing and illumination angles are at the best conditions in a day [58]. The MOD10A1 albedo substantially improved in relative accuracy from V5 to V6 [59]. A description of the MOD10A1(V6) and algorithm can be found in Riggs et al. [60]. The MOD10A1 has been widely The number of pixels available (more than 6 days of data) from 2001 to 2018 is shown in Figure 4. In some years, excessive cloud cover has reduced the number of pixels, such as 2008 with only 6 pixels available, which may cause large uncertainty in the LST analysis. Therefore, the years of 2005, 2008, 2010, 2012, 2014, and 2016 with too few pixels are excluded in the analysis of LST variations. After data screening, there are 26 pure snow/ice pixels with six days or more of data in December, January, and February together. The analysis of the LST trend in this paper is carried out on those 26 pure ice pixels (shown in Figure 5). The number of pixels available (more than 6 days of data) from 2001 to 2018 is shown in Figure  4. In some years, excessive cloud cover has reduced the number of pixels, such as 2008 with only 6 pixels available, which may cause large uncertainty in the LST analysis. Therefore, the years of 2005, 2008, 2010, 2012, 2014, and 2016 with too few pixels are excluded in the analysis of LST variations. After data screening, there are 26 pure snow/ice pixels with six days or more of data in December, January, and February together. The analysis of the LST trend in this paper is carried out on those 26 pure ice pixels (shown in Figure 5).

MODIS Albedo Products
In this study, the MODIS/Terra Snow Cover Daily L3 Global 500m SIN Grid product (MOD10A1 V6) was employed to analyze the trend of albedo; this is available from the NASA DAAC at the National Snow and Ice Data Center (NSIDC). The MOD10A1(V6) contains snow extent, snow albedo, fractional snow cover, and quality assessment data at 500 m resolution and the albedo algorithm within the MOD10A1 algorithm suite was developed by Klein et al. [57]. The MOD10A1 product is generated when the viewing and illumination angles are at the best conditions in a day [58]. The MOD10A1 albedo substantially improved in relative accuracy from V5 to V6 [59]. A description of

MODIS Albedo Products
In this study, the MODIS/Terra Snow Cover Daily L3 Global 500m SIN Grid product (MOD10A1 V6) was employed to analyze the trend of albedo; this is available from the NASA DAAC at the National Snow and Ice Data Center (NSIDC). The MOD10A1(V6) contains snow extent, snow albedo, fractional snow cover, and quality assessment data at 500 m resolution and the albedo algorithm within the MOD10A1 algorithm suite was developed by Klein et al. [57]. The MOD10A1 product is generated when the viewing and illumination angles are at the best conditions in a day [58]. The MOD10A1 albedo substantially improved in relative accuracy from V5 to V6 [59]. A description of the MOD10A1(V6) and algorithm can be found in Riggs et al. [60]. The MOD10A1 has been widely validated and is commonly used in the analysis of albedo variation for the surface of snow/ice [61,62]. The MODIS Albedo product is preprocessed with the same method as the MODIS LST data sets.

The High Asia Refined Data
The LST product of High Asia refined (HAR) data is also utilized to be compared with the MODIS LST, which is provided by the Technical University of Berlin [63]. The HAR used a daily reinitialization strategy to prevent drift from observed synoptic conditions during its total simulation period (from October 2000 to October 2014) over the Qinghai-Tibetan Plateau region. The 10-km resolution HAR Surface Skin Temperature (i.e., LST) dataset was dynamically downscaled from the global analysis data using the Weather Research Forecasting model. Although the estimation of the LST trend using HAR data may have some uncertainty, it can still enrich our understanding of temperature variations over the Purog Kangri ice field with no in situ observation.

Digital Elevation Data
The Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model (ASTER GDEM) V2.0 product (2011) with a spatial resolution of 30 m was employed to analyze the relationship between LST and altitude in order to determine if there is an EDW effect in the temperature increase of the Purog Kangri ice field. These data have been widely used for glacier height extraction and topographic analysis in High Mountain Asia [64].

Spatial Pattern in LST Trends of the Purog Kangri Ice Field
The trend in the wintertime monthly LST during clear daytime over the Purog Kangri ice field is shown in Figure 5, and the nonparametric Mann-Kendall test was used to estimate the significance of the LST trend. December showed the highest rate of warming, with 26 pixels of the ice field having an increasing LST trend of over 0.1 • C·a −1 (22 pixels are significant, p < 0.05). In comparison, although LSTs in January showed increasing trends, they were less than 0.05 • C·a −1 and were not statistically significant (p > 0.05). In February, the increasing trends in LST of most pixels ranged between 0.05 • C·a −1 and 0.07 • C·a −1 , with only one pixel showing a statistically significant increase (p < 0.05).
The trend in the wintertime monthly LST during clear daytime over the Purog Kangri ice field is shown in Figure 5, and the nonparametric Mann-Kendall test was used to estimate the significance of the LST trend. December showed the highest rate of warming, with 26 pixels of the ice field having an increasing LST trend of over 0.1 °C·a -1 (22 pixels are significant, p < 0.05). In comparison, although LSTs in January showed increasing trends, they were less than 0.05 °C·a -1 and were not statistically significant (p > 0.05). In February, the increasing trends in LST of most pixels ranged between 0.05 °C·a -1 and 0.07 °C·a -1 , with only one pixel showing a statistically significant increase (p < 0.05).

Variations in the Regional Mean LST of the Purog Kangri Ice Field
The regional mean of all the selected December LST pixels in the Purog Kangri ice field showed an increasing trend of 0.1585 ± 0.053 • C·a −1 (black dashed line in Figure 6a Figure 6a), showing a relatively robust trend. In comparison, the LST trend in January was smaller (0.0535 ± 0.106 • C·a −1 ) and was not statistically significant (black dashed line in Figure 6b). When the year with the maximum January LST (i.e., 2018) was removed, there was even a cooling trend in January (red dashed line in Figure 6b). In February, the overall LST trend of the Purog Kangri ice field was 0.0735 ± 0.077 • C·a −1 , which was not statistically significant (black dashed line in Figure 6c), and the temperature trend approached 0 • C·a −1 when the maximum February LST (i.e., 2017) was removed (red dashed line in Figure 6c). Therefore, the warming trends in January and February were sensitive to extremes and were not robust. Of note, although during the period 2001-2018 the LST maximum of all three wintertime months was observed in 2017 or 2018 ( Figure 6), a longer data set is needed to determine whether this was an anomaly or just a hint of a warmer phase of Purog Kangri ice field in winter. Figure 6 shows more available days in December than in January and February, meaning that the monthly mean LSTs of the three winter months are calculated by MODIS LST with different numbers of days, which might affect the consistency of calculated trends among different months. To investigate this issue, the mean LST of 6 random days from December, January, and February for each available pixel was used to represent the corresponding monthly mean LST (gray circles in Figure 6). In fact, there were only small differences (within ±0.2 • C) between the original and the 6-day monthly LST with the exception of some February months that had differences of approximately 0.5 • C, which indicates that using a different number of days (≥6 days) in the calculations of monthly LST has little effect on the calculated trends.

Influence of Albedo and Number of Clear Days on LST of the Purog Kangri Ice Field
LST is determined partly by surface albedo as it can change the surface energy flux by altering the amount of solar radiation absorbed by the glacier surface [13]. Previous studies (e.g., Nolin et al. [65]) on the Greenland ice sheet indicate that a slight change in ice/snow albedo can trigger substantial fluctuations in the energy flux on the ice surface. From Figure 7, the surface albedo was significantly (p < 0.05) and negatively correlated with LST in all three months of winter during 2001-2018, and the coefficients of determination were 0.53, 0.49, and 0.35 for December, January, and February, respectively. The surface albedo showed a significant decreasing trend in December at a rate of -0.0026 a -1 , corresponding with the significant rising trend in LST, whereas for January and February there were no significant trends in either albedo or LST.

Influence of Albedo and Number of Clear Days on LST of the Purog Kangri Ice Field
LST is determined partly by surface albedo as it can change the surface energy flux by altering the amount of solar radiation absorbed by the glacier surface [13]. Previous studies (e.g., Nolin et al. [65]) on the Greenland ice sheet indicate that a slight change in ice/snow albedo can trigger substantial fluctuations in the energy flux on the ice surface. From Figure 7, the surface albedo was significantly (p < 0.05) and negatively correlated with LST in all three months of winter during 2001-2018, and the coefficients of determination were 0.53, 0.49, and 0.35 for December, January, and February, respectively. The surface albedo showed a significant decreasing trend in December at a rate of -0.0026 a −1 , Remote Sens. 2020, 12, 1133 9 of 17 corresponding with the significant rising trend in LST, whereas for January and February there were no significant trends in either albedo or LST. Therefore, from 2000 to 2018 the significant increasing trend in December LST was evidently affected by the decreasing albedo and the increasing number of clear days, both of which contributed to the increase in solar radiation absorbed by the glacier surface. During the similar time period of our study, the declining trend in winter cloud cover fraction was also consistent with the increasing trend in the number of clear sky days [66,67]. The higher coefficient of determination (0.53) between albedo and LST than that of the number of clear days (0.48) indicated that the LST warming trend in December was more affected by albedo, while the increase of clear days enhanced the effect of the albedo on LST. However, it remains unclear whether the decrease of albedo is the dominant reason for the rising winter LST on glaciers, which needs further field investigation and energy/mass balance simulations. In addition, as the main energy source of the ice field, solar radiation (net shortwave radiation) was also controlled by the number of clear days during daytime. From Figure 7, there are positive correlations between the number of clear days and LST for all three months, and the coefficients of determination for December, January, and February were 0.48, 0.19, and 0.10, respectively. Only in December, with a significant LST rising trend, did the number of clear days show a significant increasing trend (about 0.32 d·a −1 ), whereas in January and February with insignificant LST change trends, the calculated trend of the number of clear days was not statistically significant (p > 0.05). Therefore, from 2000 to 2018 the significant increasing trend in December LST was evidently affected by the decreasing albedo and the increasing number of clear days, both of which contributed to the increase in solar radiation absorbed by the glacier surface. During the similar time period of our study, the declining trend in winter cloud cover fraction was also consistent with the increasing trend in the number of clear sky days [66,67]. The higher coefficient of determination (0.53) between albedo and LST than that of the number of clear days (0.48) indicated that the LST warming trend in December was more affected by albedo, while the increase of clear days enhanced the effect of the albedo on LST. However, it remains unclear whether the decrease of albedo is the dominant reason for the rising winter LST on glaciers, which needs further field investigation and energy/mass balance simulations.

Variation of December LST with Elevation on the Purog Kangri Ice Field
To investigate whether there are variations of LST with elevation in winter for the Purog Kangri ice field, December was taken as an example because it showed a significant warming trend. The altitude range of the Purog Kangri ice field (5600-6200 m) was divided into intervals of 150 m, 200 m, and 300 m so as to separately calculate the LST trend at different altitudes. Note that there are 115 temporally continuous pixels (yellow grid in Figure 8) available in December because they were extracted separately from the MODIS LST product according to the data screening rules (≥6 days). The surface temperature trend was between 0.106 ± 0.012 • C·a −1 and 0.133 ± 0.012 • C·a −1 , which decreased slightly with increase in elevation regardless of the interval used in the altitude division.

Variation of December LST with Elevation on the Purog Kangri Ice Field
To investigate whether there are variations of LST with elevation in winter for the Purog Kangri ice field, December was taken as an example because it showed a significant warming trend. The altitude range of the Purog Kangri ice field (5600-6200 m) was divided into intervals of 150 m, 200 m, and 300 m so as to separately calculate the LST trend at different altitudes. Note that there are 115 temporally continuous pixels (yellow grid in Figure 8) available in December because they were extracted separately from the MODIS LST product according to the data screening rules (≥6 days). The surface temperature trend was between 0.106 ± 0.012 °C·a -1 and 0.133 ± 0.012 °C·a -1 , which decreased slightly with increase in elevation regardless of the interval used in the altitude division.

Impact of Data Screening Rules on LST Trends
The spatial pattern in the mean number of clear days per winter month (December, January, and February) from 2001 to 2018 shows that the western part of the Purog Kangri ice field is strongly affected by clouds, with only 3-5 clear days available, while the eastern and southern parts are relatively less affected by clouds and thus have 26 pixels selected to study LST variations for all three months (Figure 9).

Impact of Data Screening Rules on LST Trends
The spatial pattern in the mean number of clear days per winter month (December, January, and February) from 2001 to 2018 shows that the western part of the Purog Kangri ice field is strongly affected by clouds, with only 3-5 clear days available, while the eastern and southern parts are relatively less affected by clouds and thus have 26 pixels selected to study LST variations for all three months (Figure 9).
In this study area, the available MODIS data in December was more than in January and February, and 115 pixels were available throughout the period from 2001 to 2018 after data screening. However, the shortage of available MODIS data in January and February meant that six years with too few pixels were removed from trend calculations, which will introduce some errors in the trend analysis when comparing all three winter months. Therefore, the December LST with more available MODIS data was used to estimate the potential influence of data screening with six years missing during the calculation of the LST trend. Comparisons of LST employing three different data sources (using 115 pixels with no years missing, 115 pixels with six years missing, and 26 pixels with six years missing as used in January and February) show almost no difference between the LST variations using 115 pixels and those using 26 pixels (blue and red dots in Figure 10). Therefore, LST using 26 pixels adequately represents the December LST of the corresponding year over the Purog Kangri ice field. However, comparisons show that the December LST trend was overestimated by 0.04 • C·a −1 (the increasing trend was 0.116 ± 0.050 • C·a −1 ) after excluding six years (black and blue dots in Figure 10). The years with less data availability were often accompanied by an increase of cloud cover. Cloud cover can affect the surface energy balance via radiation forcing, and sometimes increases solar radiation at the surface through multiple reflections between high albedo surface and clouds. Therefore, excluding some years with insufficient data might introduce an error in the LST trend. In this study area, the available MODIS data in December was more than in January and February, and 115 pixels were available throughout the period from 2001 to 2018 after data screening. However, the shortage of available MODIS data in January and February meant that six years with too few pixels were removed from trend calculations, which will introduce some errors in the trend analysis when comparing all three winter months. Therefore, the December LST with more available MODIS data was used to estimate the potential influence of data screening with six years missing during the calculation of the LST trend. Comparisons of LST employing three different data sources (using 115 pixels with no years missing, 115 pixels with six years missing, and 26 pixels with six years missing as used in January and February) show almost no difference between the LST variations using 115 pixels and those using 26 pixels (blue and red dots in Figure 10). Therefore, LST using 26 pixels adequately represents the December LST of the corresponding year over the Purog Kangri ice field. However, comparisons show that the December LST trend was overestimated by 0.04 °C·a -1 (the increasing trend was 0.116 ± 0.050 °C·a -1 ) after excluding six years (black and blue dots in Figure  10). The years with less data availability were often accompanied by an increase of cloud cover. Cloud cover can affect the surface energy balance via radiation forcing, and sometimes increases solar radiation at the surface through multiple reflections between high albedo surface and clouds. Therefore, excluding some years with insufficient data might introduce an error in the LST trend. In this study area, the available MODIS data in December was more than in January and February, and 115 pixels were available throughout the period from 2001 to 2018 after data screening. However, the shortage of available MODIS data in January and February meant that six years with too few pixels were removed from trend calculations, which will introduce some errors in the trend analysis when comparing all three winter months. Therefore, the December LST with more available MODIS data was used to estimate the potential influence of data screening with six years missing during the calculation of the LST trend. Comparisons of LST employing three different data sources (using 115 pixels with no years missing, 115 pixels with six years missing, and 26 pixels with six years missing as used in January and February) show almost no difference between the LST variations using 115 pixels and those using 26 pixels (blue and red dots in Figure 10). Therefore, LST using 26 pixels adequately represents the December LST of the corresponding year over the Purog Kangri ice field. However, comparisons show that the December LST trend was overestimated by 0.04 °C·a -1 (the increasing trend was 0.116 ± 0.050 °C·a -1 ) after excluding six years (black and blue dots in Figure  10). The years with less data availability were often accompanied by an increase of cloud cover. Cloud cover can affect the surface energy balance via radiation forcing, and sometimes increases solar radiation at the surface through multiple reflections between high albedo surface and clouds. Therefore, excluding some years with insufficient data might introduce an error in the LST trend. Figure 10. Variations in the December LST (land surface temperature) anomalies (dots) and their linear trends (dashed lines) under different conditions. Those using 115 pixels with no missing years are shown in black (a1 and a2), those using 115 pixels with six years missing are shown in red (b1 and b2), and those using 26 pixels (as used in January and February) with six years missing are shown in blue (c1 and c2).

Impact of the Cloud Detection Algorithm of the MODIS LST Product
The cloud detection algorithm in the MODIS LST product cannot identify all types of clouds. Previous research has indicated about 15% undetected cloud contamination in MODIS LST products [68]. Studies on the Arctic Svalbard Archipelago suggest that >40% of the cloud contamination was mistakenly considered as clear sky conditions in the MODIS products, resulting in a large error between the MODIS temperature value and the measured value [27]. The accuracy of the MODIS LST product is mainly limited by the errors from the cloud detection algorithm, and thus it is necessary to employ data screening rules to eliminate potential cloud-contaminated pixels before the trend calculations [25,69].
From Figure 2, it is difficult to retrieve convinced LST trends from the other three MODIS LST products with such insufficient available data, and therefore, only the MOD11A1_day could be used to analyze the LST variations in winter. The availability of MODIS LST products is closely related with cloud cover over the Qinghai-Tibetan Plateau, where there are frequent convective clouds that can blur the glacier LST retrieved from thermal infrared remote sensing. Over the Tibetan plateau, the convective clouds occur primarily at the period of 14:00-20:00, around 02:00, and around 08:00 [70,71], which are basically close to the overpass time of the other three MODIS LST products (i.e., 22:30, 13:30, and 01:30 for MOD11A1 night, MYD11A1_Day, and MYD11A1_night, respectively), while the MOD11A1_Day passes over the ice field (i.e., 10:30) in the period with relatively less convective clouds. Yu et al. [72] found that the mean daily cloud coverage exceeds 45% over the Qinghai-Tibetan Plateau, and the MODIS quality control (QC) data are related to the cloud detection calculation, while Zhang et al. [73] found that high-quality data account for 50%-70% of the total MODIS data. Even with more stringent data quality controls, such as excluding data with errors >1 K, cloud contamination still exists in MODIS LST products. In this paper, the data screening method is only a compromise to remove some of the influence of cloud contamination and will inevitably reduce the amount of available data. Therefore, eliminating the influence of cloud pollution when applying the MODIS LST product in the Qinghai-Tibetan Plateau remains a major challenge.

Comparisons Between MODIS LST and Station Observations/HAR Dataset
Although the LST trends of the present study were obtained by excluding cloud contamination effects, it is still necessary to check whether they can capture the actual LST variations over the Purog Kangri ice field. As the near-surface temperature is believed to show consistent variations with LST at regional scales [74][75][76], here we compared the 2 m mean air temperature from stations and MODIS LST data with/without data screening, to check if there is a closer relationship between them. In addition, we also added the comparisons between the HAR LST data (10 km) and MODIS LST data on the ice field. Although the estimation of LST trend using HAR data may have some uncertainty, it can still enrich our understanding of temperature variations over the Purog Kangri ice field with no in situ observations. From Figure 11, when the data screening method was utilized, there were obvious increases in correlation coefficients between MODIS LST and the near-surface temperature at both stations and HAR data, especially for February (the month with the most cloud days), which highlights the necessity of a data screening method. Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 18 Figure 11. Comparisons of anomalies in regional mean MODIS LST (land surface temperature) and HAR (High Asia refined) data over the Purog Kangri ice field as well as the near-surface temperature from two nearby stations (Nagqu and Tuotuohe) in three winter months. Left column: the regional MODIS LSTs are obtained from all the available data after quality control, a1-c1 for December, January and February; right column: the regional mean MODIS LSTs are obtained from the available data after quality control and data screening, a2-c2 for December, January and February.

Conclusions
Using the MODIS (MOD11A1) daily LST product, this study analyzed the temporal and spatial variations in winter monthly (DJF) LST over the Purog Kangri ice field from 2001 to 2018 by removing MODIS LST pixels that might be contaminated by clouds. The results showed a substantial increasing LST trend (0.116 ± 0.050 °C·a -1 ) in the ice field during December, while there was no significant LST trend in January and February. The significant reduction in surface albedo of the ice field and the significant increase in the number of clear days may have led to increased solar radiation absorption by the glacier surface, which may be the main reason for the significant warming trend over the ice field in December. This study also showed no evident changes in December LST in different altitude range (5600-6200 m) of the ice field, with only a slight decreasing trend as elevation increases. In addition, the number of pixels available for LST trend calculations was proven to have a limited impact on the trend results, while the elimination of years that were strongly affected by clouds might cause a non-negligible error in the LST trend. Therefore, when using MODIS LST products to analyze LST variations, it is necessary to pay attention to the impact of the eliminated years, and special caution is warranted for results based solely on the years with more clear days. In addition, comparisons between the MODIS LST and near-surface temperature of two nearby stations/HAR LST Figure 11. Comparisons of anomalies in regional mean MODIS LST (land surface temperature) and HAR (High Asia refined) data over the Purog Kangri ice field as well as the near-surface temperature from two nearby stations (Nagqu and Tuotuohe) in three winter months. Left column: the regional MODIS LSTs are obtained from all the available data after quality control, (a1-c1) for December, January and February; right column: the regional mean MODIS LSTs are obtained from the available data after quality control and data screening, (a2-c2) for December, January and February.

Conclusions
Using the MODIS (MOD11A1) daily LST product, this study analyzed the temporal and spatial variations in winter monthly (DJF) LST over the Purog Kangri ice field from 2001 to 2018 by removing MODIS LST pixels that might be contaminated by clouds. The results showed a substantial increasing LST trend (0.116 ± 0.050 • C·a −1 ) in the ice field during December, while there was no significant LST trend in January and February. The significant reduction in surface albedo of the ice field and the significant increase in the number of clear days may have led to increased solar radiation absorption by the glacier surface, which may be the main reason for the significant warming trend over the ice field in December. This study also showed no evident changes in December LST in different altitude range (5600-6200 m) of the ice field, with only a slight decreasing trend as elevation increases. In addition, the number of pixels available for LST trend calculations was proven to have a limited impact on the trend results, while the elimination of years that were strongly affected by clouds might cause a non-negligible error in the LST trend. Therefore, when using MODIS LST products to analyze LST variations, it is necessary to pay attention to the impact of the eliminated years, and special caution is warranted for results based solely on the years with more clear days. In addition, comparisons