3.1 Satellite images
Series of Landsat images were used for surface velocity extraction during the last three decades (1992-2019). Six images from Landsat 5 Thematic Mapper (TM) mission were processed to quantify surface velocity between 1992-94 and 2009-11. Two scenes from the Landsat 7 ETM+ (Enhanced Thematic Mapper Plus) mission were used for 2000 and 2002 (Table 1). A series of Landsat 8 OLI (Operational Land Imager) images were used to monitor glacier surface velocity variations between 2013 and 2019 (Table 1). All images were taken at the end of the ablation period (September-October) for annual velocity measurements, reducing the seasonal snow cover and shadowing effects. A complete Landsat scene (path 147 and row 37) covering the entire upper Chandra basin and adjacent regions was processed for surface velocity measurements. Glacier's outline was derived from the RGI V6 database (RGI Consortium 2017). The debris-covered outline was taken from the previous study (Herreid and Pellicciotti, 2020). Glacier outlines were manually checked and corrected based on the higher-resolution images available in Google Earth and pan-sharpened Landsat 8 images of 2020. Topographic parameters were extracted from the 30 m global ASTER GDEM.
3.2 Available image matching software
Different methods have been developed and applied to correlate optical images in glaciology, like normalized cross-correlation (Bindschadler et al. 1999; Berthier et al. 2003; Kääb, 2002), cross-correlation operated in the Fourier domain (Rolstad et al. 1997), least-squares matching (Kaufmann and Ladstädter 2003), phase correlation (Avouac et al. 2006), and orientation correlation (CCF-O) (Haug et al. 2010). However, it has been reported that CCF-O and phase correlation available in COSI-Corr are relatively more robust matching methods for global-scale mapping and monitoring of glacier velocities (Heid and Kääb 2012). A few open-source image processing software is available, e.g., CIAS, IMCORR, COSI-Corr, IMGRAFT etc. (see online resource), which utilize the above principle, but their applications have been limited due to tedious and time-consuming processing. However, the COSI-Corr tool has been performed better than others in processing time and output variables (Jawak et al. 2018). Further, this software has several pre and post-processing tools, which give additional advantages (Ayoub et al. 2009).
In this study, COSI-Corr (an add-on module in ENVI software) was used to derive surface displacement. The COSI-Corr was initially used to extract the crustal deformation field induced by an earthquake released by the tectonic observatory at Caltech in 2008 (Ayoub et al. 2009). This technique is now widely used to monitor the Earth's surface changes (Leprince et al. 2008) and glacier monitoring (Scherler et al. 2008). This software allows precise co-registration, orthorectification, and sub-pixel correlation of remote sensing images, all in one package and a more user-friendly environment. The derived output from COSI‒Corr provides three bands, i.e., east-west displacement, north-south displacement, and signal-to-noise ratio. Most importantly, the signal-to-noise band helps to assess the quality of correlation between images. During processing (correlation), the errors from different sources tend to combine and lead to a relatively higher error in the final result (Scherler et al. 2008; Ding et al. 2016). COSI-Corr includes different filtering algorithms to remove such errors (Ayoub et al., 2009). The accuracy of the final result largely depends on the precise co-registration/orthorectification, absence of cloud and shadow in the image pair.
3.3 Subpixel correlation of satellite images
The algorithm available in COSI‒Corr works on single-band greyscale images. In this study, the panchromatic band (PAN; band 8; 15 m; 0.50‒0.66 μm for OLI and 0.52‒0.9 μm for ETM+) for Landsat OLI and ETM+ and green band (band 2; 30 m; 0.52‒0.6 μm) for Landsat TM were processed. The green band of Landsat TM was chosen as its wavelength is close to the PAN band of Landsat OLI and ETM+. COSI‒Corr provides two correlation algorithms: (1) frequency correlation and (2) statistical correlation. The frequency correlation is Fourier-based and is more accurate than the statistical one (Ayoub et al. 2009). This correlation is more sensitive to noise and is recommended for optical images of good quality. The statistical correlation maximizes the absolute value of the correlation coefficient and is coarser but more robust than the frequency correlation algorithm. Its use is recommended for correlating noisy optical images that provided terrible results with the frequency correlation. In this study, the frequency correlation algorithm was used for glacier displacement measurement. The input images should be co-registered as accurately as possible as the correlation algorithm is sensitive to image misregistration errors. Besides, image noise due to stripping, sensor noise, and illumination differences caused by cloud and shadow will have a deteriorating effect on displacement results (Scherler et al., 2008). Considering the computational efficiency, the phase correlation algorithm was used in the frequency domain in COSI-Corr. In addition, the Fast Fourier Transformation (FFT) algorithm embedded in the COSI-Corr improves the traditional peak correlation methods, and its accuracy can reach up to 1/50 pixel (Avouac et al. 2006). The frequency correlation algorithm in COSI-Corr requires several initial parameter settings: (i) window size ‒ the size in pixels that will be correlated; (ii) step ‒ determines the step-in pixels between two sliding windows; (iii) robustness iteration ‒ the number of times per measurement the frequency mask should be recomputed; and (iv) mask threshold ‒ the masking of the frequency according to the amplitude of the log-cross spectrum. A detailed methodology of flowchart used in this study is presented in Fig. 2.
Here, the initial and final window size was set as 64 pixels * 32 pixels, step as 4 pixels * 4 pixels, which resulted in a ground resolution of 60 m * 60 m. In theory, the window size should be at least twice the expected displacements. Appropriate settings of these values need a priori field-based knowledge of the size of the displacements, which is not available. Based on multiple experimental studies, it has been shown that the window size of 64*32 pixels delivers promising results in terms of slowing movies debris-covered ice (Scherler et al., 2008) or crustal deformation field caused by an earthquake (Avouac et al. 2006). To reduce the effect of high-frequency noise on the accuracy of correlation results, the frequency mask threshold and robustness iteration were set as 0.9 and 2, respectively. To maintain the consistency of results, these parameter settings were used for the entire dataset.
The correlation algorithm in the COSI-Corr produces three output images. The first two images provide the 2D displacements in terms of the east-west displacement (EWD; east positive) and the north-south displacement (NSD; south positive), both expressed in meters (Fig. 2). The ground displacement (D) was measured by combining these two images by the square root of the sum of the squared displacements (Eq. 1). The surface velocity (V) was calculated by dividing the displacement with a time interval (t) between two images (Eq. 2). The signal-to-noise (SNR) image helps to assess the quality of the displacement images.
3.4 Post-processing procedures
Uncertainty in velocity measurements mainly resulted from the poor co-registration and orthorectification of images, presence of cloud cover, and shadowing effect in any optical images. In this study, terrain corrected and orthorectified (Level L1T) Landsat datasets downloaded from USGS Earth Explorer were directly used for surface velocity measurements. Few studies on the error sources of deformation/displacement field were obtained by the cross-correlation of Landsat images (Ding et al. 2016; Jeong and Howat, 2015). Ding et al. (2016) comprehensively reported that error sources in the velocity field derived from the Landsat 8 images mainly include temporal decorrelation noises, wavelength orbital error, satellite artifacts, attitude jitter distortions, and topography dependent artifacts. Rigorous post-processing was performed for the entire image pairs to remove potential errors in this study (nline resource).
3.4.1 Decorrelation noises
The radiometric properties of the Earth surface vary from features to features which produce decorrelation noises in the optical images. Invariant radiation can lead to discreet texture features, and severe radiometric changes commonly cause texture defects. Thus, it is difficult to achieve accurate correlation based on window matching. As a result, feature tracking based on pixel correlation produces inaccurate random measurements and low SNR values in corresponding regions. Erroneous measurements were masked out based on low SNR (<0.9) to minimize the effect of decorrelation noise using the discard tool available in COSI-Corr. Further, regions with cloud cover and shadow were manually checked, and resulted noises were removed based on replace/discard tool in COSI-Corr (i.e., NSD: -200 to 200 and ESD: -200 to 200). This value was chosen based on the published maximum velocity for the western Himalayan glaciers. This threshold magnitude was kept constant for the entire dataset. This resulted in data gaps in the correlated images.
3.4.2 Orbital error correction
Inaccurate geometric corrections and orthorectification can cause significant signals of linear ramp in the image deformation field. Although Landsat 8 data provided by USGS is geometrically calibrated and terrain corrected, it has some linear ramps in correlated images (Ding et al. 2016). This type of linear ramp generally increases with a short temporal timespan between two correlated images. In image preprocessing, this type of error can be eliminated by rigorous geometric correction to the reference image if the state vector of Landsat 8 images is available. In this present study, linear ramps are nearly absent due to the higher temporal span (more than one year) among image pairs. Although, a first-order polynomial fitting was performed using the detrending tool to remove such error because Landsat does not provide a state vector (Online resource). This method was also used in InSAR interferogram-derived deformation field (Moon et al. 2014; Mouginot et al. 2015).
3.4.3 Satellite stripe artifacts correction
The misalignment of Charge Coupled Device (CCD) arrays is a common problem of the push-broom scanner, producing stripe artifacts in image displacement (Leprince et al. 2008; Scherler et al. 2008; Ding et al. 2016). General orthorectification cannot eliminate this type of error. Further, satellite attitude variation produces cyclical distortions along the satellite flight direction, which are roll variations on ESD components and pitch variations on the NSD component (Leprince et al. 2008; Scherler et al. 2008). Such systematic distortions can be removed using destripe tool available in COSI-Corr (Ayoub et al. 2009). The possibility to remove these artifacts mainly depends on the fraction of visible, stable ground, i.e., off glacier area with no ice movement. Generally, the higher the amount of stable and visible ground, the better the possibilities of removing attitude effects. These corrections can be easily implemented for earthquake deformation measurements (Leprince et al. 2008; Ding et al. 2016). In this study, an entire Landsat scene (path 147 and row 37) was processed for different timespan. The scene area is mainly covered with high-rugged terrain with numerous active high-altitude glaciers. Stable terrain is manly confined to very small patches along the river channel and confluence (i.e., terraces). As a result, it is very difficult to apply destripe to the entire images. Here, destripe was performed in very small subset of correlated images which does not provided any significant corrections (Online resource). Scherler et al. (2008) also reported the similar limitations for destriping for the two high altitude regions of the Himalaya.
3.4.4 Directional filter
At the very fast step, removal of decorrelation noise based on SNR threshold and magnitude threshold of NSD and EWS mostly exclude miscorrelated patches. However, this does not exclude all miscorrelated points (i.e., opposite directional velocity field to the glacier flow direction). Studies have reported that a simple directional filter is very efficient in getting rid of this type of miscorrelation (Kääb et al. 2005; Scherler et al. 2008). This was done by defining the flow direction from flow features on the glacier surface in correlated images and allowing for some deviation, e.g., up to 20° (Scherler et al. 2008). Here, downward glacier flowlines were manually checked and digitized from Google Earth. The velocity vector fields were overlaid on the glacier flow field with the high-resolution optical image in the background to identify the misoriented pixels (Fig. 3). The miscorrelated pixels were manually checked and removed from the final velocity maps. Finally, the velocity maps were clean up and smoothed using median filter 3*3 pixels (Mouginot et al. 2012; Copland et al. 2009b; Millan et al. 2019).
3.5 Accuracy assessment
Accuracy assessments of the final velocity products were evaluated in different ways to access the quality of datasets with regard to the direction, magnitude, and gradient (Scherler et al. 2008). These include (i) a test of displacement based on the movement of stable terrain (i.e., off glacier area) and (ii) a test of the displacement direction based on streamlines along the glacier flow. The mean and standard deviation of stable terrain displacement were checked for each post-processing step (Online resource). The error of the observed velocities in the off glacier area was calculated based on the following formula (Berthier et al. 2003):
where, is the error of the off-glacier velocity, is the mean stack velocity of random points. SE is the standard error of the mean velocity, which can be measured using the following equation:
where is the standard deviation of the off-glacier area velocity. is the number of measurements in the off-glacier area. In the present study, 10000 random points were derived from the off-glacier area. Results show that mean displacement in the stable terrain is <5 m (Online resource). Flow streamlines of the Samudra Tapu glacier were manually digited from Google Earth images to check the consistency of the velocity vectors. Velocity vectors show consistent results with streamlines, indicating higher accuracy of the velocity products (Fig. 3).
3.6 Glacier classifications and topographic parameter extractions
Glaciers >5km² in size in the CBM were chosen for the present analysis. These glaciers vary significantly in terms of their morphological characteristics (debris cover, lakes), topographic settings (slope, orientation, elevation), and terminus properties (Table). Based on the amount of debris cover and lakes near the terminus, studied glaciers were categorized into three groups. Four glaciers are characterized by pro-glacial lakes near the terminus, thus categorized as lake-terminating glaciers (Table 2). The rest of the 11 glaciers are categorized based on their percentage of debris cover to total ice area (Table 2). Three glaciers are classified as debris-covered glaciers, having >30% of the area covered by debris. Eight glaciers are classified as clean glaciers (Table 2). Several topographic parameters (i.e., slope, aspect, elevation) were derived from ASTER GDEM v2. The hypsometric index was calculated based on the approach suggested by a previous study (Jiskoot et al. 2009).