General comments:

This manuscript presents a new algorithm, ‘Watta’, for automatically extracting supraglacial lake depth estimates using ICESat-2 geolocated photon heights (ATL03) which are then used to validate empirically-derived lake depths from Landsat 8, Sentinel-2, Planet Labs Skysat and PlanetScope imagery. The authors test the algorithm performance on 46 supraglacial lakes near Jakobshavn Glacier in West Greenland during an intense melt season (2019). Finally, they use this stacked dataset in combination with Operation IceBridge imagery to track volume, drainage mechanisms and ice cover evolution of two individual lakes in this region.

Supraglacial lakes form in the ablation zones of Greenland and Antarctica during the summer melt season and can impact ice sheet dynamics, making lake detection and depth retrieval important. However, lake volumes have been difficult to quantify due to a lack of in-situ measurements and uncertainties associated with image-based methods. This manuscript builds upon other recent studies by applying a novel method for lake depth extraction, which is the first application of high-resolution Planet Labs satellite imagery to calculate supraglacial lake depths in combination with other imagery sources and ICESat-2 heights. It also provides useful insight into lake dynamics and ice cover evolution, which to date have been limited by the comparatively coarse resolution of publicly-available satellite datasets (Sentinel, Landsat).

Therefore, it is my view that the findings are of broad interest to the cryospheric community and represent a promising step forward for studying supraglacial hydrology and dynamics. I look forward to seeing further development of this method and its applications elsewhere on the Greenland and Antarctic ice sheets, particularly on floating ice shelves.

In general, this is a well-written manuscript and most of my comments are relatively minor. I would like to see in places some additional detail around the discussion of lake depth retrieval methods (see specific comments attached).

Lastly, The Cryosphere’s data policy states that “Authors are required to provide a statement on how their underlying research data can be accessed. This must be placed as the section "Data availability" at the end of the manuscript.” Although the authors state at the end of the manuscript that the Matlab code will be converted and shared publicly, I would like to see this section added including statements of how Landsat, Sentinel and Planet imagery can be accessed.

Once the authors address these points and my comments below, I can therefore recommend that this manuscript is suitable for publication in The Cryosphere.

This manuscript presents a new algorithm (called “Watta”) to determine depths of supraglacial lakes in ICESat-2 satellite laser altimetry data. It adds to a growing body of literature that is using ICESat-2’s novel capability to detect shallow bathymetry features in non-turbid water. The authors demonstrate the algorithm’s ability to reliably extract meltwater depths from ATL03 photon cloud data using fifty examples of melt lakes in Greenland during the exceptionally strong melt season in summer 2019. For each of the lakes, imagery was available from Landsat 8, Sentinel-2, PlanetScope or SkySat within a few days of the ICESat-2 acquisition date. Under the assumption of exponential decay of reflectance at water depth, the authors then use the combination of imagery and Watta-derived lake depths to fit empirical relationships between water depth and reflectance in imagery, for each available combination of lake location, acquisition time and sensor type. These empirical relationships are then used to estimate the depths of each entire melt lake from all available imagery scenes, thus providing volume estimates. The authors further use the resulting depth estimates as well as the underlying altimetry and imagery data over five named lakes in combination with Operation IceBridge data to explore in detail the effects of using different bands in multispectral imagery (lake Zadie), the differences between ICESat-2 strong and weak beams (lake Cecily), the processes underlying lake filling and drainage (lakes Ayse and Julian), and ice cover evolution (lake Niels).

Greenland’s increasingly negative net mass balance is strongly driven by surface melt. Surface lakes can significantly increase surface albedo leading to furter melt, and their drainage into the subglacial environment can impact ice dynamics. It is thus crucial to better understand surface hydrological processes and their influence on mass balance to constrain future sea-level rise projections. This manuscript demonstrates that with the new ICESat-2 laser altimeter it is now possible to reliably measure supraglacial lake depths from space in an automated fashion. While this has been the topic of some other recent publications, to the best of my knowledge this manuscript is the first to utilize ICESat-2 data for satellite-derived bathymetry estimates over supraglacial lakes. The authors are also the first to combine bathymetry results from ICESat-2 with commercially available ultra high-resolution imagery from Planet Labs to demonstrate how our understanding of surface hydrology can be improved by exploiting a diverse range of satellite remote sensing products resolving complex supraglacial hydrological systems at process-relevant scales. 

I therefore believe that the methods introduced in this manuscript present a significant improvement upon the current state of research on supraglacial hydrology. Building upon the methods presented here should be a promising way forward, towards ice sheet wide generation of meltwater volumes from observations. Thus, the manuscript is an important contribution to the current literature, and its results should be of interest to the broader polar science community. 

Major concerns:

1. Lack of evidence for sub-surface ice layers: The manuscript often refers to ice cover on melt lakes as well as sub-surface ice layers without showing evidence that such features are indeed present. In fact, most of the features that are described as “sub-surface ice layer” very well match the appearance of “second return” artefacts that result from the ATLAS sensor being saturated due to specular returns from flat surfaces, and many photons returning to the instrument during its dead-time without being detected. This issue is also briefly described in Martino et al. (2020), and further detailed in the “Specular Returns” section of the known issues document for the ATL03 product. I would highly recommend that the authors consult these two sources and decide whether they are still convinced that these second returns are signals from sub-surface ice layers rather than just artefacts in the data. If so, I would like to see convincing evidence for the claimed widespread existence of such sub-surface ice layers in a revised manuscript. 
2. Insufficient information to make methods repeatable: One of the main objectives of this manuscript is to present the new Watta algorithm. Yet, the algorithm is not described in much detail, and the information provided is certainly not enough to replicate the results in this study. I think the easiest way to fix this would be to archive the already existing matlab code on a platform such as Zenodo (which was used by the authors to archive the lake depth data set). If the authors do not want to share their code, I think that they should provide some more detailed pseudo code (including a list of parameters and chosen values), or a more detailed text description at the very least. Furthermore, from the information provided in the manuscript, it does not seem possible for readers to check the underlying data themselves. This means that it is very hard to verify any claims made based on “visual inspection”. No locations are given for the lakes under consideration, and ICESat-2 data is only referred to by its track number (not the spot, so this could refer to any of the six beams) and plotted against along-track distance with the zero point seemingly arbitrarily chosen. At the very least, I think that the authors need to provide latitude and longitude coordinates for each of the lakes considered in this study and to specify the ICESat-2 spot (GT1L, GT1R, GT2L, GT2R, GT3L, GT3R) for each lake section used with Watta. This information could be included in table S1 in the supplement. Since each of of the lakes has an associated ICESat-2 overpass, this information could be automatically extracted directly from the corresponding ATL03 or ATL06 data by, for example, using the median latitude for each of the segments shown and then querying for the corresponding longitude along the ground track. (Other information useful to readers (spot, beam type, acquisition time, etc.) could easily be printed out at the same time.) 

Minor concerns / suggestions: 

1. I think some similar ICESat-2 shallow bathymetry literature should be cited. There are a few articles out there with similar methodology, just applied to satellite-derived bathymetry outside the polar regions. Examples would be Albright and Craig (2020) or Thomas et al. (2020).
2. It is unclear to me how the numbered lakes relate to the named lakes. Some sort of explanation for choosing to refer to some lakes by numbers and some by names should be included. A map with the locations of all lakes would be great. 
3. It is not clear to me from the text how geolocation and co-registration relate to each other here. Can you explain how matching the ATL03 point cloud with the GIMP-2 DEM to reduce square error will improve the co-registration between ICESat-2 and Landsat 8 data? To my understanding, GIMP-2 DEM elevations are mostly derived from WorldView stereo imagery, and if there is a significant difference in acquisition time between the image underlying the DEM and the Landsat 8 / ICESat-2 lake observations then surface topography could have changed significantly in the meantime due to ice flow or surface processes. With a geolocation accuracy of roughly 5 meters for both ICESat-2 and Landsat 8, I would assume that simply mapping both datasets to the same CRS would give better results than the intermediate use of DEM elevations. Admittedly, this might be me not fully understanding co-registration, yet it would be nice to be provided with some more detail/explanation, or to see evidence that this intermediate step using DEM elevations actually improves the method in a meaningful way. 
4. The results shown on the figures could be made somewhat more accessible to readers: I would suggest to plot any ATL03 data with latitude on the horizontal axis while also including a scale bar for along-track distance. This, along with the information about which ICESat-2 track and spot is shown on which date, would already be enough for readers to figure out where to find all the underlying data. Plotting ICESat-2 ground tracks on top of imagery or image-based depth estimated wherever applicable would help readers with visual verification of some of the claims made in the text. A graticule on some maps/imagery would help as well. 

Line-by-line comments: 

• Line 14: From just reading the abstract it is unclear what is meant by “corrected” depth. 
• Line 17: Landsat 8 
• Line 18: You are stating 46 lakes, but I am counting 45 lakes in table S1, plus the five named lakes in Figure 1 for a total of 50 lakes?
• Line 22: please spell out CAMBOT the first time you use it (Continuous Airborne Mapping By Optical Translator)
• Line 40: You mention “both ice sheets” here, yet surface melt in Antarctica has not been discussed. While this paper only uses data over Greenland, I think it would be beneficial to briefly mention surface melt in Antarctica and how it is believed to be connected to ice shelf disintegration via hydrofracture. 
• Line 54: Since you specify for Sentinel-2, can you specify what the “higher spatial and temporal resolution” is?
• Line 56: This sounds very wordy. Could simply say “ICEsat-2 now makes it possible to replace…”
• Line 60: typo (bathymmetry → bathymetry)
• Line 64: I am sceptical about the presence of ice layers under the water surface. See above in the major concerns section. 
• Line 64: Maybe here I would specify that by “the native resolution of the ATL03 photon cloud” you mean the 0.7m laser pulse frequency in along-track distance
• Line 65: typo (wen → when)
• Line 85 / Figure 1: Can you add latitude and longitude labels, or preferably a graticule in the right panel? Please also specify in the caption that RGT = ICESat-2 “Reference Ground Track”, not “repeat ground track”. The ground tracks that should be repeated (in the polar regions) are the six spots GT1L, GT1R, GT2L, GT2R, GT3L and GT3R for each numbered track. The RGT should be the point directly at the nadir, so unless ATLAS is pointing off-nadir it should fall right between GT2L and GT2R. 
• Line 92: This paper is largely about ICESat-2 so you might want to spell it out: “Ice, Cloud, and Land Elevation Satellite” and possibly ATLAS = Advanced Topographic Laser Altimeter System
• Line 98: It is unclear to me what you mean by “using x signal photons per shot”. Are you referring to the expected number of signal photons that ATLAS will detect per pulse? Are these values over land ice? Is there a citation for these values? 
• Line 101: typo (MacGruder → Magruder)
• Line 114: TOA has not been defined before → top of atmosphere
• Line 121: It is unclear to me here what the role of the DEM is in geolocation. 
• Line 124: “each area was approximately on average” makes no sense?
• Line 137: There are only 45 lakes in the supplemental table?
• Line 150: I think it might be good to point out somewhere that if the empirical estimates are “time, location and sensor specific”, then your method is currently limited to producing valid depth estimates for imagery scenes that overlap with an ICESat-2 overpass over a melt lake within that scene and a three-day window. This is a limitation that the physical models don’t have. 
• Line 166: How are “outliers” detected?
• Line 191: What you describe here sounds exactly like artefacts in the data that come from ATLAS’s dead-time when the sensor is oversaturated by a specular return. If you really believe that this is sub-surface ice in some cases, then I would need to see evidence for that to be convinced. (see major concerns section)
• Line 201: missing full stop after “lake edges”
• Lines 205-208: It is not clear from the text how matching the ATL03 point cloud with the GIMP-2 DEM to reduce square error will improve the co-registration between ICESat-2 and Landsat 8 data. GIMP-2 DEM elevations are mostly derived from WorldView stereo imagery, and if there is a significant difference in acquisition time between the image underlying the DEM and the Landsat 8 / ICESat-2 lake observations then surface topography could have changed significantly in the meantime due to ice flow or surface processes. With a geolocation accuracy of roughly 5 meters for both ICESat-2 and Landsat 8, I would assume that simply mapping both datasets to the same CRS would give a better co-registration than the intermediate use of DEM elevations. If this is not the case, it would be nice to see some sort of proof that this intermediate step using DEM elevations actually improves coregistration. 
• Line 208: Do you mean a margin of 0.2 degrees in latitude and/or longitude? 
• Line 217: two commas after NDWI_ice, missing full stop before “Boundaries”. 
• Line 224: typo (MacGruder → Magruder)
• LIne 225-226: “a line 6m in each direction perpendicular to the ICESat-2 beam” seems like a rather confusing way to describe a circle of 6 m radius around the location of the photon. 
• Line 230 / Figure 3: Please plot the ground track of the ATL03 data shown in the top left panel on top of the depth estimates shown in the bottom left panel. Please spell out “Elevation” and “Along-track distance” in the top left panel. Also, why is along-track distance going from roughly -50 m to 800 m? I think the ICESat-2 convention is that along-track distance is measured from the last equator crossing? It would probably be more helpful for the reader if elevation was plotted against latitude, with a scale bar indicating along-track distance. 
• Line 245: spell out “2” → two
• Line 259-260: If performance evaluation is done by “visual inspection”, it would be nice if the reader could also get to see a few examples of imagery with precise ICESat-2 ground tracks plotted on top, for their own visual inspection.
• Line 260: correlation coefficient between what? NDWI and Watta-derived depth?
• Line 264-265: “reference ground track (RGT) 1222, Lake 3 in Fig. S4”:  Should be referring to Fig. S3. 
• Line 265-266: “the presence of subsurface ice did not always preclude the presence of a strong bottom return” →  This suggests to me that it’s even more likely that this "subsurface" ice layer might not exist, and that it's actually the sensor saturation and dead-time effect. (see major concerns section)
• Line 266: “(e.g. Lake 7, RGT 1169, Fig. S4)”: Should be referring to Fig. S3. Also, I don’t really see anything indicative of subsurface ice in Lake 7, RGT 1169, Fig. S3. 
• Line 279: It sounds like you are using the R^2 for performance evaluation of the empirical model, but this would mean to evaluate the model on the data that was used to generate the model in the first place. So it should be made clear that the R^2 cannot be considered a performance metric for a model across an entire lake basin, and rather that it merely indicates how well you were able to fit the empirical model to the data along the given ICESat-2 ground track. However, the underlying model is rather simple and based on physics, so overfitting is probably not much of an issue here.
• Line 292: typo (there’re were → there were)
• Line 293-294: “future users would be able to select bands or combinations [...] that provide the greatest fidelity to ICESat-2 based observations”: I know what you mean by that, but the way it is phrased it sounds like a bulletproof recipe for overfitting the data...
• Line 304: Technically GT3L describes a "spot", not a beam. Two beams (one strong, one weak) will alternate in pointing at that particular spot, switching off whenever ICESat-2 performs a yaw flip. 
• Line 309: typo (lake → lakes) 
• Line 313 / Figure 5: “Sentinel-2 (l,m) and Planet SkySat (n,o)”: should be “Sentinel-2 (k,l) and Planet SkySat (m,n)”. Also, what does the red box in panel c indicate? 
• Line 328 / Figure 6: Can you show the ICESat-2 ground track on the right panels? It is very hard to see what is going on without that information. Also, it is pretty clear from context what the abbreviations Sent/LSat/SSat/PS/R/G mean here, but at least somewhere you should specify this for clarity. Also please try to stay consistent across all figures. I have noticed images with labels “Sentinel-2”, “Sentinel”, “Sent” and “S” across the figures in the paper, and they all refer to the same thing. 
• Line 344: You want to refer to Supplemental Fig. S4 here, not S5
• Line 356 / Figure 7: Can you label lake NIels and lake Julian on the left panel?
• Line 407: Ice motion should not be adjusted for in geolocation?
• Lines 409-419: What’s shown in cyan in fig 9d does not look like ice cover to me. Also none of the satellite imagery seems to show the presence of ice cover. Can you corroborate your claims about ice layers? From looking at the figure, I would guess those are specular returns from water surfaces. (see also major concerns section)
• Line 425: This is not the reference ground track for track 727. This must be GT1L (based on looking at the data myself), which is roughly 3.3 km offset from the RGT! The big stars used to show the locations A-D very much cover the actual features, which makes it hard to see any of the things discussed in the text. Can you plot the actual precise ground track 727 GT1L for this overpass on panel e as a (very) fine line, and indicate locations A-D with arrows pointing at the features without covering them. 
• Lines 458-463: This paragraph about Antarctica does not fit into the conclusions section. The information about ice shelf stability considerations, etc. would fit nicely into the introduction/background information about surface melt, where “both ice sheets” are already mentioned. Then, in the conclusion section you could just briefly mention that Watta could be used in Antarctica as well.
• Line 469-470: The goal of implementing Watta in an open-source framework is commendable and would certainly be beneficial to the scientific community. Yet, it would also be helpful to publish the already existing matlab code along with the manuscript. (also likely the easiest way to address my major concern about methods repeatability)

References:

• Albright, Andrea, and Craig Glennie. "Nearshore Bathymetry From Fusion of Sentinel-2 and ICESat-2 Observations." IEEE Geoscience and Remote Sensing Letters (2020). https://doi.org/10.1109/LGRS.2020.2987778
• Martino, Anthony, Christopher T. Field, and Luis Ramos-Izquierdo. "ICESat-2/ATLAS Instrument Linear System Impulse Response." (2020). https://doi.org/10.1002/essoar.10504651.1
• Thomas, N., Avi Putri Pertiwi, Dimosthenis Traganos, David Lagomasino, Dimitris Poursanidis, S. Moreno, and L. Fatoyinbo. "Space‐Borne Cloud‐Native Satellite‐Derived Bathymetry (SDB) Models Using ICESat‐2 and Sentinel‐2." Geophysical Research Letters (2020). https://doi.org/10.1029/2020GL092170