Linking the Recent Glacier Retreat and Depleting Streamﬂow Patterns with Land System Changes in Kashmir Himalaya, India

: This study reports the changes in glacier extent and streamﬂow similar to many Himalayan studies, but takes the unusual step of also linking these to downstream land use changes in Kashmir Valley. This study assessed changes in the area, snout, and equilibrium line altitude (ELA) of four parts of the Kolahoi Glacier using earth observation data from 1962 to 2018. Changes in the discharge of the two streams ﬂowing out from Kolahoi Glacier into the Jhelum basin were also assessed between 1972 and 2018. Additionally, satellite data was used to track the downstream land system changes concerning agriculture, orchards, and built-up areas between 1980 and 2018. This analysis suggested a cumulative deglaciation of 23.6% at a rate of 0.42% per year from 1962 to 2018. The snout of two larger glaciers, G1 and G2, retreated at a rate of 18.3 m a − 1 and 16.4 m a − 1 , respectively, from 1962 to 2018, although the rate of recession accelerated after 2000. Our analysis also suggested the upward shift of ELA by ≈ 120 m. The streamﬂows measured at ﬁve sites showed statistically signiﬁcant depleting trends that have been a factor in forcing extensive land system changes downstream. Although the area under agriculture in Lidder watershed shrunk by 39%, there was a massive expansion of 176% and 476% in orchards and built-up areas, respectively, from 1980 to 2018. The conversion of irrigation-intensive agriculture lands (rice paddy) to less water-intensive orchards is attributed to economic considerations and depleting streamﬂow.


Introduction
Himalayan glaciers are in a continuous state of retreat since the 19th century in response to climatic change and anthropogenic activities [1][2][3] except in the Karakoram region where glaciers have been reported to be in a stable phase [4][5][6]. The ever increasing temperatures have resulted in the faster melting of cryosphere reserves in the region [7,8]. Although most of the studies consider climate to be the main controlling factor in glacier recession [9][10][11], many studies identify the influence of topography [12][13][14] and debris cover on glacier retreat [15,16]. The warming trends over the region have not only accelerated the glacier melt [17][18][19] but have also changed the form of precipitation [20,21] that has resulted in persistently negative glacier mass balances in the region [22,23].
The glacier recession in the Himalayan arc has not only impacted the streamflows but also resulted in the formation of proglacial lakes [24,25], which could be potential sites for occurrence of glacial lake outburst floods. Compared to other regions of high-mountain Asia, the glaciers in the Kashmir region are retreating at an accelerated pace [26]. The shrinking of the Himalayan cryosphere has been linked Both Lidder and Sind streams are major tributaries of Jhelum River and originate from the Kolahoi Glacier. The streams are being fed many other smaller glaciers and seasonal snow-packs [46]. The area falls in a temperate climate zone with four distinct seasons: spring, summer, autumn, and winter. The mean annual precipitation recorded at Pahalgam located at an elevation of 2150 m above mean sea level (asl) is 1240 mm. A dominant amount of this precipitation falls as snow during winter and spring from westerlies, whereas summer and autumn are relatively drier. The minimum temperature during winter can drop below −18 • C, whereas the maximum temperature can rise up to 30 • C during summer [47]. The topography of the Lidder watershed is hilly, covered by snow and glaciers or rock outcrop in higher reaches (>3600 m asl). Relatively mid-altitude reaches of the watershed are covered by alpine meadows (>3000 m asl) and dense evergreen coniferous forest (1600-3600 m asl), whereas the fertile alluvial plains (1580-2200 m asl) towards the south of watershed are dominated by agriculture, orchards, and settlements [48].

Datasets and Methods
A repository of data from multiple sources comprising satellite images, digital elevation models (DEMs), topographic map, ground truth, and ancillary records have been used to complete this study. The details of the datasets are provided as Supplementary Information (Table 2).

Datasets and Methods
A repository of data from multiple sources comprising satellite images, digital elevation models (DEMs), topographic map, ground truth, and ancillary records have been used to complete this study. The details of the datasets are provided as Supplementary Information (Table 2).  This study assessed changes in the Kolahoi Glacier between 1962 and 2018. However, keeping in mind the recession of glacier and GLIMS (Global Land Ice Measurements from Space) definition that suggests that the ice divides define the hydrological boundaries of glaciers, changes in the spatial extents of four parts of the glacier (G1, G11, G111, G2) were assessed. The apparently connected units were separated in the accumulation zone using 3D information from ASTER DEM and Google Earth [49,50]. Although G1, G11, and G111 drain into Lidder, G2 drains into the Sind watershed of the Jhelum basin [51]. The changes in the overall area of Kolahoi Glacier and also its four parts were assessed using multi-date satellite data (1992-2018) and Survey of India Topographic map of 1962. Although the entire Landsat data archive is available as a georeferenced product, the topographic map was coregistered with satellite data using map-to-image georeferencing [52] in ERDAS Imagine 9.1. The ground control points, which consisted of permanent features such as bridges, mountain ridges, and prominent street crossings, were identified from both the topographic map and satellite data to ensure better geometric correction, achieving a root mean square (RMS) error of less than 1.00. This allowed co-registration of all the images and topographic map so that they could be draped over each other to precisely delineate glacier extents across timescales. Satellite data from the end of the ablation season, August to October, was used to delineate the glacier extent. Because clouds on the satellite data mask the surface features, we ensured cloud-free images for glacier delineation. Glacier boundaries were delineated using on-screen digitization at a 1:30,000 scale for 1962, 1992, 2000, 2013, and 2019 using standard false color composite band combination of 5-NIR:4-Red:3-Green for Landsat 8 OLI data and 4-NIR:3-Red:2-Green for Landsat TM Landsat ETM+ data. Additionally, 3D perspective and historical imagery from Google Earth helped in delineating areas that were influenced by terrain shadows. The snout retreat of Kolahoi Glacier was also mapped for 57 years from 1962 to 2019. Image enhancement techniques in Arc Map 10.1, histogram equalization and contrast stretch were also applied to the satellite data for better delineation of glacier extents [53]. Additionally, glacier boundaries delineated from 2013 and 2018 satellite data were also corrected using high-resolution Google Earth Imagery. The 2013 and 2018 outlines were converted into KML (Keyhole Markup Language) files and imported into Google Earth. Modifications to the KML files were made, wherever required, using high-resolution images and 3D terrain information from Google Earth.
The accuracy of glacier delineation is influenced by the spatial resolution (λ) of satellite data and error in georeferencing (ε) the satellite data [54] as where E A is the uncertainty in the glacier area and n is the number of pixels along the boundary of the glacier polygon. The uncertainties related to the changes in the glacier area (E AC ) can be expressed [55] as where E A1 and E A2 represent uncertainties in the glacier area between two time periods. The uncertainties related to the changes in glacier snout (E SC ) can be expressed as where λ 1 and λ 2 represent the spatial resolution of images between two time periods. Because the entire Landsat archive is available as an orthorectified product, the co-registration error between the images was found to be negligible and thus not influencing the uncertainties.
The changes in the equilibrium line altitude (ELA) representing the highest observed snowline at the end of the ablation season were observed using satellite data (1990 and 2018) and multi-source DEMs [56][57][58]. This involved delineation of snowline at the end of the ablation season from satellite data using manual digitization. The elevation values of those pixels were extracted from DEM of the snowline. This was completed using a module generating points at 30 m interval in Arc Map 10.1. The arithmetic mean of the extracted elevation values represented ELA.

Historical Sstreamflow Patterns
A time series of discharge data from four observation stations (1-4) along the length of the Lidder River from 1972 to 2018 were used in this analysis and were correlated with glacier changes. Additionally, discharge data from 1986 to 2017 of one station (5) having its contribution from the G2 glacier that drains into the Sindh tributary of Jhelum was also analyzed in this study ( Figure 1). The historical streamflow data were procured from the Department of Irrigation and Flood Control, Srinagar. The statistical significance of discharge trends for all the four stations was also analyzed using the Mann-Kendall test (details in [59,60]) in Trend Software (available at https://toolkit.ewater.org.au/). It is pertinent that the Mann-Kendall test has been widely used for estimating the significance of trends in time series of hydrometeorological data across the mountain regions [61][62][63][64]. The Man-Kendall statistic, denoted by S, is computed as where X represents observations of time series and n is length of time series.
For computing the Z statistic, it is imperative to assess the variance of S, V(S), as where t m is the number of ties of extent m.
Water 2020, 12, 1168 The standard test statistic Z is computed as The positive Z values indicate increasing trend and vice versa. The testing of trends is performed at a specific significance level (a), in this case at 90% (a = 0.1) and 99% (a = 0.01). When |Z| > |Z1 − a/2|, the null hypothesis is rejected, indicating that a significant trend exists in the time series. At the significance level of 90% and 99%, the null hypothesis of no trend is rejected if |Z| > 1.645 and |Z| > 2.576, respectively.
The slope of the streamflow data was estimated using the Theil-Sen approach as where b is the slope of the trend and X l is the lth observation.

Downstream Land System Changes
The land system changes downstream were mapped using satellite data of 1980, 1992, 2005, 2013, and 2019. This study primarily looked at how the area under agriculture, orchards, and built-up areas changed between 1980 and 2018. These land-use types were delineated using on-screen digitization of satellite data in Arc Map 10.1 at a 1:10,000 scale. The accuracy of land use types delineated for 1980, 1992, 2005, and 2013 was ascertained from published data [48,65]. The accuracy of 2018 generated data was obtained utilizing the field survey data about land uses. The overall accuracy and land use-wise accuracy was computed using the following formula [66]: where ρ is accuracy, n is the number of correctly classified samples, and N is the total samples taken. Additionally, informal interviews were held with 142 residents to ascertain the possible causes of land system changes in the Lidder catchment. The interview was composed of two questions: (1) "Is the area experiencing any change in climate?" and (2) "Are the changes in the land system caused by climatic or economic considerations?"

Changes in Glacier Snout
The analysis of multi-date satellite images indicated that the snouts of G1 and G2 retreated at a rate of 18.

Changes in Glacier Snout
The analysis of multi-date satellite images indicated that the snouts of G1 and G2 retreated at a rate of 18

Streamflow Patterns
Although discharge data for four gauging sites were obtained for Lidder watershed, the long-term discharge data of one downstream site was available for Sindh watershed (Figure 1). The Aru station (ID-1) located at an elevation of 2436 m asl is the highest discharge observation site located ≈8 km downstream of G1. The analysis of time series of discharge data at Aru (Figure 5a) revealed that the discharge at Aru showed a statistically significant decrease in discharge during the observation period. m, respectively, from 1992 to 2018. The results from ASTER DEM and TanDEMx indicated a mean shift of 125 m (4.8 m a −1 ) and 137 m (5.3 m a −1 ) for the glacier between 1992 and 2018. The average shift of ELA based on all the four DEMs used for the purpose indicated that the ELA of the glacier shifted by 124 m. The satellite data analysis also indicated that the ELA of the larger parts G1 and G2 shifted by 85 and 200 m, respectively, at the rate of 3.7 m a −1 and 7.7 m a −1 , respectively, from 1992 to 2018.  Although discharge data for four gauging sites were obtained for Lidder watershed, the long-term discharge data of one downstream site was available for Sindh watershed (Figure 1). The Aru station (ID-1) located at an elevation of 2436 m asl is the highest discharge observation site located ≈8 km downstream of G1. The analysis of time series of discharge data at Aru (Figure 5a) revealed that the discharge at Aru showed a statistically significant decrease in discharge during the observation period.  The annual discharge at Aru varied from 112 to 1031 cusecs with a mean annual discharge of 447 cusecs. It is pertinent to mention that 17 of the last 20 years of the observation period experienced below-average discharge at Aru. The analysis of time series of discharge data at Batakoot (ID-2), located at an elevation of 1919 m asl and 41 km downstream with respect to G1, indicated statistically insignificant depleting streamflows (Figure 5b). The annual discharge at Batakoot varied from 249 to 2879 cusecs, with an average annual discharge of 1320 cusecs. The higher streamflow at Batakoot compared to Aru is attributed to small tributaries pouring into the Lidder River. The recent streamflows indicated that 15 of the last 20 years of the observation period experienced below-average discharge at Batakoot. Discharge data of downstream station Gur (ID-3) located at an altitude of 1599 m asl and 68 km downstream of G1 snout indicated more pronounced statistically significant depleting streamflows during the observation period (Figure 5c Table 4 below. Both S-score and Z-statistic for the gauging sites confirmed the statistically significant depleting streamflows, however, Theil-Sen slope indicated higher values for Batakoot and Narayanbagh, owing to the higher degree of variability in discharge as compared to other sites.

Quantifying Land System Changes
Our analysis indicated a decrease in the area under agriculture and a pronounced increase in the area orchards and built-up areas ( Figure 6, Table 5). The area under agriculture shrunk by −39.03%   The accuracy of 2018 delineated land use types was 90%, with built-up having the highest accuracy of 92.3% followed by orchard (90.9%) and agriculture (85.7%). The accuracy of older maps of 1980, 1992, 2005, and 2013 was computed by using the data from Rashid et al. (2010;2017c). It was observed that the 1980 land cover map had an accuracy of 82%, whereas the accuracies for 1992, 2005, and 2013 maps were 85%, 87%, and 90%, respectively. The details of accuracy assessment are provided in Table 6. The location of field validation samples for agriculture, orchards, and built-up areas is provided in the Supplementary Materials section (Figures S1-3).  The accuracy of 2018 delineated land use types was 90%, with built-up having the highest accuracy of 92.3% followed by orchard (90.9%) and agriculture (85.7%). The accuracy of older maps of 1980, 1992, 2005, and 2013 was computed by using the data from Rashid et al. (2010;2017c). It was observed that the 1980 land cover map had an accuracy of 82%, whereas the accuracies for 1992, 2005, and 2013 maps were 85%, 87%, and 90%, respectively. The details of accuracy assessment are provided in Table 6. The location of field validation samples for agriculture, orchards, and built-up areas is provided in the Supplementary Materials section (Figures S1-S3).

Discussion
Kolahoi Glacier spread over 10.5 km 2 is considered to be a single glacier [47,[67][68][69]; however, a recent study [51] suggested that the glacier has two parts draining into two watersheds, Lidder and Sindh, of the larger Jhelum basin. The Randolph Glacier Inventory (https://www.glims.org/RGI/), on the other hand, indicates that the glacier has five parts. On the basis of GLIMS definition, this analysis, utilizing information from remotely sensed data and DEMs, suggests that the glacier has four parts (G1, G11, G111, and G2). The glacier G1 has fragmented into G1, G11, and G111, owing to glacier melt over the last five decades. Similar observations of glacier fragmentation have been reported across the Himalayas [70][71][72][73]. On the basis of the historical retreat rate of G11 (39 m a −1 ) and G2 (35 m a −1 ), the glaciers would vanish by 2056 and 2045, respectively. Given the exacerbated rate of melt of Kashmir Himalayan glaciers [26,74], the glaciers may perhaps vanish much earlier. A glacier cannot endure without an accumulation zone [75], and often glaciers that have an avalanche fed area can survive longer [76]. The analysis of the snout retreat of all the four parts of Kolahoi Glacier also indicated an accelerated retreat of the glacier across the analysis period. This could perhaps be attributed to the warming environment that increases ablation and [77,78] forces change in the form of precipitation [21,31,79] and potentially high ambient black carbon concentration, such as that observed near glaciers in the Kashmir region [80,81]. The accelerated rate of retreat of glaciers in northwest Himalaya has been very well documented in many recent studies [82][83][84]. The melt of Kolahoi Glacier has strongly impacted the streamflows of the two largest tributaries, Lidder and Sindh, draining into Jhelum. Our analysis of time series of streamflow data at five discharge stations from the 1970s to the present suggested depleting annual streamflows in the two watersheds of Jhelum basin, having their contribution from Kolahoi Glacier. Whereas Aru, Kirkadal, and Gur observed statistically significant depleting streamflows, Batakoot and Narayanbagh experienced a statistically insignificant decrease ( Table 4). The statistically significant depleting streamflows have also been reported by earlier works [37,85] in the Kashmir region. The depleting streamflows do not appear to have severely affected the upstream stations (Aru and Batakoot). Although statistically significant depleting flows have been observed at Aru, the trends at Batakoot are statistically insignificant. This could be due to the presence of seasonal snow in high altitude regions for the most part of the year [86], as well as contribution from other smaller glaciers for Batakoot station. However, the consistently depleting streamflows have played a role in the significant land system changes in the downstream areas of Lidder watershed. This is manifested in the satellite data analysis focusing on changes in agriculture, orchards, and built-up areas in the watershed. The analysis indicated that area under irrigation intensive rice paddies shows persistent signs of decline, whereas areas under orchards and built-up areas are increasing at unprecedented rates. The massive land system changes in the Kashmir region concerning agricultural lands, orchards, and built-up areas have been documented in many studies [39,40,87]. Although few studies suggest that the depleting streamflows are a significant factor driving the land system changes in Kashmir [88,89], the changes are also perhaps taking place due to economic considerations, especially given the fact that the production from apple orchards fetches more income compared to irrigation-intensive agriculture [46,[90][91][92][93]. The increase in the built-up area in Lidder watershed has not only resulted in shrinkage of land under agriculture but has perhaps also affected the hydrological connectivity, as a lot of smaller streams have dried due to land-filling and the associated expansion of built-up areas across Kashmir [94].
On the basis of the informal interviews with 142 people during field visits to the Lidder watershed in 2019, the majority 64% (91 people) attributed economic reasons for agriculture to orchard conversions, 30% (42 people) attributed the changes to depleting streamflows, while 6% (9 people) were clueless. There was unanimity among people that the area is experiencing climatic change. Our analysis suggests that although climatic change-induced streamflow depletion is possibly one of the causes of land system dynamics in Lidder watershed, the economic considerations perhaps play a much bigger role in forcing the massive unplanned land system changes in Lidder watershed. The depleting streamflows in the Kashmir region can massively impact water supply and demands, especially agriculture, hydropower, and domestic water use, given the changing climate prevalent over the region; however, it needs to be properly researched.

Conclusions
Kolahoi Glacier has rapidly receded due to climatic changes prevalent over the Kashmir region. The glacier has lost ≈23% area since 1962 and has fragmented into smaller parts. The snout retreat rates also suggest that the glacier has been in an imbalanced state between 1962 and 2018 and is not approaching equilibrium. The exacerbated snout retreat rates of the four parts of the glacier are driven by increasing warming that has impacted the snow accumulation on the glacier. The significant upward shift in the ELA of the glacier indicates an expansion of the ablation area that leads to a negative mass balance of the glacier. The loss of glacier area has exceeded the increase in the ablation area, translating into statistically significant depleted streamflows that have resulted in land system changes downstream. The negative glacier mass balance prevalent over the region will continue to drive the glacier recession and further decline in glacier runoff [95]. Further study can explain the differential rate of depleting streamflows among the gauging stations by streamflow partitioning using either isotopic signatures or assessing the water balance. The area under irrigation-intensive agriculture has shrunk by 39%, whereas the orchards have expanded by 177% from 1980 to 2017. On the basis of preliminary informal interviews, our analysis suggested that the land system changes are a cumulative effect of depleting streamflows and economic dividends from apple orchards. The consistently depleting streamflows will impact agricultural production, hydropower generation, and domestic water use in the region.