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 km2. We find that centre flow line velocities of lake-terminating glaciers (N = 70; umedian: 18.83 m yr−1; IQR – interquartile range – uncertainty estimate: 18.55– 19.06 m yr−1) are on average more than double those of landterminating glaciers (N = 249; umedian: 8.24 m yr−1; IQR uncertainty estimate: 8.17–8.35 m yr−1) 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.


Introduction
Himalayan glaciers provide an important baseline supply of meltwater for downstream areas 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;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 −1 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 −1 . However, within the Himalayan mountains, large intra-regional variability in glacier mass loss exists (Azam et al., 2018; Published by Copernicus Publications on behalf of the European Geosciences Union. 5578 J. B. Pronk et al.: Contrasting surface velocities between lake-and land-terminating glaciers 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 laketerminating glaciers distinctively different from their landterminating 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 laketerminating glacier flow acceleration over a similar time pe-riod (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 Shugar et al., 2020;Zhang et al., 2015Zhang et al., , 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 laketerminating 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.

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 2 , which is about 60 % of the glacierized area of the total Himalayan arc . 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).
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 2 , compared to the 5 km 2 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 2 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 glaciallake inventory of Wangchuk and Bolch (2020) and  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.   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 .
Throughout most of its mission, the multi-temporal coregistration 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  and have therefore been omitted in this study.

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.

Image processing
To reduce computational costs, a mask was applied over all non-glacierized areas and glaciers with an area below 3 km 2 . For this mask we used the glacier outlines from the Randolph 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 2 . The names of the subregions are in accordance with King et al. (2019). Country boundaries are tentative and for orientation only.  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. 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).
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 xand y-displacement pair is disregarded. Secondly, an xand ydisplacement 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) xand y-displacement field. The final 2017-2019 xand 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 featuretracking algorithm.

Uncertainty of the velocity field
Uncertainties of the median velocity field are dominated by the precision of the feature-tracking algorithm, the coregistration error, the temporal variability of glacier flow and the number of velocity estimates. For the estimation of the 95 % confidence interval (CI 95 ) of each median velocity component, we adopted the methodology of Dehecq et al. (2015) and expect the CI 95 to conform to where MAD disp (mean absolute deviation) is the dispersion at each velocity location of the N number of estimates: where T is the collection of N velocity estimates V (i, j, t) merged to obtain the median velocity V (i, j ) 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 95 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.

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 95 ), 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 5 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.

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: where u is the velocity estimate at pixel (i, j ), CI −2 95 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.

Surface elevation change, estimation of ELA and surface slope
We examined ice-thinning rates using the surface elevation change (dh dt −1 ) 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 −1 for lake-terminating glaciers (King et al., 2019), the lowermost 100 m of laketerminating 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 laketerminating 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., 2010Nick et al., , 2009Vieli 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 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) 2 where U is the vertically averaged horizontal ice velocity; ρ i = 917 kg m −3 is the density of ice; ρ w = 1000 kg m −3 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: 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 temperaturedependent rate factor and increases from a minimum of 3.5×10 −25 Pa −3 s −1 at the divide to a maximum of 1.7 × 10-24 Pa −3 s −1 at the calving front, corresponding to a depthaveraged 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: 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 concaveup 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 −1 (Dehecq et al., 2019a;Gardner et al., 2020) and found a value of A s = 2.5×10 6 Pa m −2/3 s 1/3 , which is on the low side compared to values used at marineterminating outlet glaciers of the Greenland ice sheet (i.e. Enderlin et al., 2013;Nick et al., 2009).

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 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 . 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 − ρ w ρ i D = 1), ruling out the effect of the lake on basal friction, and used a corresponding roughness parameter of A s = 9 × 10 6 Pa m −1/3 s 1/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 D = −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 −1 by tuning A s accordingly.

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 −1 (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 95 distribution is slightly less heavy tailed than MAD disp , with a median uncertainty just below 3 m (Fig. 2b).
When evaluating the CI 95 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 landterminating glaciers (Sect. 4.3). The approach of applying a Gaussian window to the velocity estimates reduced the mean CI 95 of lake-terminating and land-terminating glaciers by 24 % and 21 % respectively. In the following sections we consider velocity estimates with a CI 95 larger than 5 m yr −1 too uncertain, and these estimates are removed from further analyses.

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 −1 (Gardner et al., 2019) to well above 13 m yr −1 in our study (Fig. 4a). Velocity measurements from the three datasets agree reasonably well close to the terminus (within the 0.5 m yr −1 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-stressinduced 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.

Terminus type variability in velocity
Our velocity analysis shows (Fig. 5a) that the along-flowline mean of median centreline surface velocities (later referred to as mean velocity) of lake-terminating glaciers (median of 18.83 m yr −1 and IQR uncertainty range of 17.50 to 18.01 m yr −1 ) is substantially higher than the mean velocity of land-terminating glaciers (8.24(8.17-8.35) m yr −1 ; 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 −1 (mean velocity of 17.72(17.41-18.02) m yr −1 for lake-terminating and 3.91(3.84-3.97) m yr −1 for land-terminating glaciers) (Ta-  ble 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 laketerminating 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.

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 northwardsflowing 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).
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 2 . 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 2 for lake-terminating glaciers and 18.76(12.97-24.72) km 2 for land-terminating glaciers. Velocity outliers are particularly abundant at large (> 10 km 2 ) northwards-flowing landterminating 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 −1 . 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.

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 −1 ) and East ) m yr −1 ) (Table 4), areas with the largest proportions of clean ice. All regions show higher mean velocities for lake-terminating glaciers than for landterminating 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 −1 ) and land-terminating glaciers (6.03(5.91-6.17) m yr −1 ) are relatively modest and coincide with a high proportion of debriscovered glaciers for both terminus types. For Central 1 and 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.
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.

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 lakeland mean velocity contrast of debris-covered glaciers (11.48(11.11-11.98) m yr −1 vs. 5.96(5.90-6.05) m yr −1 ) is lower than for clean-ice glaciers (22.52(22.21-22.83) m yr −1 vs. 10.26(10.11-10.33) m yr −1 ) but indicate for both surface types a doubling in surface velocity when a lake is present. For both surface types, higher velocities for laketerminating 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 −1 ) and a less pronounced contrast for clean-ice glaciers (0.5 m yr −1 ). 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 −1 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 2 vs. 4.72(3.56-5.70) km 2 ) and longer ablation areas of cleanice, lake-terminating glaciers (3800(2910-4721) m) compared to clean-ice, land-terminating glaciers (3040(2497-3541) m).
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 2 vs. 9.42(5.92-13.02) km 2 ), 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.

Synthetic numerical experiment
Analyses of our glacier surface velocity dataset indicate substantial glacier-and regional-scale variability in terminusproximal 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.
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).

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 Figure 8. Glacier median centre flow line velocity (m yr −1 ) (a, c, e) and surface elevation change (dh dt −1 ) 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. 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.

Terminus type
Surface type Lake mean (m yr −1 ) Land mean (m yr −1 ) Both mean (m yr −1 ) 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 oc-currence 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 cleanice 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 cleanice, 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, debriscovered ablation zones of many Himalayan glaciers 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 landterminating 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 laketerminating 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 laketerminating 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.

Drivers of dynamic thinning
The acceleration towards the termini of half of the laketerminating 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, regionwide 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.
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., 2013Tsutaki et al., , 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 Figure 10. Velocity sensitivity experiment to ice thickness (a), ice thickness with varying roughness factor (b), terminus width (c) and slope (d). 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 −1 . 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. 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, laketerminating 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 , 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 landterminating 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 δH H 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 landterminating 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 laketerminating glaciers within the same area classes as landterminating 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 termi-nus 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.

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 . 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  to glacier dynamics as a fully intercoupled system.

Conclusions
In this study, we documented 2017-2019 surface velocities in the ablation zone of glaciers larger than 3 km 2 by analysing Sentinel-2 optical satellite imagery at five proglacial-lakeprevalent 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 im-proving 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 −1 vs. 8.24(8.17-8.35) m yr −1 ). 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, laketerminating 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 laketerminating 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. Figure A2. Velocity dispersion (MAD disp ). Yellow colours indicate a dispersion > 10 m yr −1 . 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 laketerminating 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. Author contributions. JBP, TB and OK designed the study. JBP performed all analyses and wrote the draft of the manuscript. JBP, TB, OK and DIB discussed the results and contributed to the writing. All authors edited and contributed to the final form of the manuscript.
Competing interests. The authors declare that they have no conflict of interest.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.