Reply on RC1

In this study, our main aim is to show that a contrast in glacier dynamics for land and laketerminating glaciers is a regional wide phenomenon in the Himalayas. This would be a valuable contribution to the topic of Himalayan lake-terminating glaciers, and laketerminating glaciers in general, and would be a solid starting point for other studies that investigate the temporal evolution and dynamic drivers of such glaciers, of which the latter is a secondary aim of this study.

The authors processed satellite data from ESA's Sentinel-2. While their processing chain is based on previous work (Gardner et al., 2020 ;Dehecq et al., 2015) and seems fairly robust, some points are not clear and needs to be more explicit (postprocessing, selection of stable areas, center line analysis, choices of repeat cycles…) and compared with existing studies that uses Sentinel-2 data to map ice flow velocity of glaciers. Specifically, if the aim of the study was to calculate a composite map, I do not understand why the author only choose to process pairs of images separated by 1 year ? Studies have shown, that using all possible repeat cycles and stacking them would largely improve (1) the signal to noise ratio and (2) the spatial and temporal coverage. This would be of great interest for the authors that are looking at glaciers with frontal velocity that rarely exceeds few tens of meters.

Response:
We thank the referee for their careful consideration of the approach we have followed to derive glacier surface velocity data. The reviewer has made a number of valid points regarding data selection which we agree is an important aspect of such a study. We feel we have largely optimised our approach to derive a high-quality velocity field. We outline our prior considerations of the points raised by the reviewer below. In line 138-139 we state: 'The maximum number of image pairs separated by one year was selected for the month of November, as this month is associated with low cloud cover and a relatively high snow line'.
Indeed, using a large number of repeat cycles improves the temporal coverage. However, after a certain number of image pairs the improvement is only marginal (Dehecq et al. (2015), and at some point, a trade-off must be made between the image matching improvement and the computing time. With an average number of 39  pairs for each image tile, we are confident that the residual error of our velocity field is low. (please see figure 9 from Dehecq et al. (2015)). Secondly, the image matching is largely hindered by cloud cover and surface cover conditions. Outside of the month November, the quality and consistency of the imagery decreases drastically because of 1) large cloud cover during the monsoonal months and 2) decreased surface contrast due to snow showers carried on westerlies which can occur across a broad elevation range. Therefore, deploying all other suitably imagery would require a tremendous amount of extra computing time but would not significantly improve the dataset. However, we agree with the reviewer that this information was lacking in the manuscript. Therefore, we will expand the text around line 138-139 to cover the complicating factors that limit the suitable time window for feature tracking.
I found the measure of the uncertainty on the ice velocity to be somewhat inconsistent throughout the paper: velocity profiles provide a measure of the median and interquartile range, the table are showing the mean and standard error of the mean, and the numbers in the text are not always clear. Something that is even more confusing, most velocity profiles are showing IQR at +/-10 m, (hence it is difficult to draw conclusion from this), while uncertainty on the velocity in the text rarely exceeds 2 m. Right now, it is hard to tell if the difference in velocity pattern between lake and land glaciers is really over the uncertainty?

Response:
We agree that the way we present the measures is somewhat confusing. Indeed, those measures represent something different: The Median and the IQR give by no means a quantification of the uncertainty, but instead provide an insight on the 'typical' (median) velocity and spread in the velocity among different glaciers. The uncertainties presented in the tables are described in section 3.5 provide an uncertainty measure of the along-flowline mean velocities of a group of glaciers, given by the mean and standard error.
As a result, these two quantities indicate something fundamentally different: One (median and IQR) tells the reader something about the characteristics of the sample group. The other shows the confidence we have in our along-glacier mean quantities. To relate those two different measures to each other in a better way, we will merge sections 3.2 and 3.5 in the revised manuscript.
Concerning the standard error of the mean given in the Table 4, how was this calculated ? The authors need to keep in mind that the ice velocities do not follow a Gaussian distribution, hence using the standard error of the mean does not apply.

Response:
The information is provided (see section 3.5 for the relevant methodology). but we agree that the writing needs to be clarified.
We will clarify this issue in the improved version of the manuscript. In short: we sample the background distribution without assuming it is gaussian. From all the resampled means we calculated the 1SE interval, as we assumed this would provide a good indication of the variation of the velocity estimates. However, those sampled means might indeed be (slightly) skewed. We therefore could present this confidence interval by presenting the 1st and 3rd quartile or add a measure of the skewness of the background population to the table with the velocity measurements, as this provides more information about the background distribution.
All of this is a bit confusing, and things needs to be more homogeneous throughout the text in order to have more confidence in the results.

Response:
We apologise that the information is not presented in a clear way. We will improve the presentation of the methods in the revised manuscript.
Additionally, the error estimation is largely based on the analysis of velocity fields on stable ground that are selected by the author. However, the selection of these regions is not clear at all throughout the text. Is it selected randomly (including valleys, mountain peaks etc.)? More details need to be provided in this regard, with potential additional figures.

Response:
We agree that more details need to be provided to be able to fully understand the approach. Our approach was the following: For each image tile we selected a square area of which we could reasonably assume the displacement to be zero. We therefore avoided glaciated terrain and high alpine terrain that might be abundant with glacial features such as rock glaciers, as we do not expect these areas to remain fixed through time. We made sure that that the stable area of each image tile was of sufficient width such that it covered multiple granules (± 25km in width).
We will improve the revised manuscript accordingly.
Moreover, despite the fact that the authors processed a large number of Sentinel-2 data, the entire study is based on the analysis of a composite map, hence completely losing the temporal variability in glacier dynamics. I was in turn, a bit surprised to see that the observation part is only based on the comparison of velocity patterns between lake and land terminating glaciers, despite assessing the real influence of lake changes on glacier dynamics.

Response:
As our glacier surface velocity dataset spans a relatively short period (2016-2019), we would suggest that this dataset is not long enough the assess the temporal variability in glacier dynamics. Also, when assessing the temporal variability of the whole region, one would likely lose again any valuable temporal signal when using a regional-wide coverage as individual glaciers will respond dynamically at different times, which would be obscured when collated or averaged out at the regional level. Therefore, a choice has to be made between focussing on the individual glacier and studying the temporal dynamics or focussing on a whole region with a composite dataset. This complicates the possibility to address the temporal variability, as our focus of the study is regional. Moreover, Dehecq et al. (2019)

discuss in the section 'Ice dynamics response to thickness change' a lag of a few years between the thickness change of a glacier and driving stress. This lag would impact the robustness of velocity change analyses over the time period covered by our velocity field.
While the dataset from ITS_LIVE is less resolved than this dataset, it would have been exciting to monitor changes in ice dynamics directly related to changes in lake heights (measured from altimetry for example) or lake area (from optical or SAR imagery), during the last 20 years. As a consequence, this would have been directly related to the second part of the paper that is modeling the influence of several parameters on glacier dynamics.

Response:
We appreciate the suggestion of the referee and agree that this would be interesting information. An investigation of lake heights and areas for a large sample over a multitemporal time frame would be something worthwhile for future studies. For now, we think that it is just as important to characterise the differences in dynamics relating to terminus type on a regional wide scale such as done in this manuscript, which has not been done before and paves the way for a multi-epoch approach such as the referee suggests. As such, we think that this suggestion is beyond the scope of the study, but needs to be addressed in the discussion, which will be done accordingly. Moreover, we would like to point out that within such a study set up, one would lose most valuable insights of the temporal variability when assessing lake-terminating glacier surface velocities on a regional wide scale. One would have to focus on specific glaciers again (Tsutaki et al., 2019), losing the regional scope of this paper.
However, we will shortly discuss the point in the discussion section.
This brings me to my final point: I found the two main sections of the paper a bit disconnected from each other's. The first one investigates the relation between several parameters (glacier orientation, area, debris cover) and velocity patterns, but there is little in common with the second section that really deal with what processes that are influencing lake-terminating glacier dynamics (where there is also no comparison with land terminating glaciers). This raises a number of interesting questions.

Response:
We agree that the two parts could have been better connected. We are convinced that the modelling section provides useful additional insights where the drivers of the observational section are investigated. A comparison to a land-terminating glacier has already been conducted by several other studies (e.g., Tsutaki et al., 2019) and would not necessarily add new insights to this field. Nevertheless, we will work on improving the revised manuscript to better connect the main sections together.
Indeed, is it right now possible to observe the influence in proglacial lake changes on changes glacier dynamics (ex: changes in lake level to be consistent with the modeling section)?

Response:
We agree with the author that this is an interesting line of investigation. Our study focusses on the diagnostic differences between land-and lake-terminating glaciers, and we therefore have no observational data on temporal changes of proglacial lakes. Although there is only limited observational evidence, in section 5.2 we mentioned observed changes of the frontal boundary condition followed up by changes in glacier dynamics: In line 576-578 we wrote: "Interestingly, an above average glacier acceleration was observed in 2006 after 5.6 m surface lake lowering as a mitigation measure in 2005(Xiao & Dai, 2011." Is there data available to see the formation of proglacial lakes, and how a change in the area of those lakes have influenced surface flow velocity?

Response:
There is certainly a robust relationship between lake area and lake depth, and consequently the area of contact between the lake and glacier terminus, which might alter the boundary conditions. However, our results indicate that the relation between lake area and surface flow velocity is ambiguous and likely lacks a direct causal relation. This should be investigated much more thoroughly, but the paucity of freely available lake bathymetry measurements across the Himalaya make this difficult to examine in detail.
Can we replicate these observations with the models presented in the last section?

Response:
Unfortunately, the accurate replication of specific examples of glacier behaviour would require a comprehensive dataset of lake bathymetry and temporal changes of the lake depth in the direct proximity of the glacier. With most of this data being yet unavailable we feel that this goal would be more suited for studies in the future.
With the large quantity of optical data to map lakes and already available velocity fields (ITS_LIVE, Golive, Dehecq et al., 2019), I think that it might be possible to assess.

Response:
We thank the referee for his suggestion. Please also see our introductory response. We feel that the quality of those datasets is not good enough to assess temporal variability on a regional scale at the terminus of glaciers < 5 km 2 . One then would have to focus on the larger glaciers, which would result in losing one of the main points of this paper, namely that special terminal glacier-lake dynamics are a regional wide phenomenon. With the improvement of satellite imagery, especially since 2016, a regional wide approach will become slowly but steadily possible.
Also, we would again like to point out that within such a study set up, one would lose most valuable insights of the temporal variability when assessing lake-terminating glacier surface velocities on a regional wide scale, as individual glaciers will respond dynamically at different times.
Finally, it is not clear to me why the author restricted their study area to the central Himalayas. It excludes an entire section of the Himalayas that is very dynamic, and with a lot more diversity in terms of glacier size, orientation, slope and velocity magnitude.

Response:
We appreciate the authors comment and will clarify these considerations properly in the revised manuscript. Please also see our introductory response. We restricted our analyses to these five regions in the central and eastern Himalayas for two reasons. Firstly, outside of the Himalayan regions of this study, the number of lake-terminating glaciers is currently limited (Nie et al., 2017). Consequently, these areas fall outside of our area of interest, as they hardly contribute to an increase in the sample size of lake-terminating glaciers. Secondly, we restricted our data to the extent of the dh/dt dataset of King et al. (2019) to be able to directly compare our velocity dataset to elevation change.

Response:
We thank the referee for this comment.
Please find specific comments below: L77. What do you mean by lake-driven changes in the velocity field ? is it related to the modeling part. Please present the objectives of this paper in the same order as it appears within the text.

Response:
We thank the referee to point this out and agree that this sentence might be confusing for the reader. In the objectives paragraph (around line 77), we write: "More specifically, we seek to investigate the attribution of lake-driven changes in the velocity field to dynamic thinning and investigate the role that debris cover plays on glacier-lake dynamics." In response to this, we will change this to simply 'changes in the velocity field' as one of the main aims is to investigate the prevalence of dynamic thinning, which is inherently lake-driven. We hope that with this, it is clear that the objective is entirely devoted to the remote-sensing part.
L102. Does that really make a difference ? Mean area of lake terminating glacier is 7 km2. How much glaciers are you adding up with these? Be more quantitative.

Response:
We thank the referee to point this out. In this study we only focus on glaciers with an area larger than 3 km2, compared to the 5 km2 previously utilised by Dehecq et al. (2019). Figure 4 shows that the quality of the dataset clearly improves with the use of 10m sentinel-2 imagery which allowed us to incorporate smaller glaciers into the dataset. Within the 3 to 5 km2 glacier area size group, we identified 14 lake-terminating glaciers, whereas the total lake-terminating dataset constitutes of 70 lake-terminating glaciers.
L105. What do you mean by "very low surface velocity" compare to what and where on the glacier ?

Response:
We agree with the referee that this has not been made clear in the manuscript. Generally, glaciers with a maximum velocity below 10 m/yr. We will make this clear in the text.
L108. Why do you restrict your study to the Central Himalayas ? By doing so, you are excluding all the glaciers in the Pamir-Karakoram, that are really diverse in terms of size, velocity magnitude, slope, debris coverage… Here you restrict yourself to glacier to mostly small size and slow-moving glaciers, which limits general conclusions that can be made.

Response:
Please see also our introductory response. We restricted our study area to the five regions in the central and eastern Himalayas for two reasons. Firstly, outside of the Himalayan regions of this study, the number of lake-terminating glaciers is limited. Consequently, these areas fall outside of our area of interest (the link between proglacial lakes and elevated glacier flow), as they do not contribute to an increase in the sample size of laketerminating glaciers. We also assembled our velocity dataset to match the extent of directly comparable dh/dt data of King et al. (2019) to be able to directly compare our velocity dataset to elevation change.

Response:
The two errors are essentially the same but are simply based on a different approach. Where Millan et al. (2019) use an average (for which we assume to be the median), we use the 95 th percentile reported by the quality report of Sentinel-2. Millan et al. (2019) finds a average error of 0.52 pixels, which roughly translates into our 95 th percentile error of 12m. We tend to follow the way of reporting the quality of Sentinel-2 by the official quality report of Sentinel-2 itself, as done in line 129.
L132. Do you mean removing the average offset calculated off glaciers? Which is mentioned in the post-processing?

Response:
Yes, this is a relatively common procedure with Sentinel-2 imagery.
L 138. Why do you restrict yourself to image pairs at 1-year interval ? Using all pairs of at least >1 month, would greatly improve your signal to noise ratio, which is really important when looking at small velocity numbers (<30 m/yr.) (cf Millan et al., 2019) Response: Please see our general response.

Response:
We processed all data between 2016 -2020 from November that was of adequate quality. If for example, data from 2016 cannot be used, the velocity data will shift to a central date more towards 2019. We thank the referee for the suggestion to use, 'central date' instead, and will do so in the revised version.
L 147-148. How is the stable ground area selected ? Is it random ? Hence including both mountain peaks (potential higher orthorectification errors) and valleys ? 300 km2 is a bit limited to see potential deviation across images and to study noise, specifically with a low number of image pairs. Furthermore, you discuss this also in section 3.1.5 right ? Please remove this part and discuss it later in the appropriate postprocessing section.

Response:
Please do also see the previous general response on the stable area selection. The selection of a stable area is not entirely random. One must be confident that most velocity estimates are zero. Very steep mountains do indeed cause potential higher orthorectification errors, and therefore such areas are generally omitted. The slopes of the glaciers do fall far out of this category, and we therefore argue that those potential higher orthorectification errors are not representative. We thank the referee for his suggestion, and we will remove this part in line 147 -148 to be clearer with our description.
L 151. Do you calculate gradient in the x and y direction ? Please specify.

Response:
We agree that this could have been specified in the text. We indeed calculate the gradient in the x and y direction. In the revised manuscript we will write (line 151): 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.
Section 3.1.4. This section is too technical and do not bring anything substantially new. The cross-correlation technique has already been widely documented in the literature. Hence, I would suggest to reduce this section and remove equations.

Response:
We thank the referee for his comment. We will condense the three processing sections and refer to existing literature in case of established methods.
Section 3.1.5. See previous comment.

Response:
See previous response.
Section 3.2. See previous comment on the error estimation. Please be consistent throughout the text between IQR, SEM, MAD….

Response:
We agree that the text has to be clearer about these intervals. Therefore, we will improve this in the revised version.
L 216. Is the use of such a large filter size limited compared to the width of the glaciers ?

Response:
The resolution of the velocity field is 80m and filtering therefore incorporates a larger surrounding area. We agree that for the very small glaciers, a filter window of 240 is still large. For this reason, indeed there is more weight on the velocity points on the very centre (see equation 4).
L 218. By how much does it increases the overall confidence ? Be more specific.

Response:
We weighted estimates with high confidence more: for example, if a neighbouring estimate next to (80m) the centreline has a much higher confidence, weighting this estimate more increases the confidence in the overall velocity estimate. From another perspective: If a centreline velocity estimate has a very low confidence, if might be worthwhile to look at the neighbouring estimate. Whether this increases the confidence drastically depends on each specific site, but it ensures that we retrieve data with the highest confidence possible.
L 256-261. How did you calculate the A value ? Does it vary spatially ? or do you take one value for the entire glacier/region ?

Response:
We feel that this is clearly written out in the text, with adequate referencing for more in depth information.
We wrote in line 256-261: "A is the temperature-dependent 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 depth-averaged ice temperature range of −10C to −2C at the ablation zone (Cuffey and Paterson, 2010), for which we follow Enderlin et al. (2013)." However, we can add more detailed information if the reviewer still thinks it would be beneficial.

Response:
We looked at the thickness of the larger, clean-ice, lake-terminating glaciers in the dataset of Farinotti et al. (2019) and used this as a rough indication of the ice-thickness in our modelling study.
L 278. How was the piezometric surface assessed ? Provide method and reference.

Response:
The piezometric surface is chosen in such a way that is starts at the lake's water table and slowly slopes upwards up-glaciers, so that it achieves a good fit with observed velocities, for which we follow Benn et al., 2007 (as referred to in the text).
L 280. What do you mean by up-glacier velocity ? Is it the ablation or accumulation zone ? Be more quantitative. Where is the speed value of 50 m/yr coming from ? Has it been taken from the measured velocity data ? Specify.

Response:
We thank the referee and agree this has to be more explicitly specified. Within our synthetic model set-up, we used 50 m/yr as a rough indication for a maximum velocity of a larger clean-ice lake-terminating glacier typically found at the end of the accumulation zone, which we based on our own velocity dataset. We will specify this in the revised version.
L 285. Ice thickness from the consensus estimate?

Response:
The referee is correct. We make this more clear by referring to H t in section 3.6.1 already, where we will write (line 275-277): "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)." L 288-289. Be more specific. What is the realistic range ?

Response:
We thank the referee for his comment. However, we think that with a clear reference of Watson et al., (2020), no range needs to be explicitly mentioned, as this only might cause confusion for the reader. In the discussion in section 5.2 we do mention these values when they are useful to be mentioned in the context of that discussion.
L 300. How do you change the ice thickness estimate ? Uniformly ? By how much ? How come you keep the maximum velocity at 50 m/yr. ? Do you still conserve mass ?

Response:
We will address these questions in the last paragraph of section 3.6.2 to make things clear. We do change the ice thickness uniformly by 50 m, as mentioned in the text. To keep the maximum velocity at 50 m/yr, we tuned the sliding parameter A s. We will mention this in the text as well.
L 305. Do you mean accuracy or precision ? I think you mean precision here.

Response:
The referee is right, but we realise that precision is also a bit misplaced here. We will rename section 4.1 to 'Algorithm Performance'.
Section 4.2. Please provide a figure illustrating the differences between each velocity dataset.

Response:
We thank the referee for his suggestion. We would like to mention that we provide a clear illustration of this in Figure 4.
Section 4.4. Considering the very large IQR, I do not find any significant differences between lake and land terminating glaciers that is above the noise. Please include Fig A3 in the main text. I think it provide more concluding evidence than Figure 6.

Response:
In the context of earlier comments about the confusion of uncertainty and the IQR, we see how these could be confused . As we now plan to merge the respective sections of the text which describe the IQR and uncertainty (3.2 and 3.5) we feel that the reader will be able to more clearly distinguish between the two in the amended manuscript.

Response:
We thank the referee for his suggestion. We think that it is clear from our previous responses that the IQR's shown in the figures are not error bars but show the variability of the background glacier population. We feel that it is important to show the spread of velocity values amongst glacier groups of different terminus type in Figure 4 and we would rather refrain from changing this figure. However, we will prepare this figure with IQR's and add it to the revised manuscript if we find that it does not obscure the original measures.
Section 4.6. Now this is a bit confusing, some much time have been dedicated to the comparison between lake and land terminating glaciers, but here we left off the land terminating ones. Why is that? Another aspect that would add even more value to the study would be to check out the influence of debris thickness on glacier velocity pattern and magnitude (check out Rounce et al., 2021 for the dataset).

Response:
This section is written in the context of both lake-and land-terminating glaciers, which we would argue is clearly illustrated in figure 8 (blue colours for lake-terminating glaciers throughout all the figures and red for land-terminating). To avoid any confusion, we will edit this section and make sure the differences are clear. We thank the referee for his suggestion on debris thickness, though we feel this would be outside the scope of this paper. Also, it is well documented that the majority of debris covered glacier area across the central Himalaya are stagnant or flowing below the level of detection of feature tracking algorithms. We do not feel that this needs re-emphasising.
L 504-509. The relation between the large-scale evolution of the Tibetan plateau and the formation of over deepening is not clear to me at this point. Please be more specific.

Response:
To make things clear, we will mention in the revised version that towards the Tibetan Plateau, the elevation generally slopes upwards (promoting overdeeping) (Royden et al., 2008).
L 510-514. Please provide a reference to a figure and section of your paper.

Response:
We thank the referee for his comment and will adopt this in the revised manuscript. Figure 11. Please provide a colorbar. Displaying the velocity on a log scale would enable to better observe the acceleration at the ice front which is not always clear (a, c, d). Would also be good to show for each glacier maps of surface lowering (Brun et al., 2017 for example).

Response:
We will provide colorbars in the revised version. If a log scale would indeed improve the referee's suggestion, we will do this as well. Directly comparing the glacier velocity might be interesting but can be problematic as you need a dataset with the same temporal coverage to get matching results for individual glaciers. For this reason, we think that this would be outside the scope of this paper.
L 517-520. It is not clear to me where this is going. Split this sentence into one or two difference sentences to make your argument clearer.

Response:
We thank the referee for his comment and will follow his suggestion to create more structure.
L 530. What about the lake temperature ? Could we potentially imagine that rising up the lake temperature would increase the melt at the front of the glacier and triggers an acceleration, as it is observed in Greenland and Antarctica ?

Response:
In the revised version, we will mention the potential importance of the lake-temperature as potential ultimate driver of changes in the boundary condition.