Glacier Surface Velocities in the Chandrabhaga Massif, Western Himalaya (India) Derived Using COSI-Corr From Landsat Images

Himalayan glaciers act as a storehouse of freshwater outside the polar region and sustain the �ow of several Asia’s major river systems, such as the Ganges, Indus, and Brahmaputra. Glaciers in High Mountain Asia are likely to be affected by climate change, posing a threat to the future water supply. Frequent mapping and monitoring of these glaciers are in urgent need to assess the future water storage and glacier-climate interaction. In this context, mapping and monitoring of surface velocity can help to infer the health of a glacier. However, systematic assessment of glacier surface velocity is limited to fewer basins or individual glaciers in the Himalaya. Here, I have characterized the spatial and temporal velocity variations of glaciers in the Chandrabhaga massif (CBM) using Landsat time-series data spanning nearly the last three decades (1992-2019). Rigorous post-processing was performed to improve the accuracy of remote sensing derived velocity products. Glaciers showed heterogeneous spatial and temporal velocity patterns based on morphological characteristics and topographical settings. Heavily debris-covered tongues showed nearly stagnant �ow (<10 m/yr) while maximum velocity (>60 m/yr) was observed for clean glaciers with steep gradients and crevasses. Surface velocity near the terminus of the lake terminating glaciers was almost double than land terminating debris-covered tongues (32.5 m/yr vs. 12.6 m/yr). An increasing trend in surface velocity is attributed to the rising trend in air temperature in western Himalaya.


Introduction
Hindukush-Karakorum-Himalaya (HKH) region contains thousands of glaciers outside the polar region.Millions of livelihoods in south and southeast Asia are dependent on the meltwater from high Asian mountain glaciers.In recent decades, Himalayan glaciers are in the main focus of debate as these glaciers are more sensitive to climate change and potential water availability for the future.In this context, frequent mapping and monitoring of multiple glaciological records (i.e., surface area, elevation, velocity change) are in urgent need to evaluate the future water availability and sustenance of glaciers on the regional scale.In these contexts, surface velocity can help to infer the health status of a glacier.
Although eld-based monitoring of surface velocity of individuals' glaciers in the Himalaya is nearly impossible considering the remoteness and highly rugged topography of the region.In this scenario, the analysis of vast freely available remote sensing datasets could help estimate the surface velocity at a broader regional scale.
Currently, three remote sensing-based methods are commonly employed to derive glacier-surface velocities: interferometry of synthetic aperture radar (SAR) imagery (Prats et al. 2009;Yu et al. 2010), SAR feature tracking techniques (Moon et al. 2012;Vijay et al. 2019), and cross-correlation of optical satellite images (Scherler et al. 2008;Dehecq et al. 2015;Fahnestock et al. 2016;Gardner et al. 2018).SAR interferometry (InSAR) is a widely used technique for deformation and velocity mapping.Although, it has several limitations in highly rugged terrain like the Himalayan region (Luckman et al. 2007).This technique may achieve high accuracies but require good coherence between the images, which sometimes is very di cult to nd due to the modi cation (i.e., melting, snowfall, avalanche) of glacier surface features (Joughin et al. 2002).Several studies have shown that image coherence is typically constrained to periods of 1, 3, or 6 days (Joughin et al. 2002).Thus, interferometry-derived velocity may be representative for the short observation period, and extrapolation of annual velocity is di cult.Offset tracking of SAR imagery is similar to the cross-correlation of the optical imagery.This can be done either using the intensity or coherence of the complex radar images (Moon et al. 2014;Mouginot et al. 2015).Compared to interferometry, feature tracking techniques using SAR images are more useful for measuring ow velocities over extended periods (Moon and Joughin 2008).The visibility of the target features is affected due to oblique viewing SAR images in such rugged terrain conditions.Further, a high incidence angle requires precise Digital Elevation Models (DEMs) to orthorectify images (Joughin et al. 2002).
Optical image correlation is another promising technique to infer displacement eld, which has been widely used in geology and earth science in recent decades.The principle involved in this technique is that two images acquired at different times are correlated to nd out the shift in the position of any moving object.Surface velocity elds of glaciers and other moving ice bodies using optical satellite images have been studied since the mid-1980s using manual tracking of features.The detail and accuracy of the measurements are limited mainly by the ground resolution of the sensor, and by the ability to precisely co-register images acquired at different dates (Leprince et al. 2008;Scherler et al. 2008).Further errors may arise from changes in the satellite attitude during scanning of the images (Scherler et al. 2008;Ding et al. 2016), and an inaccurate DEM during orthorecti cation using a rigorous model.A principal drawback of optical imagery is the dependency on cloud-free conditions.Several international initiatives are set up to monitor the surface speed of glacial based on the cross-correlation of optical images at global scale.The global land ice velocity extractions from Landsat (GoLIVE) and the intermission time series of land ice velocity and elevation (ITSLIVE) are the two main examples of such initiative.The GoLIVE and ITSLIVE datasets are primarily generated for monitoring ice sheets and large glaciers in the pole with minimum glaciers size threshold of 5 km² and a spatial resolution of 300 m (Fahnestock et al. 2016).Monitoring of heavily debris-covered small Himalayan glaciers from these datasets are di cult as reported same for the European alps (Millan et al. 2019).In the present study, an attempt has been made to map the surface velocity of glaciers at higher spatial resolution in the CBM in western Himalaya during last three decades (1992-2019).2021) also measured ice velocity for four glaciers (Samudra Tapu, Sutri Dhaka, Batal, and Gepang Gath) based on the eld mapping using stakes position.Shukla and Garg (2019) estimated the surface velocity of Panchi Nala glaciers in the Bhaga basin using Landsat images.Further, the potential in uence of several glacial morphological and topographic factors was also evaluated on surface velocity in the upper Chandra basin.Existing studies are limited by the shorter temporal scale and selected glaciers.However, no studies have been published on surface velocity in the Chandrabhaga massif addressing spatial ow patterns in association with pro-glacial lakes, debris cover, and other topographic parameters.Thus, the main goals of this study are (i) estimate the surface velocity based on the cross-correlation of Landsat images during the last three decades (1992-2019), and (ii) evaluate the spatial ow patterns of lake terminating, heavily debris-covered, and clean glaciers separately in the CBM, western Himalaya.

Study Area
The CBM lies in the greater Himalaya range of western Himalaya (Fig. 1).Part of the Bhaga and upper Chandra basin combinedly form the Chandrabhaga massif.Mulkila is the highest peak in this massif with an altitude of 6517 m a.s.l.The CBM consists of some of the largest glaciers in the Chenab basin (Samudra Tapu, Mulkila, Milang).Glaciers of the CBM are characterized by heavy debris-covered ablation zone and pro-glacial lakes near the terminus (Fig. 1).Samudra Tapu is the largest glacier in the CBM, with a surface area of ∼80 km².The Pir-Panjal range surrounds the CBM in the south and trans-Himalaya in the north.Mid-latitude westerlies in uence the climate of this region in winter and the South Asian monsoon in the summer months.Therefore, the typical geographical setting of this region makes it ideal for studying the response of glaciers to climate change.
3 Data And Methods

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.

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 eld 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, orthorecti cation, 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 nal result (Scherler et al. 2008;Ding et al. 2016).COSI-Corr includes different ltering algorithms to remove such errors (Ayoub et al., 2009).The accuracy of the nal result largely depends on the precise co-registration/orthorecti cation, absence of cloud and shadow in the image pair.

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 coe cient 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 coregistered 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 e ciency, 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 sizethe size in pixels that will be correlated; (ii) stepdetermines the step-in pixels between two sliding windows; (iii) robustness iterationthe number of times per measurement the frequency mask should be recomputed; and (iv) mask thresholdthe masking of the frequency according to the amplitude of the log-cross spectrum.A detailed methodology of owchart used in this study is presented in Fig. 2.
Here, the initial and nal 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 eld-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 debriscovered ice (Scherler et al., 2008) or crustal deformation eld 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 rst 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.

Post-processing procedures
Uncertainty in velocity measurements mainly resulted from the poor co-registration and orthorecti cation of images, presence of cloud cover, and shadowing effect in any optical images.In this study, terrain corrected and orthorecti ed (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 eld 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 eld 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).

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 di cult 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.

Orbital error correction
Inaccurate geometric corrections and orthorecti cation can cause signi cant signals of linear ramp in the image deformation eld.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 rst-order polynomial tting 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 eld (Moon et al. 2014;Mouginot et al. 2015).

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 orthorecti cation cannot eliminate this type of error.Further, satellite attitude variation produces cyclical distortions along the satellite ight 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 con ned to very small patches along the river channel and con uence (i.e., terraces).As a result, it is very di cult to apply destripe to the entire images.Here, destripe was performed in very small subset of correlated images which does not provided any signi cant corrections (Online resource).Scherler et al. ( 2008) also reported the similar limitations for destriping for the two high altitude regions of the Himalaya.

Directional lter
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 eld to the glacier ow direction).Studies have reported that a simple directional lter is very e cient in getting rid of this type of miscorrelation (Kääb et al. 2005;Scherler et al. 2008).This was done by de ning the ow direction from ow 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 owlines were manually checked and digitized from Google Earth.The velocity vector elds were overlaid on the glacier ow eld 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 nal velocity maps.Finally, the velocity maps were clean up and smoothed using median lter

Accuracy assessment
Accuracy assessments of the nal 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 ow.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).

Glacier classi cations and topographic parameter extractions
Glaciers >5km² in size in the CBM were chosen for the present analysis.These glaciers vary signi cantly 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 classi ed as debris-covered glaciers, having >30% of the area covered by debris.Eight glaciers are classi ed 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).

Results
Signi cant heterogeneous ow patterns were observed for the studied glaciers in the CBM during 1992-2019 (Fig. 4).The estimated velocity of each image pairs was analyzed in terms of spatial and temporal scale during the observation periods.Spatial velocity variations were analyzed along pro les ranging from thick debris-covered ice near terminus to clean ice in up-glacier (Fig. 4).Annual velocity products are primarily for the post-melting seasons (i.e., between August and September).Surface velocities were analyzed separately for three types of glaciers to identify the

Debris-covered glaciers
Three glaciers (Milnag, Sonapani, and G06) were classi ed as debris-covered glaciers.The surface velocity of the debris-covered tongue of the Sonapani glacier is 11 m/yr, while debris-free up-glacier showed a mean surface speed of 23.12 m/yr (Table 3; Fig. 5a).The surface velocity of the Sonapani glaciers increased at a signi cant rate during the observation period (Fig. 8b).The surface velocity of the Milang glacier varies systematically along the length of the ow (Fig. 5b).The mean velocity of the Milang glacier is 17.87 m/yr, while the heavily debris-covered tongue (up to 7.5 km from terminus) is nearly stagnant with a mean velocity of 12.87 m/yr (Table 3).Debris free up glacier zone showed a higher velocity (25.65 m/yr).During 1992-2019, the Milang glacier showed a temporal uctuation, although the overall velocity trend is nearly stationary (Fig. 8d).The mean velocity of the G06 glacier is 16.29 m/yr, while overall velocity is decreasing during 1992-2019 (Table 3; Fig. 8f).

Clean glaciers
Eight glaciers in the CBM were categorized as clean glaciers, and their ow patterns were analyzed separately (Table 3; Fig. 6).The surface velocity of the clean glaciers was found to be much higher than the debris-covered ones in the CBM, although heterogeneity in surface velocity was observed within glaciers (Fig. 6).The mean velocity of the G01 and G12 glaciers is ∼26.5 m/yr, while the G10 and G15 glaciers showed a higher mean velocity of ∼34 m/yr (Table 3; Fig. 6).Signi cant uctuations are observed in the up-glaciers for clean glaciers (Fig. 6).Temporal uctuations of clean glaciers are presented in Fig. 8.Most of the clean glaciers showed an increasing trend in surface velocity except for the G12 and G14 glaciers (Table 3: Fig. 8).

Lake-terminating glaciers
Four glaciers (Gepang Gath, G07, Panchi Nala, and Samudra Tapu) in the CBM are terminating at lakes.Except for Samudra Tapu, the ablation zone of the rest of the glaciers is covered by heavy debris.The mean velocity of the Gepang Gath glacier is 20.3 m/yr, while the thick debris-covered ablation zone showed a lower velocity of 13.16 m/yr (Table 3; Fig. 7a).Although, the frontal part of the Gepang Gath glaciers showed a maximum velocity of ∼32.6 m/yr.The Gepang Gath showed an insigni cantly increasing velocity trend during 1992-2019 (Fig. 8c).The debris-covered ablation zone (up to 2.5 km from terminus) of the G07 glacier showed lower velocity (∼6.9 m/yr) as compared to clean ice (13.3 m/yr; Fig. 7b).The debris-free accumulation zone of the Panchi Nala glacier showed nearly three-time higher velocity (22.4 m/yr) than the thick debris-covered ablation zone (∼7.7 m/yr) with a mean velocity of 17.6 m/yr (Table 3; Fig. 7c).The overall velocity of the Panchi Nala glacier showed an insigni cant decreasing trend (Fig. 8h).The mean velocity of the Samudra Tapu glacier is 42.3 m/yr, while in the lower ablation zone (up to 3 km from terminus), velocity is observed to be much lower (15.6 m/yr) than debris-free upglacier zone between 3 and 15 km along the central ow line (∼45.7 m/yr; Table 3; Fig. 7d).Although, signi cant uctuation is observed along the ow line (Fig. 7d).In general, the surface velocity of the Samudra Tapu glacier is increasing a signi cant rate during 1992-2019 (Fig. 8m).

Comparison of surface velocity with previous studies in Chandra basin
In the present study, spatial and temporal velocity variations of 15 glaciers in the CBM were estimated during the last three decades.Few studies have published on surface velocity, mainly limited to the individual glaciers in the upper Chandra basin (Tiwari et 2021) measured surface speed for Samudra Tapu (64.3 m/yr), Sutri Dhaka (52.6 m/yr), Batal (6.2 m/yr), and Gepenag Gath (26.5 m/yr) glaciers in the Chandra basin.In the present study, lower velocity was observed for the Samudra Tapu (42.5 m/yr), Sutri Dhaka (34.2 m/yr), and Gepang Gath (20.3 m/yr) glaciers.This could be attributed to the difference in observation periods, average velocity from long-term datasets, the spatial distribution of point-based velocity extractions.

Heterogeneous glacier motion
Glaciers ow down the valley due to the gravitational force; therefore, it mostly depends on ice thickness and surface gradient.Few studies have examined spatial velocity patterns and their controlling factors in western Himalaya (Scherler et al. 2011;Quincey et al. 2009;Bhushan et al. 2018).Studies have reported that catchment topography, hypsometric ice distribution, and surface gradient play an essential role in controlling ice movement (Quincey et al. 2009).In this study, an attempt has been made to evaluate the variable spatial ow patterns of lake-terminating, debris-covered, and debris-free glaciers.Possible reasons for heterogeneous spatial and temporal ow patterns are discussed below.
In the CBM, the surface velocity of debris-covered ice (< 10 m/yr) is found to be much lower than clean ice (> 35 m/yr).For example, the thick debris-covered ablation zone of the Milang glacier is nearly stagnant (< 12.5 m/yr), while nearly debris-free Samudra Tapu showed a maximum velocity of > 80 m/yr near lower accumulation zone.Glaciers covered with debris in the ablation zone in the CBM are found to have peak velocity shifted to the upper part, leading to a more asymmetrical velocity pro le (Fig. 5).Thus, the surface gradient should be greatest in the upper part compared to the lowered debris-covered zone.This pattern is prominent for the clean glaciers having top-heavy hypsometric ice distribution (Fig. 6).Previous studies have reported that a thin debris layer enhances the underneath ice temperature while thick debris cover act as a barrier of heat transfer, thus acts as an insulator of temperature (Dobhal et al. 2013).Earlier studies in the western Himalaya have reported that the area change rate of the debris-covered glaciers is less than the clean ice (Dobhal et al. 2013).However, debris-covered ice showed much higher down wasting than clean ice in the upper Chandra basin (Vijay and Braun 2016;Bhushan et al. 2018).Thus, the debris-covered glaciers in the CBM are prone to stagnancy, enhancing the formation of supraglacial ponds and further down wasting, while fronts remain relatively stable.Moreover, surface ice morphology (i.e., presence of crevasses, break of slope) controls the velocity patterns in the CBM.In the present study, the zone of crevasses with a break of slope showed maximum velocity (Fig. 8).To validate this pattern, the zone of crevasses with the break of the slope were identi ed for the studied glaciers using the high-resolution 3-D image in Google Earth (Fig. 9).An example of higher velocity associated with crevasses and break of slope observed in Mulkila glacier is presented in Fig. 9. Velocity near the break of the slope was observed at > 100 m/yr, indicating that sudden changes in surface gradient cause internal deformation on ice, thus resulted in crevasses and higher motion.In contrast, the gentle sloping debris-covered ablation zone showed a nearly homogeneous ow pattern (Fig. 9).
Previous studies have investigated the evolution and dynamics of pro-glacial lakes in the CBM based on remote sensing (Shukla and Garg 2019; Kaushik et al. 2020).Four lake terminating glaciers in the CBM showed distinct ow patterns.Out of four lake terminating glaciers, three of them (Gepang Gath, G07, and Panchi Nala glaciers) are characterized by debris-covered ablation zone (Fig. 1).Samudra Tapu glacier is nearly debris-free except for some patches of debris cover at valley edges (Fig. 1).The signi cant movement was observed at the transition zone between lake and ice for the Gepang Gath glacier (Fig. 10).This is probably because of the higher rate of frontal retreat due to the calving of ice into lakes, thus resulted in higher velocity.Further, future growth and upward shift of lake level may enhance the penetration of water into underneath bedrock, thus increase in basal sliding and surface movement.
Similar results have been reported for lakes terminating glaciers across the Himalaya (Pronk et al. 2021).Samudra Tapu glacier showed nearly stagnancy at terminus despite having a pro-glacial lake (Fig. 10).This is probably because of less contact of ice with lake water, resulting in lower calving of the front (Fig. 10).Previous studies have reported that the rate of growth of the lake in the Gepang Gath glacier is higher than in Samudra Tapu (Kaushik et al. 2020).Furthur, bottom-heavy hypsometric ice distribution of the Gepang Gath glacier may bring more ice to a lower altitude than the top-heavy Samudra Tapu glacier.
G07 and Panchi Nala glaciers also showed velocity patterns similar to Samudra Tapu.Thus the heterogeneity in surface velocity for lake and land terminating glaciers in the CBM could be attributed to the calving effect of ice, surface gradient, hypsometric ice distribution, and the morphology of the terminus.
Most of the studied glaciers (11 out of 15) showed an increasing trend in temporal velocity patterns during 1993-2019 (Fig. 8).This could be attributed to the interplay between the increasing trend in temperature in western Himalaya (Das and Sharma 2019; Mukherjee et al. 2018) and the evolution of subglacial hydrology (Satyabala 2016).Mukherjee et al. (2018) reported a signi cant upward shift of temperature after 2000 in this part of the Himalaya.An increase in temperature may enhance the meltwater discharge accompanied by a higher rate and quantity of sediment ux in glaciers exhibiting basal sliding due to changes in bed hydrology (Miles et al. 2019).The surface velocity of the Mulkila and Samudra Tapu glaciers increased at a signi cant rate (r = 0.7).Based on the analysis of grided climatic datasets, Das and Sharma (2019) reported that mean annual temperature is increasing at a nominal rate (0.0078° C per year) while winter temperature is increased signi cant (∼1.3°C) between 1948 and 2016 in this part of the western Himalaya.So, it can be noted that rising air temperature combined with several morphological and topographical factors controls the ow speed of glaciers in the CBM.

Conclusion
In the present study, I have investigated the spatial and temporal ow patterns of glaciers in the Chandrabhaga Massif, western Himalaya.Glacier velocities were derived from the cross-correlation of Landsat images spanning during 1992-2019.Glaciers in the Chandrabhaga Massif showed heterogeneous ow patterns based on their morphological and topographical characteristics.Debriscovered glaciers with a gentle slope and top-heavy hypsometric ice distribution showed nearly stagnant condition (< 10 m/yr) than clean glaciers.Clean glaciers showed maximum velocity (> 80 m/yr) near the zone of crevasses and break of slope.Surface velocity near the front of the lake terminating glaciers was nearly double than land terminating debris-covered glaciers (32.3 m/yr vs. 12.5 m/yr).Calving of ice into glacial lake and intrusion of water into bedrock increase the basal sliding thus resulted in a higher velocity near fronts of the lake terminating glaciers.Future studies should focus on the higher spatial and temporal scale remote sensing-based mapping and eld validation of surface velocity of glaciers having different morphological (i.e., debris cover, proglacial lakes) properties.

Declarations
Several studies have been published on monitoring surface area change (Das and Sharma 2019) and elevation change (Vijay and Braun 2016; Mukherjee et al. 2018) of glaciers in the upper Chandra basin.Only a few studies have been published on surface velocity estimation for the upper Chandra basin (Sahu and Gupta 2019; Patel et al. 2021).In the upper Chandra basin, most of the velocity measurements are concentrated for Bara Shigri and Chotta Shigri glaciers (Yellala et al. 2019; Tiwari et al. 2014; Patel et al. 2017).Sahu and Gupta (2019) studied the surface velocity of four glaciers (Bara Shigri, Chhota Shigri, Gepang Gath, and Hamta) during 2000-20 based on Landsat images.Another study by Patel et al. (