Glacier Surface Heatwaves Over the Tibetan Plateau

The Tibetan Plateau (TP) has warmed at a rate twice the global average and presents unique warming patterns in surface temperature changes. However, key characteristics of glacier surface heatwave duration and intensity over the TP during the present extreme warming period are still unknown. In this study, we show that surface temperatures in glacial regions of the TP (0.37 ± 0.10°C per decade) have increased faster than those in non‐glacial areas (0.29 ± 0.05°C per decade) between 2001 and 2020. Moreover, the duration (5.3 ± 3.2 days per decade) and cumulative intensity (24.9 ± 16.3 days °C per decade) of glacier surface heatwaves have increased significantly during autumn. Our results demonstrate an elevation dependence to these key warming characteristics, which we also suggest are associated with extreme glacier mass loss. Here, we highlight potential threats to the sustainability of glacier water resources and increasing risk of glacier related hazards at the “roof of the world.”

Whilst previous studies have investigated glacier surface temperatures on the TP (X. Qin, 2016;N.-L. Wang et al., 2013), observations have largely been limited in terms of their quantity, duration and/or spatial extent. In particular, previous studies have typically focused on analyzing glacier surface temperatures using either infrared thermometry (N.-L. Wang et al., 2013) or thermistor data (X. Qin, 2016) which, respectively, have only lasted a few days or months, and were also only observing a specific region (e.g., Qiumianleiketag Glacier and Laohugou No.12 Glacier). Space based observations of surface temperature derived from, for example, the AVHRR (Stroeve & Steffen, 1998), MODIS (Hall et al., 2006;Mortimer et al., 2016) and/or Landsat (Aubry-Wake et al., 2015;Lo Vecchio et al., 2018) show promise of compensating for the scarcity of in-situ observations for investigating glacier temperatures across the TP. Recent studies have analyzed the performance of glacier surface temperature retrieval from satellite thermal infrared images in this region (Wu et al., 2015). Such studies have demonstrated the accuracy of satellite-derived temperatures when compared with in-situ data, and have also suggested an increase in glacier surface temperature in recent decades in response to a warming world (Liao et al., 2020;Qie et al., 2020;Zhao et al., 2021).
While much is known about mean temperature changes on the TP's glaciers (L. Li et al., 2018;K. Yang et al., 2022), the influence of more frequent and intense extreme warm events on glacier surface temperatures has been relatively unexplored. Particularly, while atmospheric (Fischer & Schär, 2010), lake (Woolway, Anderson, & Albergel, 2021;Woolway, Jennings, et al., 2021, and marine (Frölicher et al., 2018;Laufkötter et al., 2020;Oliver et al., 2018) heatwaves have been investigated previously, the occurrence of heatwaves over glacier surfaces have been largely overlooked. Given the vulnerability of glaciers to temperature change (Bhattacharya et al., 2021;Immerzeel et al., 2010;Kraaijenbrink et al., 2017), here we explore the occurrence of glacier surface heatwaves across the TP. Moreover, we investigate the relationship between glacier surface heatwaves and mass balance throughout the region.

Glacier Surface Temperature Estimate
Glacier surface temperatures investigated in this study were processed from the MODIS Terra night-time LST product, which shows lower bias than other MODIS products over the glacierized TP (H. Zhang et al., 2018). We note that due to frequent cloud contamination in this mountainous region (H. Zhang et al., 2016), daily surface temperature retrieval is challenging. Some daily datasets have been produced by integrating night and day-time data from Terra and Aqua sensors (Chen et al., 2021;T. Zhang et al., 2022). However, glacier surface temperatures in this region fluctuate largely at different times of the same day (W. Yang et al., 2010), which may introduce unquantifiable uncertainty in trend analysis (H. Zhang et al., 2016). In this study, MODIS Terra Land Surface Temperature/Emissivity 8-Day L3 Global 1 km SIN Grid V061 product (MOD11A2) was used, which was improved by undergoing various calibration changes and polarization correction from previous versions (Wan et al., 2021). Only pixels with good data quality (without cloud contamination and with low emissivity and LST error, quality assurance flag = 64, 65, 128, 129) were used in this study. Good quality data were defined according to the MODIS Collection-6 MODIS Land Surface Temperature Products at https://lpdaac.usgs.gov/ documents/118/MOD11_User_Guide_V6.pdf. Furthermore, only pure pixels (i.e., of the MODIS LST gridded data) entirely inside the glacier boundary were selected, in order to minimize the influence of mixed pixels from other land cover types on the temperature retrieval. A total of 14,072 pure pixels are available for glacier surface temperature estimates in this study. The areal proportion of pure glacial pixels relative to the total number of pixels over the glaciated study region is 17% ( Figure S1 in Supporting Information S1).
To ensure that incomplete LST data did not influence our analysis, we excluded observations from pixels where more than 50% of the data were missing ( Figure S2 in Supporting Information S1). For the remaining pixels, we linearly interpolated the time series within each of the three 8-day windows, that is, interpolated within a 24-day period. The LST time series for each pixel was subsequently smoothed with a moving average (46 windows-length; 1-year LST record), which was used to remove the influence of the seasonal temperature signal on the trend estimation. The average LST in spring (March-May), summer (June-August), autumn (September-November), and winter (December-February) was then calculated from the gap-filled data using a 3-window moving average. The non-parametric Mann-Kendall test was subsequently used to determine the presence of a trend and the Theil-Sen's slope was used to estimate the magnitude of change in each pixel (1 km × 1 km longitude-latitude resolution). We consider a p-value of less than 0.05 (two-tailed) as statistically significant.

Glacier Surface Heatwave Duration and Intensity Estimates
In this study, we define a glacial surface heatwave as a period in which LSTs exceed a local and seasonally varying 90th percentile threshold, relative to a baseline climatological mean (the average temperature for the day/ month of year evaluated over the study period, 2001-2020). Following Hobday et al. (2016), the 90th percentile was calculated for each calendar day using daily LSTs within a 3-window size centered on the date across all years within the climatology period, and smoothed by applying a 5-window size (five 8-day LST records, 40-day) moving average. The choices of a 24 and 40-day window for the moving average were motivated to ensure sufficient sample size for percentile estimation and a smooth climatology. A seasonally varying threshold allows identification of anomalously warm events at any time of the year, rather than events only during the warmest months. When one record from MOD11A2 exceeds the threshold, we assume that one heatwave event has started ( Figure  S3 in Supporting Information S1). A specific heatwave event ends when the LST is lower than the percentile threshold. The interval between the start and end dates of the heatwave event is defined as its duration. The glacier surface heatwave intensity is defined relative to the LST anomaly, that is, the difference between the MODIS LST and the local climatology. The cumulative glacier heatwave intensity is defined as the sum of all LST anomalies during a heatwave event. For estimating the pixel-wise trends in the duration and cumulative intensity of heatwaves, we used a linear regression model, and a two-tailed F-test to assess the significance of the estimated trend. A p-value of less than 0.05 (two-tailed) is considered statistically significant.

Definition of Extreme Mass Loss
The monthly glacier surface elevation change between 2001 and 2019 based on ASTER DEM (Hugonnet et al., 2021) was utilized to analyze annual and seasonal glacier elevation changes. The glacier-area weighted elevation change rate at each mountain basin is estimated and fitted for annual and seasonal trends. For each glacier, the 90th percentile of annual elevation change during 2001-2019 was estimated as a threshold according to Vargo et al. (2020). If the annual glacier elevation change rate is below the threshold, the glacier was considered to have an extremely low mass balance in that year (extreme year) ( Figure S3 in Supporting Information S1). For glaciers with a positive mass balance, extreme mass balance corresponds to extremely low accumulation relative to the other periods. Glacier mass loss are estimated by using glacier elevation change from Hugonnet et al. (2021), area from RGI-Consortium (2017), and an ice density of 850 ± 60 kg/m 3 (Huss, 2013). A linear model was used to estimate the trends and an F-test used to quantify their significance. We only estimated the extreme mass loss for pure glacier pixels that matched with the MODIS LST data used.

Meteorological Data and the Validation of MODIS LST
Daily near-surface air temperature from China Meteorological Administration (CMA) stations between 2001 and 2020 was used in this study ( Figure 1). The daily data was averaged to 8-day intervals to match with the temporal resolution of MODIS LST 8-day product (MOD11A2). The in-situ measurements of surface temperature in 2019 for Guliya, Naimona'nyi and Dunde Glaciers were used as validation (W. Yang et al., 2021). These data are collected every 30 min, and only data coinciding with MODIS night view time at around 01:30 a.m. were used. The correlation of the smoothed temperature series between the near-surface temperature from the CMA stations and MODIS LST is strong (r = 0.78) ( Figure S4 in Supporting Information S1), indicating the robustness of the MODIS LST product for trend analysis (G. Zhang et al., 2014). The MODIS LST was compared with near-surface temperature from AWS only inside the pure pixels for Guliya, Naimona'nyi and Dunde Glaciers ( Figure S5 in Supporting Information S1). The bias of these two datasets ranged from 3.7 to 8.7°C (r = 0.85 to 0.97). The LST is lower than the near-surface temperature on three glaciers, which may be due to low incoming solar radiation and emission of longwave radiation from the snow surface during the night as suggested by Adolph et al. (2018). The colder offset with decreasing temperature was related to the solar zenith angle of the MODIS sensor (Shuman et al., 2014). In addition, the albedo derived from the Global Land Surface Satellite (GLASS) products developed by Liang et al. (2013) was also used for possible explanation of the calculated trends in glacier surface temperatures/heatwaves.

Glacial and Non-Glacial Surface Temperature Trends
The TP experienced an overall annual warming trend of 0.29 ± 0.05°C/dec between 2001 and 2020, but with clear seasonal differences (Figure 2a). The warming pattern is most pronounced in autumn, with surface temperatures increasing at a mean rate of 0.59 ± 0.37°C/dec. High TP-wide warming rates were also found during spring (0.28 ± 0.47°C/dec) with the majority of the TP, except the southern plateau (such as western-central Himalayas), experiencing considerable warming. A spatially heterogeneous, and less pronounced, warming pattern is observed during summer (0.23 ± 0.35°C/dec). During winter, a larger spatial region of the TP experienced a cooling trend, but with the northwestern plateau (e.g., eastern Hindu Kush and western Pamir) experiencing considerable warming. An overall winter average LST trend of 0.21 ± 0.76°C/dec is estimated between 2001 and 2020.
Both glacial and non-glacial regions on the TP experienced significant warming from 2001 to 2020 at annual timescales, with a rate of 0.37 ± 0.10 and 0.29 ± 0.05°C/dec, respectively (Figures 2a and 2b). The highest statistically significant warming rates were observed in autumn with the least spatial heterogeneity (higher than ∼0.5°C/dec). The warming rates of glacial areas were weak in spring (0.31 ± 0.59°C/dec) and even suggesting a slight cooling pattern (−0.08 ± 0.89°C/dec) during winter. In contrast, non-glacial areas experienced significant warming during these times of the year (higher than ∼0.2°C/dec). Interestingly, autumn surface temperatures in glacial regions are warming (1.33 ± 0.94°C/dec) twice as fast as those in non-glacial regions (0.57 ± 0.37°C/dec), but with large inter-annual variability.
At annual scales, the EDW of non-glacial areas appeared at altitudes up to 3,000 m a.s.l., and then somewhat stabilized between 3,000 and 4,000 m a.s.l. (Figure 2c). The warming rate declined with an increase in elevation from 4,000 to 6,000 m a.s.l. The EDW of glacial areas showed similar patterns with non-glacial areas at similar altitudes. The seasonal features of EDW in spring, summer, and winter for glacial and non-glacial areas showed similar patterns to those at annual timescales. However, the number of observations with a statistically significant temperature trend was very low in high-altitude regions during these seasons. In autumn, the season with the greatest number of statistically significant temperature trends, the EDW exhibited a continuous increase, particularly in glacial areas.

Glacier Surface Heatwave Duration and Intensity
Across the TP, the duration of glacier surface heatwaves showed a remarkable increase at annual timescales from 2001 to 2020 (Table S1 in Supporting Information S1), particularly in the Western Kunlun, Tibetan interior and Nyainqêntanglha Mountains (Figure 3a). Glacier surface heatwaves are more intense in regions with high surface temperature variability ( Figure S7 in Supporting Information S1). The annual glacier heatwave duration followed a statistically significant long-term trend of 7.4 ± 11.4 days/dec between 2001 and 2020, but with noticeable peaks in heatwave duration in some years (e.g., 2016) (Figure 3b). The duration of glacier surface heatwaves has also increased in many regions across the TP during autumn (5.3 ± 3.2 days/dec, p < 0.05), and to a lesser extent during summer (1.7 ± 4.0 days/dec) and winter (1.8 ± 5.2 days/dec). In contrast, a negative trend in heatwave The dashed line denotes trends with a p-value of less than 0.05 (two-tailed). The histogram shows the frequencies of glacier pixels with statistically significant trends. The observation frequencies of non-glacier pixels with significant trends are shown in Figure S6 of Supporting Information S1. The gray square indicates that there is glacier area coverage, but no pure glacier observation.
duration was calculated across the TP during spring (−0.5 ± 3.5 days/dec). The cumulative intensity of surface heatwaves followed similar temporal and spatial patterns to those in heatwave duration. At annual scales, cumulative heatwave intensity increased over much of the TP, at an average rate of 31.0 ± 61.2 days °C/dec between 2001 and 2020 (Figure 3c). An increase in heatwave cumulative intensity was also calculated in winter (9.7 ± 3.1 days °C/dec) and autumn (24.9 ± 16.3 days °C/dec), the rate of which were noticeably larger in the Western Kunlun and in the central and southern regions of the TP. A weaker long-term trend in cumulative intensity was observed during summer (5.1 ± 19.3 days °C/dec), with some regions showing a long-term decline. Finally, an overall decrease in cumulative intensity was observed during spring (−4.2 ± 18.9 days °C/dec).
Regarding the relationships between the studied heatwave metrics and surface elevation we find that, at annual timescales, trends in heatwave duration were predominantly positive at each elevation bin (Figure 3d), but followed a declining pattern from 3,500 to 4,500 m a.s.l., and an increase at elevations higher than 4,500 m a.s.l. A similar altitudinal dependence was also observed in summer. The rate of change in heatwave duration also increased with elevation during both spring and autumn, but experienced an overall decrease with elevation during winter. Generally, trends in cumulative heatwave intensity show a similar elevation dependence to those  Figures S8 and S9 of Supporting Information S1. The gray pixels indicates that there is glacier area coverage, but no pure glacier observation.
in heatwave duration. The trends of annual heatwave intensity indicated an overall non-significant decrease with increasing elevation, with similar tendencies calculated in summer and winter. In contrast, our observations suggested marginally increasing trends of heatwave intensity with elevation in spring, summer, and autumn.

Relationship Between Glacier Heatwave Duration/Intensity and Extreme Mass Loss
The calculated glacier mass loss demonstrated a considerable acceleration during the observational period, especially in the Himalayas and the Nyainqêntanglha Mountains ( Figure S10 and Table S2 in Supporting Information S1). The Karakoram and East Kunlun region show a slight negative balance. Although the average mass balance is slightly positive for West Kunlun and East Pamir, the acceleration exhibits high negative values. In this context, the extreme mass loss increased largely during the study period (2000-2019) (Figure 4a). The trends of extreme mass loss were greatest (and statistically significant) in summer and autumn. Our results demonstrate a significant increase in extreme mass loss in three seasons (spring, summer, and autumn).
The annual glacier surface heatwave duration and cumulative intensity showed a negative and statistically significant relationship with annual extreme mass loss and their variations (Figures 4b and 4d). Both annual and summer glacier surface heatwave duration/intensity showed a strong relationship with the calculated extreme mass loss. Regarding the seasonal patterns observed, the duration and intensity of glacier surface heatwaves in summer and autumn could trigger more extreme mass loss and even continue to influence the subsequent mass balance. However, for other seasons the relationship shows a weak relationship. The most extreme heatwaves correspond to the maximum variation of extreme mass loss (Figure 4c). This correlation is strong and statistically significant at annual timescales.

Pure Glacier Pixels and the Observation Period
In this study, we analyzed data only from pure glacier pixels (i.e., those that are entirely inside the defined glacier boundary) over the TP that matched with the location of the satellite derived LST from MODIS. By focusing only on pure pixels, the amount of available data over the glaciated TP was considerably lower than had we used all available observations. Critically, this criterion had a substantial impact on data availability within valley glaciers, such as the Himalaya region, the width of which is less than the 1 km resolution of the MODIS LST. Given this choice of criteria, only 17% of the total glaciated area of the TP was considered in this study. Granted that the inclusion of mixed pixels (i.e., of different land cover types) could introduce a considerable source of uncertainty in our observations, we also performed a sensitivity analysis to quantify the magnitude of change in LST had mixed pixels been considered. By including mixed pixels, the ratio of available data increased to 93%, and we estimate a noticeable increase in mean temperature over the study region ( Figure S11 in Supporting Information S1). Moreover, the calculated glacier surface temperature warming rate for all glacier pixels (0.32 ± 0.19°C/dec) was also different to that calculated using only the pure pixels (0.37 ± 0.07°C/dec, p < 0.01). While we think that it is important to highlight these differences, we also believe that using only pure pixels, as we have done in this study, is most appropriate to obtain robust estimates of glacier surface temperature changes. The use of higher resolution satellite data in the future could allow an increase in the spatial coverage of LST and should be considered in order to provide a more detailed analysis of the study region.
Furthermore, as new satellite data become available, future studies will also be able to increase the time period for which the LST climatology is calculated. In this study, the length of the satellite-derived temperature record was 20 years, which is less than the number of years recommended (30 years) by the World Meteorological Organization (WMO, 2011), to estimate a climatology. Thirty years is often considered a requirement to provide a robust baseline to define an extreme event and thus the occurrence of heatwaves. In turn, this could be considered a limitation of our investigation. However, we note that previous heatwave studies have also relied on the use of satellite Earth Observation data lasting less than 30 years (Schlegel et al., 2019;Woolway, Kraemer, et al., 2021). These studies have suggested that a 20-year period is sufficient for defining heatwave intensity (in both lake and marine environments) and duration with average differences of less than 5% when compared to those defined using a 30-year climatological period (Schlegel et al., 2019).

Seasonal and Spatial Differences in Glacier Surface Warming
The seasonal and spatial differences in glacier surface temperature trends are influenced by a combination of local driving factors as well as numerous feedback mechanisms that influence the sensitivity of the TP to climate change (Pepin et al., 2018;You, Chen, et al., 2020). For example, while the seasonal warming of glacier surfaces is significantly related to air temperature anomalies ( Figure S12 in Supporting Information S1), a positive feedback could also emerge with warmer temperatures leading to an increase in glacial melt and a subsequent decrease in albedo. In turn, this would lead to a higher absorption of short-wave radiation and thus an increase in glacier surface temperature (Cuffey & Paterson, 2010;Fujita & Ageta, 2000;Kang et al., 2020). Indeed, in this study, we demonstrated that the glacier surface warming rate is highest in autumn, mainly in 2016 and 2020 ( Figure S12 in Supporting Information S1), which could be explained by the decrease in albedo, as observed from the GLASS products ( Figure S13 in Supporting Information S1). A decrease in albedo was also observed in summer ( Figure S13 in Supporting Information S1). However, at this time, enhanced ice/snow melt simultaneously absorbed more heat (Cuffey & Paterson, 2010), which likely suppressed summer surface warming to some extent. There is an exception to this feedback in winter when the relationship between the change in glacier surface temperature and albedo is reversed. Also, other factors including cloud cover, water vapor feedback and radiative fluxes, aerosol feedback as well as their interactions likely play a role in winter temperature change on the TP (Pepin et al., 2015(Pepin et al., , 2018You, Chen, et al., 2020).
Spatially, the cooling phenomenon that appears in the southwest TP in spring, corresponds to greater snowfall in this region, which likewise increases the albedo of the glacier surface, ultimately leading to lower temperatures (Guo et al., 2019b). Moreover, owing to the atmospheric circulation patterns and the orography of the Karakoram regions, glaciers here received more solid precipitation (e.g., snow, hail, etc.) (Farinotti et al., 2020;Kapnick et al., 2014). The reduction in glacier surface albedo in the Karakorum region is smaller, and there even exists an increase for some glaciers (Y. Zhang et al., 2021), corresponding to less warming or cooling (Figure 3a). On the contrary, annual glacier surface albedo is decreasing in other regions of the TP, especially in the Tibetan interior and Nyainqêntanglha (Y. Zhang et al., 2021), corresponding to the largest warming trends (Figure 3a).

Implications of Glacier Surface Heatwaves
Our study suggests that high elevation areas have experienced, in recent years, a greater increase in heatwave duration both annually and seasonally between spring and autumn (Figure 3). Under scenarios of continued global warming, climate model projections not only suggest an acceleration of increasing temperatures in the TP , but also that EDW will climb to higher elevation (Guo et al., 2016). This means that the TP will not only experience accelerated warming this century, relative to historic trends, but that high elevation regions will warm at an even faster rate. As a result, longer lasting and more intense glacier surface heatwaves are likely to occur with higher mean temperature ( Figure S14 in Supporting Information S1). This has also recently been suggested in both marine (Oliver, 2019) and lake environments (Woolway, Jennings, et al., 2021). Future climate predictions suggest an emergence of longer lasting extreme temperature this century (L. Yin et al., 2019;Y. Zhang et al., 2006). Our study also suggests that peaks in heatwave duration coincided with anomalous El Niño years (and thus a high El Niño index (Trenberth, 2020)) during the study period (e.g., 2016). El Niño events are expected to become more frequent this century (Ying et al., 2022) which could, in turn, further increase the duration of glacier surface heatwaves in this climate-change-sensitive region.
An increase in temperature across the TP will likely lead to spatial shifts in the distribution of heatwave intensity, particularly across elevation gradients. Our study also suggests that, during the historic period, the intensity of summer and annual heatwaves is greatest at high elevation (i.e., the coldest regions) ( Figure S15 in Supporting Information S1). In the West Kunlun, Qilian Mountains, and some other areas with extremely low glacier surface temperatures (Figure 1), the duration and intensity of glacier heatwaves will likely continue to increase. In the Nyainqêntanglha Mountains, the intensity of glacier surface heatwaves may only increase at high altitudes where surface temperature is climatologically colder.
Our study identified a strong correlation between glacier surface heatwaves and extreme mass loss, suggesting that an increase in the occurrence, intensity, and duration of heatwaves will lead to higher mass loss. However, it is important to consider that this relationship may be influenced solely by an increase in mean air temperature, which are causing glaciers to thin and, at the same time, resulting in an increase in glacial surface heatwave duration and intensity. Interannual variability of precipitation may also contribute to extreme mass change, especially for glaciers with positive mass balance in Karakoram (Q. Wang et al., 2017;Yao et al., 2012), the magnitude of which may change differently to air temperature in the future. Indeed, how these interactions change in the future will have important consequences for the region.

Conclusions
In this study, we investigated glacier surface warming trends, the duration and intensity of glacier surface heatwaves, and their elevation dependence on the TP during 2001-2020. Our study demonstrated that glacial areas of the TP have experienced higher warming rates in recent decades. We also highlighted key seasonal differences in warming trends between glacial and non-glacial areas. Glacial areas showed a relatively weak warming trend in summer, whereas non-glacial areas experienced high warming rates. Our analysis also suggested that the warming trend in autumn was two times higher in glacial versus non-glacial areas as well as experiencing the least spatial heterogeneity.
The duration and cumulative intensity of glacier surface heatwaves increased between 2001 and 2020, coinciding with an increase in extreme glacier mass loss. In brief, our analysis suggests that longer lasting and more intense heatwaves are strongly correlated with extreme glacier mass loss. The duration and intensity of glacier surface heatwaves have increased considerably toward the northern TP as well as at higher elevations, which might threaten the sustainability of glacier water resources and increase the risk of glacier related hazards. We urge that additional in-situ observations and high spatial resolution satellite data of glacier surface dynamics are needed to improve our understanding of glacier responses to a rapidly warming world.

Conflict of Interest
The authors declare no conflicts of interest relevant to this study.