1. 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 to decrease in the streamflows in the Himalayan Rivers [
7,
27,
28,
29], which could affect water availability downstream [
30]. Studies suggest that the glacier recession in the Kashmir valley has already resulted in the depleted streamflows downstream [
31,
32]. The land system changes in the region [
33,
34,
35,
36] have been linked to depleting streamflows, economic considerations, and unplanned land transformation.
The use of remote sensing data for quantifying land system changes over the Kashmir region has been widely documented [
37,
38,
39,
40]. At the same time, there is substantial scientific literature detailing the glacier retreat prevalent over the region. A recent study indicated that the glaciers in the Ladakh region are retreating at 0.55% a
−1 [
41], whereas another study [
42] reported few stable and advancing glaciers in the Zanskar region of Jammu and Kashmir. Similar retreat estimates have been put forth for the Zanskar region for the glacier area changes carried between 1989 and 2007 [
43]. However, very conservative area changes (0.16% a
−1) have been reported for the neighboring Suru basin [
44] for the glacier area changes assessed for the 1977–2017 period. In another study carried out in the Lidder watershed of the Jhelum basin, the glaciers were reported to be shrinking at 0.51% a
−1 [
17]. Glacier mass change over the western Himalayas has accelerated from 0.33 to 0.5 m w.e. (water equivalent) per year between 1970–2000 and 2000–2010, respectively [
45]. There are no studies to our knowledge that link the glacier recession and streamflow patterns with downstream land system changes. In this backdrop, this study assessed the changes in the extents of Kolahoi Glacier, the largest glacier in Kashmir Valley; discharge trends and associated land system changes between 1962 and 2018 using remote sensing imageries; and in situ streamflow data at the catchment scale. We assessed area, snout, and equilibrium line altitude (ELA) in the four parts of Kolahoi using a topographic map of 1962 and remote sensing data from 1979 to 2018. The land system changes in terms of changes in agriculture, orchards, and settlements were quantified from 1980 to 2018. Additionally, historical streamflow data at five sites (four in Lidder and one in Sind) were analyzed for trends in discharge in Lidder and Sind watershed. This was aimed to ascertain the impact glacier retreat and associated hydrological changes on the land system changes downstream.
3. Results
3.1. Glacier Area Changes
This analysis suggested that the Kolahoi Glacier was spread over an area of 10.5 ± 0.84 km
2 in 2018. Analysis of 2018 satellite data suggests that the four parts of the Kolahoi Glacier are spread over 5.32 ± 0.43 km
2 (G1), 0.5 ± 0.08 km
2 (G11), 0.51 ± 0.08 km
2 (G111), and 4.17 ± 0.31 km
2, separately (
Table 3). We assessed changes in the four parts of the glacier in terms of area, snout position, and ELA from 1962 to 2018. Our analysis indicated that the glacier lost 23.6% (3.41 km
2) of its 1962 area (
Figure 2). Although the glacier in 1962 was divided into two parts (G1 and G2 draining into Lidder and Sind watersheds, respectively), the deglaciation of G1 resulted in its further fragmentation into two parts—G1 (6.41 ± 0.5 km
2) and G11 (1.85 ± 0.14 km
2)—by 1990.
Between 1990 and 2000, the G11 part further fragmented into two parts: G11 (0.69 ± 0.07 km2) and G111 (0.72 ± 0.09 km2). The areas of G1, G11, G111, and G2 were 5.63 ± 0.46 km2, 0.53 ± 0.07 km2, 0.54 ± 0.09 km2, and 4.25 ± 0.3 km2, respectively, in 2013. By 2018, the G1, G11, G111, and G2 deglaciated to 5.32 ± 0.43 km2, 0.5 ± 0.08 km2, 0.51 ± 0.08 km2, and 4.17 ± 0.31 km2, respectively. The analysis indicated that the retreat rate of Kolahoi Glacier tripled from 0.26% per year to 0.76% per year between 1962 and 2000, and 2000 and 2018, respectively. G1 and G2 lost 27.57% and 16.36%, respectively, from 1962 to 2018, of their 1962 area. Our analysis indicated that G1 retreated at a relatively higher pace (0.49% per year) compared to G2 (0.29% per year). From 1962 to 2000, G1 and G2 lost an area at a rate of 0.33% per year and 0.12% per year, respectively. The retreat rate accelerated to 0.82% per year for G1 and 0.65% per year for G2 between 2000 and 2018.
3.2. 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.3 ± 0.69 m a
−1 and 16.4 ± 0.69 m a
−1, respectively, between 1962 and 2018 (
Figure 3). However, the rate of retreat was variable across different periods. The snout of G1 and G2 retreated at 19.5 ± 1.4 m a
−1 and 19.1 ± 1.4 m a
−1, respectively, from 1962 to 1992. The snout recession decelerated during 1992–2000 to 6.1 ± 5.3 m a
−1 (G1) and 2.8 ± 5.3 m a
−1 (G2), indicating stability. The recession again picked up during 2000-2013 to 22.6 ± 3.2 m a
−1 and 20.3 ± 3.2 m a
−1 for G1 and G2, respectively. The highest rates of snout retreat at 24.4 ± 6 m a
−1 for G1 and 27.1 ± 6 m a
−1 for G2 were observed during 2013–2018. Although G11 formed in 1990, its snout retreated by 2200 ± 39 m at a rate of 39.28 ± 0.69 m a
−1 between 1962 and 2018. Similarly the fragmentation of G11 resulted in the formation of G111 in 2000, whose snout retreated by 917 ± 42.4 m at a rate of 35.27 ± 1.6 m a
−1 from 1992 to 2018.
3.3. Changes in ELA
The ELA of all the four parts of the Kolahoi Glacier showed consistent signs of an upward shift from 1990 to 2018 (
Figure 4). The ELA changes using satellite images and SRTM DEM indicated that the ELA of Kolahoi Glacier shifted upward by 120 m (4.6 m a
−1) from 4420 m asl in 1990 to 4540 m asl in 2018. The ELA of G1, G11, G111, and G2 shifted by 88, 90, 103, and 194 m, respectively, during the same period. The analysis based on ALOS DEM indicated the ELA shift of 116 m (4.5 m a
−1) from 4442 to 4539 m asl for the glacier. The ELA of G1, G11, G111, and G2 shifted by 80, 93, 102, and 190 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.
3.4. 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.
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). The annual discharge at this site varied from 59 to 366 cusecs, with an average annual discharge of 183 cusecs. Our analysis indicated that 18 of the last 20 years of the observation period at Gur experienced below-average discharge. Discharge data of another downstream station Kirkadal (ID-4) located at an altitude of 1592 m asl and 75 km downstream of G1 snout indicated statistically significant highest depleting streamflows in the Lidder watershed during the observation period (
Figure 5d). The annual discharge at this site varied from 49 to 542 cusecs, with an average annual discharge of 215 cusecs. Our analysis indicated that 17 of the last 20 years of the observation period at Kirkadal experienced below-average discharge. Discharge data of the only downstream station, Narayanbagh (ID-5), located at an altitude of 1583 m asl and 98 km downstream of G2 snout, also indicated non-significant depleting streamflows in the Sindh watershed during the observation period from 1986 to 2017 (
Figure 5e). The annual discharge at this site varied from 943 to 3528 cusecs, with an average annual discharge of 1633 cusecs. Our analysis indicated that 17 of the last 20 years of the observation period at Narayanbagh experienced below-average discharge. The details of the statistical trend analysis are provided in
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.
3.5. 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% (64.65 km
2) at a rate of 1.03% per year between 1980 and 2018. However, the rate of loss was not uniform. The agriculture declined at 0.59% per year from 1980 to 1992. The loss accelerated to 1.36% per year and 1.79% per year from 1992 to 2005 and 2005 to 2013, respectively.
The area under agriculture has remained stable since 2013 in the Lidder watershed. The area under orchards and built-up areas indicated a massive increase during the same period. Orchard area increased by 176.88% (38.33 km2) at a rate of 4.65% per year, and the rate of increase in the orchard area slowed across the analysis period from 7.7% per year from 1980 to 1992, to 4.31% per year and 2.65% per year between 1992 and 2005, and 2005 and 2013, respectively. The rate of expansion of orchards further declined to 1.43% per year between 2013 and 2018. The built-up area expanded by 476.19% (62.81 km2) at a rate of 12.53% per year between 1980 and 2018. The rate of increase in the built-up area showed a persistent increase across different periods.
The built-up area increased by 2.12% per year between 1980 and 1992, and expanded at 6.43% per year from 1992 to 2005. However, the pace of increase escalated to 23.16% per year and 36.39% per year for the 2005–2013 and 2013–2018 periods, respectively.
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).
4. 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.