Spatio-temporal variability of glacier surface area, ELA, and AAR in the Mago river basin of Arunachal Pradesh

Snow and glaciers play an important role in the global energy and hydrological cycles as they have high albedo and act as water reservoirs (Tang et al. 2019). Climate change, which is a complicated phenomenon, is mostly understood through substantial changes in the cryosphere (Kaushik et al. 2019). Outside the polar region, the largest place of snow and glaciers is found in the Himalayan region (Dixit et al. 2019). Glaciers melt and snowmelt received from the Himalayan Mountains forms the perennial river systems such as Brahmaputra, Ganga, and Indus, and the melt runoff contributes around 30–50% of the annual flow (Aggarwal et al. 1983; Jain et al. 2010; SAC 2011).

In the Himalayan Mountain glaciers, there is a considerable area and volume losses in the last 10 years due to a rise in magnitude of climate change and anthropogenic actions (Chand & Sharma 2015; Zemp et al. 2015; Singh et al. 2018; Maurer et al. 2019; Lee et al. 2021). The glacier lake outburst and retreating are very frequent phenomena that can be seen in the Indian Himalayan Region (IHR) (Bajracharya et al. 2008; Padma 2020; Dimri et al. 2021). Many studies conducted in the Indian Himalaya report that the glacier shrinkage rate is at an alarming rate (Kulkarni 2007; Mir et al. 2017; Alam & Bhardwaj 2020). The study of glacier changes is significant for monitoring and maintaining water management for ecosystem processes, irrigation practices, hydro-power generation, socio-economic development, and climate change studies due to their sensitivity towards minute changes in climate. The surface water hydrology in a basin can be considerably modified by the glaciers, even when the glacier surface coverage is as small as 5% of the basin area (Fountain & Tangborn 1985; Mcgrath et al. 2017).

The glacier area extent is key for nearly all the glaciological studies and hydrological modeling (Shukla et al. 2009). The accurate information on the glacier area extent and distribution of glaciers is necessary for water resource management, mitigation of glacial hazards, and estimation of the contributions of glaciers to sea level change in the past years and years to come (Rastner et al. 2014; Marzeion et al. 2017). The mapping of glaciers shows great significance for being a representation of global climate change (Oerlemans 2005). Changes in glacier area are primarily related to the mass balance of the glacier (Keeler et al. 2021). The annual mass excess (when net accumulation surpasses ablation) leads to glacier advancement, while a deficit causes glacier reduction (Keeler et al. 2021). Glacier mass balance (GMB) is a basic assessment to analyze the contribution of glaciers melt to the regional river systems (Kaser et al. 2010; Huss et al. 2017) and global sea level. Glaciers in Himalayan region are generally difficult to observe and assess through field measurements due to highly rugged and extremely unreachable mountainous topography (Mir et al. 2017) and thus creating a limitation in timely monitoring and collection of data through field observations (Mir et al. 2014). The use of remote sensing and GIS techniques provide an alternative method for observing and monitoring the glacier extent at regular time intervals (Mir et al. 2014). The most widely used remote sensing methods for glacier area extraction are the single band ratios (NIR/SWIR and Red/SWIR), Normalised Difference Snow Index (NDSI), and Principal Component Analysis (PCA) (Zhang et al. 2019). These techniques make use of the high reflectance of snow and ice in the visible-near infrared (VNIR) wavelength and low reflectance in the shortwave infrared (SWIR) to identify glaciers from darker areas such as rock, soil, or vegetation (Mir et al. 2014; Paul et al. 2015). These remote sensing methods tend to misclassify water features and shadowed areas as glaciers (Paul et al. 2015, 2017). In order to reduce the misclassification and to increase the accuracy of classification the Automatic Glacier Extraction Index (AGEI) is developed to reduce the misclassification and to escalate the accuracy of classification (Zhang et al. 2019).

The assessment of glacier's response to climate change is better measuring Equilibrium Line Altitude (ELA) of a glacier than glacier area or length changes (Keeler et al. 2021). The estimation of the ELA of a glacier is often required when there is no direct glacier mass balance observation (Osmaston 2005). The ELA of a glacier is defined as the line that demarcates the accumulation and ablation areas and represents the altitude at which the yearly amount of mass added exactly equals the yearly mass lost (Cuffey & Paterson 2010). The ELA is hardly detected as unbroken line at the same altitude traversing the entire extent of the glacier, due to localized topographic and climatic disparities in ablation and accumulation (Hughes 2009; Bakke & Nesje 2011). The ELA is mainly affected by precipitation and temperature (Porter 1975; Rupper & Roe 2008; Sagredo et al. 2014) due to this it is a good indicator of the glacier's response to climate change at the regional scale (González-Reyes et al. 2019). Therefore, variations in the altitude of the equilibrium line on a particular glacier can be used as an indicator of climatic fluctuation (Kuhn 1981). There are several methods that have been developed to determine ELA when there is no available direct measurements of mass balance. Some of the most extensively used methods are the accumulation area ratio (González Trueba & Serrano 2004), the toe-to-headwall altitude ratio (Meierding 1982), the altitude area (AA), the altitude area balance ratio (AABR) (Osmaston 2005), the mean glacier elevation (MGE) (Kurowski 1891). In order to reduce the laborious work of calculating the ELA, the automatic GIS tool based on AAR, AA, AABR, and MGE was developed (Pellitero et al. 2015). The different approach is known as topographic map-based method in which the ELA is estimated from the inflections of contour on topographic maps (Leonard & Fountain 2017). All of the methods have their own benefits and limitations, and selection of the method may depend on the availability of data as well as on the scale of investigations, however, remote sensing methods might be a better option at a smaller scale. (Racoviteanu et al. 2019).

Using remote sensing methods, the ELA of a glacier has been determined from snow line altitude (SLA) at the end of the melting season on the glacier surface through satellite imagery (Lei et al. 2012; Racoviteanu et al. 2019; Rastner et al. 2019; Tang et al. 2019). In the hydrological field, the snow line is recognized as the line separating snow area from snow free area (Kaur et al. 2010). So, mapping of the minimum seasonal snow cover extent on the glacier will give a proxy for ELA (Rastner et al. 2019). There is a significant correlation between snow lines obtained from satellite data and ELAs obtained from field measurement (Paul et al. 2009). The SLA on a glacier signifies that the altitude higher than SLA is the accumulation zone and the altitude lower than SLA is the ablation zone of a glacier. Separation of snow on glaciers helps to determine the accumulation area ratio (AAR) of a glacier which is defined as the ratio of the accumulation area to the total area of the glacier (Meier & Post 1962). The ELA and AAR are directly proportional to each other when the ELA shifted up the accumulation area reduces and vice-versa. ELA-AAR have been used to evaluate glacier mass balances at regional scale using remote sensing methods (Kulkarni et al. 2004; Rabatel et al. 2008). The main problem in mapping snow cover on glaciers is the similarity of spectral signatures of both glacier and snow and creating difficulties in differentiating them. The most extensively applied approach for mapping snow cover on glaciers is applying threshold value on single-band reflectance image, preferably in the NIR in order to avoid saturated values over snow with 8-bit Landsat data (Hall et al. 1987; Williams et al. 1991; Rabatel et al. 2012; Shea et al. 2013). The selection of the threshold value is the main key in separation of snow and ice.

The study of changes in glacier surface area, ELA and AAR have been carried out in the Central Himalaya, Western Himalaya, some parts of Eastern Himalaya, and in the Karakoram ranges (Wagnon et al. 2007; Dobhal et al. 2008; Ramanathan 2011; Azam et al. 2012; Vincent et al. 2013; Pratap et al. 2014; Mir et al. 2017; Rai et al. 2017; Racoviteanu et al. 2019; Alam & Bhardwaj 2020).

In the present study, a glacierized basin in the Eastern Indian Himalaya was analysed using tools developed by Pellitero et al. (2015) and algorithm developed by Rastner et al. (2019) and Zhang et al. (2019).

For mapping glaciers surface area and snow cover over the glacier area, multi-spectral satellite imagery of Landsat 8 OLI/TIRS and Landsat 5 TM were downloaded from United States Geological Survey (USGS) website (https://glovis.usgs.gov/) for the period of 1988 to 2019. The primary criteria for downloading the Landsat imagery data was based on less snow coverage, i.e., end of snow melting season and low cloud cover. The downloaded data were the best possible combinations that can be obtained in the study area. The digital elevation model (DEM) was required for determining the slope, aspect, contour lines and ELA of the glacier which was obtained from NASA website (http://search.earthdata.nasa.gov/). The downloaded DEM is of Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Digital Elevation Model (DEM) of 30 m spatial resolution. The details of downloaded Landsat data are given in Table 1.

The ELA of the glaciers greater 1 km² was determined using two approaches, namely GIS tool developed by Pellitero et al. (2015), and using NIR band of Landsat data and Otsu thresholding (Rastner et al. 2019). Both the methods were compared with high resolution Google Earth images for the selection of the satisfactory method in the study area.

The GIS tool consists of four methods: Accumulation Area Ratio (AAR) (González Trueba & Serrano 2004), Area-Altitude Balance Ratio (AABR) (Osmaston 2005), Area-Altitude (AA) and Mean Glacier Elevation (MGE) (Kurowski 1891). This toolbox simplified the process of ELA estimation and successfully worked on both single glacier and multiple glaciers. This toolbox was coded in Python and run in ArcGIS requiring only the reconstructed surface of the paleo-glacier (a DEM) as input. After providing the DEM, we need to give a ratio value for AAR and AABR. The balance ratio was obtained from Rea (2009) for global (1.75 ± 0.71) and Central Asia (1.75 ± 0.56). The accumulation ratio of 0.5 means the ablation and accumulation zone are exactly equal, so a ratio of 0.5–0.65 was used (Quesada-Román et al. 2020). After checking different ratios within the given range for both AAR and AABR, a ratio of 0.6 and 1.75 was found to be the most suitable ratio in the study area.

Snow cover on glaciers were extracted by applying a threshold value interactively generated with the Otsu algorithm (Otsu 1979). The algorithm is a widely used automatic thresholding method with the aim to maximize inter-class variance and minimise intra-class variance. Otsu automatically detected a threshold value, t, that separates snow and glacier. The Otsu thresholding was performed in ArcGIS Pro where the in-built binary threshold is based on the Otsu algorithm (Otsu 1979). From this binary image snow cover on the glacier was obtained. For determination of SLA, firstly DEM of the glacier was converted into 20 m elevation bins. The snow cover map on the glacier was then digitally intersected with the elevation bins. The snow and ice pixels in the bins were counted and the lowest bin that showed snow cover percentage more than 50% was selected as SLA (Rastner et al. 2019), this meant all the upper bins (above the selected bin) contained snow cover percentage more than 50% (Figure 2). The average value of the selected bin was taken as ELA and using the determined ELA, AAR was calculated.

In other parts of the Himalayan region, reduction in glacier surface area was observed too. In the Tirungkhad basin, Western Himalaya, a total shrinkage of 29.1 km² (26.1% reduction) was observed during 1966 to 2011 (Mir et al. 2014). In the Garhwal Himalaya, glaciers in Saraswati/Alaknanda basin and upper Bhagirathi basin lost glacier surface areas of 18.4 ± 9.0 km² (5.7 ± 2.7%) and 9.0 ± 7.7 km² (3.3 ± 2.8%), respectively, from 1968 to 2006 (Bhambri et al. 2011). In the Himachal Himalaya, glaciers in the Chenab basin retreated at about 0.53% per year, in Parbati at about 0.56% per year and in Baspa at about 0.48% per year from 1962 to 2001 and this investigation had shown the overall deglaciation of 21% during the study periods (Kulkarni 2007). The five selected glaciers (Bara Shigri, Chota Shigri, Hamtah, G4, and Parvati glacier) in Himachal Pradesh, showed reduction in glacier surface area from 154.58 to 123.39 km² (20.17% deglaciation) during the period from 1976 to 2013 (Rai et al. 2017). In Baspa basin of Western Himalaya, a shrinkage in glacier area was recorded at a rate of 1.18 ± 0.3 km² per year from 1976 to 2011 (Mir et al. 2017). The temporal variation in glacier's area analysis in Sikkim, Eastern Himalaya showed that shrinkage of glacier area of about 22 ± 6% from 1992 to 2015 (Alam & Bhardwaj 2020). Upon comparing the present findings and other results, it was observed that glaciers in Mago river basin showed similar deglaciation percentage and retreating rate with other glaciers in the Himalayan mountainous region.

The ELA of the 15 glaciers were determined using both the automatic GIS tool, and NIR band and Otsu thresholding methods. For the selection of the best method, high resolution Google Earth images were used as reference due to lack of observed field data because of the inaccessibility of the area. The Google Earth images were available for the year 2015 and 2016 in the glaciers area; therefore, the ELA was calculated first for these two years. From the GIS ELA tool, three methods, namely, MGE, AAR, and AABR were used. Among these methods, the AABR showed a feasible ELA lines for all the glaciers compared to the other two methods. Using the same automatic GIS tool, the ELA was calculated in the Shigar river basin of Karakoram Range, Pakistan, from AAR and AABR, it was found that AABR method was appropriate to estimate ELA of clean-ice covered and snow-fed glaciers (Baig & Muneeb 2021). The empirical evidence from modern glaciers showed that MGE overestimated the ELA and was unable to take the variations in morphology into account, which intensely affects the area-elevation distribution of a glacier (Bakke & Nesje 2011).

For the analysis of the ELA variations and rate of variations in ELA position, the analysis period was divided into three periods, i.e., 1988–2010, 2010–2019, and 1988–2019. From the average ELA changes, the ELA was increased by 86.7, 50.7, and 137 m in the analysis period between 1988–2010, 2010–2019, and 1988–2019. The average rate of ELA variations showed that ELA was increased at the rate of 3.94, 7.85, and 4.43 m/year in the analysis period of 1988–2010, 2010–2019, and 1988–2019, respectively. The analysis period of 2010–2019 showed highest ELA changes rate. The glaciers G36 and G41 showed same elevation of ELA at the start and end of the year in all the analysis periods, though there were some fluctuations between years. From the glacier's orientation, by taking the average of the elevation changes of the glaciers having same orientation, it was observed that in the analysis period between 1988 and 2010, SE facing glaciers had highest ELA increasing rate (7.28 m/year); in 2010–2019, South West facing glaciers had highest increasing rate (13.33 m/year) and in 1988–2019, East facing glaciers had highest ELA increasing rate (6.78 m/year).

The study of annual variations in glacier ELA in the western Alps from 1984 to 2010 for 43 glaciers from the end-of-summer snow line altitude showed the average ELA increased by 170 m and the time series showed an increasing trend of 6.4 m/year (Rabatel et al. 2013). Considering each glacier, the increasing trend of the ELA varied between less than 1 m/year and more than 13 m/year and glaciers facing east displayed a more pronounced increasing trend (Rabatel et al. 2013). Over the period studied (2000–2016), the annual ELAs fluctuated from 4,917 to 5,336 m (419 m change) in Hunza basin, Karakoram and 5,395 to 5,565 m (170 m change) in Trishuli basin, eastern Himalaya (Racoviteanu et al. 2019). So, the annual ELA variation rate was found to be 26.18 and 10.63 m/year, in Hunza and Trishuli basins, respectively.

Using the calculated ELA, the actual AAR of each glacier for all the study years was determined (Table 6). From the average value of AAR throughout the study period, it was found that in the study area, the AAR varied from 0.44 to 0.97 with an average of 0.75. The large range of AAR was observed within the study area and it helped us in understanding that the accumulation zone and ablation zone depend not only on the climatic conditions but also on the topography and the morphometry of the glaciers. Even if the glaciers were within the same catchment, it does not mean that the ratio will be same for all. So, there was disadvantage of applying a constant ratio value for all the glaciers as the accumulation and the ablation zone responds differently. The ratio is different from glacier to glacier. The analysis based on glaciers orientation was carried out, and it was observed that, by taking the average of the AAR of the glaciers having same aspect, the glaciers of South, South East, and North had AAR of around 0.8, South West and West facing glaciers showed around 0.6, and 0.74 was observed in East facing glaciers.

The study of AAR changes in Chhota Shigri glacier in Chandra basin, Western Himalaya, from 2005 to 2011 showed the lowest AAR of 0.29 and highest of 0.74 (Wagnon et al. 2007; Ramanathan 2011; Azam et al. 2012; Vincent et al. 2013) and the ELA changes during this period was found out to be 118 m (from 4,855 to 4,973 m). The AAR for Dokriani glacier in Bhagirathi basin, Central Himalaya, varied from 0.66 to 0.70 during 1993–2000 (Dobhal et al. 2008; Pratap et al. 2014) and the ELA changes during this period was 65 m (from 5,030 m to 5,095).

All the glaciers experienced shrinkage in the surface area during the analysis period from1988 to 2019. The glaciers of G1, G9, G15, G35, G41, and G42 showed significant decreasing trend of the surface area as their R² was greater than 0.6. The glacier shrinkage rate was found to be higher in recent years' analysis (2010–2019) than the old periods (1988–2010). Considering single glacier loss from 1988 to 2019, glacier G40 showed highest deglaciation percentage and G3 showed lowest deglaciation percentage. It was observed that the glaciers facing South West and West lost more surface area compared to other glaciers. The ELA determined using the NIR band of Landsat dataset and Otsu thresholding gave a satisfactory result for all the tested glaciers in comparison with other methods based on Google Earth images. The ELA of all the glaciers showed an increasing trend during the analysis period between 1988 and 2019. The ELA was shifted to the higher elevation zone. The significant increasing trend was observed in glacier G2, G15, and G40 as the linear trend R² was greater than 0.6. The analysis period of 2010–2019 showed highest ELA changes rate. The highest ELA changes rate was observed in glaciers G8 (22.22 m/year) during the 2010–2019 analysis period. The glaciers G36 and G41 showed same elevation of ELA at the starting and ending of the year in all the analysis periods, though there were some fluctuations in between years. From the glacier's orientation, it was observed that in the analysis period between 1988 and 2010, South East facing glaciers had highest ELA changes rate, and in 2010 to 2019, South West facing glaciers had highest variation rate and 1988–2019, East facing glaciers had highest ELA changes rate. The large range of AAR was observed in the study area, from the average value of AAR throughout the study period, it was found that in the study area, the AAR varied from 0.44 to 0.97 with an average of 0.75. It was also observed that the glaciers of South, South East, and North had AAR of around 0.8; South West and West facing glaciers showed around 0.6, and around 0.7 was observed in East facing glaciers.

The authors gratefully acknowledge the help, encouragement, and financial support provided by the Climate Change Programme (CCP), Strategic Programmes, Large Initiatives and Coordinated Action Enabler (SPLICE), Department of Science and Technology, Govt. of India under National Mission on Sustaining Himalayan Ecosystem (NMSHE) through Grant No. DST/CCP/NHC/154/2018 and DST/CCP/MRDP/184/2019.

ArcMap and ArcGIS Pro are used under license which needs to be procured from ESRI.

V. Nunchhani conducted data acquisition, developed the methodology, prepared the original draft. V. Mekro prepared the data and developed the methodology. Arnab Bandyopadhyay supervised the process, edited the manuscript, communicated and helped receiving the grant. Aditi Bhadra conceptualized the whole article, supervised the work, visualized the article, and helped receiving the grant.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.