Anthropogenic climate change drives melting of glaciers in the Himalaya

The Himalayan glaciers provide water to a large population in south Asia for a variety of purposes and ecosystem services. As a result, regional monitoring of glacier melting and identification of the drivers are important for understanding and predicting future cryospheric melting trends. Using multi-date satellite images from 2000 to 2020, we investigated the shrinkage, snout retreat, thickness changes, mass loss and velocity changes of 77 glaciers in the Drass basin, western Himalaya, India. During this period, the total glacier cover has shrunk by 5.31 ± 0.33 km2. The snout retreat ranged from 30 to 430 m (mean 155 ± 9.58 m). Debris cover had a significant impact on glacier melting, with clean glaciers losing ~ 5% more than debris-covered glaciers (~ 2%). The average thickness change and mass loss of glacier have been − 1.27 ± 0.37 and − 1.08 ± 0.31 m w.e.a−1, respectively. Because of the continuous melting and the consequent mass loss, average glacier velocity has reduced from 21.35 ± 3.3 m a−1 in 2000 to 16.68 ± 1.9 m a−1 by 2020. During the observation period, the concentration of greenhouse gases (GHGs), black carbon (BC) and other pollutants from vehicular traffic near the glaciers increased significantly. Increasing temperatures, caused by a significant increase in GHGs, black carbon and other pollutants in the atmosphere, are driving glacier melting in the study area. If the current trend continues in the future, the Himalayan glaciers may disappear entirely, having a significant impact on regional water supplies, hydrological processes, ecosystem services and transboundary water sharing.


Introduction
Cryosphere is an important component of the Earth's natural system, especially in the Himalayas (Lemke et al. 2007; Barry and Gan 2011). Although mountain glaciers constitute only ~ 3% of the global glacial area (Arendt et al. 2002), the need for precise areal and volumetric estimations of mountain glaciers is well established, as they serve as important freshwater reserves and contribute significantly to meeting the water demands of a large population of about 1.4 billion in the South Asia and neighbouring regions (Immerzeel et al. 2010). Because of their sensitivity to climate change, melting glaciers in the Himalayas contribute to the sea-level rise and an increase in the frequency of glacier hazards such as GLOFs (glacial lake outburst floods), avalanches, landslides and permafrost creep due to their sensitivity to climate change (Kääb et al. 2005;Snehmani et al. 2014;Thompson et al. 2020;Khan et al. 2021). Continuous monitoring of glacier recession and dynamics is therefore critical for observing the direct impacts of climate change on water security, water supplies, future sea levels and glacier-related disasters.
One of the most noticeable impacts of climate change in the Himalaya is the significant retreat and thinning of glacier (Shrestha and Aryal 2011;Ali et al. 2021;Dixit et al. 2021). The observed retreat and thinning, however, are marked by heterogeneity across the Himalaya (Abdullah et al. 2020;Farinotti et al. 2020). The influence of climate change, glacier morphology and local topography, as well as other influencing factors, has been extensively investigated in order to understand and explain the heterogeneous glacier response (Fujita and Nuimura 2011;Garg et al. 2017;Murtaza and Romshoo 2017;Abdullah et al. 2020). Several studies have estimated two-dimensional (i.e. glacier area, length or snout) glacier changes (Cogley 2011;Zemp et al. 2015;Romshoo et al. 2020a), which represents an assessment of indirect, delayed and filtered response of changing climate (Scherler et al. 2011;Gardelle et al. 2013;Bhattacharya et al. 2016).
Most of the studies have reported the general state of glacier recession over the last few decades across the Himalaya, with the exception of some areas in the Karakoram where glaciers have advanced, surged, behaved irregularly and even retreated (Bhambri et al. 2011;Bolch et al. 2012;Paul 2015;Murtaza et al. 2021;Nüsser and Schmidt 2021).
However, studying other parameters such as glacier surface ice velocity, elevation and mass balance changes, in addition to the two-dimensional changes, to assess overall dynamic response of glaciers to changing climate is now widely promoted (Gardelle et al. 2013;Vijay and Braun 2016;. The changes in glacier mass balance and thickness demonstrate the direct influence of climate change on glacier melting and the resulting changes in the streamflow. Similarly, glacier velocity can be used to assess changes in glacier ice thicknesses, mass changes, glacier advance and retreat and glacier ice flow regimes (Tiwari et al. 2014;Singh et al. 2021). Unfortunately, only a few mass balance, glacier thickness and glacier velocity studies at a basin or at range scale have been conducted in the Himalayas (Dehecq et al. 2019;Gardelle et al. 2013) limiting our understanding of glacier response to climate change in the region. The insufficient and sometimes lack of observations of the driving factors of glacier melting related to climatic, hydrological and atmospheric processes complicates the understanding of cryospheric melting in the data-scarce Himalaya (Nüsser and Schmidt 2021).
Despite the fact that the role of climatic and topographical factors in glacier retreat and thinning has been scientifically well studied and understood, however, the effects of other key factors such as greenhouse gases (GHGs), black carbon concentration (BC) and other pollutants from vehicular traffic that get deposited in glacier environments are still in their infancy particularly in the Himalayas. Preliminary findings of some studies aimed at understanding the impacts of these parameters on glacier melting (Luo et al. 2004;Ramanathan and Feng 2009; suggest that an increase in the concentration of these parameters causes warming of the environment, which in turn has strong influence on glacier melting. Glaciers in the study have receded in accordance with the global glacier recession trend since the Late Holocene (Koul et al. 2016;Rashid et al. 2021). Although a few studies have estimated the recession of glaciers in this region in terms of area loss, snout retreat and glacier length changes, with a particular focus on the Machoi Glacier (Koul et al. 2016;Rashid et al. 2021;Pall et al. 2019), the overall glacier response and impacts of climate change in the study area have not been thoroughly investigated, and thus, the role of controlling factors, particularly impacts of GHGs, BC and other vehicular pollutants on the observed glacier recession, remains largely unknown. In order to better understand the influence of controlling factors driving glacier recession, this study compared the observed glacier shrinkage, snout retreat, changes in glacier thickness, mass loss and glacier velocity, based on Landsat images, topographic maps, SRTM DEM and TerraSAR-X/TanDEM-X data, with the observed changes in the greenhouse gas emissions, BC concentration, pollutants, temperature and precipitation trends obtained from various sources.

Study area
Drass, the gateway to the cold desert valley of Ladakh and the world's second coldest inhabited region, is located in the north-western Himalayas, almost 140 km from Srinagar, Kashmir (Fig. 1). The Drass Valley, which stretches from 34.10° to 34.45°N and 75.10° to 76.05°E with an elevation range of 2705-5942 m, begins at the base of the Zojila Pass, on the way to Ladakh. The Greater Himalaya and the Drass Mountains encircle the Drass Valley, which has cold semiarid type of climate receiving precipitation in the form of snowfall during the winter from the western disturbances and rainfall during the summer from both the westerlies and southern jet stream movements (Koul et al. 2016;Romshoo et al. 2020b). Despite the fact that glaciers are distributed throughout the region, the majority of them are concentrated in the southern parts. One of the biggest glaciers in the study area is the Machoi Glacier (G35), which is easily accessible from the Srinagar-Ladakh national highway. In the study area, many glacial landforms such as U-shaped valleys, arêtes, erratic boulders, lateral moraines and terminal moraines have been well preserved. Moraines are the most common landforms that have survived in the region. The abundance of erosional and depositional glacial landforms in the study area (Pall et al. 2019) indicates that the area has been significantly glaciated in the recent past.

Satellite data
In this study, various snow-and cloud-free satellite data from different sensors was employed for determining glacier area changes, snout changes, velocity measurements, glacier thickness and mass change and greenhouse gas changes (Table 1). Landsat ETM + and OLI satellite images, acquired at the end of the ablation season, were chosen for glacier area and snout-retreat mapping, because this period is usually cloud-free and there are almost no snow packs in the valleys adjacent to glaciers (Kääb 2002;Paul et al. 2002). For glacier velocity estimates, panchromatic images of ETM + (1999)(2000) and OLI (2019-2020) were used for glacier velocity estimations. Various glacier morphological and topographical parameters were generated from the Shuttle Radar Topography Mission (SRTM) DEM. SRTM DEM and TanDEM-X data were employed for estimating glacier thickness changes. Furthermore, satellite data from the Atmospheric Infra-Red Sounder (AIRS), Greenhouse

Glacier shrinkage and retreat mapping
For mapping and monitoring spatial and temporal changes of glaciers, remote sensing methods are helpful (Racoviteanu et al. 2008;Bajracharya et al. 2014). Although field-based monitoring is required to study the glacier parameters such as the extent, snout position, glacier depth and mass balance, however, it is extremely challenging, labour-intensive and time-consuming, especially in mountainous regions like the Himalaya, where harsh climatic conditions, rugged terrain, security challenges and remoteness are all present (Gao and Liu 2001;Murtaza and Romshoo 2014). Glacier outlines were mapped using visual image interpretation and digital image processing of the satellite data (Paul et al. 2002;Murtaza et al. 2021). Major changes which usually occur in the ablation area of a glacier are more common than in the accumulation area (Paul et al. 2017;); hence, we limited the glacier-area mapping and change detection analysis to the ablation zone. The presence of some distinguishing features, such as ice wall, proglacial lakes and the emergence of streams from snout, helped to delineate glacier snout (Murtaza and Romshoo 2017). Wherever possible, Google Earth images were used to double-check the accuracy of glacier outlines.

Topographical characterization
This study looked at the impact of local morphology and topographic characteristics on glacier recession. To understand the heterogeneous behaviour of glaciers, influencing factors such glacier size, altitude and slope were studied. Various studies have reported that local topographic parameters have significant impact on glacier recession (Bolch et al. 2012;Murtaza and Romshoo 2017;Romshoo et al. 2020a).
Based on their size, the glaciers were categorized into three classes, namely, small (0.27-1.00 km 2 ), medium (1-5 km 2 ) and large (5-14 km 2 ). The glaciers were categorized into three categories based on their altitude: glaciers between 4600 and 5000 m, glaciers between 5001 and 5400 m and glaciers between 5401 and 5800 m. Similarly, the glaciers were classified into two slope categories based on their slope: glaciers with a gentle slope (< 15°) and glaciers with a steep slope (> 15°).

Debris-cover mapping
The debris-covered parts of the glaciers were visually identified from Landsat images, based on thermal characteristics and were delineated using onscreen manual digitization technique in GIS environment (Dobhal et al. 2013). Only areas with a somewhat continuous debris-cover extent over the glaciers were considered and small isolated patches of debris that appeared sporadically over the glaciers were excluded from the debris-cover mapping. The supra-glacier debriscover changes in space and time have a significant impact on glacier recession. As a result, the glaciers in the study area were categorized into clean glaciers (CGs) with less than 5% debris cover and debris-covered glaciers (DCGs) with more than 5% debris cover in order to better understand the role of debris cover in glacier recession.

Surface ice velocity (SIV) measurement
Because small glacial morphological features are not identifiable from Landsat images, 39 glaciers with surface area more than 1 km 2 were selected for SIV estimations . Feature-tracking method using remotely sensed data has proved useful in monitoring glacier velocity (Scherler et al. 2008;Wu et al. 2020;Zhao et al. 2020). Co-registration of Optically Sensed Images and Correlation (Cosi-Corr) approach, that is based on the repeat use of orthorectified snow-and cloud-free Landsat ETM + and OLI images, was acquired at the end of ablation season between 1999/2000 and 2019/2020, for SIV estimation. Cosi-Corr is a Fourier-based highly sophisticated matching algorithm for the detection of co-seismic displacements (Leprince et al. 2007), which offers sub-pixel accuracy and has been widely used by researchers to measure glacier displacements, primarily for push-broom sensors (Tiwari et al. 2014; Sahu and Gupta 2019; Shukla and Garg 2020). Scherler et al. (2008) provide the SIV details of the Cosi-Corr algorithm, which employs three stages to calculate SIV: orthorectification, coregistration and correlation. The displacement images were checked and corrected for several probable errors and several filters were applied before calculating SIV. Low signal to noise ratio (SNR) values less than (0.9) were initially filtered out to remove the poorly correlated pixels. The SIV values > 200 m a −1 were also removed to exclude the incorrect SIV values caused by poor correlation (Bhattacharya et al. 2016;Bhushan et al. 2018).
Landsat images offer good spatial and temporal coverage (Ding et al. 2016; Shukla and Garg 2020) but they also feature subpixel noise caused by attitude changes (Heid and Kääb 2012;Bhattacharya et al. 2016). The average image-to-image registration precision for ETM + and TM sensors is ~ 5 and ~ 6 m, respectively (Storey and Choate 2004;Bhattacharya et al. 2016), which is difficult to model due to their whisk-broom nature. However, this precision is justified for glaciological studies (Heid and Kääb 2012).

Glacier mass changes
The elevation differences in the SRTM DEM (2000) and TanDEM-X (2012TanDEM-X ( -2015 obtained ~ 12 year apart were used to assess glacier thickness change. The universal co-registration algorithm was used to remove the horizontal and vertical offset between the two DEMs (Nuth and Kääb, 2011). The elevation difference (dH/dT) map over the glaciated terrain was generated by comparing the co-registered DEMs. Before generating the DEM difference map at pixel level, the DEMs were corrected for radar penetration bias. The volume changes were then converted into glacier mass changes using an ice density value of 850 kg m −3 . The detailed procedure of glacier elevation changes and mass thickness is described in Abdullah et al. (2020).

Vehicular traffic and BC data analysis
To assess the direct effects of the vehicular movement on glacier melting, the glaciers were grouped into two categories based on their proximity to the national highway (NH), by creating a buffer of 6 km 2 on either side of the highway, i.e. proximal glaciers situated within the buffer zone and distant glaciers situated outside the buffer zone.
The emission quantities of GHGs and particulate pollution were calculated from vehicle data, obtained from traffic department, Ladakh government, and also, the light and heavy motor vehicles plying daily on the national highway were analysed. The primary fuel used by vehicles is petrol and diesel. Emission factors (EFs) were calculated using appropriate approaches recommended by various researchers for different gases and particulate matter (Richardson 1999;Baidya and Borken-Kleefeld 2009;Ramachandra et al. 2015;Prakash and Habib 2018) and the estimation is based on the vehicle age, engine technology and fuel type. Emissions from road traffic were estimated based on the number of vehicles and distance travelled by the vehicles in a year (Ramachandra et al. 2015).
During one of the field expeditions conducted from 29 Sep to 03 Oct 2014, black carbon (BC), the main constitutent of vehicular pollution, was measured in the study area on the ablation zone of one of the glaciers, namely, Machoi Glacier, for about 9-10 h daily using the portable AE-42 Athalometer (Magee Sci., Inc., Berkeley, CA, USA) (Bhat et al. 2017). However, to determine the long-term changes in BC in the study area, the area averaged BC surface mass concentration for a longer period from 1980 to 2020, available from the Modern-Era Retrospective Analysis for Research and Applications (MERRA) at a spatial resolution of 0.5° × 0.625°, was also used.

Climate change data analysis
Because the meteorological observations over the area are very scarce due to the challenging rugged topography, remoteness and harsh climatic conditions, a time series of mean annual temperature and precipitation derived from Climatic Research Unit (CRU) time series (TS) v-4.04 grids of dimension 0.5° × 0.5° was used (Harris et al. 2020) to ascertain the long-term climate trends in the study area. The non-parametric Mann-Kendall test was used to determine the significance of the trends in the time-series data from 1901 to 2019 (Mann 1945;Kendall 1975). Because of the presence of three distinct periods showing specific climate trends, the time-series climate data was grouped into three bins: 1901 to 1940, 1941 to 1980 and 1981 to 2019. Climate data trends were compared to the observed glacier recession in the study area to assess the impact of climate change on glacier recession.

Greenhouse gas data analysis
We used a time series of satellite data of different periods from the GOSAT, AIRS and OMI satellites (Table 1) to evaulate the temporal patterns of CO 2 , O 3 , NO 2 , CH 4 and atmospheric water vapour in order to determine the trend of greenhouse gases (GHGs) and their impact on glacier retreat in the study area (Table 1). The gridded AIRS version 6, monthly weighted means, level 3 product, was referred to as AIRX3STM v006 (Chahine et al. 2006). GOSAT satellite data was used for retrieving O 3 data from January 2003 to December 2020. From January 2009 to December 2019, GOSAT satellite data was used to calculate CO 2 column density. The data has a relative accuracy of 1%, which equates to around 4 ppm for CO 2 and accuracy of 2% at a spatial scale of 1° (Kuze et al. 2009). The gridded AIRS version 6, monthly weighted means, level 3 product, referred to as AIRX3STM v006 was used to retrieve CH 4 data from January 2003 to December 2020. CH 4 is retrieved using AIRS channels around 7.6 µm, which are based on the atmospheric temperature and water moisture profiles, as well as surface temperature and emissivity (Kramer et al. 2021). The total atmospheric column NO 2 , for all atmospheric conditions, was retrieved from OMI (Ozone Monitoring Instrument) level 3 daily global-gridded (0.25° × 0.25°) (OMNO2d) data. Visible channel was employed for NO 2 retrieval in sky conditions where cloud fraction was less than 30% (Krotkov et al. 2016). From January 2003 to December 2020, columnar water vapour was retrieved from the AIRS version 6-product referred to as AIRS3STD V006. The AIRS water vapour product provides a standard physical retrieval that contains both AIRS and AMSU radiances, which are used to remove clouds (Susskind et al. 2003).

Glacier recession
The study involved extracting various glacial parameters utilizing multi-temporal and multi-source satellite data; as a result, it is important to quantify the uncertainties associated with mapping glacier outlines using satellite data. The positional and mapping errors, caused by the usage of multi-temporal satellite data, are the sources of error. The degree of uncertainty in glacier mapping is determined by the resolution of the imagery used and the meteorological conditions at the time of image acquisition. It is possible to obtain an accuracy of less than half a pixel under ideal conditions (Bolch et al. 2010;Murtaza et al. 2021). The error caused due to the use of different resolution satellite data was eliminated using identical resolution Landsat satellite images. For the assessment of the uncertainties in the glacier-area mapping, the buffer approach proposed by Bolch et al. (2010) and Granshaw and Fountain (2006) was employed. For the uncertainty assessment, a buffer width of half the pixel resolution of the source image was used. The glacier terminus change uncertainty (U) was estimated (Hall et al. 2003) as follows: where a and b are the spatial resolution of the two satellite images used for glacier mapping and σ represents the error in image registration.

Glacier velocity
The lack of field-observed glacier velocity measurements makes it difficult to assess the accuracy of the remotely sensed SIV measurements. The quality of the remote sensing-derived SIV is mainly affected by the image quality (snow and cloud cover), pre-processing (ortho-rectification, co-registration), the matching method itself and the surface transformation between the two image acquisitions (Scherler et al. 2008;Heid and Kääb 2012;Bhushan et al. 2017). To reduce errors due to the poor image quality, the images used for SIV were acquired towards the end of ablation season, with minimal snow and cloud cover . The uncertainty due to the pre-processing and matching method is usually the easiest to correct (Scherler and Strecker 2012). Orthorectification corrects the horizontal shifts that are caused by inaccuracies in the base DEM that was used in the procedure, rather than the displacement. The U = √ a 2 + b 2 + uncertainties in the SIV from the DEM can be reduced, if the correlated images are acquired from the same path and row (Heid and Kääb 2012). Since the Landsat images are already co-registered, therefore, the error is negligible and falls below the acceptable limits for glaciological studies (Storey and Choate 2004;Heid and Kääb 2012;Bhattacharya et al. 2016). The accuracy of image correlation is in the range of 1/20-1/10 of the pixel size, implying an error of about ± 1.5 m a −1 (Tiwari et al. 2014). The error caused due to surface transformation between the two image acquisitions is difficult to measure because of the deformation, melting, opening or closing of crevasses and spatial distribution of debris cover on glacier and, hence, persists as a residual error . The total uncertainty in the SIV estimates was calculated using Scherler and Strecker (2012) and Garg et al. (2017) approaches, with the assumption that the SIV over the non-glaciated terrain is zero. The mean standard deviation (σ stable area ) and mean velocity over the non-glaciated stable area (M stable area ) were summed up and divided by the interval between the correlated images.

Glacier area changes
Using Landsat images from 2000 and 2020, 77 glaciers in the Drass region were identified and mapped (Fig. 2). The glaciers range in size from 0.27 to 14.01 km 2 with an average of 2.30 km 2 . The glacier area in the region has decreased from 176.77 ± 10.70 km 2 in 2000 to 171.46 ± 10.37 km 2 in 2020, indicating a loss of 5.31 ± 0.33 km 2 (3%) during the last two decades. The pace of glacial recession, on the other hand, varies greatly among glaciers ranging from 0.24% (G57) to 15% (G15), with an average of ~ 5% (Table S1), controlled primarily by the glacier size and topographic parameters such as elevation, slope and aspect of the glaciers. Overall, glacier cover in the study area has reduced by 0.26 km 2 yr −1 during the last 20 years. Despite being an important resource of water, glaciers in the northwestern Himalaya have been studied far less than those in the neighbouring Himalayan regions. Only a few studies regarding glacier recession have been caried out till date in the study area (Koul et al. 2016;Rashid et al. 2021). Koul et al. (2016) investigated 81 glaciers in the region using SOI topographic map of 1976 and LISS-III imageries of 2001 and 2013. They claimed that the glaciers had shrunk from 187.9 km 2 in 1965 to 156.65 km 2 in 2013, owing to morphological, topographical and climatic changes. The Machoi Glacier in the Drass was recently studied by Rashid et al. (2021) and found that it has shrunk by ~ 29% (0.62% yr −1 ) betweeen 1972 and 2019. The results of glacier recession observed in this study were compared to previous studies (Chudley et al. 2017;Romshoo et al. 2020a;, which revealed that the recession rate of Drass glaciers is in accordance with the glacier recession rate of other adjacent glaciated regions with some minor variations (Table 2), owing to differences in the observation period, datasets used and number of glaciers investigated.

Influence of topographic and morphological factors on glacier recession
The observed glacial recession indicates the declining health of glaciers in the study area during the observation period. However, the rate of the observed glacial recession varies primarily because of the variable climatic and non-climatic parameters in the study area. A number of studies have highlighted the significance of topographic and morphological   (Jiskoot et al. 2009;Scherler et al. 2011;Bolch et al. 2012;Murtaza and Romshoo 2017). Therefore, the impact of topographic and morphological parameters such as size, altitude and slope on glacial recession has been evaluated and discussed below.

Impact of altitude on glacier recession
Glaciers in the study area range in elevation from 3600 to 6700 m, with the majority of them situated between 4600 and 5800 m. Although glaciers are melting at all elevations, the rate of retreat varies depending on the elevation. Based on their altitude, the glaciers were categorized into three categories: glaciers between 4600 and 5000 m, glaciers between 5001 and 5400 m and glaciers between 5401 and 5800 m ( Table 4). Analysis of the glacier recession w.r.t. elevation suggested that glaciers at lower elevations have receded by 4.10%, whereas glaciers at mid and higher elevations have receded by 3.23% and 1.46%, respectively, over the observation period. The general trend indicates that glaciers at higher elevation have lower temperature and receive more snow precipitation, resulting in less ablation or even mass gain. The results are consistent with the normal behaviour of recession in relation to elevation. Several studies have found significant altitudinal influence on glacier recession, with glaciers at lower elevations receding faster than those at higher elevations, and viceversa (Fujita and Nuimura 2011;Venkatesh et al. 2012;Kumar et al. 2021).

Impact of slope on recession
The impact of topographic slope on glacier recession has been highlighted in a number of studies (Venkatesh et al. 2012;Nainwal et al. 2016). Slope affects ice velocity, mass flux and snow accumulation rates through its effect on avalanche transportation of snow over the glacier surface (Hoelzle et al. 2003). The slopes of glaciers in the study area range from 11 to 28.5° with a mean slope of 18°. To evaluate the influence of slope on glacier recession, glaciers were categorized into two slope categories (Table 4): glaciers with gentle slope (< 15°) and those with steeper slopes (> 15°). It was observed that glaciers with gentle slope experienced a lower recession of ~ 2.30% than glaciers with steeper slope, which experienced a recession of ~ 3.60%. These findings however contradict those of a number of other studies (Venkatesh et al. 2012;Patel et al. 2018;. The reason for the higher recession rate shown by the glaciers on steeper slopes could be attributed to their lower mean elevation, as glaciers at lower elevations are more susceptible to higher ablation rates (Fujita and Nuimura 2011;Pandey and Venkataraman 2013;Murtaza et al. 2021).

Snout recession
Snouts of all the glaciers in the study area have retreated according to data analysed from 2000 to 2020; however, the retreat rates varies grately among the glaciers. The snout retreat ranged from 30 to 430 m with an average of 155 ± 9.58 m during the study period (5.75 ± 0.47 m a −1 ). Furthermore, it was observed that that the larger glaciers (5-14 km 2 ) have a higher annual snout retreat rate (10.55 ± 0.47 m a −1 ) than the smaller glaciers (6.81 ± 0.47 m a −1 ), which is also reported by Romshoo et al. (2020a) elsewhere in the Himalaya.

Debris cover changes
The extent of debris cover on the glacier in the study area varies from 0.10 to 2.12 km 2 (Fig. 3). The variation in the debris cover of the glaciers situated within the same geographic setting is caused either by the differences in the surrounding lithology or varying geomorphic settings (Scherler et al. 2011). The temporal analysis showed that the debris cover increased from 11.10 km 2 in 2000 to 14.20 km 2 in 2020, a 28% increase during the period (0.15 km 2 yr −1 ). The findings are consistent with the past studies conducted in other regions of the Himalaya, which suggest that glacier debris cover is increasing (Bolch et al. 2008;Bhambri et al. 2011;. The analysis also found that the 54 CGs retreated faster (5.11%) than the 23 DCGs (1.79%) implying that the debris cover considerably modifies the glacier surface thermal regime, which in turn effects the ablation process beneath the debris layer (Table 3). The findings are consistent with other recent studies conducted on the impact of debris cover on glacier melting in the Himalaya, which have suggested that debris cover has a major impact on glacier recession (Gardelle et al. 2013;Pratap et al. 2015;Ali et al. 2017;Debnath et al. 2019).  Table S2 summarizes the respective glaciers minimum, maximum and average SIV. SIV averaged across all the 39 glaciers ranges between 21.35 ± 3.3 m a −1 in 1999/2000 and 16.68 ± 1.9 m a −1 in 2019/2020. Analysis indicates a substantial decrease in SIV from 2000 to 2020 (Table S2), with significant variation among the glaciers. The highest SIV is shown by G13 glacier (46.41 m a −1 ) and slowest by G36 glacier (8.0 m a −1 ). The highest decrease in SIV is noticed for the G34 glacier (31.57%) and the lowest for the G48 glacier (1.26%) during the observation period (Table S2). Between 1999/2000 and 2019/2020, SIV of all glaciers decreased, on an average, by 23.27% (1.16% a −1 ) with an annual average decrease of 0.23 m a −1 . The displacement observed over the non-glaciated terrain is considered an uncertainty in SIV estimation which ranges from 1.9 to 3.3 m (Table S3) in this study and is in accordance with the published studies (Scherler et al. 2011;Garg et al. 2017;Bhushan et al. 2017;. The decreasing observed SIV trend is in line with the various other studies reported in other parts of the Himalaya. A recent study by , based on 18 glaciers in the central Himalaya, reported an average SIV of 22.63 ± 5.8 m a −1 in 1993-94, which decreased to 17.32 ± 3.1 m a −1 in 2000-01 and further slowed to 11.50 ± 1.7 m a −1 in 2015-16. Bhattacharya et al. (2016) reported SIV of ~ 46 ± 7.5 m a −1 during 1993-94, ~ 50 ± 7.2 m a −1 in 1998-99 and 43 ± 5.1 m a −1 during 2013-14 for the Chorabari Glacier in the central Himalaya. Similarly, Bhushan et al. (2018) reported an average SIV of 42.42 ± 5.6 m a −1 for glaciers in the Zanskar Himalaya during 1999-2000, which decreased to 35.44 ± 5.6 m a −1 by 2013-14. The main reasons for the reduction in SIV are retreat, thinning and mass loss of glaciers and the consequent decrease in the glacier gradient and basal sliding (Bhushan et al. 2017). All these studies suggested that variations in SIV reflect changes in the evolution of ice flux in the middle and upper parts of glacier, consequently causing dynamic feedbacks in mass budgeting of a glacier.

Glacier mass changes
Accurate assessment of glacier mass loss is critical for understanding glacier sensitivity to climate change. There is a spatial variation in the observed glacier thickness changes in the study area. Thinning of glaciers suggests mass loss. Table S4 shows specific mass balance (SMB) rate of the selected glaciers. It was observed that the SMB rate of different glaciers ranges from − 1.55 ± 0.37 to − 0.36 ± 0.33 m w.e.a −1 . The highest negative SMB rate is shown by G17 glacier (− 1.55 ± 0.37 m w.e.a −1 ), whereas G18 glacier experienced the lowest mass loss rate (− 0.36 ± 0.33 m w.e.a −1 ) among the 77 glaciers investigated in the study area. The mean SMB rate of all the glaciers is − 1.08 ± 0.31 m w.e.a −1 . The mass loss is often regarded as the most credible indicator of glaciers responding to climate change.

Vehicular pollution and black carbon concentration
Several studies have found that an increase in BC concentration in the atmosphere contributes significantly to atmospheric warming as BC is highly heat absorbing atmospheric constituent, thus enhancing glacier ablation (Ramanathan et al. 2007;Xu et al. 2009). From the analysis of black carbon observations in the study area, it was observed that BC concentration ranged from 287 to 3726 ng/m 3 with an average of 1518 ng/m 3 (Fig. 5), which is markedly higher compared to the BC concentration reported from other high altitude locations in the Hindu Kush Himalaya (Table S5) (Zeb et al. 2020;Nair et al. 2013). Thakur et al. (2021) observed an average BC concentration of 127 ng g −1 and 131 ng g −1 from snow pits of Khardung and Phuche Glaciers, Ladakh. Furthermore, trend analysis of the time series of reanalysis BC data from 1980 to 2020 (Fig. 6) showed that BC concentration has increased significantly from ~ 338 ng/m 2 in 1984 to ~ 634 ng/m 2 in 2020. It is thus inferred that the increasing BC concentration, due to the proximity to the national highway (NH), has significantly affected the glacier health  as is evident from the observed glacier area shrinkage, snout retreat and mass loss in the study area. The proximal glaciers, totalling 17, showed higher glacier shrinkage (4.11%) and snout retreat (209 m) than the distant glaciers, numbering 60, situated away from the NH which showed glacier shrinkage of 2.82% and snout retreat of 148 m during the observation period (Table 5).

Vehicular pollution
Transport sector forms the bulk of CO 2 and other greenhouse gas emissions from fossil fuel combustion. About 48% of the total freight plying on the national highway passing through the study area is composed of the heavy motor vehicles, which includes buses, trucks, minibuses, heavy-duty trucks, sumo, jeeps and trailers, which run on diesel fuel. Light motor vehicles, such as two-wheelers, three-wheelers, cars and others, which operate on gasoline, account for the remaining 52% of the traffic. Figure 7 and Table 6 show the annual CO 2 , NOx, BC, CH 4 and PM emissions for various vehicle categories from 2011 to 2020. GHGs and particulate matter emissions increased from 1728, 116, 22.34, 5.66, 0.64 and 1.28 tonnes in 2011 to 4413, 56.36, 13.20, 1.43 and 3.05 tonnes in 2020 indicating an annual growth rate of 11.8%, 11.07%, 11.4%, 11.4% and 11.6% for CO 2 , NOx, BC, CH 4 and PM, respectively, due to the increase in the traffic load. Overall, emissions of CO 2 , NOx, BC, CH 4 and PM have increased by around 39% during the study period. It was observed that the heavy motor vehicles are the primary sources of CO 2 emissions accounting for 90% of total CO 2 emissions. Light motor vehicles that run on gasoline contributed about 56% of total CH 4 emissions primarily because of their large number travelling on the NH. For NO X emission, heavy motor vehicles account for 88% of the total NO X emissions. Buses, goods-carrying vehicles and other heavy vehicles are the primary sources of PM pollution   GHG emissions are increasing, which contributes to global warming and enhances the melting of glaciers in the region. Heating by black carbon (BC), on the other hand, warms the atmosphere higher elevations from 2 to 6 km (Ramanathan and Feng 2009), where the majority of glaciers in the study area are located, amplifying the effect of GHGs on the retreat of snow packs and glaciers in the region.
Greenhouse gas data analysis CO 2 , which accounts for 9-26% of the greenhouse gas effect on earth, is followed by water vapour, which accounts for 36-70%; O 3 , which accounts for 3-7% and CH 4 , which accounts for 4-9% (Griggs and Noguer 2002). The average monthly values of the GHGs were used in this study to determine GHG trends over the study area (Tables 7, S6). Figure 8 depicts the monthly distribution of GHG data from AIRS, GOSAT and OMI. The highest CH 4 was found in the month of August, September and October (Table S6). The time-series data showed an increasing trend with 18.3 ppb year −1 , with CH 4 values ranging between 1812 and 944 ppb and the highest concentration occurring during the spring and summer months (Fig. 8I). The increasing trend of CH 4 is primarily due to the global increase in emissions from natural and anthropogenic sources; however, the semi-arid climate of the study area with almost no vegetation exacerbates the problem (Kavitha et al. 2018).
The highest levels of O 3 concentration are seen in the month of April, May and June (Tables 7, S6). The timeseries data shows an increasing trend of 0.25 ppb year −1 with values ranging from 43 to 69 ppb (Fig. 8II). CO 2 sequestration is essentially non-existent in the study area, due to its cold desert nature, which is characterized by semi-arid climate, and a lack of vegetation cover (Juyal 2014). Figure 8III shows the monthly daytime/night-time retrievals of water vapour distribution and the total columnar water vapour values ranging from 7 to 47 kg/m 2 . The time-series water vapour data shows no overall trend. The monthly mean values, however, indicate clear seasonal variation, with the highest (average of 39 kg/m 2 ) and lowest (average of 10.5 kg/m 2 ) during the summer and winter seasons, respectively (Table S6).
Mean monthly CO 2 data from January 2009 to December 2019 calculated from GOSAT measurements are presented in Fig. 8IV. CO 2 peaks in March, April, May and June almost every year, with troughs in August and September (Table S6). CO 2 levels range from 379 to 413 ppm (Fig. 8IV), with winters showing the highest and autumn the lowest (Tables 7, S6). The overall time-series data shows an upward trend (1.5 ppm year −1 ). CO 2 emissions, transport and military operations in the region all have the potential to increase local CO 2 levels (Al-Bayati and Al-Salihi 2019). Figure 8V shows the monthly variation of NO 2 in the troposphere column over the study area, with an upward trend detected from 2005 to 2015. The annual increase of 0.72% is observed in NO 2 concentration, which is 1.2 × 10 13 molecules cm 2 year −1 . The average value of NO 2 over the study area was 1.3×10 15 molecules cm 2 . The highest concentration was observed in the summer months of May and June and the lowest in the winter month of February (Tables 7, S6). Several factors are responsible for the observed seasonal changes and the increasing trend of the GHG (Lamsal et al. 2014) in the study area but the close proximity to the national highway and the heavy traffic plying in the area is one of the most important factors contributing to the higher values of the NO 2 . Figure 9 depicts a trend analysis of mean annual temperature data from 1901 to 2019. Mann-Kendall analysis of the mean annual temperature demonstrated a statistically significant increasing trend at < 0.01 with a Z statistic value of > 2.57, indicating high significance of the observed trend. It was also observed that the mean annual temperature during 1901-1940 period (1.13 °C) was lower than that observed during 1941-1980 period (1.43 °C). The mean annual temperature further increased sharply to 1.69 °C during 1981-2019, indicating a significant increase in the mean temperature in the region during the last few decades. This significant rise in the mean temperature, due to the increasing GHGs, PM and BC is one of the likely reasons for glacier recession observed in the region. Several studies have found a significant negative correlation between glacier recession and temperature (Zhang et al 2014;Mandal et al. 2016).

Climate data analysis
Similarly, CRU precipitation data analysis revealed that an average precipitation of 635 mm was observed between 1901 and 2019 (Fig. 10). The trend analysis shows that there is inter-annual variability in overall precipitation, with a decreasing but non-significant trend from 1901 to 1940, increasing but non-significant trend from 1941  to 1980 and almost flat average precipitation from 1981 to 2019, indicating that precipitation over the region has remained almost constant over the last few decades, which corresponds to the observed glacier recession. This shows that the region has had a warmer climate for several decades, which could be the cause for the observed glacier recession in the study area under the no-change precipitation scenario.

Conclusions
A bidecadal basin-wide assessment of a variety of glacier parameters such as glacier shrinkage, snout retreat, glacier thickness change, mass loss and glacier velocity revealed a strong correlation with increasing concentrations of the GHGs, BC and other pollutants caused by anthropogenic activities. The resulting climate change primarily drives the glacier melting in the basin, with an average glacier mass loss of − 1.08 ± 0.31 m w.e.a −1 . The increasing BC concentration, in the vicinity of glaciers as a result of their proximity to the NH, has had a significant impact on glacier health, as evidenced by the glacier recession observed in the study area. Given the varying responses of the Himalayan glaciers to climate change, and the fact that the observation of climate data, BC, GHGs and other pollutants is few and far between in the inaccessible regions, the use of frequent, synoptic and high-resolution remotely sensed observations of the glacier parameters, GHGs, pollutants and other driving factors provided a reliable alternative. Local emissions of black carbon and other short-lived pollutants from the combustion of fossil fuels and woody biomass exacerbate effects of climate change on glacier melting and must be considered when assessing the impact of climate change on glacier melting. Given the limited observations of BC and other pollutants in the basin, it is important to strengthen the observation network in the data-scarce Himalaya in order to fill knowledge gaps about the impacts of local anthropogenic drivers and global climate change on cryosphere melting.
It is feared that if the observed trends of the climate change continue in the future due to increased greenhouse gas emissions and increase in BC concentration and other anthropogenic pollutants, glaciers in the Himalaya may disappear entirely, having a significant impact on regional water supplies, hydrological processes, ecosystem services and transboundary water sharing.