The West Kunlun Glacier Anomaly and Its Response to Climate Forcing during 2002–2020

: Research into glacial mass change in West Kunlun (WK) has been sufﬁcient, but most of the existing studies were based on geodetic methods, which are not suitable for speciﬁc health state anal-yses of each glacier. In this paper, we utilize Advanced Spaceborne Thermal Emission and Reﬂection Radiometer (ASTER) imagery, applying the continuity equation to obtain altitudinal speciﬁc mass balance (SMB) for 615 glaciers (>2 km 2 ) during 2002–2011, 2011–2020, and 2002–2020 to research glacial health and its response to climatic forcing. The results show dissimilar glacier SMB patterns between 2002–2011 (0.10 ± 0.14 m w.e. a − 1 ), 2011–2020 (–0.12 ± 0.14 m w.e. a − 1 ) and 2002–2020 ( − 0.01 ± 0.07 m w.e. a − 1 ). Additionally, the glacier equilibrium line altitude (ELA) in WK was 5788 m, 5744 m, and 5786 m, respectively, and the corresponding accumulation area ratios (AARs) were 0.59, 0.62, and 0.58, during 2002–2011, 2011–2020, and 2002–2020, respectively. Regarding glacier response, compared with the ordinary-least-square (OLS) model, the artiﬁcial neural network (ANN) model revealed a respectively less and more sensitive glacier SMB response to extreme negative and positive summer skin temperatures. In addition, the ANN model indicated that the glacier ELA was less sensitive when the integrated water vapor transport (IVT) change exceeded 0.7 kg m − 1 s − 1 . Moreover, compared with IVT ( − 121.57 m/kg m − 1 s − 1 ), glacier ELA shifts were chieﬂy dominated by summer skin temperature (+154.66 m/°C) in the last two decades. From 2002–2011 and 2011–2020, glacier SMB was more susceptible to summer skin temperature ( − 0.38 m w.e./°C and − 0.16 m w.e./°C, respectively), while during 2002–2020, it was more inﬂuenced by IVT (0.45 m w.e./kg m − 1 s − 1 ). In contrast with eastern WK, glaciers in western WK were healthier, although mitigation measures are still needed to safeguard glacier health and prevent possible natural hazards in this region. Finally, we believe that the inconsistent change between glacier SMB and ELAs from 2002–2020 was connected with ice rheology and that the combined effects of skin temperature and IVT can explain the WK glacier anomaly.


Introduction
In the past two decades, glaciers worldwide have experienced accelerated ice mass loss due to global warming [1,2], which will result in more than one billion people facing water shortages and food insecurity in the upcoming three decades [3,4]. High Mountain Asia (HMA), regarded as "The Third Pole" [5], contains approximately 7.0 ± 1.8 × 10 3 km 3 ice [6] and includes 95,536 glaciers with a total area of approximately 97,605 ± 7935 km 2 , based on the Randolph Glacier Inventory version 6.0 (RGI v6.0) [7,8]. Accordingly, it is crucial for glaciologists to research glacier change in HMA [9]. Previous studies have utilized the ERA5.1 reanalysis dataset [48] to yield the mean annual skin temperature and integrated water vapor transport (IVT) trends during 2002-2020, and estimate the nonlinear response to these of glacier SMB and ELA in WK in the last two decades. Overall, the contribution of this work can be summarized into three points: (A) obtaining via continuity equation the multi-temporal glacier SMB in WK and its two subregions during 2002-2011, 2011-2020, and 2002-2020; (B) estimating via Dice coefficient the multitemporal glacier ELA and AAR variation in WK in the past two decades; (C) researching nonlinear glacier SMB and ELA response to climate change in WK during 2002-2020.

Study Area
WK (35 ∘ -37 ∘ N, 76 ∘ -83 ∘ E) is widely covered by glaciers (see Figure 1). We divided WK into two subregions, western WK and eastern WK. Based on RGI v6.0, WK includes 4819 glaciers, and the total area of glaciers (>0.01 km ) is 7855 km , making it one of the most extensively glaciated areas in China. In western WK and eastern WK, there are 2958 glaciers (2899 km ) and 1861 glaciers (4956 km ), respectively. These glaciers are characterized as cirque, valley, and ice-cap types. Herreid and Pellicciotti [34] found that the area of debris-covered glaciers in WK was 245 km . The Guliya Ice Cap is the largest ice cap in mainland Asia (~370 km ) [49], and the Duofeng glacier, the largest mountain glacier in the WK, has an area of ~240 km . The main ridge of WK separates the steep northern slope, descending to the Tarim Basin, and the southern slope, which has a lower gradient towards the Tibetan Plateau [50,51]. The local climate is frigid and dry and is primarily controlled by midlatitude westerlies. The annual precipitation and mean temperature near the ELA (~5930 m a.s.l.) are approximately 300 mm (mainly concentrated in May-September) and −13.9 ℃, respectively [52]. Summer accumulation is dominant for glaciers in WK [53].

Data
This study calculated glacier SMB by the continuity equation, leveraging Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) DEMs, ITS_LIVE ice velocity, and ice thickness during 2002-2011, 2011, and 2002. In addition, we used the ERA 5.1 reanalysis dataset to obtain the annual mean skin temperature and IVT trends in WK from 2002-2020 (see Table 1).

Glacial Elevation Change Generation
To provide full coverage of glaciers in WK in 2002, we used ASTER DEMs from 2001, 2003, 2004, and 2005 to compensate for the dearth of ASTER image coverage in 2002. This situation also occurred in 2011 and 2020. We chose 2002, 2011, and 2020 as representatives because the greatest quantities of ASTER imagery were available for these three time points. We calculated the elevation change for each glacier within its own specific time offset (Supplementary Information Discussion). In addition, we chose predominantly winter and spring ASTER images because they could minimize seasonal offset when acquiring glacier elevation change (hereafter referred to as "seasonal bias") [54]. To account for the seasonal bias, we used a winter accumulation rate of +0.15 m per month [55,56].
We used MicMac ASTER software to obtain ASTER DEMs. We chose this software because it could provide higher-accuracy DEMs with sensor motion ("jitter") correction, compared with silcAst software [57]. Next, we applied the methods proposed by Nuth and Kääb [58] to coregister DEMs for 2002, 2011, and 2020 (excluding glacier area using RGI v6.0), to obtain glacier elevation change results for the three different periods. Here, we used ASTER DEMs from 2020 as a primary map to coregister ASTER DEMs from 2002 and 2011. In addition, during the coregistration procedure, we applied fifth-order polynomial correction for curvature bias, which has been described and applied previously [11,59]. One example of the coregistration results can be seen in Figure 2. Here, we used Median and Normalized Median Absolute Deviation (NMAD) as statistical tools to evaluate the effects of DEM coregistration [60].
Furthermore, after we obtained glacier elevation change in WK, we used ASTER DEMs for the three different periods to acquire the corresponding slope map, and then excluded from the slope map glacier elevation change pixels > 45 • to eliminate erroneous pixels in the glacier elevation change maps [61].

Continuity Equation
In this study, we used the continuity Equation (1) to calculate SMB for each glacier i WK [62][63][64]: is the annual rate of elevation change on the glacier surface and b is th annual SMB, including the surface, internal, and basal mass balance combined; her however, we refer to b as the surface mass change, assuming that frontal ablatio englacial and subglacial components are negligible [65,66]. q ∇ ⋅  is the annual ice flu divergence. All these equation elements are in the unit of meters of ice equivalent per yea

Continuity Equation
In this study, we used the continuity Equation (1) to calculate SMB for each glacier in WK [62][63][64]: where ∂H ∂t is the annual rate of elevation change on the glacier surface and b is the annual SMB, including the surface, internal, and basal mass balance combined; here, however, we refer to b as the surface mass change, assuming that frontal ablation englacial and subglacial components are negligible [65,66]. ∇ · → q is the annual ice flux divergence. All these equation elements are in the unit of meters of ice equivalent per year (m, i.e., a −1 ). To transfer each unit into meters of water equivalent per year (m w.e. a −1 ), we can modify and reorganize Equation (1) into Equation (2): where b w refers to SMB in meters of water equivalent per year [23]. ρ ∂H and ρ ∇· → q are ice densities on the glacier surface, ρ H 2 O is the water density (999.972 kg m −3 ).

Ice Flux and Flux Divergence
To acquire glacier SMB, it was necessary to obtain ice flux divergence. We calculated the ice flux vector at each pixel according to Equation (3), with the ITS_LIVE multiyear mean ice velocity dataset for three periods (2002-2011, 2011-2020, and 2002-2020) and the modelled glacier thickness dataset. Surging glaciers were not excluded, because Yasuda and Furuya [51] concluded that there were only three surging glaciers in the West Kunlun Shan from 2000 to 2014, while Dehecq et al. [47] reported that there were only seven glaciers in a surging state.
where → q represents the ice flux vector, H is the ice thickness in meters, γ → V s is the columnaveraged glacier velocity, and the constant γ represents the relative importance of basal motion and vertical ice shear deformation. To determine γ for each glacier, we applied a Monte Carlo analysis to estimate the depth-integrated velocity per pixel [23].
The ice flux divergence represents the vertical component of glacier surface velocities, resulting in submergence in divergent flow zones and emergence in convergent flow zones. Based on this, we calculated ice flux divergence ∇ · → q using Equation (4) [23]: where ∂ → q x ∂x and ∂ → q y ∂y represent ice flux divergence values in the x and y directions, respectively.

Density Correction
In this study, we used transformed Equation (2) to obtain SMB as the unit of meters of water equivalent per year (m w.e. a −1 ). Thus, it was necessary to consider glacier surface density variations. We did not consider temporal mean ice density variation, because this led to an uncertainty of only 2% [67]. Nevertheless, we had to consider the ice density variation between the accumulation and ablation zones of each glacier, and the dissimilarity of snow and ice densities. In general, snow and firn usually have a lower density than clean ice [68].
We took three steps to correct glacier density. First, we assumed that clean ice was the dominant type spreading over the glacier surface; thus, we gave each glacier an ice density of 900 kg m −3 . After that, we combined the elevation change signal and flux divergence to classify the accumulation and ablation zones for each glacier. If elevation change and flux divergence were both positive, we regarded this area as an accumulation zone and assigned a density of 600 kg m −3 . If they both represented negative values, we assumed this corresponded to the glacier ablation zone, and a density of 900 kg m −3 was assigned. Finally, for the rest of the glacier area, we assumed a density of 850 kg m −3 . We assumed that glacier density uncertainty was 60 kg m −3 [23].

Uncertainty of Glacial Elevation Change
We suggest that elevation change uncertainty involves three factors: residual bias related to DEM coregistration, glacier area change, and voids in the DEM difference map. Here, we calculated the glacier elevation change uncertainty by Equation (5): where E ∆h represents the glacier elevation change uncertainty and E coreg is the uncertainty associated with the residual error in the stable area of the DEM difference maps after DEM coregistration. E area represents the uncertainty caused by glacier area changes in the last 20 years. E voids refers to the void-filling method [69], leading to offsetting of glacier elevation change.
We assert that in the same slope range, the elevation changes in the stable area and glaciated area have the same elevation change offset patterns, which has been described previously [70,71]. To acquire stable regions, we excluded lakes according to the Hydro-LAKES v1.0 dataset [72]. We also tried to remove glacial lake impacts [73] from the DEM difference maps, but we found very few glacial lakes in WK (11 glacial lakes with a total area of~5.28 km 2 ). Hence, we ignored glacial lake influences when estimating E coreg .
In this study, we generated a slope map for each glacier and stable area. Then, we classified slope with a bin interval of 1 • . We took the DEM difference of each slope bin in the stable region as a reference to estimate E coreg , which can be expressed by Equation (6): where E coreg represents the uncertainty of each slope band for each glacier, n is the number of slope bin, dh i is the mean elevation difference at each slope bin. Furthermore, we assumed that the uncertainty associated with glacier area change during 2002-2020 was mainly focused on glacier termini; thus, we chose a 250 m altitude band above each glacier terminus as the ablation area. We assumed that 250 m was enough for glacier recession, since previous studies found that glaciers in WK mainly show advance [74]. Afterwards, we extended this ablation area (250 m above glacier termini) with an outward buffer distance of 500 m to compensate for the glacier mass outside RGI v6.0, and the RGI v6.0 boundary uncertainty. This determined the maps of two regions: the ablation zone and the buffered zone. To account for the contribution of ice mass loss to ice velocity, which flows into the buffered region and the uncertain RGI v6.0 boundary, we assumed that approximately 5% of buffered regions are covered by ice mass. Therefore, we used this ratio (5%) to multiply the mean difference of buffered region for each glacier to estimate E area , which can be expressed by Equation (7): where dh i represents the elevation difference per pixel in the 500 m-buffered region of each glacier, n is the number of pixels in the 500-m buffered region, A ratio is 5%. Unavoidably, we needed to exclude some outliers after differencing DEM maps for different periods, which led to some gaps in glaciated areas. For the period 2002-2020, we excluded pixel values outside the range of −17 to +14 m. For 2002-2011 and 2011-2020, we judged that the reasonable glacier elevation change rates were −13∼33 m and −50∼10 m, respectively. To estimate the uncertainty caused by gaps in glaciated areas, we divided each glacier into 50 m interval elevation bands and then excluded outliers at each elevation band (>3σ). Then, we used the mean glacier elevation change for each bin to fill the gaps [75].
Although we filled the gaps, glacier elevation change offsets could still occur; thus, we needed to account for the uncertainty E voids caused by the void-filling approach, which can be calculated by Equation (8): where N is the pixel count in the glacier voided area. C is a factor that is determined by Equation (9). Here, we selected dh = −0.0440 m and σ dh = 6.1212 m [75].

Uncertainty of SMB Results
Van et al. [64] assumed that the main source of SMB uncertainty E ∆M can be primarily attributed to the inaccuracy of surface elevation change and ice flux divergence. Therefore, in this study, we calculated E ∆M by Equation (10): where E ∆h is the glacier elevation change uncertainty, E ∇q is the ice flux divergence uncertainty, ρ ∆h and ρ ∇q represent glacier surface density variations (see Section 2.4.3, E ρ∆h and E ρ∇q represent glacier surface density uncertainty, ∆h is the annual glacier elevation change, |∇q| is the magnitude of ice flux divergence. To calculate E ∇q , Miles et al. [23] assumed that we can use the 68th centile value of reported error in glacier surface velocity magnitude, and normalized ice thickness standard errors for each glacier, which can be expressed by Equation (11) where V s is the mean glacier surface velocity, σ v is the standard inaccuracy of the glacier surface velocity magnitude, H is the mean ice thickness for each glacier, σ H is the standard error of ice thickness for each glacier.

ELAs and AARs Calculation
ELA is the altitude where the mean glacier mass balance is zero. Hence, it divides glaciers into ablation and accumulation zones, which is significant for glaciologists' investigations of glacier health [76]. AAR represents the accumulation area ratio of each glacier and is determined by the ELA. Here, we used three stages to retrieve ELAs across the whole WK. First, on the basis of glacier SMB results and ASTER DEM, we used possible ELA to classify pixels within each glacier into ablation and accumulation zones. The resultant possible ELA divided each glacier into predicated accumulation and ablation regions. Then, we moved the possible ELA across the whole elevation range of each glacier, at a 25 m altitudinal distance with each movement, and determined the optimal classification outcomes by using the best-performing Dice coefficient results [77] according to Equation (12): where DC represents the Dice coefficient, N SMB is the total pixel count of the SMB map, N DEM is the number of total pixels in the DEM map, N c is the true prediction of number of pixels in the glacier accumulation zone, which can be calculated by the number of pixels that represent positive values in both the DEM and SMB maps. N a is the true prediction of the number of pixels in the glacier ablation zone, which can be described as the negative value pixel count in both the DEM and SMB maps. Finally, for those glaciers with optimal ELAs at either end of the glacier elevation range (glaciers in a complete ablation or accumulation state), we fitted a linear trend between the SMB area ratio and each altitude bin (25 m), to determine the theoretical climatic ELA. Here, for each glacier, the SMB area ratio was defined as the mean SMB at each elevation band multiplied by its area ratio. We applied the SMB area ratio technique for all glacier ELA generation, because of the unpredictable impacts of high ice flow and glacier terminal thickening during glacier ELA acquisition. After all the ELA-retrieving procedures were completed, AARs could be straightforwardly calculated for each glacier.

Glacial Health Index (GH Index)
In this study, we estimated the GH index for 2002-2020 for each glacier in WK, using Equation (13). We divided glacier health status into two groups and four subgroups. For the two groups, we classified glaciers in a healthy or unhealthy state as those with GH index > 0 or GH index < 0, respectively. With respect to the four subgroups, we first selected 0 < GH index < 10 or −10 < GH index < 0, and labelled these glaciers as low-quality healthy or sub-healthy. For those glaciers with GH index < −100 or GH index > 100, we called these superiorly healthy and unhealthy glaciers, respectively.
where ∆h represents glacier surface elevation change, ∆t is the time shift, AAR latter is the AAR of each glacier in the later time span (2011-2020), AAR f ormer is the AAR of each glacier in the earlier time span (2002-2011). OLS modelling is a common machine learning technique for estimating linear regressions equations with minimum squares error. In this work, we applied this method to investigate linear glacier SMB and ELA shift response to climatic forcing (summer skin temperature and IVT), as in previous studies [46,78]. Using the OLS model, calculation of glacier SMB response to climate change was carried out with Equation (14), and glacier ELA shift sensitivity to meteorological factors was calculated using Equation (15): where ∆b ∆t is glacier SMB (m w.e. a −1 ), ∆T is June-July-August (JJA) skin temperature annual trend (°C a −1 ), ∆P is IVT annual trend (kg m −1 s −1 a −1 ), db dT is glacier SMB sensitivity to JJA skin temperature (m w.e./°C), db dP is glacier SMB sensitivity to IVT (m w.e./kg m −1 s −1 ).

ANN Model
The ANN model utilizes deep learning techniques to research nonlinear glacier SMB and ELA response to climatic forcing (summer skin temperature and IVT), chiefly to observe how glacier SMB and ELA respond to extremely positive and negative climate change, and further to assist projecting glacier response to multiple future climate forcings [43].
Unlike Bolibar et al. [46], we applied only four layers (two hidden layers) since we did not consider glacial topography impact on SMB and ELA changes. For the two hidden layers, we assigned 20 neurons to each. For cross-validation, we utilized the LeavePOut method (p = 10) to validate the training results. For the training and test datasets, we use mean squared error as loss function and determination coefficient R 2 as metrics.

Varied SMB Spatiotemporal Patterns
In this study, we estimated glacier SMB (>2 km 2 ) in WK (see Figure 3).  We found that during 2002-2011, glacier SMB in eastern WK was 0.05 ± 0.15 m w.e. a , while in western WK, it was 0.20 ± 0.11 m w.e. a , which showed that glaciers were thickening at a higher rate in western WK. From 2011 to 2020, the glacier SMB in western WK was −0.04 ± 0.08 m w.e. a , and in eastern WK it was more negative (−0.16 ± 0.17 m w.e. a ). During 2002-2020, glacier SMB in eastern WK was −0.06 ± 0.07 m w.e. a , while in western WK, it was 0.12 ± 0.06 m w.e. a (see Table 2).

Periods
Western WK Eastern WK WK Figure 3. Glacier SMB spatial pattern differences in western WK (in the left three subplots) and eastern WK (in the right three subplots) in the past two decades.
We found that during 2002-2011, glacier SMB in eastern WK was 0.05 ± 0.15 m w.e. a −1 , while in western WK, it was 0.20 ± 0.11 m w.e. a −1 , which showed that glaciers were thickening at a higher rate in western WK. From 2011 to 2020, the glacier SMB in western WK was −0.04 ± 0.08 m w.e. a −1 , and in eastern WK it was more negative (−0.16 ± 0.17 m w.e. a −1 ). During 2002-2020, glacier SMB in eastern WK was −0.06 ± 0.07 m w.e. a −1 , while in western WK, it was 0.12 ± 0.06 m w.e. a −1 (see Table 2). Moreover, according to RGI v6.0, in WK there are 471 and 84 glaciers on the northern and southern slopes, respectively. In western WK, 238 glaciers are located on the northern slope, and 26 glaciers on the southern slope. In eastern WK, there are 233 glaciers on the northern slope and 58 glaciers on the southern slope. We found that in the past two decades, glaciers on the northern slope in eastern WK displayed regular trends of increased gain or reduced loss of ice mass, which is consistent with outcomes obtained by Zhang et al. [79]. However, this trend was reversed in Western WK, except during 2002-2011 (see Table 3). Overall, glaciers showed different SMB variations between western and eastern WK. Glaciers in western WK predominantly underwent greater gain or reduced loss of ice mass compared with eastern WK. Additionally, we found that there were dissimilar slope effects on SMB distribution in western and eastern WK. The northern slope in western WK showed the greatest loss of mass in the last two decades, while in eastern WK, the southern slope lost more ice mass.
Gardner et al. [80] estimated that the average rate of elevation change was 0.17 ± 0.15 m a −1 in WK during 2003-2009, similar to our results (0.18 ± 0.10 m a −1 ). Zhang et al. [79] found that the glacier elevation change was 0.002 ± 0.003 m a −1 in eastern WK during 2000-2018, which also fits our results (−0.01 ± 0.03 m a −1 ). Shean et al. [81] calculated that glacier mass balance in WK was slightly positive (0.01 ± 0.14 m a −1 ) during 2000-2018; the result in ours study was 0.06 ± 0.04 m a −1 during 2002-2020. We also compared glacier elevation change outcomes during 2002-2011 and 2011-2020 with those of Hugonnet et al. [2], which showed only a minor difference, with a maximum dissimilarity of 0.04 m a −1 during 2011-2020. Brun et al. [11] estimated that the glacier mass balance in WK was 0.14 ± 0.08 m w.e. a −1 during 2000-2016, which represents a large difference compared with our results but does overlap within uncertainty. There are two reasons for this difference. First, our time span was 2002-2020, which differed from that studied by Brun et al. [11]. Changes in glacial elevation described by Hugonnet et al. [2] confirmed that glaciers can experience remarkable mass gain or loss with a five-year period. Second, we used the continuity equation rather than the rationally estimated ice density value (850 kg m −3 ) (see Table 4).

ELAs and AARs Outcomes
Based on the aforementioned SMB results for three different time spans, we obtained glacier ELAs and AARs across the whole WK. We found that the area-weighted mean and median ELAs in WK from 2002-2020 were 5786 m and 5535 m, respectively, and their respective AARs were 0.58 and 0.78. In HMA, the area-weighted mean and median ELAs and AARs were 5283 m and 5349 m, and 0.51 and 0.44, respectively [23]. This showed that glacier ELAs and AARs (area-weighted mean and medium) in WK were larger than the average levels in HMA (see Figure 5).

ELAs and AARs Outcomes
Based on the aforementioned SMB results for three different time spans, we obtained glacier ELAs and AARs across the whole WK. We found that the area-weighted mean and median ELAs in WK from 2002-2020 were 5786 m and 5535 m, respectively, and their respective AARs were 0.58 and 0.78. In HMA, the area-weighted mean and median ELAs and AARs were 5283 m and 5349 m, and 0.51 and 0.44, respectively [23]. This showed that glacier ELAs and AARs (area-weighted mean and medium) in WK were larger than the average levels in HMA (see Figure 5).  respectively. We found that ELAs underwent a slight decline and AARs saw a slight increase from the earlier subperiod to the latter, which was inconsistent with the corresponding SMB variations (see Table 5). Owing to the lack of in situ measurements for WK, we were reliant on previous studies to verify our ELA and AAR outcomes [23]. We found that dissimilarities between our results and others were slightly substantial, at 322 m for ELAs and 0.19 for AARs. We assumed that these results were mainly due to different time spans, spatial coverage dissimilarities, glacier samples, and slight differences in the ELA retrieval methods.

Skin Temperature and IVT Results
To assess climatic forcing in WK, we extracted the mean annual summer skin temperature (JJA skin temperature), winter skin temperature (December-January-February skin temperature, DJF skin temperature), mean annual IVT, and corresponding summer IVT (JJA IVT) and winter IVT (DJF IVT) trends, for 2002-2011, 2011-2020, and 2002-2020 within the glaciated area [82,83] (the significance and correlation per pixel of climatic factors are described in the Supplementary Information Discussion). Sakai and Fujita [78] assumed that the annual range of monthly temperature (annual temperature) and the ratio of summer IVT have close interplay with glacier mass balance in HMA (Supplementary Information Discussion). However, in this study, we found that the summer IVT ratio and annual skin temperature were less significant than the annual IVT and summer skin temperature when depicting glacier responses to meteorological variations.
Generally speaking, in WK, JJA skin temperature has been increasing, while DJF skin temperature decreased during 2002-2011, 2011-2020, and 2002-2020. In the two subperiods, 2002-2011 and 2011-2020, decreasing JJA skin temperature trends could be observed (from 0.04°C a −1 to 0.02°C a −1 ), and the DJF skin temperature showed the same tendency (from −0.15°C a −1 to −0.25°C a −1 ) (see Table 6). Table 6. June-July-August (JJA) and December-January-February (DJF) skin temperature (unit:°C a −1 ) trends in western WK, eastern WK and WK overall, over three time periods. Additionally, during the full period (2002-2020), we found that the JJA skin temperature increased by 0.05°C a −1 in western WK, while in eastern WK it increased more obviously (0.08°C a −1 ). From 2002-2011, the JJA skin temperature increased by 0.01°C a −1 in western WK, while in eastern WK it increased to 0.07°C a −1 . During 2011-2020, the JJA skin temperature in western WK maintained the same increasing pace (0.01°C a −1 ) as the earlier subperiod, while in eastern WK its rate of increase declined (0.04°C a −1 ) (see Figure 6).

Remote Sens. 2022, 14, x FOR PEER REVIEW 16 of 27
a ) as the earlier subperiod, while in eastern WK its rate of increase declined (0.04 ℃ a ) (see Figure 6). Moreover, during 2002-2020, the DJF skin temperature in western WK dropped by 0.09 ℃ a , and in eastern WK it declined by 0.01 ℃ a . During the earlier subperiod, the DJF skin temperature in western WK decreased by 0.33 ℃ a , while in eastern WK Moreover, during 2002-2020, the DJF skin temperature in western WK dropped by 0.09°C a −1 , and in eastern WK it declined by 0.01°C a −1 . During the earlier subperiod, the DJF skin temperature in western WK decreased by 0.33 • C a −1 , while in eastern WK it increased by 0.02°C a −1 . During the later subperiod, the DJF skin temperature decreased by 0.38°C a −1 in western WK, while in eastern WK the decrease was −0.15°C a −1 , showing a remarkable difference between western and eastern WK (see Figure 7).  Overall, JJA and DJF skin temperatures in western and eastern WK had similar trends during 2011-2020 and 2002-2020, except 2002-2011. Furthermore, across the whole WK, the JJA skin temperature was rising, and the DJF skin temperature was declining.
Additionally, we extracted the annual mean IVT inclination in WK for the last two decades (see Figure 8). The results showed that although IVT inclination in WK continued to rise, the rate of increase was declining. IVT increased at 0. 19  In two subregions, western WK and eastern WK, there was also an indication that the rate of increase in the annual IVT trend declined from the earlier subperiod (2002-2011) to the later subperiod. In addition, there were almost no differences in the annual IVT trends between western and eastern WK in the three periods (2002-2011, 2011-2020, In two subregions, western WK and eastern WK, there was also an indication that the rate of increase in the annual IVT trend declined from the earlier subperiod (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011) to the later subperiod. In addition, there were almost no differences in the annual IVT trends between western and eastern WK in the three periods (2002-2011, 2011-2020, and 2002-2020) (see Table 7). Furthermore, we also researched the JJA and DJF IVT trends in WK for the past two decades (see Table 8) (correlation coefficient R and significance p value per pixel are described in Supplementary Information Figures). Overall, we found no obvious spatial differences between western and eastern WK for 2002-2011, 2011-2020, and 2002-2020. Nevertheless, we discovered obvious differences in the JJA and DJF IVT trends between the three periods (2002-2011, 2011-2020, and 2002-2020). Specifically, we found a remarkable increasing JJA IVT trend (0.55 kg m −1 s −1 a −1 ) during 2011-2020, which was apparently higher than the periods 2002-2011 and 2002-2020. In addition, we found that the DJF IVT was decreasing, except during 2011-2020 when it rose by 0.02 kg m −1 s −1 a −1 .

Glacial Health in WK
ELAs or AARs can be used to denote the accumulation zone over the glacier surface; thus, they can be used to estimate glacier health in WK. Miles et al. [23] demonstrated that glaciers in WK exhibit much healthier states compared with other HMA regions. Nonetheless, due to limited glacier sample size and the lack of multitemporal comparisons, evidence for glacier well-being in WK remains scarce.
In western WK, we estimate that 69% of glaciers are in a healthy state (regarding glacier healthy state classification here and below, see Section 2.6.2). The remaining unhealthy glaciers show no remarkable thinning pattern, indicating that they are unable to supply much downstream water. In eastern WK, 62% of glaciers are in a healthy state, which demonstrates that glaciers in western WK are healthier and more stable than those in eastern WK. Furthermore, low-quality healthy glaciers account for 84.15% of total glaciers in western WK, while in eastern WK, this figure decreases to 74.72%. This result shows that, when accounting for global warming, western WK exceeds eastern WK in its stores of potential downstream water. We also researched sub-healthy glaciers in WK. The outcomes show that there are 179 glaciers in a state of subhealth in WK. Nonetheless, in contrast to the varied spatial distribution pattern of low-quality healthy glaciers between western WK and eastern WK, the allocation of sub-healthy glaciers is dispersed relatively evenly between the two subregions. The results show that 85 sub-healthy glaciers are located in western WK and 94 glaciers are located in eastern WK.
Additionally, 34 glaciers in WK currently show a superiorly healthy state. The distribution of these glaciers is unevenly dispersed between western and eastern WK; 10 glaciers are located in western WK, and 24 glaciers are located in eastern WK. Furthermore, 15 glaciers are in a supremely unhealthy state in WK, and the number of worst-state glaciers is 14 in eastern WK and one in western WK. Overall, glaciers in eastern WK generally have a more detrimental status and have been the primary downstream water providers for mountain dwellers in the last two decades. However, it will be necessary to pay careful attention to glacier variations in western WK in the upcoming decades, as they hold a substantial quantity of potential meltwater and are more likely to be affected by global warming and become natural hazards [72,[84][85][86].
Brun et al. [87] concluded that mean glacier elevation can explain glacier SMB variability for most regions in HMA to some extent (>8%). However, we found no remarkable altitude characteristics in WK that can explain why some glaciers are healthy and others are not. However, we recognize that glacier aspect has some connections with dispersion of glacier well-being in our study area. The results show that 80% of healthy glaciers are on the northern slope, while for unhealthy glaciers, this ratio decreases to 70%. This may contribute to topographic shading effects, altering the potential clear-sky solar radiation received by glaciers [88,89].

SMB Response to Climate Forcing
Che et al. [40] concluded that the ratios of glacier mass balance to summer skin temperature and precipitation are nonlinear. Wang et al. [90] assumed that sensitivity of glacier mass balance to summer skin temperature is nonlinear, and that the relationship between glacier mass balance and precipitation is linear. Sakai and Fujita [78] assumed a linear mass balance response to both of these climatic factors. However, via deep learning technology, Bolibar et al. [43] suggested that the relationship between glacier mass change and meteorological elements is nonlinear. We compared linear and nonlinear relationships between glacier SMB and climate forcing, using an OLS multilinear regression model and an ANN model, respectively (see Supplementary Information Discussion). Conversely, for extreme positive summer skin temperature changes, the ANN model demonstrated larger glacier SMB sensitivity than the OLS model. For extreme negative summer skin temperature changes, glacier SMB sensitivity was smaller indicated by the ANN model. The results we revealed for the periods 2011-2020 and 2002-2020 match well with Bolibar et al. [43], except for the period 2002-2011.
Generally, both the OLS and ANN models depicted that WK glaciers are more sensitive to IVT than to summer skin temperature in the long run, which may be attributed to the high altitudes of WK glaciers [53,83,90,91]. However, in the two subperiods, glaciers were more sensitive to summer skin temperature.
During 2002-2011, glacier SMB sensitivity to summer skin temperature was −0.38 m w.e./°C; with a 1°C increase in summer skin temperature, glacier SMB would decrease 0.38 m w.e. During the period 2011-2020, glacier SMB sensitivity to summer skin temperature was larger (−0.16 m w.e./°C). However, during the full period (2002-2020), glacier SMB sensitivity to summer skin temperature was −0.13 m w.e./°C, which was less than in the period 2002-2011. We found that glacier SMB sensitivity to summer skin temperature during 2011-2020 and 2002-2020 was substantially smaller than in other HMA regions [92][93][94][95][96][97]. We believe our results are feasible, because previous studies have shown that glacier SMB sensitivity to summer skin temperature varies in different regions and, generally speaking, continental glaciers are the least sensitive to climate change [90]. Specifically, Ebrahimi and Marshall [98] found that glacier mass balance sensitivity to summer skin temperature was −0.236 m w.e./°C in the Canadian Rocky Mountains, which is considerably different from glacier sensitivity to summer skin temperature in Switzerland (−0.63∼−0.72 m w.e./°C) [99]. This shows that glaciers in different regions can display different climatic sensitivities, as they are surrounded by varied environmental conditions. Furthermore, glaciers in WK are located at comparatively higher altitudes than other HMA regions, and they are considered continental glaciers; thus, it is reasonable to suggest that their sensitivity to summer skin temperature is smaller than that of glaciers in other HMA regions.
Furthermore, glacier SMB sensitivity to IVT was found to be similar to the summer skin temperature sensitivity pattern. From 2002-2011, glacier SMB sensitivity to IVT was +0.22 m w.e./kg m −1 s −1 . Combined with the annual mean IVT during this period, this is equivalent to 0.02 m w.e. (50%) −1 . This result was smaller than that for its neighbouring region, Karakoram (0.02 m w.e. (10%) −1 ), and smaller than subcontinental glacier SMB sensitivity to precipitation (0.09 m w.e. (10%) −1 ) [90]. We believe that a primary reason for this apparent discrepancy is the dissimilar estimation of glacier SMB (Supplementary Information Discussion), as well as the differences between IVT and precipitation [100]. During 2011-2020, glacier SMB sensitivity to IVT was 0.07 m w.e./kg m −1 s −1 , equal to approximately 0.01 m w.e. (50%) −1 . This outcome demonstrated that during the later subperiod, glaciers in WK were less sensitive to IVT than in the earlier subperiod. Overall, during the full time period (2002-2020), glacier SMB sensitivity to IVT was +0.45 m w.e./kg m −1 s −1 , which is approximately 0.02 m w.e. (50%) −1 .

ELAs Response to Climate Forcing
Sagredo et al. [27] showed that glacier ELA shifts are linearly related to summer skin temperature and precipitation variations. Conversely, motivated by the nonlinear glacier SMB response [43], we researched the nonlinear relationship between glacier ELA shifts and climate forcing (see Supplementary Materials). The result showed that glacier ELA was almost linear with climate change. However, the ANN model captured a significant glacier ELA response to IVT. With the increment of IVT exceeding 0.7 kg m −1 s −1 , the glacier ELA response was less sensitive, and its shift was smaller than the value displayed by the OLS model, indicating that increasing glacier IVT cannot effectively mitigate the impacts of increasing summer skin temperature on variations in glacier ELA. For the following analysis, we exclusively utilized the linear glacier ELA response.
We found that glacier ELA sensitivity was mainly driven by summer skin temperature (+154.66 m/°C) rather than IVT (−121.57 m/kg m −1 s −1 ). Zhang et al. [101] found that HMA glaciers are more sensitive to summer skin temperature than to precipitation. Wang et al. [90] concluded that continental glacier ELAs are relatively insensitive to skin temperature change, in contrast with maritime glaciers. The ELA sensitivity of maritime glaciers in HMA is generally +200∼300 m/°C, and ELA can sometimes be lifted above the maximum elevation with a merely 1 K increase of summer skin temperature. Furthermore, glacier ELA sensitivity to summer skin temperature varies globally. Noël et al. [102] found that the ELA sensitivity of Svalbard glaciers to summer skin temperature was +129.66 m/°C. Sagredo et al. [27] found that Andes glacier ELAs sensitivity to summer skin temperature ranged from +140∼230 m/°C. Réveillet et al. [103] found that the ELA sensitivity of the Zongo glacier to summer skin temperature would respond by +120∼180 m/°C under different climatic model scenarios. Based on our results, with a 1°C increase in summer skin temperature, glacier ELAs in WK would increase 143.76 m, which is reasonable and rational.
Furthermore, Bolibar et al. [43,46] have confirmed that glacier mass change is connected with glacier topography and climatic change. Thus, we can conclude that glacier ELA shift is also connected with these two parameters. However, in the long run, glacier topographic change is comparatively small compared with climatic change, thus, in the current study, we have not accounted for its impact on glacier ELA shift. Based on our retrieved glacier SMB and ELA results, we found that in eastern WK average glacier SMB changed from positive to negative from 2002-2011 to 2011-2020, while glacier ELA decreased in the same period, which is an abnormal phenomenon. We believe that this may be connected with ice flow on the glacier surface. Specifically, if we presume that glaciers are stationary and we ignore the snow drift effect, then the ELA would rise with increasing summer skin temperature. However, the glacier surface velocity in WK is remarkable [47,51,104]; therefore, the ice mass in the accumulation zone would be transferred to lower elevations, which would decrease the ELA. We found that in WK there are 99 glaciers showing negative SMBs and decreasing ELAs. Meanwhile, ice velocity was apparent, and ice flux was positive close to ELA in the accumulation zones of these glaciers. This evidence suggests that the ice mass in the accumulation zone, driven by glacier dynamics, results in the abnormal phenomenon of decreasing ELA. Additionally, this could explain why glacier SMB was negative in the later subperiod (2011-2020), while its corresponding ELA decreased.

Conclusions
In this study, we calculated glacier SMB in WK via the continuity equation, for the periods 2002-2011, 2011-2020, and 2002-2020. The results show that in the former time span, glaciers were in a state of mass gain (0.10 ± 0.14 m w.e. a −1 ), while during the later subperiod, glaciers were in a state of mass loss (−0.12 ± 0.14 m w.e. a −1 ). Overall, in the last two decades, glacier SMB in WK has been slightly negative (−0.01 ± 0.07 m w.e. a −1 ). We estimated glacier ELAs and AARs in WK and its two subregions, western and eastern WK, for 2002-2020. We found that ELAs slightly declined from the earlier time span (5788 m) to the later (5744 m), in contrast to SMB variations. We believe that this abnormal phenomenon was due to increasing glacier surface velocities, especially during 2011-2020. We found that glacier surface velocity was apparent in the vicinity of ELAs, and ice flux in the neighboring ELA zones was positive, which could lead to snow or ice slipping over the ELA and amassing in the upper position of the ablation zone. Hence, ELA would decrease, even with increasing summer skin temperature. Furthermore, ELAs in western WK were lower than those in eastern WK, because of the higher median glacier altitude in eastern WK.
Additionally, we found that glaciers in western WK are healthier than those in eastern WK, and that necessary monitoring measurements are still indispensable to safeguard substantially vulnerable low-quality healthy glaciers. Eastern WK glaciers are in a worse state of well-being than western WK glaciers. Furthermore, we established a quantitative relationship between glacier SMB and meteorological elements, summer skin temperature and IVT. The results showed that SMB was mainly affected by summer skin temperature during the two subperiods (2002-2011 and 2011-2020). However, in the full period (2002-2020), glacier SMB was more prone to IVT, which can explain why glacier SMB in WK is less negative than other HMA regions in the last two decades. Overall, glacier SMB was more sensitive to extreme positive summer skin temperatures and less sensitive to extreme negative summer skin temperatures.
Finally, we also observed a quantitative relationship between glacier ELA shifts and the climatic factors of summer skin temperature and IVT. Glacial ELA variations were found to be more sensitive to summer skin temperature (+154.66 m/°C) than to IVT (−121.57 m/kg m −1 s −1 ). The ELA shift sensitivity equation has relatively large uncertainty, which could be associated with glacier dynamics. To conclude, the WK glacier ELA is less sensitive to an extremely positive IVT.   Table S1: The number of ASTER images for different years from 2001-2005. Glacier counts and area (>2 km 2 ) were shielded by ASTER images in each year during this period, Table S2: The number of ASTER imagers for different years from 2009-2013. Glacier counts and area (>2 km 2 ) were shielded by ASTER images in each year during this period, Table S3: The number of ASTER images for different years from 2018-2021. Glacier counts and area (>2 km 2 ) were shielded by ASTER images in each year during this period, Table S4: Correlation coefficients (R) and significances (p-value) between glacier mass balance and climatic factors (JJA surface temperature and annual IVT) for three periods (2002-2011, 2011-2020, and 2002-2020), Table S5: Correlation coefficient (R) and significance (p-value) between glacier ELAs shifts and climatic factors (JJA surface temperature and annual IVT) during the period 2002-2020. Data Availability Statement: The ASTER imagery was obtained from Earthdata (https://earthdata. nasa.gov/, accessed on 23 April 2022). The glacier velocity dataset was provided by NASA MEaSUREs ITS_LIVE project (https://its-live.jpl.nasa.gov/, accessed on 23 April 2022). the glacier thickness dataset is available at https://doi.org/10.3929/ethz-b-000315707, accessed on 23 April 2022. The RGI v6.0 dataset is available at https://www.glims.org/RGI/rgi60_dl.html, accessed on 23 April 2022. The supraglacial debris extents are available at https://zenodo.org/record/3866466#.YYOuu21By00, accessed on 23 April 2022. The HydroLAKES v1.0 dataset is available at https://www.hydrosheds. org/page/hydrolakes, accessed on 23 April 2022. The glacial lakes dataset is available at https: //nsidc.org/data/HMA_GLI/versions/1, accessed on 23 April 2022. The ERA 5.1 reanalysis datasets were provided by EVMWF and can be downloaded at this website: https://www.ecmwf.int/en/ forecasts/datasets/reanalysis-datasets/era5, accessed on 23 April 2022. The DEM co-registration code is from Github: https://github.com/dshean/demcoreg, accessed on 23 April 2022.