1. Introduction
The majority of mountain glaciers worldwide are undergoing rapid mass loss in response to climate warming (Dehecq and others, Reference Dehecq2019; Hugonnet and others, Reference Hugonnet2021; Rounce and others, Reference Rounce2023). Glaciers in High Mountain Asia (HMA) show substantial variability in mass change between regions, with those in the Himalaya losing mass at accelerating rates, while those in the Karakoram and Pamir have remained relatively stable, at least until recently (Bolch and others, Reference Bolch2012; Dehecq and others, Reference Dehecq2019; Hugonnet and others, Reference Hugonnet2021; Nie and others, Reference Nie2021; Millan and others, Reference Millan, Mouginot, Rabatel and Morlighem2022). Understanding intra-regional pattern in glacier mass loss is important because glacial meltwater is an essential resource for hydropower, agriculture and sanitation for 20% of the global population (Immerzeel and others, Reference Immerzeel, Van Beek and Bierkens2010; Pritchard, Reference Pritchard2019; Rounce and others, Reference Rounce2023). Reducing the uncertainties associated with regional projections of glacier mass change during the 21st century require improved representation of the processes, feedbacks and tipping points that are affecting glacier evolution (Nie and others, Reference Nie2021).
Intra-regional variability in glacier mass loss across HMA suggests that observed changes in patterns of ice flow do not result from climate change alone (King and others, Reference King, Dehecq, Quincey and Carrivick2018) but are modulated by glaciological processes including the feedbacks between supraglacial debris and ice flow, lakes in contact with glacier margins (hereafter; ice-contact lakes) and surge-type behaviour (Quincey and others, Reference Quincey2011, Reference Quincey, Glasser, Cook and Luckman2015; Rowan and others, Reference Rowan, Egholm, Quincey and Glasser2015; Brun and others, Reference Brun2019). Ice-contact lakes influence glacier behaviour through two main mechanisms; (1) subaqueous melt and (2) calving of the glacier terminus, which together promote faster ice flow and dynamic thinning through the lower part of the glacier compared to climatically equivalent land-terminating glaciers (Carrivick and Tweed, Reference Carrivick and Tweed2013; King and others, Reference King, Quincey, Carrivick and Rowan2017, Reference King, Dehecq, Quincey and Carrivick2018; Carrivick and others, Reference Carrivick, Tweed, Sutherland and Mallalieu2020; Zhang and others, Reference Zhang2023). Ice-contact lakes that are shallow and/or in the early stages of evolution may have limited impact on ice flow, but can enhance mass loss through frontal ablation (Carrivick and Tweed, Reference Carrivick and Tweed2013). When ice-contact lakes reach a sufficient depth relative to the ice thickness, the effective pressure at the glacier bed (i.e. the difference between subglacial water pressure and ice overburden pressure; Harper and others, Reference Harper, Humphrey, Pfeffer and Lazar2007) is reduced, ice flow increases and terminus flotation may occur (Benn and others, Reference Benn, Warren and Mottram2007). These mechanisms can increase lake depth and/or cause a glacier to recede into deeper water, setting up a positive feedback that increases the proportion of the ice margin that is in contact with the lake (Benn and others, Reference Benn, Warren and Mottram2007; King and others, Reference King, Dehecq, Quincey and Carrivick2018; Pronk and others, Reference Pronk, Bolch, King, Wouters and Benn2021).
The area occupied by glacial lakes in HMA has expanded over the last 30 years, with an increase of over 45% between 1990 and 2018 (Shugar and others, Reference Shugar2020). This trend is expected to continue, or accelerate, as glaciers continue to lose mass leaving behind moraine dams and newly exposed subglacial overdeepenings (Furian and others, Reference Furian, Maussion and Schneider2022). The effects of this lake expansion on glacier mass balance has been characterised by previous remote-sensing studies, which have focussed on the Eastern Himalaya and Central Himalaya, Sikkim, the Everest region and along the wider Himalayan arc (Basnett and others, Reference Basnett, Kulkarni and Bolch2013; King and others, Reference King, Quincey, Carrivick and Rowan2017; Brun and others, Reference Brun2019; Tsutaki and others, Reference Tsutaki2019; Liu and others, Reference Liu2020; Lee and others, Reference Lee2021). King and others (Reference King, Quincey, Carrivick and Rowan2017) observed 32% more negative mass balance for lake-terminating glaciers compared to land-terminating glaciers in the Everest region between 2000 and 2015. Across HMA, rates of mass loss for lake-terminating glaciers were 18–97% more negative than regional means (Brun and others, Reference Brun2019). In terms of ice velocity, King and others (Reference King, Dehecq, Quincey and Carrivick2018) found that lake-terminating glaciers in the Everest region demonstrated contrasting trends over time, with one group accelerating between 2000 and 2015, and the other decelerating. Pronk and others (Reference Pronk, Bolch, King, Wouters and Benn2021) found that, between 2017 and 2019, glaciers in contact with lakes flowed at around twice the rate of glaciers without lakes, but that debris cover affected the magnitude of these differences.
This paper seeks to build on these previous studies and address intra-regional patterns in glacier behaviour in further detail. Specifically, we integrate regional datasets of ice-contact lake characteristics, glacier surface elevation change and glacier velocity to investigate; (i) the extent to which the evolution of individual glaciers across the Himalaya differs depending on their terminus environment and surface characteristics, and (ii) the extent to which the influence of ice-contact lakes propagates up-glacier, in terms of changes in surface elevation and ice velocity.
2. Datasets and methods
2.1 Selection of glaciers
We selected 352 representative lake-terminating and land-terminating glaciers >1 km2 across the Himalaya within three sub-regions based on the Global Terrestrial Network for Glaciers (GTN-G) glacier sub-regions (Fig. 1). The Western Himalaya sub-region is in the transition zone between the Westerly and monsoon-influenced parts of the range, and the Central Himalaya sub-region and Eastern Himalaya sub-region are in the monsoon-dominated part of the range (Bookhagen and others, Reference Bookhagen, Thiede and Strecker2006). The glaciers within each sub-region differ in characteristics such as surface area, aspect and geometry. Since lake-terminating glacier density varies spatially, the population size is inconsistent between sub-regions; we selected 78 glaciers in Western Himalaya (22% of total sample), 89 in Central Himalaya (25% of total sample) and 185 in Eastern Himalaya (53% of total sample). The population size is inconsistent between sub-regions, but closely reflects the relative distribution of glacial lakes throughout Western (26%), Central (19%) and Eastern Himalaya (54%) (Shugar and others, Reference Shugar2020).
Map of Western, Central and Eastern Himalaya with glacier sample numbers of lake- and land-terminating glacier types.
Glaciers were classified by terminus and surface cover types as either lake-terminating debris-covered, lake-terminating clean-ice, land-terminating debris-covered or land-terminating clean-ice glacier types following Lee and others (Reference Lee2021). Glaciers were classified as lake-terminating when contact between the glacier terminus and an ice-contact lake was evident in both the most recent optical imagery (Landsat 8 or Sentinel-2) and in imagery prior to 2000 (Landsat 5) broadly following the approach of King and others (Reference King, Dehecq, Quincey and Carrivick2018). Glaciers were classified as debris-covered when the extent of supraglacial debris was >7% of total glacier area following the approach of Herreid and Pellicciotti (Reference Herreid and Pellicciotti2020).
2.2 Analysis across normalised glacier elevations
To account for topographic differences between each of the three study sub-regions, we generated ten normalised elevation bins per glacier using outlines from the Randolph Glacier Inventory (RGI Consortium, 2017) and the ASTER DEM v.3 (NASA and others, 2019) (Fig. S2). Glacier mass change and velocity were analysed across these bins following similar approaches to King and others (Reference King, Bhattacharya, Bhambri and Bolch2019) and Hugonnet and others (Reference Hugonnet2021). Bins were numbered from 0.1 to 1.0 according to their distance from glacier termini. The normalised elevation bins were not altered to account for glacier area change to maintain consistency between time steps in our data analysis.
2.3 Surface elevation change (data)
We extracted surface elevation change from the dataset of Hugonnet and others (Reference Hugonnet2021). This dataset is derived using a stacked DEM approach, where a Gaussian process regression is applied to all available observations at each pixel through time, rather than simply subtracting one DEM from another. Hugonnet and others (Reference Hugonnet2021) provide gridded data at 100 m spatial resolution, calculated over 5-year observation periods between 2000 and 2019 (2000–2004, 2005–2009, 2010–2014, 2015–2019). Glacier surface elevation change was calculated over both the entire glacierised area and also within the ten normalised elevation bands per glacier (section 2.2). Data points were removed from the analysis where the 19-year elevation change was larger than five times the normalised median absolute deviation within each elevation bin following the approach of Hugonnet and others (Reference Hugonnet2021).
2.4 Glacier velocity (data)
Annual velocity data were extracted for all glaciers for the period 2000–2018 from NASA's Inter-Mission Time Series of Land Velocity and Elevation (ITS_LIVE) product, derived using the auto-RIFT processing scheme applied to all Landsat 4, 5, 7 and 8 images within the time period (Gardner and others, Reference Gardner, Fahnestock and Scambos2019). These annual composite velocities are at 240 m resolution and were created by taking the error-weighted average of all image pairs that have a time span of <546 days and a centre-date that falls within the relevant calendar year between 2000 and 2018 (Dehecq and others, Reference Dehecq2019; Gardner and others, Reference Gardner, Fahnestock and Scambos2019). Velocity data were filtered to remove velocity pixels that had errors >5 m a–1 following the approach of Dehecq and others (Reference Dehecq2019). However, velocity filtering in this study deviates from that implemented by Dehecq and others (Reference Dehecq2019), where they additionally removed pixels with velocity lower than 5 m a–1. We found that implementing this filtering step removed over 78% of velocity pixels within glacier outlines, many of which were likely robust (with a relative error <50%), which precluded the analysis of velocity across normalised glacier elevation (section 2.2).
We calculated the annual ablation area velocity anomaly for each glacier type to normalise velocity change between glaciers and to allow trends to be identified (cf. Dehecq and others, Reference Dehecq2019). Specifically, we used linear regressions of glacier ablation area velocity with time to identify patterns of speed-up or slow-down between glacier types. The glacier ablation area was defined as elevation bins 0.1–0.6 (section 2.2), where glacier median elevation was located within bin 0.5 (51%) or bin 0.6 (49%) for the glacier sample (Dehecq and others, Reference Dehecq2019). The velocity anomaly was defined as the difference between the annual ablation area velocity of an individual glacier and the annual mean velocity of all glaciers in our sample (n = 352), where a positive value indicates ice flow faster than the sample mean and a negative value indicates ice flow slower than the sample mean (cf. Dehecq and others, Reference Dehecq2019).
2.5 Ice-contact lake area (data)
All available ice-contact lake outlines were collated for our sample glaciers, providing outlines for periods 1990–1999, 2000–2004, 2005–2009, 2010–2014 and 2015–2018 (Shugar and others, Reference Shugar2020). Repeat lake outlines for 1990 and 2018 from Wang and others (Reference Wang2020) were used where lake outlines from Shugar and others (Reference Shugar2020) were not available. Lake area change was quantified as a percentage based on the earliest and latest observations available (i.e. 1990 and 2018).
3. Results
Our results are presented as the median of each glacier type (unless otherwise specified) to identify differences in glacier change over time related to the characteristics of the glacier terminal and or surface environment. Rates of change in velocity and surface elevation change are given as median values for individual glaciers. The median value for each type of glacier (e.g. lake-terminating glaciers, debris-covered glaciers) is termed the ‘glacier type average’.
3.1 Glacier surface elevation change
During the 19-year observation period (2000–2019), only glaciers in the Eastern Himalaya (n = 185, p < 0.001) experienced statistically different changes in glacier-wide surface elevation change between lake-terminating glaciers (−0.68 ± 0.05 m a−1) and land-terminating glaciers (−0.54 ± 0.04 m a−1) (Fig. 2). There were no statistically significant differences in glacier-wide surface elevation change between lake-terminating and land-terminating in the Western or Central Himalaya. In the Western Himalaya (n = 78, p = 0.42), lake-terminating glaciers exhibited a mean surface elevation change rate of −0.54 ± 0.09 m a−1, while land-terminating glaciers exhibited a mean surface elevation change rate of −0.49 ± 0.07 m a−1. In the Central Himalaya, surface elevation change rates of lake-terminating (−0.72 ± 0.07 m a−1) and land-terminating glaciers (−0.70 ± 0.07 m a−1) were also comparable (n = 89, p = 0.66). Furthermore, lake-terminating and land-terminating glaciers with different surface characteristics (debris-covered and clean-ice) showed no statistically significant difference in surface elevation change in Western and Central Himalaya (Table S1), although there is visual evidence of the lake effect within the Eastern Himalaya dataset (Fig. 2).
Boxplots summarising rate of mean surface elevation change between 2000 and 2019 for glaciers by terminus type in (a) Western Himalaya (n = 38 and 40), (b) Central Himalaya (n = 45 and 44) and (c) Eastern Himalaya (n = 93 and 92) and glaciers by terminus and surface cover type for (d) Western Himalaya, (e) Central Himalaya and (f) Eastern Himalaya.
In all sub-regions, the mean surface elevation change was negative even at the highest elevations where glaciers may be expected to accumulate mass (Fig. 3). For example, <17% of glaciers exhibited positive mean surface elevation change in the highest elevation bin.
Mean rate of surface elevation change of land-terminating and lake-terminating glaciers between 2000 and 2019 across normalised glacier elevation for (a) and (b) Western Himalaya (n = 38, n = 40), (c) and (d) Central Himalaya (n = 45, n = 44) and (e) and (f) Eastern Himalaya (n = 93, n = 92).
In the Western and Central Himalaya, all glacier groups irrespective of terminus or surface cover type exhibited increasingly negative surface elevation change rates with distance downglacier (Fig. S3). This was also the case in the Eastern Himalaya, with the exception of land-terminating debris-covered glaciers where surface elevation change became less negative in the lowest elevation bin. In the Central Himalaya, the mean surface elevation change of lake-terminating and land-terminating glaciers was not statistically different at any elevation (Fig. 3, Table S3). Lake- and land-terminating glaciers in the Western Himalaya exhibit statistically significant differences in mean surface elevation change within the lowest elevation bin where lake-terminating glaciers had higher rates of mean surface elevation change (−1.52 ± 0.26 m a−1) than land-terminating glaciers (−1.15 ± 0.21 m a−1), a difference of 0.37 ± 0.05 m a−1 (p = 0.03) (Fig. 3, Table S2).
In the Eastern Himalaya, the mean surface elevation change of lake- and land-terminating glaciers was statistically different in the three lowest elevation bins (0.1–0.3) (Fig. 3, Table S4). In the lowermost elevation bin, lake-terminating glaciers experienced surface elevation change at a rate of −1.80 ± 0.18 m a−1, which was 0.82 ± 0.05 m a−1 more negative than land-terminating glaciers at this elevation (p < 0.001). In the second (0.2) and third (0.3) elevation bins, lake-terminating glaciers exhibited greater surface elevation change rates with differences of 0.39 ± 0.03 m a−1 (p < 0.001) and 0.16 ± 0.01 m a−1 (p = 0.01).
3.2 Glacier velocity and velocity anomaly
In all three sub-regions, lake- and land-terminating glaciers exhibited median velocities of <5 m a−1 between 2000 and 2018, although individual glaciers flowed up to speeds of 44 m a−1. In the Eastern Himalaya, lake-terminating glaciers showed a significantly higher median velocity (3.26 ± 0.56 m a–1) in comparison to land-terminating glaciers (2.22 ± 0.34 m a–1) in elevation bins between 0.1 and 0.8 (Fig. 4). In the Western and Central Himalaya, lake-terminating glaciers exhibited statistically higher median velocities than land-terminating glaciers only in the lowest two (0.1–0.2) and three elevation bins (0.1–0.3) (Fig. 4). In the Western Himalaya, lake-terminating glaciers exhibited median velocities of 1.88 ± 0.55 m a–1, while land-terminating glaciers exhibited median velocities of 1.60 ± 0.35 m a–1 in the lowest two elevation bins. While in the Central Himalaya, lake-terminating glaciers exhibited median velocities of 2.18 ± 0.80 m a–1, while land-terminating glaciers exhibited median velocities of 1.86 ± 0.39 m a–1, in the lowest three elevation bins.
Median velocity across normalised glacier elevation between 2013 and 2018 for lake-terminating and land-terminating glaciers in (a) Western Himalaya, (b) Central Himalaya and (c) Eastern Himalaya, where 1 is the maximum normalised glacier elevation. The interquartile range is indicated by the shading for each line.
In the Western and Central Himalaya, negative trends of median velocity anomaly were observed across all glacier types, regardless of terminus environment. In the Western Himalaya, lake-terminating and land-terminating glaciers exhibited the most negative trends of median velocity anomaly of −0.32 ± 0.04 m a–1 decade–1 (R 2 = −0.38, p = 0.1) and −0.52 ± 0.04 m a–1 decade–1 (R 2 = −0.52, p = 0.02) between 2000 and 2018 (Fig. 5). Lake- and land-terminating glaciers in the Central Himalaya exhibited less negative trends of median velocity anomaly of −0.04 ± 0.02 m a–1 decade–1 (R 2 = −0.09, p = 0.7) and −0.11 ± 0.03 m a–1 decade–1 (R 2 = −0.2, p = 0.42), respectively, although these were not statistically significant over time. In the Eastern Himalaya, positive trends of median velocity anomaly were observed in both lake-terminating (0.14 ± 0.03 m a–1 decade–1) (R 2 = 0.27, p = 0.27) and land-terminating (0.12 ± 0.03 m a–1 decade–1) (R 2 = 0.19, p = 0.45) glaciers, although similarly these were not statistically significant over time. Lake-terminating glaciers exhibited statistically significant positive trends of median velocity anomaly of 0.28 ± 0.02 m a–1 decade–1 (R 2 = 0.67, p = 0.002) and 0.30 ± 0.02 m a–1 decade–1 (R 2 = 0.54, p = 0.02) in the Central and Eastern Himalaya, respectively.
Median velocity anomaly, mean velocity anomaly and its interquartile range for lake-terminating and land-terminating glaciers in the Western Himalaya (a) and (b) (n = 40, n = 38), Central Himalaya (c) and (d) (n = 44, n = 45) and Eastern Himalaya (e) and (f) (n = 92, n = 93) between 2000 and 2018. Velocity anomaly is the difference between the annual velocity of an individual glacier and the mean velocity of the total glacier sample (n = 352).
3.3 Lake area change
Ice-contact lakes in the Western and Central Himalaya had a median surface area of 0.21 ± 0.07 and 0.27 ± 0.12 km2, respectively, whereas those in the Eastern Himalaya had a larger median surface area of 0.37 ± 0.19 km2 between 2015 and 2018 (Fig. 6). Lake area change during 2015–2018 relative to the earliest observation of each individual lake was comparable in Western, Central and Eastern Himalaya, although reducing with longitude, measuring 38 ± 39, 30 ± 43 and 27 ± 53% (median values), respectively. Furthermore, lake-terminating glaciers in the Central and Eastern Himalaya hosted a larger proportion of lakes with earliest observations that occurred before 2000, with 82 and 85%, respectively. In comparison, Western Himalaya lake-terminating glaciers hosted a smaller proportion of 70%, implying lake age of the sample increases moving eastwards across the range.
Distribution of ice-contact lake area (in 2018) (upper histogram), distribution of normalised lake area change (2018 relative to earliest observation) (lower histogram) and the number of ice-contact lakes whose earliest observation that fall within that year (pie chart) for (a) Western Himalaya (n = 38), (b) Central Himalaya (n = 44) and (c) Eastern Himalaya (n = 91).
4. Discussion
4.1 Limitations of our analyses
Although our analyses draw on the latest available datasets, there are some inherent assumptions and uncertainties associated with our approach. We sampled an approximately equal number of lake-terminating glaciers and land-terminating glaciers, where land-terminating glaciers were selected based on their proximity and similarities in attributes to the corresponding lake-terminating glacier (e.g. surface area, aspect). However, the relative proportions of our sample population are not representative of the relative proportions of these glacier types across the Himalaya. Therefore, these results should only be scaled for regional analysis after considering the relative proportions of each glacier type within the total population; for example, lake-terminating glaciers equate to about 50% of our sample, whereas they equate to only about 5% of the entire population (Lee and others, Reference Lee2021). There are uncertainties associated with the datasets used in this study; Dehecq and others (Reference Dehecq2019) and Hugonnet and others (Reference Hugonnet2021) quantified the uncertainties associated with glaciers velocity (median uncertainty of 2.0 m a–1 for the Central and Eastern Himalaya) and surface elevation change (median uncertainty of 0.5 m a−1 for glaciers with areas >1 km2).
4.2 Influence of ice-contact lakes and supraglacial debris on glacier change
The spatially heterogeneous impact of ice-contact lakes on glacier change is evident from the differences in glacier-wide surface elevation change of lake-terminating and land-terminating glaciers, which becomes more pronounced from the west to the east. Specifically, the influence of ice-contact lakes on glacier-wide surface elevation change is not evident in the Western and Central Himalaya, whereas lake-terminating glaciers in the Eastern Himalaya exhibited more negative glacier-wide surface elevation change rates than their land-terminating counterparts (Fig. 2). These sub-regional differences become more pronounced upon analysis within normalised elevation bins, which reveal that the influence of ice-contact lakes on glacier surface elevation change extends to the lowermost 30% of the normalised glacier elevation in the Eastern Himalaya. In contrast, no statistical differences were evident at any elevation between lake-terminating and land-terminating glaciers in the Central Himalaya, and only in the lowest 10% of normalised glacier elevation in the Western Himalaya (Fig. 3).
In all sub-regions, we find in line with previous studies that the presence of extensive supraglacial debris has little influence on glacier-wide surface elevation change (Kääb and others, Reference Kääb, Berthier, Nuth, Gardelle and Arnaud2012; Brun and others, Reference Brun2019). However, our analysis across normalised glacier elevations highlights the insulating effect of supraglacial debris whereby rates of surface elevation change were less negative in the lowest 10% of normalised glacier elevation for land-terminating debris-covered glaciers (cf. Rowan and others, Reference Rowan, Egholm, Quincey and Glasser2015). A similar effect was not observed for lake-terminating debris-covered glaciers where the presence of an ice-contact lake is the primary control on surface elevation change close to the terminus. While this result demonstrates that supraglacial debris is a secondary control on glacier surface elevation change in comparison with the presence of an ice-contact lake, velocity anomalies in the Eastern Himalaya suggest that supraglacial debris still exerts an influence on glacier behaviour, as lake-terminating clean-ice glaciers showed the greatest increases in velocity anomaly in both Central and Eastern Himalaya (Figs S5, S6). Pronk and others (Reference Pronk, Bolch, King, Wouters and Benn2021) reported that lake-terminating clean-ice glaciers in these regions exhibited the fastest ice velocities of the groups observed over a 2-year period (2017–2019). Our analyses show a similar result, in that velocity anomaly for lake-terminating clean-ice glaciers in Central and Eastern Himalaya showed the greatest increase over the period 2000–2018 (Figs S4–S6), which would then explain the high values observed at the end of our study period by Pronk and others (Reference Pronk, Bolch, King, Wouters and Benn2021).
4.2.1 Influence of ice-contact lakes on glacier mass balance
The lack of any statistically robust differences in glacier-wide mass balance between glacier terminus types in Western and Central Himalaya is in contrast to findings from previous studies (e.g. Brun and others, Reference Brun, Berthier, Wagnon, Kääb and Treichler2017; King and others, Reference King, Bhattacharya, Bhambri and Bolch2019). This is largely driven by more negative mass loss rates for land-terminating glaciers reported in this study compared to previous work. This can be partly attributed to differences in methodology whereby King and others (Reference King, Dehecq, Quincey and Carrivick2018) used discrete epochs for DEM differencing rather than a stacked-DEM approach, and differences in reporting period whereby Brun and others (Reference Brun, Berthier, Wagnon, Kääb and Treichler2017) used a time series up to the end of 2016 whereas we include more recent observations. However, given that the differences we have identified are subtle, and that the anomalous behaviour of glaciers within the western part of the range appears to have ceased in the most recent years, we are confident that the longer time series used by Hugonnet and others (Reference Hugonnet2021) provides the best available dataset for a contemporary analysis of these effects. Despite these differences, our absolute values of lake-terminating glacier mass balance are comparable to those previously published. We find that lake-terminating glaciers exhibit median mass balances between 2000 and 2019 of –0.43 ± 0.08 m w.e. a–1 (n = 38) in the Western Himalaya, −0.63 ± 0.06 m w.e. a–1 (n = 44) in the Central Himalaya and 0.57 ± 0.05 m w.e. a–1 (n = 91) in the Eastern Himalaya. King and others (Reference King, Dehecq, Quincey and Carrivick2018) reported similar values for lake-terminating glaciers in their Western (−0.49 ± 0.08 m w.e. a−1), Central (−0.67 ± 0.10 m w.e. a−1) and Eastern (−0.59 ± 0.12 m w.e. a−1) Himalayan regions.
4.2.2 Influence of ice-contact lakes on glacier velocity
Our analysis of glacier velocity demonstrates the spatially heterogeneous influence of ice-contact lakes on both a sub-regional and glacier scales. A clear pattern is evident as lake-terminating glaciers in the Eastern Himalaya show greater mean increases in ice velocity between 2000 and 2018, in addition to the influence of ice-contact lakes propagating further up-glacier in this region than in the Western Himalaya. While our focus is on velocity change over time, the median ice velocities (2000–2018) presented here are generally lower than those reported in previous studies (e.g. Dehecq and others, Reference Dehecq2019; Pronk and others, Reference Pronk, Bolch, King, Wouters and Benn2021). This study and Dehecq and others (Reference Dehecq2019) both used the ITS LIVE dataset (Gardner and others, Reference Gardner, Fahnestock and Scambos2019). Dehecq and others (Reference Dehecq2019) aggregated velocity pixels to quantify velocity anomaly regionally while this study used an individual glacier approach to resolve sub-regional differences between lake-terminating and land-terminating glaciers. The differences in median ice velocities are primarily attributed to differing data processing techniques, as Dehecq and others (Reference Dehecq2019) removed all pixels with values lower than 5 m a−1 and with error values >5 m a−1 while we filtered these data only based on the associated error grids by removing pixels with error values >5 m a−1 (section 2.4). Analysis of median velocity anomaly of all glaciers in this study (n = 352) agrees well with the negative velocity anomalies reported by Dehecq and others (Reference Dehecq2019) for their corresponding regions (West Nepal, East Nepal and Bhutan).
4.3 Temporal evolution of lake-terminating glaciers across the Himalaya
Our sub-regional analysis of glacier surface elevation change and glacier velocity, particularly across normalised glacier elevation, yields new evidence of the mechanisms by which ice-contact lakes can impact glacier behaviour (King and others, Reference King, Quincey, Carrivick and Rowan2017, Reference King, Bhattacharya, Bhambri and Bolch2019; Tsutaki and others, Reference Tsutaki2019; Pronk and others, Reference Pronk, Bolch, King, Wouters and Benn2021; Carrivick and others, Reference Carrivick2022a; Sato and others, Reference Sato, Fujita, Inoue and Sakai2022). We propose that the differences in ice-contact lake-induced surface elevation change and velocity anomaly observed between 2000 and 2019 are indicative of the typical evolutionary stages of ice-contact lakes in these regions. Based on the findings of this study and sub-regional differences in climatic conditions, glacier mass loss and the number and size of glacial lakes reported in previous studies, we conclude that the evolutionary stage of ice-contact lakes in terms of their coupling with glacier dynamics are more advanced in the Eastern Himalaya than the Central and Western Himalaya, and this is reflected in the differing behaviour of lake-terminating glaciers in these regions.
4.3.1 Lake-terminating glacier evolution in the Western and Central Himalaya
Analysis across normalised elevation bins on glaciers in the Western and Central Himalaya revealed subtle, and often statistically non-significant, differences between glacier termini types in both surface elevation change and velocity anomaly. Despite the lack of statistical differences, surface elevation change rates were more negative in the lowest 10% of the normalised glacier elevation, suggesting that surface elevation change at lower elevations has switched from mostly by surface ablation to mostly by frontal ablation (Truffer and Motyka, Reference Truffer and Motyka2016). These findings, coupled with relatively low incoming ice flux at lake-terminating glacier termini (2.3 m a–1 for Western Himalaya and 2.2 m a–1 for Central Himalaya), suggest that subaqueous melt as opposed to glacier calving is the dominant component of frontal ablation for glaciers with ice-contact lakes in these regions (Röhl, Reference Röhl2006; Truffer and Motyka, Reference Truffer and Motyka2016; Carrivick and others, Reference Carrivick, Tweed, Sutherland and Mallalieu2020). This is often the case for ice-contact lakes at earlier evolutionary stages, where they have not reached sufficient fetch (>80 m; Sakai and others, Reference Sakai, Nishimura, Kadota and Takeuchi2009) or become sufficiently deep relative to the ice thickness of their host glaciers to establish a hydraulic connection with the subglacial drainage system (Carrivick and others, Reference Carrivick, Tweed, Sutherland and Mallalieu2020). This behaviour is exemplified in the Western and Central Himalaya where surface elevation change rates between lake-terminating and land-terminating glaciers do not differ significantly above the lowest 10% of the glacier elevation.
4.3.2 Lake-terminating glacier evolution in the Eastern Himalaya
Lake-terminating glaciers in the Eastern Himalaya display much clearer evidence of lake-induced changes within the lowermost 30% of normalised elevation for glacier surface elevation change and the lowermost 80% of normalised elevation for velocity. This behaviour results from the development of a positive feedback whereby decreasing ice thickness and increased longitudinal strain promote deeper crevasses, which trigger calving events and greater terminus recession (Benn and others, Reference Benn, Warren and Mottram2007; Sugiyama and others, Reference Sugiyama2011; King and others, Reference King, Dehecq, Quincey and Carrivick2018). The pronounced influence of ice-contact lakes on glaciers in this sub-region indicates that the ice-contact lakes have likely reached a sufficient size in terms of their depth or ice-contact lake boundary length to induce a dynamic response in their host glaciers (Carrivick and others, Reference Carrivick2022a, Reference Carrivick2022b), and are therefore at a more mature evolutionary stage than those in Central and Western Himalaya.
Ice-contact lakes in the Eastern Himalaya generally formed earlier than those in the Western Himalaya, evidenced by a greater proportion of ice-contact lakes that formed prior to 2000 in the Eastern Himalaya (Fig. 6). This pattern is replicated in the larger population; Shugar and others (Reference Shugar2020) showed that the areal extent of glacial lakes in the Eastern Himalaya (55.1 km2) was 110% greater than in the Western Himalaya (26.2 km2) in 1990–1999. Differing sensitivity of Western Himalaya sub-continental glaciers and Eastern Himalayan monsoon-influenced glaciers coupled with rising air temperatures across the Himalaya and non-uniform changes in precipitation patterns has driven spatially heterogeneous glacier mass loss rates (Yao and others, Reference Yao2012; Sakai and Fujita, Reference Sakai and Fujita2017; Wang and others, Reference Wang, Liu, Shangguan, Radić and Zhang2019; Farinotti and others, Reference Farinotti, Immerzeel, de Kok, Quincey and Dehecq2020; Rounce and others, Reference Rounce, Hock and Shean2020; Nie and others, Reference Nie2021). The current rapid rates of glacier mass loss across the Himalaya are a key factor in ice-contact lake formation and expansion, where glacier recession and ice surface lowering promotes the exposure of glacially carved bedrock overdeepenings that fill with meltwater and supraglacial pond formation and coalescence (Carrivick and Tweed, Reference Carrivick and Tweed2013; Carrivick and others, Reference Carrivick, Tweed, Sutherland and Mallalieu2020; King and others, Reference King, Turner, Quincey and Carrivick2020). This glacier terminus environment evolution is exemplified in the Eastern Himalaya where the greatest increases in both number and areal extent of glacial lakes are coincident with highest rates of glacier mass loss between 2000 and 2020 (Shean and others, Reference Shean2020; Shugar and others, Reference Shugar2020).
There is strong evidence for the stage of lake development being an important control of glacier behaviour, and explaining the observed differences in lake-terminating glacier behaviour across the Himalaya. This is a subtle yet important advance in understanding from merely considering glacier-wide surface elevation change or terminus-wide velocity as indicative of how glaciers are changing. We show that ice-contact lakes in the Western and Central Himalaya have generally formed more recently than those in the Eastern Himalaya (Fig. 6) and that these lakes are less likely to be of sufficient depth or ice-contact boundary length to exert a significant control on near-terminus glacier dynamics (Carrivick and others, Reference Carrivick2022a, Reference Carrivick2022b). However, we do not have sufficient data relating to exact dates of lake formation, evolving lake volumes and ice thickness to establish a precise relationship between lake evolutionary stage and the expected glacier response, and it is clear from previous work that such a relationship would be complex and likely glacier specific. For example, King and others (Reference King, Dehecq, Quincey and Carrivick2018) observed four Central Himalayan glaciers with substantial ice-contact lakes (0.48–1.38 km2) that showed no evidence of dynamic influence on their host glaciers, perhaps indicating that the lake effect begins to diminish beyond a threshold, for example, as the glacier recedes out of an overdeepening.
4.4 Impact of ice-contact lakes on future glacier change in the Himalaya
The distribution of ice-contact lakes will determine how lake-terminating glaciers may behave in the future at a regional scale. While it is difficult to project future lake development through supraglacial pond formation and coalescence, prediction of the locations and extent of exposed subglacial overdeepenings due to glacier recession is possible in alpine settings (Linsbauer and others, Reference Linsbauer2016; Farinotti and others, Reference Farinotti, Round, Huss, Compagno and Zekollari2019; Carrivick and others, Reference Carrivick2022b). In the Western Himalaya, the number of ice-contact lakes may increase by 65% more than in the Eastern Himalaya until 2100 (Furian and others, Reference Furian, Maussion and Schneider2022), suggesting that the influence of glacial lakes on ice recession may become more pronounced within the Western Himalaya and less pronounced in the Eastern Himalaya over time. However, the increase in the number of lake-terminating glaciers may be matched with many existing glacier–lake systems reaching the latter stages of lake development. The influence of the ice-contact lakes may then continue, however it is more likely that the importance of most present-day ice-contact lakes for glacier change will decrease as their host glaciers recede out of the overdeepening that the lake currently occupies (Truffer and Motyka, Reference Truffer and Motyka2016). Therefore, it is difficult to use the findings of studies such as those presented here to extrapolate into the future, which further supports the need for a process-based, time-evolving understanding of glacier–lake interaction. Empirical observations of the impact of ice-contact lakes are thus critical to be able to project glacier response to future climate scenarios. Simplified process representation can overcome challenges of computational expense, particularly when simulations spin up from, and extend to, long timescales (e.g. Sutherland and others, Reference Sutherland2020).
5. Conclusions
Analysis of freely available datasets of glacier surface elevation change and velocity across the Himalaya between 2000 and 2019 revealed a spatially heterogeneous influence of ice-contact lakes on glacier behaviour related to the stage of lake development. Lake-terminating glaciers exhibit statistically greater surface elevation change rates than land-terminating glaciers in the Eastern Himalaya (a difference of 0.14 m a−1) but are similar in the Western and Central Himalaya. We conclude that ice-contact lakes affect the temporal evolution of glacier surface elevation change and velocity only when they have reached an advanced evolutionary stage, as seen in the Eastern Himalaya, which as a region has experienced greater mass loss rates to date in comparison to that in the Western and Central Himalaya driving the formation and expansion of a greater number of glacial lakes. The distance to which the influence of ice-contact lakes propagates up-glacier from the glacier–lake interface in cases of advanced lake evolution is important in the lowermost 30% of the glacier elevation. Therefore, projections of glacier evolution need to incorporate the spatial variability in surface elevation change within glaciers that terminate in lakes.
Our results demonstrate that the response of a glacier to the presence of an ice-contact lake evolves over decadal time scales. However, process-based understanding of glacier–lake interaction in the Himalaya is still limited by a lack of in situ observations. In particular, understanding is limited of the physical characteristics of ice-contact lakes (e.g. their thermal regime and bathymetry), their influence on frontal ablation and how to best represent these processes in broad-scale projections. These issues should be prioritised if the influence of ice-contact lakes is represented in projections of glacier change in this region and more widely in other glacierised mountain ranges.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2024.9.
Data
This study used open access datasets and tools; velocity data for 2000–2018 were generated using auto-RIFT (Gardner and others, Reference Gardner2018) and provided by the NASA MEaSUREs ITS_LIVE project (Gardner and others, Reference Gardner, Fahnestock and Scambos2019). Mean surface elevation change data for 2000–2019 were from Hugonnet and others (Reference Hugonnet2021). Ice-contact lake area data for 1990–2018 were from Shugar and others (Reference Shugar2020) and Wang and others (Reference Wang2020).
Acknowledgements
A. C. S. is a recipient of a NERC Panorama PhD studentship (NE/S007458/1). A. V. R. was supported by a Royal Society Dorothy Hodgkin Research Fellowship (DHF\R1\201113). We would like to express our thanks to the editors, Fanny Brun and Hester Jiskoot, and two anonymous reviewers for their constructive comments and suggestions which ultimately led to an improved manuscript. The University of Leeds is acknowledged for funding the article publication.
Open access
