4 | Results and Discussions

4.1 | Ice Thickness Measurements

We selected six typical glaciers of the UIB region to conduct GPR measurements, their geographical location and information was presented in Figure 3, they have different shapes and are suitable for GPR investigation because of spatial representation and location accessibility. Batura, Pasu, and Barpu glacier belong to Hunza sub-basin, the 59.2 km length makes Batura ranked 4thon the eight longest glaciers in the middle-low latitude region, it has the maximum value of vertical height difference and was classified as a mustagh-type glacier for multiple branches nourished by large avalanches (Boerst, Winiger & Bookhagen, 2013; Hewitt, 2011). Pasu, is adjacent to Batura, has retreated about 1.5 km since the end of the 19th century, it is characterized by small end moraines that consist mostly of till in the form of hammocks and small parallel ridges (Owen & Leicester, 1989). Barpu is lying in the southwestern part of the Karakoram range. Sachen and Chungphar glacier belong to the Astore sub-basin, Sachen is located on the eastern slope of Nanga Parbat massif nourished mostly by ice-fall avalanches, in contrast, Chungphar is situated on the southern foothills, and one tributary has shown a low retreat rate, which is characterized by extremely steep topography and vertical gradients. Gharko glacier belongs to UIB_D sub-basin. Figure 3 (b)-(c) exhibited the geomorphological feature on Barpu and Pasu glacier, we examined these glaciers were debris mantled as a consequence of abundant supply by rock avalanching process, with low ice velocities and no efficient system of debris removal (Shroder, Bishop, Copland & Sloan, 2000).
GPR work was from August to October during 2016 to 2018, the surveying data were recorded with a distance ranging from 30 to 50 m on the glacier surface along with the transverse profiles, the coincided GPS with GPR measured locations and routes were depicted in Figure 3 (d)-(i) subplots respectively. We completed four transverse and three longitudinal profiles from the terminus to the upper reaches on Sachen glacier along with the central flowline and conducted six and five profiles on Barpu and Gharko glacier, respectively. Three transverse profiles in the tongue area of Batura, four profiles in the fissure zone of Pasu, and three profiles within a small elevation below 3000 m on Chungphar were finished. Totally, 26 GPR topographic profiles and 187 points, the ice thickness of the glacier boundary was assumed to be zero, including rock outcrops.
FIGURE 3 The GPR observation overview, (a) Location of six surveyed glaciers in the UIB sub-basins. (b)-(c) field working photography that displayed geomorphological features on Barpu and Pasu glacier. (d)-(i) The GPR surveyed routes and points, the solid blue dots and black lines denote surveyed points and profiles, respectively
The discrimination of ice thickness mainly relies on the quality of the GPR data, Figure 4 has shown some representative raw radargrams. Hyperbola presented the englacial object signature, they have a complete and a sharp radar wave that we can distinguish the interface between ice and bedrock and correct the dt -value on the GPR screen in the case of the noise. Any change in the dielectric constant of ice and rock shows a change in reflection, sometimes, the non-bed reflector occurred because of the entrained debris layer, vertical discontinuity, and meltwater in summer. The whole conducted points and calculated results were posted in the appendix information Table S1.
FIGURE 4 The representative GPR radargrams of six selected glaciers

4.2 | Optimized parametric schemes of GlabTop2

GlabTop2 estimates were based on SRTM-C DEM that released in 2011, while GPR measurements were carried out from 2016 to 2018. The estimated ice thickness and surveyed ones could not be compared directly because the ice surface changed during the past decades, nevertheless, ice bed erosion occurred over the geological time scale, and it was usually calculated on a centennial-scale even though in rapid tectonic regions. The glacier bed erosion rate generally ranges from 1 to over 10 mm yr-1, the ice bed position was assumed to be relatively much more stable than the ice thickness calculated on the modern time scale (Koppes & Montgomery, 2009).
The GlabTop2-modeled result is intensely dependent on the parameters ofτ and f , especially for the upper limit onτ max. It has been revealed the optimalτ- value of five glaciers in northwestern China was between 50 to 175 kPa, and the fluctuation interval is large (Li, Ng, Li, Qin & Cheng, 2012), τ = 150 kPa was commonly used in the Karakoram range (Frey et al., 2014).. Regarding the f- value, it represents the friction between ice body and valley sides, depending on the aspect ratio of glacier width and midpoints’ depth along the flow line. Namely, when the measured half-width becomes wider or the midpoints’ depth turns shallower, the f -value will turn larger, vice versa. It was verified that the measured glaciers in HKKH ranges have narrow shape, the values of half-width were relatively small, accompanied with thicker ice depth on central flowline.
As presented in Figure 5, we compared the GPR-surveyed ice bed topography with GlabTop2-modeled results under nine parametric schemes, theτ max-value was processed ranging from 100 to 150 kPa that increased by 20 kPa and 30 kPa gradually, coupled with thef -value ranging from 0.7 to 0.9 that increased by 0.1 grades. To assess the reliability of the GlabTop2 intuitively, three commonly used evaluation indexes were quoted, including the Root Mean Square Error (RMSE), Mean Deviation (MD), and Nash Sutcliffe Efficiency (NSE):
\(\text{RMSE}=\sqrt{\frac{1}{n}\sum_{i=1}^{n}\left(h_{o}-h_{m}\right)^{2}}\)(5)
\(MD=\frac{1}{n}\sum_{i=1}^{n}\left(h_{o}-h_{m}\right)\) (6)
\(NSE=1-\frac{\sum_{i=1}^{n}\left(h_{o}^{n}-h_{m}^{n}\right)^{2}}{\sum_{i=1}^{n}\left(h_{o}^{n}-\overset{\overline{}}{h_{o}}\right)^{2}}\)(7)
where ho refers to the observed ice bed elevation; hm refers to the modeled one;\(\overset{\overline{}}{h_{m}}\) refers to the mean value of observed results; and n means the number of measured points, these variables in equations (5) to (7) are identical. GlabTop2 performs high quality and better credibility when NSE-value is close to 1, together with the smaller RMSE and MD-values. There is a concept that ice bed elevation was deduced by ice surface elevation minus ice thickness, the thicker of ice depth, the lower of the ice bed elevation, implying that these two variables made an inversely proportional relationship on magnitude. GlabTop2 overestimated the actual ice bed elevation compared to GPR-measured results when τ- value became smaller and f -value turned higher, that means we underestimated the actual ice thickness. In contrast, GlabTop2 overestimated the actual ice thickness compared to the surveyed results when the τ- value became higher and the f- value turned smaller.
From the perspective of individual glaciers, we processed GlabTop2 on three kinds of glaciers (small, medium, and large-sized). Fig 5 (a)-(c) presented the comparison results of Batura, Barpu, and Sachen glaciers that possessed considerable GPR- measured samples. Batura (243.5 km2) belongs to the large-sized glacier, its RMSE, MD, and NSE-values of the GlabTop2-modeled ice bed elevation versus the GPR-measured ones were coincided best when τ = 150 kPa and f = 0.8, Barpu (90.6 km2) belongs to the medium-sized glacier, the values of these referral indicators performed best when τ = 120 kPa andf = 0.8, while Sachen (9.5 km2) belongs to a small-sized glacier, GlabTop2 experiment was not very ideal for the NSE-values were close to 0.5 even extremely negative, displaying a divergent condition, the simulated ice bed elevation versus surveyed results fit better when τ = 150 kPa and f = 0.8, even the MD-value reached to -1.64 m. The number of observed GPR samples on Pasu, Gharko, and Chhungpar was relatively inadequate, thence they were not explained separately, they have been discussed entirely in Figure 5 (d).
If we choose an identical parametric scheme to apply on all glaciers that various in geometry and acreage, it is a high probability that all glaciers would not perform best under one certain scheme, a reasonable solution is to select the most comprehensive one. From the perspective of total GPR-surveyed points of six glaciers, the NSE-values of these nine parametric schemes were all greater than 90%, the MD-values of all validation points were between negative and positive 50 m, and the RMSE-values were below 80 m. Amazing consistency appeared on GlabTop2-estimated ice bed elevation versus GPR-surveyed results, demonstrating that GlabTop2 performed effectively because it kept the small gap with our observed data. Among nine parametric schemes, we choose four scenarios that raced better with smaller RMSE (≤ 65 m), MD-values (≤ 25 m), and reasonable NSE-values (≥ 0.967), the first two,τ = 150 kPa and f = 0.8 & 0.9; the last two, τ = 120 kPa and f = 0.7 & 0.8. As demonstrated in Figure 6, by comparing the profiles of the GPR-surveyed ice bed elevation versus GlabTop2-reconstructed results along with the transverse and vertical lines, it was discovered that two schemes had a relatively large gap when τ = 150 kPa and f= 0.8 & 0.9. Prominently, the simulated ice bed landform was much shallower when the f -value turns higher. However, these two simulated results were extremely close when τ = 120 kPa andf = 0.7 & 0.8, the f -value had substantially less effect on estimated results than τ -value for the HKKH region. In summary, the surveyed ice bed morphology and estimated results were fit better on BaG_01, BaG_03, BpG_01, BpG_03, PaG_03, PaG_04, SaG_03, SaG_04, SaG_06, GhG_02, GhG_05, and ChG_03 profiles when τ = 120 kPa, it has covered three kinds of sized glaciers.
FIGURE 5 A comparison between the GPR-measured ice bed elevation with the GlabTop2-estimated results on (a)-(c) various sized glaciers and (d) total samples under nine parametric schemes
FIGURE 6 A comparison between the GPR-surveyed ice bed elevation with the GlabTop2-simulated results along with the measuring profiles under four selected parametric schemes
In order to select the uniform parameters scheme, Figure 7 exhibited the ice reserves and average ice thickness reconstructed by GlabTop2, GlabTop, Volume and Topography Automation model (Volta), and IDW interpolated methods, the IDW results of Sachen and Barpu glacier relying on a considerable amount of GPR points. These three SIA models have the identical principle that just DEM and glacier outline needed. Among them, ice thickness is obtained along several manually digitized branch lines by GlabTop, and Volta is a python script installed in ArcGIS that generated centerline network automatically based on the GIS hydrological tool, distance, and geographic analysis. For small-sized glacier, Sachen, the ice volume and average ice thickness reconstructed by GlabTop2, GlabTop, and Volta were almost kept in the order of magnitude. Taking the IDW-calculated ice reserves (0.8 km3) and mean ice thickness (89.9 m) as the most precise reference due to its valuably, GlabTop2 preserved the most similar ice volume quantity (1.0 km3), mean ice thickness (105.0 m), and spatial distribution feature when τ = 120 kPa and f = 0.8, the estimated results from GlabTop and Volta deviated far from the IDW-calculated results. It was demonstrated that the actual condition was agreed with the GlabTop2 simulation effect whenτ = 120 kPa. However, the model performed best whenτ = 150 kPa in the previous comparison of ice bed topography plots, indicating that the τ- value impacted the final calculated results weakly on a small-sized glacier, even can be ignored. The trustworthiness of the selected parameter scheme (τ = 120 kPa) was proved, despite slight differences (0.2 km3). Moreover, the glacier branch lines generated by Volta were much more straightforward than GlabTop that extracted from the contour lines manually. For medium-sized glacier, Barpu, its IDW-calculated ice volume (6.6 km3) was less than GlabTop2 estimates (7.3 km3) slightly when τ = 120 kPa and f = 0.8, but far away from the GlabTop simulated result (7.8 km3), clarifying that this parametric scheme was also suitable for the medium-scale glacier. As for large-sized glacier, Batura, the Volta modeled ice volume (23.1 km3) was close to GlabTop2 simulated result (25.7 km3) when τ = 120 kPa and f = 0.8.
It was worth noting that the IDW results usually generated shallower mean ice thickness and less ice reserves than physical models, overestimated phenomena occurred on models’ application. The ice thickness distribution reconstructed by GlabTop2 looks more authentic because it has taken surface morphology of per glacial grid into account. GlabTop and Volta emerged relative ideal status because the thickest parts always posited on the centerlines. Moreover, the vast majority of glaciers in the HKKH chains are cramped and slender, a moderate f -value conforms to the reality. Ultimately, it is more suitable to choose the parametric scheme (τ = 120 kPa andf = 0.8) applied to the whole UIB region.
FIGURE 7 The comparison of the estimated ice thickness among GlabTop2 (τ = 120 kPa andf = 0.8), GlabTop, Volta models, and IDW interpolated results on three kinds of sized glaciers, the solid black points and lines represented the GPR points and flowlines

4.3 | Spatial Distribution of the Simulated Ice Reserves

We processed GlabTop2 depends on the verified optimal parametric scheme, three consistent boundary files were prepared, including glacier mask, SRTM-C DEM, and generated slope. The estimated ice thickness ranges distribution of the UIB sub-basins and the statistical glacier information were summarized in Figure 8 & 9, containing the total number, glacier area, basin area, average ice thickness, as well as ice reserves.
GlabTop2-estimated ice depth was varied from 0 to 488.1 m, with an average value of 74.4 m. The average ice thickness in the Karakoram range (81.5 m) was remarkably thicker than the Western Himalaya (60.3 m) and Hindukush range (45.0 m), it was lower than previous research slightly (Frey et al., 2014). The most abundant ice resources were focused on the Karakoram range that a semi dispersive glacier area concentrated massively, Hunza (224.1 km3), Shigar (229.4 km3), and Shyok (585.6 km3) high-altitude sub-catchments dominated an absolute advantage (81.8%). Less ice volume was stored in the Western Himalaya range, consisting of Astore (15 km3), Shiquanhe (7.7 km3) and Kharmong (120.9 km3) sub-basins that occupied 11.3% proportion, while the least quantity was found in the Hindukush range, it is composed of Gilgit (50.3 km3) and UIB_D (37.0 km3) sub-watershed that just accounted for a small proportion (6.9%). Undoubtedly, from the viewpoint of geological structural units, significant spatial heterogeneity of glacier resources distribution exists in the UIB sub-basins, mainly was decided by the terrain and climatic condition. Ice thickness will become much shallower along the steep terrain at higher altitudes above 4000 m, on the contrary, the thicker ice depth was distributed along the low-lying terrain below 4000 m, most importantly, orographical precipitation provides abundant material supply. Therefore, the favorable environment for glacier development contained flat topography of the high-altitude region, integrated with sufficient materials accumulated, the terrain and climatic factors in the Karakoram range is much more beneficial for glacier emerging and preserving.
FIGURE 8 The GlabTop2-estimated ice thickness ranges distribution of the UIB sub-basins
FIGURE 9 The GlabTop2-estimated ice thickness ranges distribution and statistic glacier information of the UIB region
The heterogeneity distribution of water resources converted from ice volume (assuming a mean glacier density of 900 kg m-3), corresponding water level equivalents (WLE), and its available consumption time in the UIB sub-catchments were revealed in Figure 10 (a)-(c) successively. When discussing the converted water resources, its ranking and proportion feature was same as the ice reserves. Shyok (527.1 km3), Shigar (206.4 km3), and Hunza (201.6 km3) were in three top-ranked, can be regarded as the adequate region that possessed ample water resources more than 200 km3. The distribution of normal regions was dispersed relatively, including Kharmong (108.8 km3), Gilgit (45.3 km3), and UIB_D (33.3 km3) for they maintained the resemble quantity that less than 120 km3, in contrast, Astore (13.3 km3) and Shiquanhe (6.9 km3) were classified as the poor district that owned a few water resources (≤20 km3).
We took sub-watersheds’ area as the referenced divisor when assessing the WLE, its distribution feature and ranking gradation were different from the ice reserves. The maximum value of water resources per square kilometers emerged in Shigar (29.3 m) because of the minimal basin area, Shyok (15.9 m) and Hunza (14.7 m) ranked behind inseparable, nevertheless, the WLE-value of Gilgit (3.6 m), Astore (3.3 m), Kharmong (2.4 m), UIB_D (1.2 m), and Shiquanhe (0.3 m) are below 5 m, the most extensive drainage area of Kharmong river causing the water scarcity occurred. Kharmong, overlapped with the Kashmir, is confronted with the dilemma of water resources, resulting in the territorial and resource dispute between Pakistan and India. Shiquanhe, located in the northwest TP, sorted bottom in all assessed units, it illustrated that short supply of water resources exists in the headwater region of the UIB mostly, guaranteeing water supply for the Shiquanhe sub-basin has become an urgent task between Pakistan and China. Not only upper reaches are basically in deficit condition, but also the water resources distribution of the downstream area is pessimistic. The UIB relies on water resources in the middle reaches that was comprised of Shigar, Hunza, and Shyok sub-catchments because they are keeping in a self-sufficient even surplus state.

4.4 | Spatial Homogeneous of the Hydrological Effect

Ice reserves that stored in the UIB region is a significant-high-quality natural substantial reservoir for the downstream area of the Indus river, the importance of its hydrological effect is never overemphasized. According to the discharge data of hydrological stations in the UIB sub-basins that acquired from the Pakistan Water & Power Development Authority (WAPDA), we analyzed the multi-year average discharge and interannual change trend by statistical unary linear regression as exhibited in Figure 10 (d)-(k). The annual discharge of Kharmong (456.3 m3s-1), Shyok (363.3 m3s-1), Hunza (298.9 m3s-1), Gilgit (284.6 m3s-1), and Shigar (207.3 m3s-1) were quite higher than Astore (143.3 m3 s-1) and Shiquanhe (76.3 m3 s-1) sub-basins.
FIGURE 10 The heterogeneity feature of (a) converted water resources, (b) water level equivalents, and (c) available consumption time of water resources of the UIB sub-basins without considering glacier process and other variables. (d)-(k) demonstrated the interannual change trends of discharge based on hydrological stations
The available consumption time of glacier resources in the UIB sub-basins various in magnitude noticeably, it was inconsistent with the ice reserves. An available consumption time of Shyok is 46.0 y, that is to say, the ice resources will be depleted within 46.0 y under ideal status, Shigar (31.6 y) and Hunza (21.4 y) will be depleted within 32.0 y. On the contrary, the available consumption time of Kharmong (7.6 y) and Gilgit (5.0 y) sub-basins were less than 8.0 y, even the depleted time of Astore (3.0 y) and Shiquanhe (3.0 y) sub-catchments were less than 3.0 y, they have the high expending speed of water reserves. Outstandingly, Kharmong owns a small WLE-value, but it can supply 7.6 y for the outlet of the tributary, indicating that its dependence on glacier melting is relatively lower than Gilgit, Astore, and Shiquanhe sub-basins. Precipitation variations influence the Astore river predominantly, a slight decrease in annual and summer precipitation was observed in the last 30 years (Farhan, Zhang, Ma, Guo & Ma, 2015). Accompanied by the annual discharge increased slightly from 1980 to 2015 as revealed in Figure 10 (i), it was confirmed that Astore sub-watershed was facing the lack of precipitation replenishment and accelerated glacier melting meanwhile. A total of 1270 km3 ice volume, corresponding to 1143.7 km3 water resources, there was 2249.2 m3 s-1 discharge at the Besham hydrological station near the Taberla Dam annually, the converted water resources will be exhausted within 16.1 y as runoff for the UIB outlet on average. The above explanation proves that the UIB will be up against water shortage in the future evenly, especially in the sub-watersheds that located in the upper reaches and downstream area.

4.5 | Uncertainty Analysis

We analyzed what extent uncertainties in the inputs of GlabTop2 by incrementally changing parameters (τ and f ). Figure 11 (a)-(d) described the boxplots of estimated ice volume and average ice thickness under nine parametric schemes, Table 1 can be referred to as detailed counted results. Discussed by glacier scale, the largest error of ice volume emerged on the large-sized glacier, Batura, maximum ice volume was twice than the minimum result, and the median value was about 25.0 km3 that close to the result of the optimal parametric scheme aforementioned (τ = 120 kPa and f = 0.8). For medium-sized glacier, the median values of Barpu and Pasu glacier were not different from the optimal modeled results strongly, it should be noted that for a small-sized glacier, the parameters’ influence on the final simulated result was minimal, almost negligible. Discussed by sub-basins, the largest uncertainty exists in Shyok, Shigar, Hunza, the median value of Shyok was less than 600 km3 approximately that meet the final estimated result from the optimal parametric scheme selected previously, testifying that our estimation has no excessive deviation. As for ice thickness, various sized glaciers and sub-basins were basically below 150 m, the largest value occurred on Sachen glacier and Shigar sub-catchment, the IDW interpolated results of Barpu and Sachen were also fitted with the median values reconstructed by GlabTop2.
We defined ”sensitivity” (ΔV ) as various in estimated ice thickness divided by unperturbed value, expressed as a percentage, the unperturbed value is the mean ice thickness deduced by the optimal parametric scheme (τ = 120 kPa and f = 0.8). Figure 11 (e)-(f) plots (әh/әf & әh/әτ )/h for an input perturbation of τ = 50, 100, 120, 150 kPa and f = 0.5, 0.6, 0.7, 0.8, 0.9, 1. In absolute terms, while a decrease of 20 kPa in τ- value can reduce the estimated ice thickness by 20%, and a 30 kPa plus inτ- value increased simulates about 25%. By comparison, a 0.1 rising the f -value reduced ice thickness by about 10%, and a 0.1 dropped f- value increased 15% estimate, it was proved that the final estimates were less sensitive to the f -value calculated by half-width.
FIGURE 11 Boxplots of the GlabTop2-estimated (a)-(b) ice volume and (c)-(d) average ice thickness on six selected glaciers and UIB sub-basins, the green triangle donates the median value. The sensitivity analysis of estimates for an increase on the (e) τ -value and (f)f -value