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