Influence of Topographic Shading on the Mass Balance of the High Mountain Asia Glaciers

Abstract

Most studies attribute the glacier mass balance within High Mountain Asia (HMA) to climate change, ignoring the influence of its complex terrain. Knowledge of the influence of this complex terrain is crucial for understanding the spatial variability in its mass balance. However, there is a lack of any systematic assessment of this influence across HMA. Therefore, in this study, we used the glacier outlines and raster data (SRTM DEM, slope and aspect) to calculate the topographic shading of all 97,965 glaciers within HMA during the ablation period, which is regarded as a major index of the influence of complex terrain on the mass balance. The results showed that 27.19% of HMA glacier area was subjected to topographic shading, and regional differences were significant with respect to both their altitudinal and spatial distributions. The topographic shading contributed to the protection of the smallest glaciers from solar illumination. Furthermore, we found a significant correlation between the topographic shading and mass balance in these small north-facing glaciers. However, these small glaciers were most prevalent in the north-facing orientation, especially in West Kunlun, East Kunlun, Inner Tibet Plateau and Qilian Shan, where shading was found to increase with decreases in the glacier area. This indicates that complex terrain can affect the spatial distribution of the mass balance by altering the solar illumination pattern.

1. Introduction

High Mountain Asia (HMA) is known as the Asian water tower [1], as it constitutes the largest concentration of glaciers and ice caps outside the polar regions, accounting for approximately 30% of the total area of all mountain glaciers on Earth [2,3]. The total volume of the glaciers in HMA is 6.14 × 10³ km³, providing meltwater supply to rivers that meet the requirements of nearly a quarter of a billion people [4], as well as agriculture on the Indo-Gangetic plain [5]. Due to ongoing climate change, most of the glaciers in HMA have suffered remarkable mass loss, becoming vulnerable and unsustainable due to their higher ablation than accumulation in these regions [6]. The mass balance spatial heterogeneity loss has been obvious across HMA in recent decades [7,8,9,10,11,12], with extremely negative mass balance leading to marked glacier thinning in Tian Shan and Hengduan Shan. The exception is the Karakoram, with a slightly positive mass balance having been observed [13]. Nevertheless, the current climate factors do not fully explain the mass balance differences, such as the Karakoram Glacier Anomaly [6], as numerous gaps remain in our knowledge of regional glacier behavior and the interactions of glaciers with complex topography. To improve our understanding of the relationship between mass balance and topography, it is crucial to carry out mass balance simulations and water resource estimates.

Although the mass balance within HMA is a combined climate response derived from numerous influencing factors [14,15,16,17,18,19], climate change is still a primary driver [14]. Spatial heterogeneity of the mass balance within HMA is driven by different atmospheric circulation patterns, which transform the annual distribution and source of precipitation [8]. Moreover, glaciers under these atmospheric circulation patterns can be classified as maritime, continental and subcontinental glacier types, with each having differential climate sensitivities [20,21]. However, the problem still exists that climate factors are not able to explain the difference in mass balance under the same climate conditions at a smaller scale. It is evident that the mass balance for north-facing glaciers is less than that of those facing south [19], and the same is true for the glacier area change [22]. It has been pointed out that the currently stable regional mass balance of Karakoram Anomaly glaciers might be attributable to topographic factors [6]. Within HMA, which is often referred to as the Third Pole [8], glaciers exist where cold climate, high altitude and complex terrain are the main environmental characteristics, all of which promote the glaciers’ evolution.

The effect of complex topography on glacier mass balance reflects the influence of topographic shading. By obstructing solar illumination, topographic shading exerts a direct impact on solar radiation reaching the glacier surface and has a significant effect on glacier mass balance [23]. Consequently, most studies take topographic shading into account in the distributed mass balance to improve the simulation accuracy [24,25,26,27]. As well as providing topographic shading data for models, these data also help to interpret glacier changes. A previous study found that the change in area of some small glaciers was not significant over their monitoring period, due to the topographic shading’s obstruction of direct solar radiation [28]. York (2013) found obvious differences in the mass balance for south- and north-facing glaciers [29], which were likely attributable to differences in the glacier aspects and solar illumination, with south-facing, higher solar illumination glaciers losing the most mass and vice versa for north-facing glaciers. Radiation and topographic shading are more important than wind in affecting glacier aspects and particularly affect the mass balance in accumulation areas [30]. It can be seen that topographic shading is not only important for a distributed mass balance model but also for understanding glacier changes.

Several topographic shading methods have been proposed. Mcdonnell et al. (1998) carried out the hill shading calculations in ArcGIS by determining the hill shading value of the center cell of a moving window using five surrounding grid units [31], ignoring the influence of the surrounding topographic grid cells that were further away. Some studies have considered the distant terrain obscuring effect and applied the ray tracing method to calculate topographic shading, which requires a digital elevation model (DEM) covering the glacier to be large enough so that the ray can trace the edge of the DEM to improve the calculation accuracy [24,32,33,34]. This approach is widely used to constrain the solar radiation over rugged terrain; hence, it is used in distributed mass balance models [26,27,35,36,37]. Topographic shading is particularly important for understanding glacier differences and model applications. However, systematic estimates of the relationship between topographic shading and mass balance remain lacking within HMA.

Here, we used the Shuttle Radar Topography Mission (SRTM) DEM data with a resolution of 90 m that covered all the glaciers within HMA, together with data on the glacier slope and aspect. We calculated the grid topographic shading by ray tracing according to the location between the grid and solar beam. The main aim of this study was to calculate the topographic shading of 97,965 glaciers included in the Randolph Glacier Inventory (RGI) version 6.0 on each pixel scale. Our goal was to provide a systematic estimate of the relationship between mass balance and topographic shading across HMA.

2. Materials and Methods

2.1. Study Area

HMA includes 16 mountain ranges and covers many countries in Asia, spanning from the Altai and Tian Shan Mountains to the Tibetan Plateau (TP) (Figure 1). Due to its high elevation and cold climate, this region holds the largest area of glaciers outside the polar regions [8]. According to the global Randolph Glacier Inventory version 6.0 [38], there are 97,965 glaciers within HMA, covering a total area of 98,752.33 km² (

Table 1. Summary of the glacier characteristics statistics for different regions within HMA.

RGI Region	Nglacier	Anorth (km2)	Aeast (km2)	Asouth (km2)	West (km2)	A (km2)
01 Altay and Sayan	2431	551.58	325.48	132.75	151.77	1161.57
02 East Tianshan	5227	1373.07	665.27	292.97	521.88	2853.18
03 West Tianshan	9739	3869.71	2048.93	1555.03	2055.02	9528.69
04 Hissar Alay	3151	804.03	480.33	217.39	343.69	1845.43
05 Pamir	10,233	4625.88	2520.35	1307.45	1781.13	10,234.80
06 West Kunlun	5397	3379.15	1946.58	1064.66	1759.46	8149.85
07 East Kunlun	3519	1448.55	717.18	507.55	577.28	3250.56
08 Qilian Shan	2730	869.49	337.21	144.58	285.83	1637.10
09 Karakoram	13,759	8285.10	5807.44	4308.07	4457.47	22,858.08
10 Inner Tibet Plateau	9365	3228.15	2122.98	1198.37	1382.54	7932.05
11 South and East Tibet	5063	1283.68	1083.66	730.77	772.95	3871.05
12 Hindu Kush	4401	1257.87	702.68	415.83	563.01	2939.39
13 West Himalaya	9829	3092.48	1844.76	1160.46	1627.41	7725.10
14 Central Himalaya	4530	1786.11	1372.67	1165.31	1168.30	5492.39
15 East Himalaya	4237	1399.29	1277.93	1165.56	1045.61	4888.39
16 Hengduan Shan	4354	1380.20	1274.30	884.10	846.10	4384.70
Total	97,965	38,634.33	24,527.73	16,250.83	19,339.45	98,752.33

The number of glaciers (N) and their corresponding areas in each aspect (A). Orientations are classified as north (315–360°, 0–45°), east (45–135°), south (135–225°) and west (225–315°), with respective areas Anorth, Aeast, Asouth and Awest, and A is the total area of all orientations.

). The Karakoram has the largest area covered by glaciers (22,858.08 km²), accounting for 23% of the total HMA glacier area. The Altay and Sayan have the smallest glacier area, 1161.57 km², accounting for approximately 1% of the total area. The glacier area depends strongly on orientation in HMA. Nearly 40% of the glacier area is orientated north, of which 53% is in Qilian Shan and 28% is in the East Himalayas. The north-facing aspects are most suitable for glacier development, while glaciers oriented south have the smallest areas, because the stronger solar illumination hinders glacier development. South-facing glaciers account for only 15% of the total glacier area. There is no marked difference in glacier areas between east and west orientations, which account for 24% and 19%, respectively.

The rugged landscape affects the glacier development and distribution within HMA and also influences climate models applied to this region. Glaciers in HMA are primarily affected by atmospheric circulation patterns, including westerlies, the Indian monsoon and the East Asian monsoon [8]. The weather systems associated with these climatic patterns significantly affect glaciers by influencing the annual distribution of precipitation and temperature. Glaciers in different regions show different patterns of summer versus winter accumulation types [20,21]. The westerlies-dominated region receives abundant winter precipitation, such that glaciers in Western HMA (for example, the Tian Shan and Karakoram) are winter accumulation-type glaciers [39]. The Indian and East Asian monsoons strongly influence the Southwestern HMA, and summer precipitation is dominant in this region. Mountains block the transfer of most moisture to the TP, with precipitation decreasing sharply away from the plateau margin in both the monsoon and westerlies regions [18]. Consequently, most glaciers are east-facing in this region. Overall, both topography and climate affect the glacier distribution.

2.2. Data

The glacier outlines within HMA were taken from the RGI 6.0 [38]. Data for most glaciers in China were extracted from the second Chinese Glacier Inventory (CGI), covering the period 2006–2010 [40]. Glaciers outside China were based on RGI 6.0, which is based on satellite imagery for 1999–2010 [41]. Raster data for each glacier were established using the DEMs, aspect, slope and glacier outlines.

Glacier surface topography was extracted from the SRTM DEM version 4 [42] for the region 65–110° E and 25–55° N, fully covering all the glaciers in HMA. The DEM was transformed from geographical coordinates to Asia North Lambert Conformal Conic projection coordinates with a spatial resolution of 90 m. These transformed data were then used to calculate the aspect and slope.

The glacier-wide geodetic mass balance data were calculated from the rate of elevation change extracted from a linear fit of multitemporal DEMs obtained from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) [7], which covers most glaciers within HMA, except Altay and Sayan. The individual glacier mass balances have a median uncertainty of ±0.22 m water equivalent (w.e.)/a, ranging from a minimum of ±0.14 to a maximum of ±0.89 m w.e./a [7]. Only glaciers larger than 2 km² were included to limit the error for the glacier mass balance, which decreases with the increasing glacier area [43]. Since small glaciers have a larger uncertainty for mass balance, we calculated the correlation between topographic shading and mass balance to maximize the information available for the mass balances of small glaciers. The glacier mass balances were calculated using the RGI 6.0 glacier mask to match the topographic shading raster.

2.3. Topographic Shading Calculation

For the topographic shading calculations, we first determined the solar ray location of each glacier grid cell. The solar ray location mainly involves calculating the solar altitude angle and solar azimuth of each pixel for individual glaciers within HMA using the following equations [44,45,46]:

$$\begin{array}{c}\alpha =asin(sin\left(\Phi \ast \frac{\pi }{180}\right)\ast sin\left(\delta \ast \frac{\pi }{180}\right)+cos\left(\Phi \ast \frac{\pi }{180}\right)\ast cos\left(\delta \ast \frac{\pi }{180}\right)\ast \\ cos\left(h\ast \frac{\pi }{180}\right))\end{array}$$

(1)

$$\beta =acos\left(\frac{sin\left(\Phi \ast \frac{\pi }{180}\right)\ast sin\left(\alpha \right)-sin\left(\delta \ast \frac{\pi }{180}\right)}{cos\left(\Phi \ast \frac{\pi }{180}\right)\ast cos\left(\alpha \right)}\right)$$

(2)

where α denotes the solar altitude angle, β is the solar azimuth, Φ is the geographical latitude (positive value in the Northern Hemisphere and negative in the Southern Hemisphere), δ is the solar declination and h is the hour angle.

2.3.1. Solar Declination

The solar declination (δ) is the angle between the sun’s rays and the equatorial plane, which is a function of the day of the year:

$$\delta =0.3948+{\delta }_1+{\delta }_2+{\delta }_3$$

(3)

$${\delta }_1=-23.2559\ast cos\left(2\pi \ast \frac{t_j}{365}+0.1582\right)$$

(4)

$${\delta }_2=-0.3915\ast cos\left(4\pi \ast \frac{t_j}{365}+0.0934 \right)$$

(5)

$${\delta }_3=-0.1764\ast cos\left(6\pi \ast \frac{t_j}{365}+0.4539 \right)$$

(6)

where t_j is the Julian date (number of day of the year).

2.3.2. Hour Angle

The hour angle is the angle through which the Earth must turn to bring the meridian of the site directly under the sun. This angle is a function of the time of day:

$$h=15\ast \left(t-12\right)$$

(7)

where t is the local apparent solar time.

In astronomy, time is divided into reference solar time and local apparent solar time. For glaciers within HMA, we referred to Beijing time (longitude 120°), which is not the local apparent solar time. When the solar ray location is calculated, the time must be the local apparent solar time. Therefore, we translated Beijing time to local apparent solar time using the following equations:

$$t=MT+\frac{lon-lo{n}_{ref}}{15}+\Delta T$$

(8)

where MT indicates the reference time and is assumed to be the Beijing time, lon is the grid pixel longitude, lon_ref is 120° and ΔT is expressed by:

$$\mathsf{\Delta }T=-0.128\ast \sin \left(\left(\theta -2.8\right)\ast \frac{\pi }{180}\right)-0.165\ast sin\left(\left(2\ast \theta +19.7\right)\ast \frac{\pi }{180}\right)$$

(9)

$$\theta =\frac{t_j}{365.25}\ast 360$$

(10)

2.3.3. Topographic Shading of Glaciers

Glacier surface topographic shading occurs when the sun is obscured by topography during the day. Here, the solar illumination time is regarded as the period when the sun rises above the horizon in the morning and sinks below the horizon in the evening. This period changes with space and time. In this study, we calculated the solar altitude angle and solar azimuth at every hour during the ablation period (from June to September), as this is when the shading has a more significant influence on glacier melting compared to any other time period. We also restricted our focus to this period because glacier shading calculations are computationally expensive when considering a spatial resolution of 90 m for each hour from June to September within HMA.

Glacier surface topographic shading was calculated following the work flow shown in Figure 2 using the DEM and glacier outlines dataset. This process can be divided into three key steps as follows:

• Data preprocessing

To obtain the glacier slope and aspect data, we transformed the DEM data from geographical coordinates to Asia North Lambert Conformal Conic projection coordinates, which is the projection data used to calculate the slope and aspect of glaciers within HMA. Since this process changed the raster data resolution, we translated all those data (slope, aspect and DEM) from the projection coordinates to geographical coordinates when considering the solar beam tracing for each grid in the third step, thereby ensuring a consistent resolution for the slope, aspect and DEM. These datasets were clipped by the extent of each glacier within HMA to form the basic dataset before the next step.

2.

Determine whether a grid is self-shading

For any glacier grid, we must consider whether the grid itself is obscuring the solar beam. If the slope of the glacier grid is greater than the solar altitude angle, when the solar azimuth is not within ±90° of the grid aspect range, then the glacier grid is obscured by its own grid. If not, the third step is required.

3.

Determine whether the solar beam is obscured by another DEM grid.

When tracing the glacier grid, it is necessary to determine a maximum tracking distance. In most cases, this maximum distance is equal to the maximum elevation of the region where the glacier is located divided by the tan(α). If this distance is greater than 50 km, 50 km should be chosen as the maximum tracking distance. Hock and Holmgren (2005) set zenith angles greater than 78° to 78°, allowing for refraction and other effects to be considered [25]. We initially designated all grid cells in the DEM as being illuminated. For each grid cell of the glacier DEM, the computational algorithm was applied to the start cell along the path of the solar beam and then to each subsequent cell until the maximum distance was reached. If any of the DEM elevations along the solar beam path were higher than the solar beam in that grid, the glacier grid was obscured. Otherwise, the glacier grid was solar-illuminated.

Finally, we obtained hourly glacier topographic shading data for all glaciers within HMA during the ablation period following the above steps. Taking Koxkar glacier (RGI60-13.43232) as an example, Figure 3 shows the topographic shading of that glacier at different time points. The topographic shading at noon was less extensive than that at 08:00 and 18:00 in June and varied greatly with time. The spatial frequency of shading differed significantly during the ablation period (shading frequency refers to the fractional duration of shading during the whole ablation period). The topographic shading was concentrated on the accumulation zone above the median elevation. Some areas were permanently shaded and had a shading frequency equal to 1. Most low-frequency areas below median elevation received more solar radiation, leading to glacier melting.

3. Results

3.1. The Topographic Shading Area Differences between Each Region of HMA

To quantify the topographic shading area extent within each region of HMA, we calculated the regional percentage of the shading area relative to the regional glacier area using the DEM pixels. As shown in Figure 3, the shading area varied greatly with the time. Therefore, we averaged the shading area for individual glaciers during the ablation period. Our results showed apparent differences in the topographic shading areas between each region due to the complex terrain influence, as shown in Figure 4. The regional distribution of the shading area percentage varied from 17.76% to 30.43% during the ablation period. The larger shading area percentages were distributed in the Karakoram, Pamir and West Tianshan, where the shading area percentages were 30.43%, 30.35% and 30.20%, respectively (Table 1), meaning that more than 30% of the regional glacier area was shaded. This can be attributed to the high mountains in these regions. The East Kunlun and Inner Tibet Plateau had smaller shading percentages (19.73% and 17.76%, respectively). The glacier inventory showed that glaciers in the Inner Tibet Plateau have a smaller slope than those in other regions [40]. This indicates that the terrain is relatively flat, which resulted in the small shading area percentages in this study. The overall shading area percentage of all glaciers within HMA was found to be 27.19%.

The glaciers in HMA are located in the middle latitudes of the Northern Hemisphere, and the sun mainly illuminates their south-facing slopes. Therefore, the north-facing slopes of the glaciers are often shaded. Although glaciers are influenced by the orientation of the mountains, the topographic shading is still strongest for the northern aspect. The average shading area percentages for the north, east, south and west aspects for all the glaciers in HMA were found to be 10.97%, 6.72%, 4.36% and 4.31% (

Table 2. Glacier shading area percentages in the regions within HMA.

RGI Region	Pnorth (%)	Peast (%)	Psouth (%)	Pwest (%)	Ptotal (%)
01 Altay and Sayan	14.83	7.56	3.64	3.10	29.13
02 East Tianshan	14.04	6.58	2.83	4.13	27.57
03 West Tianshan	12.23	7.05	5.47	5.45	30.20
04 Hissar Alay	12.79	9.09	4.11	3.80	29.79
05 Pamir	13.56	8.36	4.42	4.02	30.35
06 West Kunlun	9.61	5.21	2.30	4.02	21.13
07 East Kunlun	9.77	4.04	2.82	3.11	19.73
08 Qilian Shan	12.03	3.79	2.03	3.77	21.62
09 Karakoram	11.02	7.97	6.16	5.28	30.43
10 Inner Tibet Plateau	7.77	4.62	2.66	2.71	17.76
11 South and East Tibet	8.86	7.53	5.72	5.04	27.15
12 Hindu Kush	12.89	8.00	4.53	4.20	29.62
13 West Himalaya	10.89	6.85	4.36	4.85	26.96
14 Central Himalaya	9.29	7.24	6.85	5.69	29.07
15 East Himalaya	8.00	7.00	6.72	5.50	27.23
16 Hengduan Shan	8.02	6.79	5.24	4.36	24.40

Pnorth, Peast, Psouth and Pwest are the shading area percentages relative to the glacier areas in the north, east, south and west aspects, respectively. Ptotal is the shading area percentages relative to all the glacier area in HMA.

), respectively. This finding also implies that glaciers on the north-facing slopes are less prone to melting, which is especially worth considering in the simulations of glacier area changes and glacier mass balances in the northern regions.

For the regions in HMA (e.g., Altay and Sayan, East Tianshan, Pamir and Qilian Shan), the shading area percentages of the northern aspects were found to be far greater than those of the other aspects. This phenomenon is not evident for glaciers in the Himalayas and Hengduan Shan, where the shading areas vary less with the aspect, perhaps due to the south–north orientation of the mountains. We also found that East Tianshan, Qilian Shan and West Kunlun glaciers had greater shading percentages of the west aspects than the south aspects. This shows that the regional differences in shading areas is complex due to the influence of the terrain.

3.2. Altitudinal Distribution of Glacier Shading Areas

The shading areas vary with the altitude, as well as horizontally. To evaluate how shading depends on elevation, we calculated a pixel-wise altitudinal distribution of the shading areas for different aspects in each region of HMA (Figure 5). The shading areas were found to increase with the decreasing elevation, but the greatest shading area was typically around the median altitude. However, there were obvious regional differences in altitudinal shading areas for the different aspects. For the glaciers in Altay and Sayan, West Tianshan, East Kunlun and Qilian Shan, the shading areas on the west-, south-, east- and north-facing aspects successively increased from high altitude to low altitude; in addition, high mountain areas were found to obscure illumination of the west- and south-facing slopes in the high elevation zone, while they obscured the east- and north-facing slopes in the low elevation zone. In general, areas in high elevations were found to be in the glacier accumulation zone, while those at low elevations were in the ablation zone. This distribution of shading favors glacier development for the east and north orientations.

In contrast, the shading areas for the west, south, east and north aspects increased from high elevation to low elevation in the West Himalayas, which indicates that the larger mountains cause shading at low elevations in the west and south and at high elevations in the east and north. This pattern promotes glacier development for the west and south aspects, while glaciers for the east and north aspects undergo more intense ablation due to the reduced shading.

3.3. The Relationship between Glacier Area and Shading Area Percentage

Topographic shading protects glaciers from solar radiation. The strength of this protection changes with the glacier shape and area due to variations in the area covered by shading but is greatest for two cases: one in which the shading area percentage decreases with the decreasing glacier area (Figure 6), as commonly seen in most regions, such as glaciers in West Tianshan, Hindu Kush and the East Himalayas (Figure 7); the other pattern is that the shading area percentage increases with the decreasing glacier area, which is more commonly seen for glaciers in West Kunlun, East Kunlun, Qilian Shan and the Inner Tibet Plateau.

These two patterns of shading area are not related to the number of glaciers but are dependent on the glaciers’ orientations (Figure 8). Glaciers in West Kunlun, East Kunlun, Qilian Shan and the Inner Tibet Plateau tend to be distributed on the north-facing aspects, which is remarkably different from the other regional glacier distributions. Shading is more likely to protect small glaciers, because smaller glaciers tend to have a greater shading area percentage. We infer that this is because glaciers in areas with little shading are already absent, while those in well-shaded areas remain, entering a state of self-protection. In contrast, other regions (such as the Himalayas) still have extensive glaciers on south-facing slopes, even if the ablation and net mass loss are intense, as they receive abundant precipitation from the Indian and East Asian monsoons [7,18]. These south-facing glaciers will receive relatively less shading as their areas decrease, making them more susceptible to future retreat.

3.4. Glacier Shading Ratios Beneath the Median Elevation

The glaciers within HMA show significant mass losses, except for those in the Karakoram region [7]. Glacier melting mainly occurs in the ablation zone at low altitude but can be reduced if shading blocks the solar radiation. Therefore, shading at low altitudes has an important influence on glacier melting. To evaluate the shading area ratios at low altitude for each glacier, we calculated the shading area ratio (SAR) at elevations lower than the median elevation. The SAR refers to the ratios of shading area lower than the median elevation to glacier areas lower than the median elevation.

Our results found that the SAR values varied widely between glaciers across HMA, with these values ranging from 0 to 0.97 (Figure 9). Low values indicate glaciers with little shading area under the median elevation, while high values indicate that most of the ablation area of a given glacier is shaded. Over 50% of glaciers in HMA were found to have SAR values greater than 0.25 during the ablation period, with less than 50% of glaciers thus being found to have SAR values less than 0.25. For the glaciers with SAR values less than 0.25, the Inner Tibet Plateau glaciers accounted for 7.5%, while Altay and Sayan glaciers only accounted for 0.6%. Compared with the glaciers that recorded high SARs values greater than 0.25, these glaciers were distributed in Karakoram, Pamir, West Tianshan and the West Himalayas, accounting for 9.5%, 7.3%, 6.5% and 5.3%, respectively. Particularly for the glaciers in Karakoram, the shading area covered most of the glaciers in the low elevation zone, which provides protection from melting and contributes to their positive mass balance state. This supports the hypothesis that topographic factors might contribute to the currently stable regional mass balance of Karakoram glaciers [6]. Similarly, studies also found a positive mass balance for glaciers in Pamir [47,48], which is likely related to the shading area.

Few glaciers recorded high SAR values. For example, glaciers with SAR values greater than 0.34 accounted for 18% of all the glaciers in HMA, of which 3.7% were in Karakoram, 2.5% in Pamir and 2.2% in West Tianshan. This pattern supports the suggestion that the shading area is important for glacier melting in these regions. Glaciers in the Himalayas and Hengduan Shan were found to have an extremely negative mass balance, but less than 1% of the glaciers had SAR values greater than 0.34.

3.5. The Effect of Topographic Shading on the Glacier Mass Balance

Topographic shading protects glaciers from solar illumination and thus affects the glacier mass balance. Figure 10 shows the spatial distribution of the mass balance, topographic shading and DEM for one glacier (RGI60-14.07524). We found that both the mass balance and shading area frequency show broadly consistent spatial distribution, in which the lower shading frequency areas had a more negative mass balance and vice versa, except for the low-elevation ablation zone, where the influence of higher temperatures became stronger (Figure 10c). This indicates that the shading area frequency affects the spatial distribution of the mass balance.

To evaluate whether this relationship held for all the other glaciers, we searched for correlations between the mass balance and shading frequency for all glacier-covered grid points. We found that approximately 20% of the glaciers showed significant correlations (p < 0.05). However, the remaining 80% of the glaciers did not pass this significance test, perhaps because the mass balance is influenced by other climate factors, which masks the influence of the topography. The significant correlation between the spatial distribution of the mass balance and shading for 2000–2016 is shown in Figure 11. The mass balance heterogeneity was obvious for individual glaciers, ranging from −1.82 to 1.55 m w.e. The glaciers with the most positive mass balances appeared in West Tianshan, Pamir, West Kunlun and other regions, illustrating that topographic shading has a striking protective effect on these glaciers.

Table 3. Summary of the regional glaciers with mass balances that were significantly correlated with topographic shading during 2000–2016.

RGI Region	Nglacier	Nnorth	Neast	Nsouth	Nwest	MB (m w.e.)	Area (km2)
01 Altay and Sayan	−	−	−	−	−	−	−
02 East Tianshan	957	570	166	59	162	−0.47	1.254
03 West Tianshan	1792	1084	252	189	267	−0.27	2.747
04 Hissar Alay	479	277	94	38	70	−0.09	1.531
05 Pamir	2112	1275	416	163	258	−0.10	2.59
06 West Kunlun	890	588	111	60	131	0.10	4.422
07 East Kunlun	344	222	37	44	41	−0.09	3.41
08 Qilian Shan	328	203	43	27	55	−0.35	1.642
09 Karakoram	2375	1360	384	250	381	−0.09	5.285
10 Inner Tibet Plateau	1189	705	230	106	148	−0.26	2.712
11 South and East Tibet	615	213	152	109	141	−0.68	3.058
12 Hindu Kush	811	389	172	92	158	−0.17	1.8
13 West Himalaya	1891	976	364	218	333	−0.49	2.215
14 Central Himalaya	1099	382	252	211	254	−0.38	2.732
15 East Himalaya	981	298	193	237	253	−0.48	3.018
16 Hengduan Shan	850	324	248	136	142	−0.72	2.662

Nnorth, Neast, Nsouth and Nwest are the numbers of glaciers with north, east, south and west aspects, respectively. MB is the regional average mass balance, and Area is the regional average glacier area.

summarizes how many regional glaciers had a mass balance that was significantly correlated with topographic shading, showing that most of these glaciers were located on the northern aspects, with relatively few being found on the southern aspects. The proportions of these glaciers on the eastern and western aspects varied from region to region. Overall, shading was mainly found to affect glaciers on the north-facing aspects.

4. Discussion

4.1. Factors Influencing Glacier Mass Balance

Among all the influencing factors, climate forcing is a major external control on the mass change of a glacier. The general patterns of mass balance within HMA follow atmospheric circulation patterns [8]. Recent studies have found that the Indian monsoon is weakening [49] and that the westerlies are strengthening [50], which can influence changes in precipitation patterns. The glacier responses to climate change have exhibited significant spatial heterogeneity in recent decades [10]. Temperature as a main indicator of climate change has accelerated glacial shrinkage within HMA [8]. These previous studies showed that climatic differences between low and high altitudes involve important issues relating to warming. An increasing warming trend at higher elevations has been observed within HMA, and the warming rate has been found to increase with the elevation before becoming quite stable, with a slight decline near the highest elevations observed [51,52]. Moreover, the mass balance is more sensitive to rises in the temperature, particularly for marine glaciers such as those located in the Himalayas and Hengduan Shan [20]. Increases in the precipitation cannot balance the glacier mass loss due to the increasing temperatures. Therefore, the glaciers are in a continuous state of mass loss [53]. This is the reason for the accelerating trend of mass loss.

Solar radiation is the dominant component in the surface energy balance, and it controls glacier ablation [54,55]. It accounts for 60–90% of the total energy budget for glaciers in HMA during the ablation period [17,56,57,58]. However, the spatial distribution of solar radiation is likely influenced by topographic shading and clouds, particularly for rugged terrains. The top-of-atmosphere shortwave radiation of Laohugou glacier No. 12 was reduced by ~30% on average owing to the effects of clouds and the atmosphere [59]. Topographic shading directly impacts solar radiation, altering the radiation intensity at the glacier surface by obstructing solar illumination. Our results confirm that the shading percentages are different for glaciers of different aspects. One of the major limitations of our work is the use of the significant correlation method to statistically analyze the relationship between mass balance and topographic shading, the deficiency of which being that the relationship between them cannot be quantified. To determine the strength of the impact, we resorted to a modeling approach through which we could accurately determine and quantify the impact of topographic shading on the mass balance for each glacier.

We did not consider the effect of debris on the mass balance when we carried out this significant correlation analysis. Debris-covered glaciers are common in the rugged Central Himalayas, but they are almost absent in the subdued landscapes of the TP, where the retreat rates are higher. In contrast, more than 50% of observed glaciers in the westerlies-influenced Karakoram region in the Northwestern Himalayas are advancing or stable [17,60]. The accelerating or retarding effect of debris cover on glacier melting depends on the debris thickness [15]. Debris exceeding a certain thickness will retard glacier melting [61]. Otherwise, it will accelerate glacier melting [62]. This could explain why these are the glaciers with the lowest correlation coefficients between mass balance and topographic shading. In addition, cloud cover and albedo also indirectly have an important impact on the glacier mass balance. Quantifying the impact of all factors on the glacier mass balance will be a key point of future research.

4.2. Comparison to Topographic Shading Estimates

Studies of topographic shading throughout HMA are quite scarce. There is also a lack of comparisons between previous results. Past studies have used topographic shading data to drive mass balance models for glaciers in Swedish Lapland [32], Switzerland [24] and elsewhere in the world [25,26,37]. The hill shading tool in ArcGIS can be used to calculate topographic shading [28,31]. Although the hill shading value of the center cell of a moving window can be calculated in relation to neighboring cells, this method only considers the influence of the local grid without considering high mountains at greater distances. To compare the ArcGIS results with our results, we calculated the hill shading used in ArcGIS, setting the solar azimuth to 225° and the solar altitude angle to 55°. We also calculated the topographic shading with the same azimuth and altitude parameters as those used in the ArcGis calculations. We found that the topographic shading area calculated in ArcGis was far smaller than that calculated by our method under the same conditions (Figure 12). Mountain glaciers (such as the Koxkar glacier) in HMA are distributed in the mountain valley. The surrounding high mountains obscure solar illumination, leading to high proportions of shading on the region’s glaciers. Since the ArcGIS calculations ignored this distant mountain influence, there remains great uncertainty in its results. However, we considered this distant mountain influence by tracing the solar rays for a distance of 50 km to more accurately determine whether the glacier surface was shaded. This approach considerably improved the shading area calculation accuracy.

4.3. Uncertainty in Topographic Shading and Mass Balance Data

This study used topography data extracted from SRTM DEM version 4 with a spatial resolution of 90 m [42], which was the basis for the topographic shading calculations. However, when clipping the raster data (such as DEM) using the glacier vector boundaries, both of the nonoverlapping boundaries will inevitably cause differences in the area. According to the glacier inventories used, the area uncertainty for each glacier was assumed to be 3.2–10% [40,63]. All these factors contributed to the uncertainty in our results. Although the uncertainties would decrease if the resolution improved, this approach would increase the computational costs. Consequently, studies of glacier thickness estimates use DEMs with a low resolution, typically 100 m or more for larger glaciers [2], especially when considering all HMA glaciers or the global scale. In addition, the DEM data represent the average elevation within the 90-m grid, which does not reflect the sub-grid scale topographic features. The real topography may obscure the sun’s radiation, even if the average DEM grid elevation used to calculate the topographic shading is lower than the solar altitude angle. This could particularly influence the results in rugged terrain.

The topographic shading in this study focused on the ablation period and was calculated without considering the specific year. Although the annual topographic shading changes over long time scales, such changes are extremely slow, as they are related to the varying distances from the Earth to the Sun caused by the variations in the Earth’s orbit (orbit scale: 10–100 kyr) [64]. Consequently, the orbital influences on topographic shading were ignored. Our results are valid for energy balance studies over long time periods, but it is worth noting that the glacier inventories used in this study were based on data for 1999–2010 [40,41]. Topographic shading is closely related to the inventory data, because glaciers are in a state of change (e.g., due to temperature change), which alters the shading area of each glacier.

Alternative mass balance data resources are also subject to uncertainty with regards to their data quality. Individual mass balance data have a relatively high level of uncertainty, with a median uncertainty of 0.22 m w.e./a having been recorded (in the range from 0.14 to 0.89 m w.e./a); This depended on the glacier area, on the proportion of the glacier surface surveyed and the number of DEMs available to extract a reliable rate of elevation change [19]. The median uncertainty was slightly higher than the root mean square error of 0.17 m w.e./a found when comparing the ASTER mass balance with higher-resolution geodetic mass balances for 60 glaciers [7].

5. Conclusions

Our goal was to provide a systematic estimate of the relationship between mass balance and topographic shading across HMA. To achieve this, we used the glacier outlines and raster data (DEM, slope and aspect) to systematically estimate for topographic shading across HMA. We found a clear correlation between mass balance and topographic shading, which was most obvious for glaciers of the north-facing aspects. The main conclusions are as follows:

Overall, 27.19% of the HMA glacier area was found to be affected by topographic shading. The north, east, south and west aspects accounted for 10.97%, 6.72%, 4.36% and 4.31% of this total shading, respectively. The differences in the topographic shading area were apparent between each region, and the regional distribution of the shading area fraction varied from 17.76% to 30.43% during the ablation period. Glaciers with the greatest shading percentages were found to be distributed in the Karakoram, Pamir and West Tianshan, recording shading area percentages of 30.43%, 30.35% and 30.20%, respectively. Moreover, regional differences were found to exist in the altitudinal distribution and the SAR below the median elevation. The glacier shading areas under the median elevation in the Karakoram, Pamir, West Tianshan and West Himalayas were found to be mainly distributed at low elevations, which accounted for 9.5%, 7.3%, 6.5% and 5.3% of the total HMA glacier area. This was especially true for the glaciers in Karakoram, where shading was found to cover most glaciers in the low-elevation zone, thus protecting them from radiative melting and contributing to their positive mass balance state.

Topographic shading was found to contribute to the protection of the smallest glaciers from solar illumination. These small glaciers are located on the north-facing aspects, especially in West Kunlun, East Kunlun, the Inner Tibet Plateau and Qilian Shan, where the shading area percentages increased with the decreasing glacier area. In contrast, in other regions such as the Himalayas, glaciers on the south-facing aspects maintained a large glacier area, but the shading percentages decreased as the glacier area decreased. These glaciers are at risk of disappearing in the future.

Approximately 20% of the glaciers with mass balance data showed a significant correlation between topographic shading and mass balance, passing the significance test at p < 0.05. These glaciers are mainly located at the northern aspect, indicating that complex terrain can affect the spatial distribution of the mass balance by altering the solar illumination pattern.