Contrasting surface velocities between lake- and land-terminating glaciers in the Himalayan region

Meltwater from Himalayan glaciers sustains the flow of rivers such as the Ganges and Brahmaputra on which over half a billion people depend for day-to-day needs. Upstream areas are likely to be affected substantially by climate change, and changes in the magnitude and timing of meltwater supply are expected to occur in coming decades. About 10 % of the Himalayan glacier population terminates into proglacial lakes, and such lake-terminating glaciers are known to exhibit higher-than-average total mass losses. However, relatively little is known about the mechanisms driving exacerbated ice loss from lake-terminating glaciers in the Himalaya. Here we examine a composite (2017–2019) glacier surface velocity dataset, derived from Sentinel 2 imagery, covering central and eastern Himalayan glaciers larger than 3 km². We find that centre flow line velocities of lake-terminating glaciers (N = 70; u_median: 18.83 m yr⁻¹; IQR – interquartile range – uncertainty estimate: 18.55–19.06 m yr⁻¹) are on average more than double those of land-terminating glaciers (N = 249; u_median: 8.24 m yr⁻¹; IQR uncertainty estimate: 8.17–8.35 m yr⁻¹) and show substantially more heterogeneity than land-terminating glaciers around glacier termini. We attribute this large heterogeneity to the varying influence of lakes on glacier dynamics, resulting in differential rates of dynamic thinning, which causes about half of the lake-terminating glacier population to accelerate towards the glacier termini. Numerical ice-flow model experiments show that changes in the force balance at the glacier termini are likely to play a key role in accelerating the glacier flow at the front, with variations in basal friction only being of modest importance. The expansion of current glacial lakes and the formation of new meltwater bodies will influence the dynamics of an increasing number of Himalayan glaciers in the future, and these factors should be carefully considered in regional projections.

1 Introduction

Himalayan glaciers provide an important baseline supply of meltwater for downstream areas (Bolch, 2017; Immerzeel et al., 2010; Pritchard, 2019; Viviroli et al., 2007). A large decrease in runoff from the rivers that drain this mountain range will have major implications for downstream water security, particularly in the populous catchments of the Ganges, Indus and Brahmaputra rivers. Although a drastic reduction in glacier area and mass is projected in the Himalaya over the 21st century (Kraaijenbrink et al., 2017; Rounce et al., 2020), large uncertainties in the pace of the loss exist (Lutz et al., 2013). Hence, there are also large uncertainties in future meltwater supply, and an improved understanding of the evolution of Himalayan glaciers is needed.

Himalayan glaciers have been retreating and losing mass since the mid-19th century, and rates of mass loss have been increasing over at least the last 4 decades (Azam et al., 2018; Bhattacharya et al., 2021; Bolch et al., 2012, 2019; King et al., 2019; Maurer et al., 2019). Various studies report Himalayan-averaged glacier mass losses of around −0.40 ± 0.10 m w.e. yr⁻¹ since the beginning of this century (Brun et al., 2017; Shean et al., 2020), which roughly translates into a total mass loss of 7.5 Gt yr⁻¹. However, within the Himalayan mountains, large intra-regional variability in glacier mass loss exists (Azam et al., 2018; Bolch et al., 2012; Brun et al., 2017; King et al., 2019; Maurer et al., 2019), which indicates that there are factors capable of exacerbating – or reducing – glacial mass losses that are at least partially decoupled from climate.

The development of proglacial lakes in direct contact with the glacier terminus has been linked with enhanced glacier mass loss in the Himalayan region (King et al., 2019; Maurer et al., 2019). This contrast in mass loss with land-terminating glaciers manifests itself in two ways, namely by elevated terminal retreat rates and by enhanced surface lowering towards the terminus of the lake-terminating glaciers (King et al., 2019). The latter indicates that proglacial lakes can influence the flow characteristics of their host glacier through dynamic thinning.

A factor that further complicates the dynamics and mass loss rate of lake-terminating glaciers in the Himalaya is the presence of a thick layer of debris, which is widespread on Himalayan glaciers (Herreid and Pellicciotti, 2020). The low-gradient, debris-covered portions of many Himalayan glaciers are preconditioned for meltwater ponding and eventually proglacial-lake development, which often results from a deepening and coalescence of supraglacial lakes (Benn et al., 2012; Quincey et al., 2007) which are bounded by a stagnant, ice-cored moraine dam. The combination of the morphology, insulating characteristics of debris and lake development may cause a response to climate forcing that is strongly non-linear (Benn et al., 2012), though only little is known about how such a transition develops.

Two key factors can be identified which make lake-terminating glaciers distinctively different from their land-terminating counterparts, namely the stresses at the bed and the terminus of the glacier. Firstly, a body of water exerts a buoyancy force on the host glacier, reducing the effective pressure and consequently reducing the basal resistance, which ultimately can result in faster glacier flow (Benn et al., 2007b). Secondly, dynamical changes can result from a force imbalance that act at the terminus and trigger a retreat and reduce along-flow resistive stresses (Nick et al., 2009). This can be especially important in rapidly evolving environments (Benn et al., 2007b), such as the Himalayan region where the number and area of proglacial lakes have increased substantially over a few decades. In alpine settings, the transition from a land-terminating glacier to a lake-terminating glacier could therefore change the dynamic regime of the glacier, and such a transition might be partially decoupled from climate (Benn et al., 2012).

Several recent remote-sensing-based studies on glacier surface velocity indicated the divergent evolution of the dynamics of lake- and land-terminating glaciers. Dehecq et al. (2019b) documented widespread land-terminating glacier slowdown since the start of the 21st century across High Mountain Asia (HMA) in response to diminished driving stress caused by long-term ice thinning. In contrast, more localized studies have shown several examples of lake-terminating glacier flow acceleration over a similar time period (King et al., 2018; Liu et al., 2020; Song et al., 2017; Tsutaki et al., 2019). The number and total area of proglacial lakes in the Himalayan region has increased (Nie et al., 2017; Shugar et al., 2020; Zhang et al., 2015, 2019), a trend which is likely to continue in the near future, as many glacier beds are characterized by overdeepenings (Linsbauer et al., 2016). Therefore, a robust understanding of the behaviour of lake-terminating glaciers is crucial. However, a spatially comprehensive assessment of the contrasting dynamics between land-terminating and lake-terminating glaciers has yet to be undertaken.

Numerical ice-flow models have been utilized to investigate the dynamic thinning of marine-terminating outlet glaciers (Benn et al., 2007a; Enderlin et al., 2013; Nick and Oerlemans, 2006; Vieli et al., 2001; Vieli and Nick, 2011) and more recently of lake-terminating glaciers in alpine regions (Sutherland et al., 2020; Tsutaki et al., 2019). Tsutaki et al. (2019) employed a diagnostic 2D model setup to show that the transition from a land- to a lake-terminating boundary condition will significantly increase the surface flow velocities near the calving front. Sutherland et al. (2020) used a numerical transient model setup to show that a proglacial lake was a dominant control on the ice velocity during times of glacial-lake growth after the Last Glacial Maximum in New Zealand. Ice thickness data and a quantification of the hydrological characteristics of a glacier's proglacial lake are important components of such model setups (Carrivick et al., 2020), but these data are very limited in the Himalaya, making a realistic model setup problematic.

The main aim of this study is to examine the influence of proglacial lakes on Himalayan glacier dynamics, in order to improve the current understanding of the large subregional heterogeneity of glacier behaviour. More specifically, we seek to investigate the attribution of changes in the velocity field to dynamic thinning and investigate the role that debris cover plays on glacier–lake dynamics. To achieve this, we used Sentinel-2 satellite imagery to derive a large-scale contemporary Himalayan glacier velocity dataset at an improved resolution compared to studies to date. We compare the velocity dataset against surface elevation change data (King et al., 2019) to discuss the role of proglacial lakes and debris cover on glacier dynamics. Finally, we employ a numerical flow model to investigate the factors that dominate the lake-terminating glacier dynamics and explore their potential to accelerate current and future mass losses.

2 Study area

Our study area covers five subregions within the central and eastern (CE) Himalaya (Fig. 1) in which glacial lakes have been widespread over an extended time period. Glaciers in the CE Himalaya cover an area of ∼ 13 900 km², which is about 60 % of the glacierized area of the total Himalayan arc (Bolch et al., 2012). The location of the Himalaya is around the southern rim of the Tibetan Plateau (TP), with the CE Himalaya being the source of two major trans-boundary rivers, namely the Ganges and the Brahmaputra. The extreme Himalayan topography exerts a strong influence on north–south contrasting precipitation patterns by forming an orographic barrier and depleting the monsoonal air of its moisture on southern windward slopes, resulting in relatively dry slopes on the TP (Ageta and Higuchi, 1984).

Figure 1 Map showing the regional subdivisions (red rectangles) with the associated glacier characteristics, including terminus type and surface cover as a fraction of the total subregional glacierized area. The excluded fraction represents the glacierized area from glaciers smaller than 3 km². The names of the subregions are in accordance with King et al. (2019). Country boundaries are tentative and for orientation only.

Related to this stark contrast in north–south relief is the distribution of clean-ice and debris-covered glaciers (Scherler et al., 2011). Glaciers in low-relief areas, sloping northwards and facing the TP, generally show little or no debris cover and have extensive accumulation areas. In contrast, glaciers surrounded by much steeper topography on the southern side of the main orographic divide receive a large proportion of their accumulation by snow avalanching from steep hillslopes. Steep hillslopes supply large fluxes of rocky material to the glacier, and, as a result, glaciers in such settings often have an extensive debris cover, which can range from a few centimetres to several metres in thickness (Scherler et al., 2011).

In this study we focus on glaciers with an area larger than 3 km², compared to the 5 km² previously utilized by Dehecq et al. (2019b), which enables us to add substantially more glaciers to our dataset (Fig. 6b). Glaciers smaller than 3 km² are often located at a high elevation and do not typically host a proglacial lake and thus fall beyond the scope of our study. Also, glacier volume scales exponentially with glacier surface area, which increases the representativeness of our study onto potential ice volume losses. For the glacier outlines and corresponding surface area we used the dataset from the Randolph Glacier Inventory (RGI 6.0) (RGI Consortium, 2017).

To study glacier–lake dynamics, we select five subregions within the CE Himalaya with a high density of proglacial lakes, according to the lake inventory of Zhang et al. (2015) (Fig. 1; Table 1). We classify glaciers as lake-terminating when the glacier shows a clear terminal ice cliff in direct contact with the lake and base our classification on the glacial-lake inventory of Wangchuk and Bolch (2020) and Zhang et al. (2015) using multiple sources of satellite imagery. The classification of debris cover is binary (debris covered or clean ice), and for this we follow the criteria defined by Brun et al. (2019), classifying glaciers as debris covered where more than 19 % of their area was mantled by debris.

Table 1 Regional distribution of debris-covered and clean-ice, lake- and land-terminating glaciers.

^* The percentage of glacierized area covered (coverage) over the whole region is relatively low, since it also incorporates all CE Himalayan glaciers outside the five subregions.

3 Methods

3.1 Surface ice velocity

3.1.1 Sentinel-2 satellite imagery

The European Copernicus Sentinel-2 series consists of two satellites: Sentinel-2a, launched in June 2015, and Sentinel-2b, launched in March 2017. The two satellites have a combined revisit time of 5 d, and their orthorectified image products (Level-1C) are freely available at https://earthexplorer.usgs.gov/ (last access: 5 June 2020) and https://scihub.copernicus.eu/ (last access: 24 April 2020). In this study the 10 m near-infrared band (band 8) was used to exploit the contrasting spectral properties of fresh snow, firn and clean ice at this wavelength, an approach which has proven to work well for feature tracking on a variety of glacier surfaces (Kääb et al., 2016).

Throughout most of its mission, the multi-temporal co-registration accuracies of Sentinel-2 products from the same orbit have been below 12 m (95.5 % confidence interval), which is reported monthly by the European Space Agency (ESA) (Clerc et al., 2020). When co-registering two Level-1C images from the same relative orbit (repeat orbit), DEM (digital elevation model) effects will be present but have the same pattern in both datasets (because they have a similar off-nadir cross-track look angle) so that they can be eliminated by calculating the average offset field obtained from correlating the two images. Products acquired from neighbouring overlapping swaths translate into additional offsets of up to 5.9 m (Kääb et al., 2016) and have therefore been omitted in this study.

3.1.2 Image pair selection

We selected multiple image pairs of the same locality to maximize the velocity field coverage, which is often limited by shadows, cloud cover, low visual contrast or sensor saturation. The selection of multiple image pairs over the same location also increases the overall confidence in the velocity estimates. The final velocity field is then an average of all the valid velocity estimates, a strategy explored by several studies (i.e. Dehecq et al., 2015; Scherler et al., 2011; Willis et al., 2012). The maximum number of image pairs separated by 1 year was selected for the month of November. We limit the image selection to this month because (1) from April until October, the optical satellite imagery is largely obscured by monsoonal cloud cover and (2) from December until April, the glacier surface contrast is generally low due to low-altitude, westerly induced snow cover, causing the image-matching algorithm to perform poorly. Also, after a certain number of image pairs, the reduction of the residual error in the averaged velocity field is only marginal (Dehecq et al., 2015). This approach produced a dataset of 149 images and 427 image pairs with a central date of 24 August 2018 (Table 2). Note that due to the lower Sentinel-2 repeat cycle in 2016, this central date is centred towards the end of the November 2016–November 2019 interval.

3.1.3 Image processing

To reduce computational costs, a mask was applied over all non-glacierized areas and glaciers with an area below 3 km². For this mask we used the glacier outlines from the Randolph Glacier Inventory (RGI 6.0) (RGI Consortium, 2017). Heid and Kääb (2012) evaluated several feature-tracking methods and showed that orientation correlation performed best under most circumstances, which we employ in this study. This method, developed by Fitch et al. (2002), creates two complex orientation images from the original image pairs. Each pixel represents the orientation of intensity gradient in the x and y direction at that pixel, making the method invariant to illumination change, which is a desired property for feature-tracking algorithms.

Table 2 Number of November Sentinel-2 images and number of image pairs. The date range of all tiles are from November 2016 to November 2019.

From the two orientation images a search window (i.e. a squared collection of pixels) and a reference window centred around the same location (x, y) were extracted and matched using correlation computed in the frequency domain with fast Fourier transforms (FFTs), according to the convolution theorem (McClellan et al., 1999).

We matched the orientation of the intensity gradient that is contained in the phase of the orientation image. After this initial estimate we then refined the maximum estimation by upsampling the product of the two orientation images in the frequency domain in a small neighbourhood of the initial maximum (Guizar-Sicairos et al., 2008). This matching process was repeated over all the glacierized areas with steps equalling half of the search window size, for all image pairs. This leaves us with N pairs of velocity data matrices (x and y displacement) with a resolution of half times the search window size for each given satellite scene (Table 3).

Table 3 Parameters used in image matching using Sentinel-2 10 m pixel resolution imagery.

To remove matching blunders present in the derived velocity fields, we largely adopted a strategy proposed by Gardner et al. (2020): we used a disparity filter that consists of two components. First, the filter checks for the uniqueness of each x and y displacement by comparing each element with their surrounding neighbours that are co-located in a 5 × 5 kernel. If less than 9 of the 25 co-located are within a 25 % range of the search limit of the algorithm, the x- and y-displacement pair is disregarded. Secondly, an x- and y-displacement pair is also considered to be a blunder when one of the two elements deviates more than 3 times the interquartile range from the median of all co-located pixels.

To calculate the median offset in the x and y direction, we selected for each 100 km × 100 km tile a large stable area of which we can reasonably assume the displacement to be zero. We therefore avoided extremely high alpine terrain that might be abundant with periglacial features such as rock glaciers or solifluction lobes, as we do not expect these surfaces to remain motionless through time. However, topography that characterizes the general environment of mountain glaciers is required, and we therefore also avoided the selection of flat terrain. The stable area of each image tile is of sufficient width such that it covered multiple granules (±25 km in width). To reduce the noise in the velocity data, we subtracted this median offset from the whole (glacierized and stable) x- and y-displacement field. The final 2017–2019 x- and y-displacement field was then created by taking the median of all the image pairs for both velocity components. Finally, we used these displacement results in these stable areas to evaluate the precision and uncertainty of the feature-tracking algorithm.

3.2 Uncertainty of the velocity field

Uncertainties of the median velocity field are dominated by the precision of the feature-tracking algorithm, the co-registration error, the temporal variability of glacier flow and the number of velocity estimates. For the estimation of the 95 % confidence interval (CI₉₅) of each median velocity component, we adopted the methodology of Dehecq et al. (2015) and expect the CI₉₅ to conform to

$$\begin{array}{}\text{(1)}& {\mathrm{CI}}_{\mathrm{95}}=\mathit{\kappa }{\displaystyle \frac{{\mathrm{MAD}}_{\mathrm{disp}}}{N^{\mathit{\alpha }}}},\end{array}$$

where MAD_disp (mean absolute deviation) is the dispersion at each velocity location of the N number of estimates:

$$\begin{array}{}\text{(2)}& {\mathrm{MAD}}_{\mathrm{disp}}\left(i,j\right)=\mathrm{1.483}\times {\mathrm{median}}_{t\in T}\left\{\left|V\left(i,j,t\right)-\stackrel{\mathrm{‾}}{V}\left(i,j\right)\right|\right\},\end{array}$$

where T is the collection of N velocity estimates $V\left(i,j,t\right)$ merged to obtain the median velocity $\stackrel{\mathrm{‾}}{V}\left(i,j\right)$ at pixel (i,j). The number 1.483 is a scale factor that relates the MAD to the standard deviation. Parameters κ and α determine the width and the thickness of the tail of the distribution and have yet to be estimated. Equation (2) leaves us with three unknowns but can be solved at stable areas where CI₉₅ can be determined as a function of MAD_disp and N for each tile location with its corresponding stable area, providing an uncertainty estimation for areas with actively flowing ice.

3.3 Glacier group uncertainty

The glacier group uncertainty (i.e. the uncertainty in the median velocity of a collection of glaciers along the normalized flow line) depends on the uncertainty of the individual glacier velocity measurements (CI₉₅), the spread between the velocity points u_c among the sample group and the number of velocity points N_u. We estimated this uncertainty by applying a Monte Carlo simulation which draws 200 random points from the uncertainty distribution of each individual velocity point u_c in the region of interest, where u_c is the centreline velocity (Sect. 3.4). Then for each sample round, following the bootstrap method, we drew N_u samples with replacement to calculate the median and repeated this 500 times. This results in 10⁵ estimates of the median, from which we determined the interquartile range (IQR), and we used this as a primary estimator of our regional mean median velocity uncertainty.

3.4 Glacier centre flow line analysis

We analysed surface velocity, elevation change and slope (Sect. 3.5) along the main glacier centre flow line (hereafter, centreline), an approach adopted by several earlier studies (i.e. Liu et al., 2020; Nagler et al., 2015; Scherler et al., 2011). Centrelines were produced with the Open Global Glacier Model (OGGM; Maussion et al., 2019) using a slightly adapted algorithm from Kienholz et al. (2014) and glacier outlines from RGI 6.0 using the SRTM (Shuttle Radar Topography Mission) DEM. All centrelines were manually adapted using 2019 Sentinel-2 satellite data and velocity data from this study to ensure that the centrelines end at the 2019 terminus position and that they follow the main flow tributary.

To extract centreline velocity data, we conducted a nearest-neighbour sampling every 80 m and averaged the velocity estimate by using a 3 × 3 (240 m × 240 m) Gaussian window:

$$\begin{array}{}\text{(3)}& u_c={\sum }_{j,i=-\mathrm{1}}^{\mathrm{1}}{\displaystyle \frac{u_{i,j}}{{\left({\mathrm{CI}}_{{\mathrm{95}}_{i,j}}\right)}^{\mathrm{2}}}}{e}^{-\frac{\mathrm{1}}{\mathrm{2}}\frac{i^{\mathrm{2}}+j^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}},\end{array}$$

where u is the velocity estimate at pixel (i,j), ${\mathrm{CI}}_{\mathrm{95}}^{-\mathrm{2}}$ is the weighting factor and σ=0.7 is the standard deviation of the Gaussian window. This approach increases the overall confidence of our median velocity estimates. The Gaussian window also prevents pixels further away from the centre flow line skewing the averaged data, which may result in an underestimation of the velocity values.

Then, to compare the velocity profiles for multiple glaciers at the ablation zone, we normalized the glacier centrelines horizontally along their ablation zone length. To achieve this, we first selected all discrete centreline velocity data points starting at the equilibrium-line altitude (ELA) and upsampled all centrelines with an ablation area length below 4000 m and downsampled the rest. We took 4000 m as the most representative length, as the results showed that this approaches the overall median length of the ablation zone for the whole glacier population. Our choice to analyse the centreline data along the normalized length of the ablation zone provides information on the dynamic influence of terminus type and surface cover but limits our ability to evaluate the effect of climate on surface elevation change, such as done by King et al. (2019) and Maurer et al. (2019). This also restricts the possibility to quantitatively attribute contrasting surface elevation change rates of lake-terminating and land-termination glaciers to dynamic thinning, which is especially true for clean-ice glaciers whose thinning rates appear to be highly dependent on climate, hence elevation (Scherler et al., 2011). Therefore, this study will be restricted to a qualitative analysis when evaluating the velocity data in the context of surface elevation change rates.

3.5 Surface elevation change, estimation of ELA and surface slope

We examined ice-thinning rates using the surface elevation change (dh dt⁻¹) dataset from King et al. (2019). The elevation change field is a mean estimate derived from 499 DEMs generated from Worldview and GeoEye optical stereo pairs spanning the period 2012–2016, with a central date around mid-2015 (Shean, 2017), and the SRTM DEM from 2000. The central dates of our elevation change and surface velocity datasets are therefore separated by ∼3 years. Considering an average retreat rate of 26.8 ± 1.4 m yr⁻¹ for lake-terminating glaciers (King et al., 2019), the lowermost 100 m of lake-terminating glaciers are likely to be devoid of surface velocity data when the two datasets are compared. However, this could be considerably more for quickly retreating glaciers, which must be considered when interpreting our results.

In this study we focused our analyses on the ablation zone of glaciers, for which we need an estimation of the ELA. We followed the approach of Braithwaite and Raper (2009) and considered the median altitude of each glacier as a proxy for the ELA, defined by the elevation of the 50th percentile in glacier area. This estimate of the ELA is available within the Randolph Glacier Inventory (RGI 6.0) by the RGI Consortium (2017) (Pfeffer et al., 2014).

To calculate the slope of the ablation zone we used the Advanced Land Observing Satellite (ALOS) World 3D DEM (Tadono et al., 2014), which is available at a 30 m resolution and is based on DEMs generated from stereo image pairs collected over the period 2006–2011.

3.6 Numerical flow line model

3.6.1 Model description

Substantial variability in the behaviour of calving glaciers is common even when they are located in similar environmental and therefore climatic settings (Enderlin et al., 2013; Truffer and Motyka, 2016). Such variability is often attributed to factors such as glacier shape (Enderlin et al., 2013), the thermal regime of the proglacial water body (Truffer and Motyka, 2016) and glacier bed topography (Enderlin et al., 2013). Similar variability can be expected for lake-terminating glaciers in the Himalaya, where glacier morphology and glacial-lake geometry is diverse, but little work has been done to examine why the frontal dynamics of the lake-terminating glaciers may vary. In combination with our analyses of remotely sensed glacier surface velocities, we carried out a synthetic diagnostic numerical experiment to examine the response of terminus-proximal flow rates of lake-terminating glaciers to adjustments in glacier geometry (thickness, slope and width) and lake depth. We use a flow line model that was previously utilized on predominantly marine-terminating outlet glaciers (Enderlin et al., 2013; Nick et al., 2010, 2009; Vieli and Payne, 2005; Vieli and Nick, 2011). The depth integration of the model implicitly employs the shallow-shelf approximation (SSA), which is not fully appropriate for the entire model domain. However, the model results in the ablation zone, where the surface slope is generally low and where we assume sliding to dominate the glacier flow (Liu et al., 2020; Tsutaki et al., 2019), should still be adequate to provide us useful insight into the relevant physical processes in operation (Le Meur et al., 2004). The governing force-balance equation determined through conservation of momentum is (Nick et al., 2010)

$$\begin{array}{}\text{(4)}& \begin{array}{rl}\mathrm{2}{\displaystyle \frac{\partial }{\partial x}}\left(Hv{\displaystyle \frac{\partial U}{\partial x}}\right)& +A_{\mathrm{s}}{\left[\left(H-{\displaystyle \frac{{\mathit{\rho }}_w}{{\mathit{\rho }}_i}}D\right)U\right]}^{\frac{\mathrm{1}}{m}}\\ & +{\displaystyle \frac{H}{W}}{\left({\displaystyle \frac{\mathrm{5}U}{\mathrm{2}AW}}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}}={\mathit{\rho }}_{i}gH{\displaystyle \frac{\partial h}{\partial x}},\end{array}\end{array}$$

where U is the vertically averaged horizontal ice velocity; ρ_i= 917 kg m⁻³ is the density of ice; ρ_w= 1000 kg m⁻³ is the density of fresh water; H is the ice thickness; A_s is the sliding parameter; D is the ice thickness submerged under the lake level; W is the glacier width; g is gravitational acceleration; h is the ice surface elevation; m is the bed friction exponent; and v is the depth-averaged effective viscosity, which is defined as follows:

$$\begin{array}{}\text{(5)}& v=A^{-\mathrm{1}/\mathrm{3}}{\left|{\displaystyle \frac{\partial U}{\partial x}}\right|}^{-\frac{\mathrm{2}}{\mathrm{3}}}.\end{array}$$

The right-hand side of Eq. (5) is the gravitational driving stress, which is balanced by gradients in longitudinal stress (first term on the left-hand side), basal resistance (second term) and lateral resistance (third term). A is the temperature-dependent rate factor and increases from a minimum of $\mathrm{3.5}\times {\mathrm{10}}^{-\mathrm{25}}$ Pa⁻³ s⁻¹ at the divide to a maximum of 1.7 × 10–24 Pa⁻³ s⁻¹ at the calving front, corresponding to a depth-averaged ice temperature range of −10 to −2 ^∘C at the ablation zone (Cuffey and Paterson, 2010), for which we follow Enderlin et al. (2013).

The assumption is made that basal drag depends on sliding velocity and effective basal pressure non-linearly (Bindschadler, 1983; van der Veen and Whillans, 1996; Vieli and Payne, 2005), for which we choose m = 3. Resistance from drag along the lateral margins is estimated by integrating the force-balance equation over the width of the glacier assuming a constant ice thickness such that lateral drag supports the same fraction of driving stress along a transect across the glacier. Consequently, the model assumes a flow through a rectangular basin, with lateral support that is independent from effective pressure. The up-glacier boundary is the upper bound of the glacier (U = 0) and at the calving front the longitudinal stress is balanced with the difference in hydrostatic pressure between the ice and lake water, which results in the following depth-averaged stretching rate:

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\partial U}{\partial x}}=A{\left[{\displaystyle \frac{{\mathit{\rho }}_{i}g}{\mathrm{4}}}\left(H-{\displaystyle \frac{{\mathit{\rho }}_w}{{\mathit{\rho }}_i}}{\displaystyle \frac{D^{\mathrm{2}}}{H}}\right)\right]}^{\mathrm{3}}.\end{array}$$

The synthetic model domain extends 8000 m horizontally, 1000 m vertically and resembles the main characteristics of relatively large clean-ice, lake-terminating glacier of which numerous examples can be found flowing northwards onto the TP (Fig. A1 in the Appendix). We assumed a concave-up profile resulting in a slope of about 4.5^∘ within 2 km of the terminus, and our interpretation in the results will solely be focused on this part of the glacier. For glacier width (W), we used a value of 600 m at the terminus which increases to 2.5 km up-glacier to reduce the influence of lateral margins. We used a maximum ice thickness (H) of 230 m and an ice thickness of 120 m at the terminus (H_t), values in line with ice thickness estimates of the larger Himalayan glaciers (Farinotti et al., 2019). Subglacial water pressure is assumed to follow a piezometric surface rising up-glacier (Benn et al., 2007b), and we estimated the piezometric surface to be equivalent to 60 % of the ice thickness away from the calving front, accounting for a simplified basal hydrology. We then tuned A_s such that the maximum velocity near the ELA of the larger clean-ice, lake-terminating glaciers reaches a typical value of 50 m yr⁻¹ (Dehecq et al., 2019a; Gardner et al., 2020) and found a value of A_s=  2.5×10⁶ Pa m${}^{-\mathrm{2}/\mathrm{3}}$ s${}^{\mathrm{1}/\mathrm{3}}$, which is on the low side compared to values used at marine-terminating outlet glaciers of the Greenland ice sheet (i.e. Enderlin et al., 2013; Nick et al., 2009).

3.6.2 Experimental design

To examine the importance of the frontal boundary condition, we varied the height of the terminal ice cliff above buoyancy by lowering the lake level by

$$\begin{array}{}\text{(7)}& D={\displaystyle \frac{{\mathit{\rho }}_i}{{\mathit{\rho }}_w}}{H}_{\mathrm{t}}+\mathrm{\Delta }\mathrm{D},\end{array}$$

where D was the ice thickness submerged under the lake level (i.e. the lake depth at the terminus), H_t is the ice thickness at the terminus and ΔD is the lake surface level change. We examined the glacier dynamics when the calving front is exactly buoyant (i.e. ΔD=0 m) and for an increased ice cliff height resulting from a lake surface level change of −10 and −15 m, keeping the terminus position and ice thickness constant. The varying height of the terminal ice cliff above buoyancy are chosen to fall within a realistic range based on the limited observational evidence on terminal ice cliffs (Watson et al., 2020).

For each ice cliff height (i.e. lake surface level) configuration, we also ran the numerical model by keeping the basal friction independent from the effective pressure $(H-\frac{{\mathit{\rho }}_{\mathrm{w}}}{{\mathit{\rho }}_i}D=\mathrm{1}$), ruling out the effect of the lake on basal friction, and used a corresponding roughness parameter of $A_{\mathrm{s}}=\mathrm{9}\times {\mathrm{10}}^{\mathrm{6}}$ Pa m${}^{-\mathrm{1}/\mathrm{3}}$ s${}^{\mathrm{1}/\mathrm{3}}$. Note that lowering the lake surface is just one way of manipulating the frontal boundary condition and that an instantaneous glacier retreat could be an alternative way to alter this balance. Also note that this experimental design investigates lake-terminating glaciers under varying frontal and basal conditions and prohibits the analysis of a land-terminating glacier, which would ask for a different glacier terminus morphology.

We performed a basic sensitivity analysis in the situation of $\mathrm{\Delta }\mathrm{D}=-\mathrm{10}$ m by varying the glacier width by ±300 m, ice thickness uniformly by ±50 m and surface slope by ±1.5^∘. We performed additional analyses to investigate the model sensitivity to the large current uncertainty of ice thickness estimates for Himalayan glaciers (Farinotti et al., 2019) by varying the ice thickness again but keeping the maximum velocity at 50 m yr⁻¹ by tuning A_s accordingly.

4 Results

4.1 Algorithm performance

The co-registration error calculated on non-glacierized, stable areas enabled us to reduce the dispersion (MAD_disp) by 56 %, resulting in a MAD_disp distribution with a median at 4.15 m yr⁻¹ (Fig. 2a). The distribution is heavy tailed, with the largest uncertainties found over accumulation zones where the algorithm was unable to remove all mismatches (Fig. A2 in the Appendix). Another large source of uncertainty is the interannual variability in glacier flow, resulting in high dispersion in areas with an overall high flow velocity. The CI₉₅ distribution is slightly less heavy tailed than MAD_disp, with a median uncertainty just below 3 m (Fig. 2b).

Figure 2 Dispersion (MAD_disp) (a) and 95 % confidence interval (CI₉₅) (b) of the velocity estimates at glacierized areas. (a) The red distribution represents the MAD_disp before subtracting V_off from each image pair, which realized a 56 % reduction of the median dispersion, resulting in the grey dispersion distribution. (b) CI₉₅ resulting from the estimated distribution of median displacement vector as a function of the number of velocity estimates and MAD_disp (Eq. 2).

When evaluating the CI₉₅ along the centrelines, a consistent trend is apparent (Fig. 3). The uncertainty decreases from the ELA moving further into the ablation zone due to the enhanced pixel contrast. Close to the terminus however, the uncertainty rises again, which can be related to relatively large interannual changes in surface properties, resulting in reduced algorithm performance. Interestingly, lake-terminating glaciers have consistently higher uncertainty along the ablation zone, which likely results from the large velocity differences between lake-terminating and land-terminating glaciers (Sect. 4.3). The approach of applying a Gaussian window to the velocity estimates reduced the mean CI₉₅ of lake-terminating and land-terminating glaciers by 24 % and 21 % respectively. In the following sections we consider velocity estimates with a CI₉₅ larger than 5 m yr⁻¹ too uncertain, and these estimates are removed from further analyses.

Figure 3 Median CI₉₅ (m yr⁻¹) of the velocity estimates for lake-terminating and land-terminating glaciers along the normalized glacier centre flow line at the ablation zone, with the terminus positioned at the right end of the figure. The spread of CI₉₅ along the centre flow line among the glacier population is represented by the interquartile range (IQR). Velocity estimates above the CI₉₅ threshold (5 m yr⁻¹) are removed from the dataset.

4.2 Comparison to other glacier velocity datasets

The lack of ground-truth velocity measurements hinders simple evaluation of remotely sensed measurements in most cases (Scherler et al., 2008). To assess the quality of our measurements, we compare them with two region-wide velocity datasets, both processed with predominantly Landsat 8 imagery. Dehecq et al. (2019a) produced a composite glacier surface velocity field for the Pamir–Karakoram–Himalaya for the years 2013–2015. Velocity fields are available at 120 m resolution and produced using a 240 m reference window (Dehecq et al., 2015). Another region-wide dataset was generated using autoRIFT (autonomous Repeat Image Feature Tracking; Gardner et al., 2020) and provided by the NASA MEaSUREs ITS_LIVE project (Making Earth System Data Records for Use in Research Environments Inter-Mission Time Series of Land Ice Velocity and Elevation; Gardner et al., 2020). This velocity field spans from 1985 to 2020, but we compare our data to an ITS_LIVE velocity field with a central date of around spring 2018, which is available at 120 m resolution and again is computed using a 240 m reference window.

Substantial differences exist in the region-wide median centreline surface velocity (later referred to as velocity) between the three datasets, with maxima ranging from just above 5 m yr⁻¹ (Gardner et al., 2019) to well above 13 m yr⁻¹ in our study (Fig. 4a). Velocity measurements from the three datasets agree reasonably well close to the terminus (within the 0.5 m yr⁻¹ range) where velocities are expected to be close to stagnant, which indicates that differences between flow fields are proportional to the magnitude of the regional median centreline velocity. Dehecq et al. (2019a) observed a slowdown of glacier velocities for all our subregions, ranging from −14.5 ± 1.3 % to −21.0 ± 2.3 % per decade, in response to climate-induced changes in slope and ice thickness. These reductions in surface velocity only partly explain the differences in velocity between Dehecq et al. (2019a) and Gardner et al. (2019).

The limited width of some Himalayan valley glaciers, which can be as little as 300 m, may explain discrepancies in the velocity fields derived using different methods. Narrow valley glaciers are subject to considerable lateral-stress-induced transverse velocity gradients. A large reference window size might consequently result in an underestimation of the centreline velocity, as it is simply unable to resolve this velocity gradient. The usage of a variable reference window size by Dehecq et al. (2019a) and Gardner et al. (2020) up to 4 times the original reference window size could also potentially explain these differences, although we are unable to verify this, as the spatial variability of the effective reference window size is not documented in these datasets. Selecting only large glaciers, often with a wider ablation zone, largely diminishes this discrepancy (Fig. 4b), which supports our hypothesis, indicating that the employment of the Sentinel-2 satellites improved the resolution and therefore the analytical potential of the glacier centreline velocity data.

Figure 4 Regional median centre flow line surface velocity (m yr⁻¹) comparison between Dehecq et al. (2019b), Gardner et al. (2019) and this study of glaciers greater than 3 km² in area (a) and glaciers greater than 10 km² in area (b).

4.3 Terminus type variability in velocity

Our velocity analysis shows (Fig. 5a) that the along-flow-line mean of median centreline surface velocities (later referred to as mean velocity) of lake-terminating glaciers (median of 18.83 m yr⁻¹ and IQR uncertainty range of 17.50 to 18.01 m yr⁻¹) is substantially higher than the mean velocity of land-terminating glaciers (8.24(8.17–8.35) m yr⁻¹; please note that the median is followed by the IQR uncertainty estimate in brackets) (Table 4). Differences are negligible at the ELA but become steadily larger throughout the ablation zone. Over the lower ablation zone, the differences in surface velocity reach 13.8 m yr⁻¹ (mean velocity of 17.72(17.41–18.02) m yr⁻¹ for lake-terminating and 3.91(3.84–3.97) m yr⁻¹ for land-terminating glaciers) (Table A1). Land-terminating glaciers show a stagnant terminus with only little spread around the median velocity among the glacier population. On the contrary, the median velocity of lake-terminating glaciers decreases only slightly but shows a very large spread, indicating a large heterogeneity in lake-terminating dynamical behaviour.

Overall, lake-terminating glaciers cover a larger surface area and show a slightly higher mean surface slope over the ablation zone (Table 5), which might partially contribute to the overall contrast in the mean velocities (Bahr et al., 1997; Scherler et al., 2011). However, this does not explain the large contrast in velocity or heterogeneity at the glacier terminus. Interestingly, when only focusing on the lowermost portion of the ablation zone, where the velocity contrast between land-terminating and lake-terminating glaciers is greatest, the mean surface slope of lake-terminating glaciers (7.2(9.7–4.9) ^∘) is within the range of the slope of land-terminating glaciers (8.2(11.4–5.5) ^∘) (Table A1), suggesting that factors other than slope are responsible for the velocity contrast close to the glacier terminus.

Figure 5 Median centre flow line surface velocity (m yr⁻¹) (a) and coverage (%) of the velocity estimates (b). (a) The spread of the velocity along the centre flow line among the glacier population is represented by the IQR. (b) The coverage is defined by the percentage of valid velocity estimates (CI₉₅ < 5 m yr⁻¹) at a given position along the centre flow line.

Table 4 The median of the along-flow-line mean regional centre flow line velocities of lake-terminating and land-terminating glaciers. Uncertainty estimates are represented by the IQR of the sample median estimates, calculated according to Sect. 3.3. For the location of the subregions, see Fig. 1.

Table 5 Key characteristics of lake-terminating and land-terminating glaciers. Values refer to the median, whereas the spread among these parameters is represented by the IQR.

4.4 Velocity dependence on orientation and surface area

We noted that glaciers flowing north onto the TP typically have larger accumulation zones and less debris cover compared to glaciers located in catchments draining to the south of the main Himalayan orographic divide, as also reported elsewhere (Scherler et al., 2011). Visual inspection indicated that the highest velocities are found at those northwards-flowing glaciers, especially in Central West 1, Central 1 and East Himalaya, implying a positive correlation between glacier orientation and mean velocity. Concurrently, a large fraction of the total number of lake-terminating glaciers are orientated northwards, which might falsify the apparent relationship between surface velocity and terminus type proposed in the previous section.

To investigate the link between the dynamics and orientation, we therefore subdivided our dataset depending on the orientation of the glacier ablation zone (Fig. 6a). The results show a large heterogeneity for lake-terminating glaciers, with the highest velocities shown for glaciers with their ablation zone orientated to the north. Notwithstanding, for all orientations, lake-terminating glaciers show a higher mean velocity than land-terminating glaciers, although the contrast is minor for glaciers flowing east- or southwards. When only considering the lower half of the ablation zone however, the contrast between land-terminating and lake-terminating glaciers becomes substantial for all orientations (Fig. A3a in the Appendix).

Figure 6 Boxplot showing the mean velocity contrast between lake-terminating and land-terminating glaciers depending on the orientation of the ablation zone (a) and surface area (b). The IQR (boxes) represents the spread within the mean sample group. Points outside of the third quartile plus 1.5 times the IQR range are plotted explicitly.

We examined the relationship between glacier surface area and mean glacier velocity by separating glaciers into three area bins with equal sample sizes. Higher mean velocities are apparent for lake-terminating glaciers for each glacier area bin (Fig. 6b), with the largest contrast for glaciers greater than 10 km². Note that in the largest size bin glaciers are not bounded by an upper area limit but nevertheless show a comparable median area of 19.19(14.01–24.07) km² for lake-terminating glaciers and 18.76(12.97–24.72) km² for land-terminating glaciers. Velocity outliers are particularly abundant at large (> 10 km²) northwards-flowing land-terminating glaciers, such as the clean-ice Zeng glacier (28^∘14^′ N, 90^∘14^′ E) in East Himalaya, which shows a mean velocity of about 93 m yr⁻¹. Again, contrasting velocities between lake-terminating and land-terminating glaciers increase when solely considering the lowermost portion of the ablation zone (Fig. A3b in the Appendix), indicating that regardless of orientation and size, substantial contrast in glacier surface velocity is related to terminus type and increases towards the glacier tongue.

4.5 Regional variability

We find large variability in mean velocities between different regions (Fig. 7), with the highest mean velocities in Central West 1 (13.02(12.78–13.36) m yr⁻¹) and East Himalaya (13.07(12.86–13.36) m yr⁻¹) (Table 4), areas with the largest proportions of clean ice. All regions show higher mean velocities for lake-terminating glaciers than for land-terminating glaciers, though large variability between regions is apparent. In Central 2 mean velocity differences between lake-terminating (11.78(10.25–13.16) m yr⁻¹) and land-terminating glaciers (6.03(5.91–6.17) m yr⁻¹) are relatively modest and coincide with a high proportion of debris-covered glaciers for both terminus types. For Central 1 and East Himalaya, a substantial part of the velocity contrast can be attributed to the relatively high abundance of lakes at large clean-ice, northwards-flowing glaciers, explaining the large velocity contrast which is already substantial at the ELA. Finally, in the regions Central West 1 and Central 1 we observe an increase in velocity towards the terminus, indicating that most of the glaciers accelerate towards the ice–water interface. Trends in the velocity of lake-terminating glaciers in some of these regions should the treated with caution however, where the population of lake-terminating glaciers is very limited (6 glaciers in Central West 1 and 10 glaciers in Central 2). Nonetheless, all regions show a large contrast in heterogeneity close to the terminus, suggesting that the influence of proglacial lakes on glacier dynamics is a region-wide phenomenon.

Figure 7 Subregional glacier median centre flow line velocity estimates and their location along the CE Himalaya (red rectangles). The spread of the velocity along the centre flow line among the glacier population is represented by the IQR, with the number of glaciers shown in the legend between the brackets.

4.6 Impact of surface cover on glacier–lake dynamics

To examine the role of debris cover on glacier–lake dynamics, we subdivided our dataset into glaciers with > 19 % debris cover and those with < 19 % debris cover, which we classify as clean-ice glaciers (Fig. 8a, c, e). We measured substantially higher velocities for lake-terminating glaciers, both debris covered and clean ice (Fig. 8c, e; Table 6), although large differences are apparent depending on surface type. Most notably, the absolute lake–land mean velocity contrast of debris-covered glaciers (11.48(11.11–11.98) m yr⁻¹ vs. 5.96(5.90–6.05) m yr⁻¹) is lower than for clean-ice glaciers (22.52(22.21–22.83) m yr⁻¹ vs. 10.26(10.11–10.33) m yr⁻¹) but indicate for both surface types a doubling in surface velocity when a lake is present.

Figure 8 Glacier median centre flow line velocity (m yr⁻¹) (a, c, e) and surface elevation change (dh dt⁻¹) estimates (after King et al., 2019) (b, d, e) for lake-terminating and land-terminating glaciers. A further subdivision is made between debris-covered (c, d) and clean-ice glaciers (e, f). The spread of the velocity along the centre flow line among the glacier population is represented by the IQR.

For both surface types, higher velocities for lake-terminating glaciers are coincident with elevated surface lowering over coincident portions of glacier ablation zones (Fig. 8b, d, f). Again, large differences exist depending on surface type, with a very large contrast in surface lowering close to the termini of debris-covered, lake- and land-terminating glaciers (about 1.5 m yr⁻¹) and a less pronounced contrast for clean-ice glaciers (0.5 m yr⁻¹). We relate this to the distinct differences in surface mass balance properties between clean-ice and debris-covered glaciers, which becomes clearly visible for land-terminating glaciers.

For clean-ice, lake-terminating glaciers, the velocity profile remains close to constant towards the terminus and shows a very large spread among the glacier population (Fig. 8e). Enhanced surface lowering at clean-ice, lake-terminating glaciers steadily grows to −0.5 m yr⁻¹ towards the terminus and is less pronounced than the surface lowering contrast at debris-covered glaciers. We find no substantial differences in altitudinal distribution between lake-terminating and land-terminating glaciers that could partly explain this difference. The large velocity contrast coincides with a larger surface area of clean-ice, lake-terminating glaciers compared to clean-ice, land-terminating glaciers (8.45(5.01–11.92) km² vs. 4.72(3.56–5.70) km²) and longer ablation areas of clean-ice, lake-terminating glaciers (3800(2910–4721) m) compared to clean-ice, land-terminating glaciers (3040(2497–3541) m).

Table 6 The median of the along-flow-line mean regional centre flow line velocities of lake-terminating and land-terminating glaciers, subdivided by surface type. Uncertainty estimates are represented by the IQR of the sample median estimates, calculated according to Sect. 3.3.

Debris-covered, lake-terminating glaciers do not have the same concave-up velocity profile which characterizes their land-terminating counterparts. Their velocity profiles show a larger spread than debris-covered, land-terminating glaciers close to their termini. Debris-covered, lake-terminating glaciers show distinctly enhanced surface lowering very close to their termini, with rates of surface lowering exceeding those on the ablation zone of land-terminating glaciers. Notably, these lake-terminating glaciers are generally smaller in surface area (6.78(4.23–9.11) km² vs. 9.42(5.92–13.02) km²), and the ablation zones are much shorter (2720(1834–3680) m vs. 5680(4003–7218) m) than clean-ice, lake-terminating glaciers. Many of the proglacial lakes in front of these small debris-covered, lake-terminating glaciers are between 1 and 4 km in length, which likely explains the difference in debris-covered glacier length depending on terminus type.

4.7 Synthetic numerical experiment

Analyses of our glacier surface velocity dataset indicate substantial glacier- and regional-scale variability in terminus-proximal ice flow. With the synthetic numerical experiment, we can gain useful insights into the potential drivers of this glacier-scale variability in ice flow. The results indicate that both changes at the frontal boundary condition and in basal friction can alter lake-terminating glacier dynamics (Fig. 9), with velocities at the terminus increasing in response to a reduction in effective pressure and in response to lake surface lowering. When the glacier front reaches flotation (ΔD=0 m), the velocity within ∼1 km proximity to the terminus increases significantly if the glacier bed friction is dependent on the effective pressure. An acceleration towards the terminus is shown when the ice cliff height increased through a sufficient lowering of the surface lake level (ΔD≥10 m), resulting in a larger frontal imbalance. We find that a smaller cliff height (ΔD<10 m) has only a limited effect on the glacier dynamics, as the basal drag quickly increases when moving away from a buoyant situation, which is related to the non-linear sliding law we adopted in this study (m = 3). Beyond a certain increase of the cliff height (ΔD>10 m), a further increase of the cliff height above buoyancy through the lowering of the lake level rapidly increases the surface velocity.

Figure 9 Velocity results from the numerical experiment for three varying lake levels (ΔD), for both effective-pressure-dependent roughness parameter and a constant roughness parameter (A_s). Blue dotted line indicates the piezometric surface for when ΔD = 0 m.

Our sensitivity experiments show a high dependence on ice thickness (Fig. 10a, b) for both a constant roughness parameter A_s and for a varying A_s to account for the ice thickness uncertainty. This is a direct result of the frontal dynamic boundary condition (Eq. 7), which is a function of the ice thickness. Also, the velocity sensitivity at the terminus to ice thickness might explain the limited dynamic impact at debris-covered glaciers, which are often relatively thin at their termini. Glacier width or surface slope also modify the velocity field substantially but do not heavily influence the relative change in velocity towards the glacier terminus (Fig. 10c, d).

Figure 10 Velocity sensitivity experiment to ice thickness (a), ice thickness with varying roughness factor (b), terminus width (c) and slope (d).

5 Discussion

5.1 Dynamic heterogeneity of lake-terminating glaciers

The contrasting morphological attributes of lake-terminating glaciers on either side of the main orographic divide likely substantially influences their flow regime and tendency to retreat. Clean-ice, lake-terminating glaciers, mainly found on the north facing slopes of the TP, are larger in size and have a longer ablation zone than clean-ice, land-terminating glaciers. Debris-covered, lake-terminating glaciers, predominantly found on the southern side of the main orographic divide, are generally much smaller and have shorter ablation zones than their land-terminating counterparts. This contrast in glacier dimensions illustrates how the evolution and occurrence of lake-terminating glaciers must be put in context of the surface cover properties of lake-terminating glaciers, which is related to the morphological settings in which clean-ice and debris-covered glaciers are prone to develop.

Clean-ice, land-terminating glaciers with extensive accumulation zones, often flowing onto the TP, possessed enough erosive power to form overdeepenings and large terminal moraines (Scherler et al., 2011). Overdeepenings beneath glaciers flowing onto the TP are not only promoted by terminal moraines but also inherent features of the reversed slope of the TP itself (Royden et al., 2008), making these localities a hotspot of proglacial-lake development. Himalayan clean-ice, lake-terminating glaciers are often found at such localities, are large and show consequently higher velocities over the entire ablation length than their land-terminating counterparts.

Unlike for clean-ice glaciers, our results suggest that there is no clear glacier surface-area-related preference for proglacial-lake development on debris-covered glaciers. Akin to previous studies, we find that low-gradient, debris-covered ablation zones of many Himalayan glaciers (Steiner et al., 2019; Wijngaard et al., 2019) are the result of reversed surface mass balance gradients (Bisset et al., 2020). Such ablation zones act as a sweet spot for proglacial-lake development (Benn et al., 2012; Quincey et al., 2007), which become bounded by a stagnant, ice-cored moraine dam. The development of a proglacial lake leads to a transformation which is associated with a drastic increase in retreat rates (Basnett et al., 2013; King et al., 2019; Watson et al., 2020) and results in lake-terminating glaciers of shorter length than land-terminating glaciers.

As a result, debris-covered, lake-terminating glaciers develop from the glaciers whose area is close to the median of the land-terminating glacier population, whereas clean-ice, lake-terminating glaciers predominantly evolve from larger land-terminating glaciers. This, together with the overrepresentation of clean-ice glaciers in the total lake-terminating glacier population (50 out of 70), explains a large part of the velocity contrast between land-terminating and lake-terminating glaciers (Fig. 8a). Also, this overrepresentation of clean-ice glaciers in the lake-terminating glacier population explains a significant part of the contrasting thinning observed in Fig. 8b, which makes it erroneous to attribute this contrast entirely to dynamic thinning. Notwithstanding, the large heterogeneity in velocity at the terminus of lake-terminating glaciers, with an accelerating velocity for almost half of the population (Fig. A4 in the Appendix), and elevated surface lowering for both debris-covered, lake-terminating and clean-ice, lake-terminating glaciers clearly show that dynamic thinning is a process that must be considered.

5.2 Drivers of dynamic thinning

The acceleration towards the termini of half of the lake-terminating glaciers in our study area (Fig. A4) strongly indicates an influence of lakes on glacier dynamics. A visual inspection of Fig. 11 underlines the regional extent over which dynamic thinning can be observed, although again, a large heterogeneity in the acceleration at the terminus is visible. Our numerical experiment indicates that changes in the force balance at the glacier termini, ice thickness and to a lesser degree basal friction likely play a profound role in accelerating the glacier terminus. As a consequence, large heterogeneity in the acceleration might be at least partly attributed to the varying influence of these factors. However, region-wide measurements in the Himalaya of the ice cliff heights or ice thickness are absent or highly uncertain, restricting us to a mainly qualitative evaluation of processes that drive dynamic thinning.

Figure 11 Examples of lake-terminating glaciers accelerating towards their terminus with glacier attributes within 2 km of the terminus. White numbers represent the glacier surface velocity in m yr⁻¹. Colour scale of plotted velocity data is indicative and varies among glaciers. Velocity data and RGB (red–green–blue) images are retrieved from Sentinel-2 Multispectral Instrument (MSI) data.

A reduction in basal friction caused by the development of a proglacial lake is often designated as one of the main drivers of dynamic thinning (Carrivick et al., 2020; King et al., 2018; Liu et al., 2020; Sutherland et al., 2020; Tsutaki et al., 2013, 2019), though several remarks must be made when evaluating the importance of this resisting force. Before the development of a proglacial lake, a hydrological base level is already present due to impermeable moraine at the front. The hydraulic potential at the glacier bed in these overdeepening localities cannot be less than that at the base level, and the presence of a proglacial lake will therefore be of no influence on the effective pressure if the glacier surface had no time to dynamically adjust, as described by Tsutaki et al. (2019). In that case, the force balance at the glacier front entirely determines the dynamic influence of a glacial lake and results from a reduction in buttressing caused by a hypothetical instantaneous expansion of the proglacial lake at the expense of the glacier front.

However, such a situation might be unrealistic, as we have just concluded that a glacier will dynamically respond and thin, bringing the calving front closer to flotation. The reduction in basal stress then largely depends on the sliding law, for which we took a non-linear Weertman type that includes effective pressure dependency (m = 3). Here, basal stress dramatically decreases when the glacier front approaches flotation, which in turn leads to further acceleration of glacier flow (Benn et al., 2007b). Nevertheless, as our model results suggest, the reduction in the force imbalance at the front through thinning might dominate in the dynamic response, resulting in an overall decelerating flow when the glacier reaches flotation. This is indicative of a dynamic regime where imbalances at the frontal boundary, through for example glacier retreat rates or a lake level lowering, are balanced by enhanced dynamic thinning rates, resulting in a reduced acceleration of the flow, which is conceptually in line with Nick et al. (2009).

At an early stage of lake development, however, when a clear calving front is yet absent, Tsutaki et al. (2013) showed that a reduction in basal friction can have a key control on the velocity evolution of the glacier through dynamic thinning and subsequently promoting the development of transverse crevasses and a disintegration over the expanded area. Also, when the glacier only thins locally at the calving front, the increase in surface slope might promote an acceleration of the flow again (Benn et al., 2007a), something this study has not been able to test without violating the assumptions of the SSA model and which requires the use of a higher-order model for further investigation.

The scarcity of local data on glacier ice thickness or lake depth at the calving front makes it difficult to evaluate the representativeness of the frontal configuration we used in our experimental setup. Watson et al. (2020) found a mean calving front height of three Nepalese debris-covered, lake-terminating glaciers varying between 27 and 41 m, which indicates that we can assume with a relatively high confidence that these glaciers are at least 10 m above flotation, considering the upper-range ice thickness estimates from Farinotti et al. (2019). Also, for all three glaciers, a surface acceleration in the proximity of the glacier front was observed (Watson et al., 2020), indicating that the force balance at the boundary condition might dominate over the importance of an in situ reduction in basal drag. Remarkably, these lake-terminating glaciers were heavily debris covered, indicating that also for these glaciers similar processes are relevant.

The frontal ice thickness itself is a variable that needs more consideration when evaluating drivers of frontal ice velocity. Evidently, lake-terminating glaciers are thicker near the terminus than land-terminating glaciers, since they end at a calving cliff rather than at a front that thins to zero. As ice thickness drives ice flow (see the right-hand side of Eq. 5), a substantial part of the terminal velocity contrast between lake-terminating and land-terminating glaciers could then be attributed to this difference in ice thickness. Indeed, comparison of the median ice thickness of our glacier sample group (Fig. A5), using the ice thickness estimates from Farinotti et al. (2019), indicates that lake-terminating glaciers are on average thicker over the ablation zone (∼ 110 m) than land-terminating glaciers (∼ 100 m). This difference in thickness becomes especially substantial near the terminus.

Using the approach of Dehecq et al. (2019b), where changes in surface velocity are related to changes in driving stress, a first-order approximation can be made to estimate the importance of the ice thickness onto glacier dynamics, assuming factors such as surface slope and glacier width remain constant and a proglacial lake plays no role in controlling the surface velocity. This approximation only holds in the limit of $\frac{\mathit{\delta }H}{H}\ll \mathrm{1}$, and therefore only the ice thickness contrast over the entire ablation zone is considered. Here, we find that 28 % to 64 % of the observed velocity contrast can be explained by this approximation, indicating that factors other than ice thickness play a role in controlling the ice dynamics.

This indicates that ice thickness data are unable to explain the whole velocity contrast at the glacier terminus. Also, Fig. A5 shows a clear decrease in ice thickness for both land-terminating and lake-terminating glaciers towards the terminus. At the same time, the lake-terminating glacier velocity does not show a deceleration towards the terminus and even accelerates for half of the glacier sample group (Fig. 5). Similarly, we measured substantially faster flow over lake-terminating glaciers within the same area classes as land-terminating glaciers (Fig. 6b), and as ice thickness generally correlates with glacier area (Bahr et al., 2015), it is not unreasonable to expect glaciers within these same area classes to be of a comparable thickness. It is noteworthy that uncertainties in ice thickness estimates are significant and that errors could be systematic depending on surface type. Therefore, these results of analyses which incorporate ice thickness estimates should be treated with caution until a sufficient number of direct measurements of terminus ice thickness are available.

Lake-terminating glacier dynamics can only be understood as inseparable parts of an intricately coupled system with frontal ablation (Benn et al., 2007b), consisting of mechanical calving and subaqueous melt (Carrivick and Tweed, 2013). If frontal ablation rates are considerably high, dynamic thinning rates are expected also to be high, in order to bring the ice front closer to a balanced state. This is in line with annual velocity and glacier retreat observations at the Longbasaba Glacier, where Liu et al. (2020) found the glacier acceleration to be driven by glacier retreat through calving since the mid-1990s. Interestingly, an above-average glacier acceleration was observed in 2006 after 5.6 m surface lake lowering as a mitigation measure in 2005 (Xiao and Dai, 2011). Again, a local increase in the surface slope as a result of dynamic thinning could play a key role here by promoting calving fluxes under close to buoyant conditions (Benn et al., 2007b).

Another important factor is the reverse gradient of the bed slope, which is an inherent property of proglacial lakes (i.e. Somos-Valenzuela et al., 2014), resulting in an instability as the force imbalance at the front will increase with a glacier retreating to deeper water with greater ice thickness. This phenomenon is known for marine ice sheets as marine ice sheet instability (MISI) (Katz and Worster, 2010; Weertman, 1974) and is shown to be also an important mass loss feedback mechanism in lake-terminating settings (Sutherland et al., 2020). With such a feedback mechanism, the longitudinal-stress gradient plays a large role by distributing the dynamic response triggered at the calving front from the reduction of buttressing (Benn et al., 2007b). This is in line with findings on Alaskan lake-terminating glaciers, where the specific geometry and topographic setting of proglacial lakes are suggested to play a key control on glacier evolution (Field et al., 2021).

Finally, it needs to be considered that glacier width and surface slope are further relevant factors in the dynamic evolution of lake-terminating glaciers, as local variations of glacier width or high slopes in close proximity to the terminus are known to be imperative in the context of the transient evolution of these glaciers (Benn et al., 2007b). However, in this study we merely focus on direct processes that accelerate the glacier terminus, and drawing conclusions about these factors in a transient context would be outside the scope of this study.

5.3 Implications for future evolution of Himalayan glaciers

Regional ice loss through lake-terminating dynamics will remain important in the near future, given the sustained expansion of proglacial lakes across the Himalayan region (Chen et al., 2021; Nie et al., 2017; Zhang et al., 2015) and the susceptibility of many debris-covered glaciers for proglacial-lake development. Our results also emphasize the importance of clean-ice, lake-terminating glaciers terminating in overdeepenings (Linsbauer et al., 2016), and their future contribution to regional ice loss might be disproportionately large, considering their active flow and propensity to thin dynamically. Many of these clean-ice glaciers drain northwards into the tributaries of the Brahmaputra, a river of which the meltwater supply is of high importance during the dry season (Pritchard, 2019), and changes in ice mass loss projections are of essential importance for millions of people in downstream regions (Immerzeel et al., 2020).

In order to better understand the impact of proglacial lakes on glacier dynamics and to find out whether the contribution of lake-terminating glaciers to Himalayan ice mass loss may increase further, spatially resolved, multi-temporal analyses of glacier–lake dynamics are needed, such as those by Liu et al. (2020) and Watson et al. (2020). In this context, measurements of lake depth, ice cliff height, ice thickness and surface slope in the proximity of the calving front are essential and will help to constrain factors controlling the dynamics of these lake-terminating glaciers. There is also a need for more detailed modelling studies on glacier–lake dynamics with a particular focus on basal friction during the transient evolution of lake-terminating glaciers in alpine settings with respect to the importance of the force balance at the glacier front. Finally, considerable progress still can be made by linking frontal mass loss processes that are found to be characteristic for lake-terminating settings (Watson et al., 2020) to glacier dynamics as a fully intercoupled system.

6 Conclusions

In this study, we documented 2017–2019 surface velocities in the ablation zone of glaciers larger than 3 km² by analysing Sentinel-2 optical satellite imagery at five proglacial-lake-prevalent subregions in the Himalaya. Our results show that the enhanced resolution of Sentinel-2 (10 m) with respect to Landsat 8 (15 m) improves image matching and yields a better resolved velocity field over small glaciers, thereby improving the potential for the analysis of glacier centreline velocities.

Analysis of the centreline velocity profiles revealed that lake-terminating glaciers display substantially higher flow velocities than land-terminating glaciers (18.83(18.55–19.06) m yr⁻¹ vs. 8.24(8.17–8.35) m yr⁻¹). This finding is consistent regardless of the orientation, glacier size and subregion of the glacier population. The velocity contrast between clean-ice, lake- and land-terminating glaciers is much greater than for debris-covered glaciers, and therefore a major contribution of the mean velocity difference can be attributed to the overrepresentation of large clean-ice glaciers in the lake-terminating population.

Notwithstanding, both clean-ice and debris-covered, lake-terminating glaciers show heterogeneous behaviour at the glacier terminus, remain dynamically active along the entire flow and show an accelerating trend for almost half of the glacier population, revealing that dynamic thinning is prevalent in the Himalayan region. In line with this, we found that a positive correlation between high terminal velocities and elevated surface lowering is evident for both surface types.

Our synthetic numerical ice-flow model experiment revealed that the surface flow velocity is most sensitive to changes in the boundary condition at the terminus of a lake-terminating glacier, with variations in basal friction playing a less prominent role in our model setup. Rapid changes in terminus position or proglacial-lake level could therefore contribute greatly to the dynamic evolution of the glacier front itself. Further analyses emphasized the importance of ice thickness in the glacier–lake dynamics, which might explain the limited dynamic impact at the, often relatively thin, termini of debris-covered glaciers.

The contribution to ice mass loss and, hence, runoff from lake-terminating glaciers is unlikely to diminish in the near future, but the exact contribution to downriver meltwater supply in the next decades is still highly uncertain. An improved understanding of lake-terminating glacier dynamics, both by field observations and numerical studies, is therefore imperative for the future of those people that depend on a year-round meltwater supply from the major Himalayan rivers.

Appendix A

Table A1 The median of the along-flow-line mean regional centre flow line velocities and median slope of lake-terminating and land-terminating glaciers for the second half of the ablation zone. Uncertainty estimates for the velocities are represented by the IQR of the sample median estimates, calculated according to Sect. 3.3. The spread among the slope is also represented by the IQR.

Figure A1 Examples of lake-terminating glaciers with glacier attributes within 2 km of the terminus. RGB images are retrieved from Sentinel-2.

Figure A2 Velocity dispersion (MAD_disp). Yellow colours indicate a dispersion > 10 m yr⁻¹. The hillshade is produced using ALOS World 3D DEM.

Figure A3 Boxplot showing the mean velocity contrast between lake-terminating and land-terminating glaciers depending on the orientation of the second half of the ablation zone (a) and surface area (b). The IQR (boxes) represents the spread within the mean sample group. Points outside of the third quartile plus 1.5 times the IQR range are plotted explicitly.

Figure A4 Separation of accelerating and decelerating lake-terminating glaciers. A glacier is considered to be accelerating if within 500 m of the terminus it shows higher surface velocities than the 500–1500 m proximity range. The spread of the velocity along the centre flow line among the glacier population is represented by the IQR.

Figure A5 Ice thickness estimates for land-terminating and lake-terminating glaciers, using the dataset of Farinotti et al. (2019). The spread of the ice thickness along the centre flow line among the glacier population is represented by the IQR.