3 | DATA AND METHODOLOGY

3.1 | Data (Glacier inventory and digital elevation model)

Spatial data included the digital elevation model (DEM) and PGI, SUPARCO released the available latest PGI database in cooperation with ITPCAS, previous studies rely on the old version of the glacier database mostly, therefore the accumulated glacier area was too large. PGI was originated from Landsat-8 (30 m) imagery that acquired from 2013 to 2015. The potential error in the glacial boundary estimated measurement (Eb ), in case of the debris-covered is 2.6%, whereas for the clean glacier area is 1.9%, with the calculated adjustment error is below 0.35%. Similarly, to assess the average error (θ ), the PGI was compared with the boundaries delineated through some high-resolution multispectral images of Pleiades and SPOT-5, 6 & 7. In total, 32 glaciers were randomly selected for this comparison, the result shows less than 9% difference in the debris-covered outlines, on the other hand, the θ -value in clean glacier outline is found less than 5% (SUPARCO & ITPCAS, 2015). We counted smaller glacier area because of the improved mapping accuracy, and PGI excluded glaciers with an area that less than 0.01 km2. Meanwhile, we used the Chinese Glacier Inventory version 2.0 (CGI2.0) (http://data.tpdc.ac.cn/) and GLIMS Randolph Glacier Inventory version 6.0 (RGI6.0) (http://www.glims.org/RGI/) as the supplemental database to determine glacier outline of a small part of Shyok and Kharmong sub-basins that PGI uncovered (Guo et al., 2015).
The quality of GlabTop2-modeled ice thickness largely depends on DEM, and we derived the topographic inputs from the free Shuttle Radar Topography Mission Digital Elevation Model (SRTM-C DEM) (http://dds.cr.usgs.gov/srtm/version2_1/SRTM1), which has a high spatial resolution of 1 arc-second (approximately 30 m). SRTM-C DEM was released in 2011 by NASA and NGA and acquired using a radar interferometry technique. We can collect the most accurate glacier surface elevation because it can penetrate many meters into snow firm efficiently and has been used successfully over 80% of the earth’s land surface between 60°N and 56°S latitude interval. In the subsequent study, we mosaiced the SRTM-C DEM, and projected it to Universal Transverse Mercator Projection (UTM43N) and World Geodetic System 1984 ellipsoidal elevation (WGS84), all GPS data were kept the unified projections.

3.2 | Methodology

3.2.1 | Ground-penetrating radar and GPS field investigation

In the 1980s, the Lanzhou Institute of Glaciology and Geocryology, Chinese Academy of Science conducted the B-1 GPR work on the Antarctic Nelson, Urumqi glacier No. 1, and Hailuogou glacier (Su, Ding & Liu, 1984; Zhang, Zhu, Qian, Chen & Shen, 1985). A direct comparison between the GPR-surveyed ice depth with the results derived from the hot water drilling technique indicated that radar results were usually within ± 5% deviation (Liu, 2012), it was verified by the precise ice core drilling data down to the bedrock in the ablation area of the Hailuogou glacier, the B-1 GPR were quite satisfied with an error below 5 m (Wang, Li, Shen & Huang, 1996). In this study, we used an impulse B-1 modified GPR system with a separate transmitter and receiver antenna separated by a fixed distance of 5 or 10 m, the low frequency (approximately 2–220 MHz) is much more suitable for probing mountain glaciers, this GPR has a standard offset geometry with point measuring mode with a 5-MHz resistively loaded dipole antenna length of 10 m and wavelength of 33.8 nm. On the two-dimensional radar graphic, we derived ice thickness (h ) from the vertical axis radar wave and calculated the two-way travel time by the following equation:
\(h=\frac{\sqrt{v^{2}t^{2}-x^{2}}}{2}\) (1)
where t refers to the two-way travel time of radar wave, xrefers to the distance between transmitter and receiver antenna, andv refers to the velocity of radar signal in the glacier. The speed of electromagnetic wave propagation in the ice was assessed to be 0.169-0.171 m ns-1, we set it to be 0.169 ± 0.002 m ns-1, ensuring the relative error was within the accuracy demands of glaciology research.
During the actual operation process, firstly, we adjusted a start point of echo wave to check the whole wave shape was displayed on the GPR screen, then measured the t -value between the direct wave arrival through the air and the reflections from the ice bedrock, finally, we determined the h -value of measured points by identifying the ice and bedrock interface in the radar graphic. Glacier surface elevation of GPR points was surveyed simultaneously by a differential global positioning system (DGPS) device (SF-3050 GNSS, NavCom Technology, Inc., accuracy < 0.05 m) and a portable GPS (Shtech, accuracy≤ ± 1 m). The accuracy of ice thickness was decided by two factors: one is the measuring system, another is the property of ice and bedrock, the accuracy of travel time was determined to be ± 10 ns (1.6 m) from the oscilloscope trace. In our study, Batura glacier, the maximum error was ± 2.4 m (the highest thickness measurement 201 m), within the quoted error at ± 5%.

3.2.2 | Description of GlabTop2 Model

GlabTop2 is a grid-based and slope-dependent ice thickness estimate model, Linsbauer et al. (2009) established and developed it based on the flow mechanics of the relationship between average basal shear stress (τ ) with shape factor (f ) (Nye, 1952; Paterson, 1970), the detailed formula explanation of GlabTop2 were depicted as follows:
\(h=\frac{\tau}{\text{fρg\ }\sin(\alpha)}\) (2)
\begin{equation} H\leq 1.6\ km,\tau=0.005+1.598H-0.435{H}^{2}\nonumber \\ \end{equation}
\(H>1.6\ \text{km},\ \ \tau=150\ kPa\) (3)
\(f=\frac{w}{h_{m}}\) (4)
where the quantity of τ is under a perfect-plasticity assumption for the flow rheology. τ is parameterized with the glacier vertical elevation range (ΔH ) (Haeberli & Hölzle, 1995). The f -values depends on the aspect ratio of half-width (w ) to the midpoints’ depth on central flowline (hm ) along the glacier cross-section, it is usually set to 0.8 for all glaciers (Paterson, 1994). ρ refers to the ice density (900 kg m-3), refers to the gravitational acceleration (9.81 m s-2), andα refers to the surface slope. A detailed working flow chart for processing data, GPR diagram, and GlabTop2 illustration schematic were presented in Figure 2. GlabTop2 calculated ice thickness for selection of a set of randomly picked DEM cells within the glacierized area, then interpolated to the entire mask, which was fully automated and required only glacier mask, DEM, and surface slope as inputs. It is indispensable for solving the problem of the non-measured area (Farinotti, 2017; Frey et al., 2014), and the prospect is considerable by using GPR-measured ice thickness to modify the GlabTop2 results.
FIGURE 2 (a) Working flow chart for processing data, (b) GPR (Liu, 2012), and (c) GlabTop2 illustration schematic (Frey et al., 2014)