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

Abstract

This study reports the changes in glacier extent and streamflow 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 flowing 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⁻¹ and 16.4 m a⁻¹, 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 streamflows measured at five sites showed statistically significant 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 streamflow.

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⁻¹ [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⁻¹) 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⁻¹ [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.

2. Materials and Methods

2.1. Study Area Description

This study was carried out on Kolahoi Glacier, spread between Lidder and Sind watersheds, in the greater Himalayan mountain range of Kashmir valley, India (

Table 1. Characteristics of the Lidder and Sind watersheds.

Characteristic	Lidder	Sind
Latitude	33.73°–34.26° N	34.12°–34.45° N
Longitude	75.08°–75.54° E	74.60°–75.6° E
Minimum elevation	1585 m asl	1570 m asl
Mean elevation	3060 m asl	3200 m asl
Maximum elevation	5300 m asl	5400 m asl
Area	≈1230 km2	≈1560 km2

, Figure 1).

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].

2.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. Details of the datasets used in this study.

Dataset	Path/Row	Acquisition/Release Date	Spatial Resolution
A. Satellite Data
Landsat Multi-Spectral Scanner (MSS)	160/36	18 September 1980	60 m
Landsat Thematic Mapper (TM)	149/36	7 August 1990	30 m
	149/36	24 October 1992	30 m
Linear Imaging Self-Scanning Sensor	92/46	19 October 2005	23.5 m
	92/47	19 October 2005	23.5 m
	93/46	24 October 2005	23.5 m
Landsat Enhanced Thematic Mapper+ (ETM+)	149/36	27 August 2000	30 m
Landsat 8 Operational Land Imager	149/36	25 October 2013	30 m
	148/36	29 September 2018	30 m
High-resolution Basemap imagery	ESRI Arc Map 10.1	27 September 2018	≈1 m
B. Topographic data
Survey of India Topographic map	43N/8	1962	1:50,000 scale
Shuttle Radar Topography Mission (SRTM) DEM	N34E075	2000	30 m
Advanced Land Observing Satellite ALOS) DEM	N34E075	2006–2011	30 m
Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) GDEM v2	N34E075	2012	30 m
TerraSAR-X add-on for Digital Elevation Measurement (TanDEMx)	N34E075	2011–2014	30 m
C. Discharge data	-	1972–2018	Point data
D. Ground Truth data (GPS validation)	-	2017–2018	Point data
E. Interviews	-	2017–2018	Point data

).

2.2.1. Quantifying Glacier Changes

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

$$E_A=n\times \lambda \times \sqrt{{(\lambda /2)}^2+{\epsilon }^2}$$

(1)

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

$$E_{AC}=\sqrt{E_{A1}{}^2+E_{A2}{}^2}$$

(2)

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

$$E_{SC}=\sqrt{{\lambda }_1{}^2+{\lambda }_2{}^2}+\epsilon$$

(3)

where λ₁ and λ₂ 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.

2.2.2. 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

$$S={\displaystyle \sum _{i=1}^{n-1}{\displaystyle \sum _{j=i+1}^n\mathrm{sgn}(X_j-X_i)}}$$

(4)

where X represents observations of time series and n is length of time series.

$$\mathrm{sgn}(X_j-X_i)=\{\begin{array}{c}1\\ 0\\ -1\end{array}\begin{array}{c}if(X_j-X_i)>0\\ if(X_j-X_i)=0\\ if(X_j-X_i)<0\end{array}$$

(5)

For computing the Z statistic, it is imperative to assess the variance of S, V(S), as

$$V(S)=\frac{n(n-1)(2n+5)-{\displaystyle \sum _{m=1}^n{t}_{m}m(m-1)(2m+5)}}{18}$$

(6)

where t_m is the number of ties of extent m.

The standard test statistic Z is computed as

$$Z=\{\begin{array}{cc}\frac{S-1}{\sqrt{V(S)}}& S>0\\ 0& S=0\\ \frac{S+1}{\sqrt{V(S)}}& S<0\end{array}$$

(7)

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

$$b=Median\left(\frac{X_j-X_l}{j-l}\right)\forall l<j$$

(8)

where b is the slope of the trend and X_l is the lth observation.

2.2.3. 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]:

$$\rho =\left(\frac{n}{N}\right)\times 100$$

(9)

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?”

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² in 2018. Analysis of 2018 satellite data suggests that the four parts of the Kolahoi Glacier are spread over 5.32 ± 0.43 km² (G1), 0.5 ± 0.08 km² (G11), 0.51 ± 0.08 km² (G111), and 4.17 ± 0.31 km², separately (

Table 3. Changes in the glacier-wise and cumulative area of Kolahoi Glacier from 1962 to 2018.

	Glacier Area (km2)
Year	G1	G2	G11	G111	Total
1962	8.74	4.99	NA	NA	13.73
1992	6.411	4.89	1.85	NA	13.16
2000	6.211	4.76	0.69	0.72	12.38
2013	5.633	4.25	0.53	0.54	10.96
2018	5.316	4.17	0.5	0.51	10.49

). 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²) 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²) and G11 (1.85 ± 0.14 km²)—by 1990.

Between 1990 and 2000, the G11 part further fragmented into two parts: G11 (0.69 ± 0.07 km²) and G111 (0.72 ± 0.09 km²). The areas of G1, G11, G111, and G2 were 5.63 ± 0.46 km², 0.53 ± 0.07 km², 0.54 ± 0.09 km², and 4.25 ± 0.3 km², respectively, in 2013. By 2018, the G1, G11, G111, and G2 deglaciated to 5.32 ± 0.43 km², 0.5 ± 0.08 km², 0.51 ± 0.08 km², and 4.17 ± 0.31 km², 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⁻¹ and 16.4 ± 0.69 m a⁻¹, 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⁻¹ and 19.1 ± 1.4 m a⁻¹, respectively, from 1962 to 1992. The snout recession decelerated during 1992–2000 to 6.1 ± 5.3 m a⁻¹ (G1) and 2.8 ± 5.3 m a⁻¹ (G2), indicating stability. The recession again picked up during 2000-2013 to 22.6 ± 3.2 m a⁻¹ and 20.3 ± 3.2 m a⁻¹ for G1 and G2, respectively. The highest rates of snout retreat at 24.4 ± 6 m a⁻¹ for G1 and 27.1 ± 6 m a⁻¹ 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⁻¹ 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⁻¹ 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⁻¹) 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⁻¹) 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⁻¹) and 137 m (5.3 m a⁻¹) 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⁻¹ and 7.7 m a⁻¹, 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. Significance of trends in streamflows in Lidder and Sind watersheds.

Gauging Station	S-Score	Z-Statistic	Trend	Theil–Sen Slope
Aru	−191	−1.80	Significant (at a < 0.1)	−3.93
Batakoot	−91	−0.85	Insignificant (at a = 0.1)	−6.19
Gur	−291	−2.75	Significant (at a < 0.01)	−2.39
Kirkadal	−389	−3.67	Significant (at a < 0.01)	−4.36
Narayanbagh	−92	−1.48	Insignificant (at a = 0.1)	−12.66

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. Area changes in land use between 1980 and 2018.

Land Use Type	Area (km2)
	1980	1992	2005	2013	2018	Area Change	% Change
Agriculture	165.65	153.9	124.66	101	101	−64.65	−39.03
Orchard	21.67	41.7	53.85	58.45	60	38.33	176.88
Built-up	13.19	16.54	27.56	52	76	62.81	476.19

). The area under agriculture shrunk by −39.03% (64.65 km²) 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 km²) 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 km²) 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. Accuracy assessment of land use maps of 1980, 1992, 2005, 2013, and 2018.

Year		Land Use Types			Overall Accuracy
		AG	BU	OR
2018	n	24	36	30	90%
	N	28	39	33
	ρ	85.71	92.3	90.91
2013	n	18	17	18	83.33%
	N	20	20	20
	ρ	90	85	90
2005	n	17	18	19	90%
	N	20	20	20
	ρ	85	90	95
1992	n	18	18	16	86.66%
	N	20	20	20
	ρ	90	90	80
1980	n	16	17	17	83.33%
	N	20	20	20
	ρ	80	85	85

. 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² 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⁻¹) and G2 (35 m a⁻¹), 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.

5. 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.