My review is attached.

Review of “Mapping the aerodynamic roughness of the Greenland ice sheet surface using ICESat-2: Evaluation over the K-transect” by van Tiggelen et al.

General comments

This is strong manuscript that demonstrates impressive proficiency with many different sources of data (AWS, UAV, ICESat-2, modeling). The methods are generally well-described. The results section is very interesting and the development of spatially extensive aerodynamic roughness lengths for the K-Transect from ICESat-2 is commendable.

However, I do recommend some revisions.In its current form, the introduction is poor. Some of the terminology is vague, references are lacking and the overall research is poorly motivated. I encourage the authors to revise it thoroughly and have provided some ideas for doing so below.

While it is useful to know that the commonly used method for deriving z0m from ICESat-2 (i.e. the standard deviation of ATL03 heights) tends to overestimate z0m, the new measure is slightly unsatisfactory if it underestimates z0m by a factor two. Without looking at the data, it is difficult to discern why. It could be due to the slightly arbitrary choice of filtering (qlow = 1 and qhigh = 2) to remove photons above and below the median. It could due to the choice gaussian covariance function, window size or assumed wavelength. Given that this is one of the first papers to investigate roughness lengths using ICESat-2 and availability of ground-truth data, it would be useful if the authors could develop a more unbiased method. I would encourage the authors to perform some sensitivity tests with these choices to see if they would reduce bias in their ICESat-2 z0m products.

Specific comments

L16: Please consider capitalizing “ice sheet”. It’s the Amazon River, the Tibetan Plateau and should be the Greenland Ice Sheet. Indeed the Nature paper that you cite (Shepherd et al., 2020) has it this way.

L19: If you define an acronym, it is usually appropriate to use it here and elsewhere (e.g. L50, L146).

L18-21: Please provide some references for these two statements. A lot of work has been done on these topics and it is negligent to overlook it.

L20: “can be” is poor rationale for studying something. Please revise with something stronger, perhaps relative to radiative heat fluxes.

L22-26: Again, please provide references to backup these statements. A paragraph in the introduction without any references indicates that the research is poorly motivated or that the authors have a complete lack of respect for previous research on this topic. Please revise.

L37: What do you mean by “confined accessible areas”? Please provide some examples.

L39: Consider replacing “unmanned” with an ungendered term.

L40: What do you mean by “limited”. Please be more specific.

L41: I am not aware of a satellite altimeter that maps the surface roughness of entire glaciers. The ground sampling distance is not small enough. This sentence also makes it sound like UAVs are completely unnecessary. Please revise and be more specific.

L42-44: This sentence about sea ice does not fit here in a paragraph about glaciers and ice sheets, please move somewhere else.

L99-100: Presumably Fig. 1b could be referenced here?

L145: missing an “of” between transect and AWS.

L226: I thought you just said that this approach did not require interpolation to 1 m profile?

L252-259: This text would be more useful in the introduction.

L260-274: Some more references to Fig. 6 in this paragraph would be useful to the reader.

L285: Please clarify what is mean by “satellite backscatter”. I presume you are referring to a satellite radar instrument since ICESat-2 does not measure backscatter.

L288: Fig. 6? This figure does not show an elevation profile.

Consider swapping Sections 4.1 to 4.2 and Fig. 5 and Fig. 6. I think it would make more logical to move from small to large scale.

L396-397: It would be useful to briefly state again why Lettau (1969) is not recommended. Some people may only read the abstract and conclusions.

L399-402: I’m not sure I follow this logic. How do you know that ICESat-2 does not capture snow sastrugi or ice hummocks > 1000 m a.s.l. when your UAV surveys are constrained to < 600 m a.s.l.?

L405: It’s a bit of stretch to say ICESat-2 cannot map z0m above 1000 m when this study presents no UAV surveys above > 1000 m.

Figure 1. Most of panel (a) is irrelevant, given that data from S9 are not used in this study. It makes it difficult to see how the ICESat-2 tracks intersect the UAV survey grids (A and B). Please consider removing the picture of S9 and providing a zoomed version of the UAV survey grids around the margins of the ice sheet. In the caption please specify if these are the ICESat-2 reference ground tracks or from an actual ICESat-2 beam (e.g. 1r).

Figure 2: What is the rationale for these wind directions? Prevailing wind direction from AWS? Please clarify.

Figure 6: There is no reason for such large x and y axis limits on this figure which makes it difficult to determine the correspondence between the SEC and VPEC dots and modeled lines. Please provide a zoomed version of this figure.

General Comments

----------------

This manuscript addresses retrieval of surface roughness length on ice sheets using ICESat-2 data profiles. Empirically based retrieval of surface roughness length from satellite observations is an enormously important task; the parameter modulates energy fluxes between the atmosphere and the cryosphere, changes in both space and time, and is poorly known. Current methods of retrieval generalize single point measurements to large expanses of the ice sheet; not only do we not have spatially resolved estimates of this parameter, we lack comprehensive understanding of the variance, range, and uncertainty of the parameter. Thus, the present work is extremely timely and important to the community at large. That said, there are several shortcomings with this work; the applicability is limited to a narrow range of surface types and elevations that form a minority of the ice sheet area, the measurements themselves have large uncertainties and are resolved for specific wind directions that do not match prevailing katabatic patterns, and the validation strategy and data are marginally matched to the task. This study is undeniably useful as it forms a basis for future work to build on; the problem under study is a hard task, and incremental progress should be recognized and iterated with new, separate publications that extend to the rest of the ice sheet. In short, this work is worthy of publication following revisions-- there are specific tasks and issues that should be addressed in the revision, and other issues that can be deferred as 'out of scope' and addressed in distinct publications rather than in the current work.

Specific Comments

-----------------

The spatial and temporal mismatch of the validation data from S5 is the largest issue with the current submission. The S5 site is the only location that ties together all three components of data used in this study-- structure from motion DEMs, ICESat-2 tracks, and empirical measurements of surface roughness length by in situ measurements. Although both the 'A' and 'B' structure for motion boxes are bigger, and overlap with ICESat-2 tracks, the lack of weather station data forces S5 to be the primary validation loci, despite the smaller than desired area/fetch coverage. The large temporal gap (September DEM and in situ measurements, March ICESat-2 data) isn't ideal, and needs to be fixed (preferably), or explained and justified in greater detail. ICESat-2 was operating in September of 2019-- why isn't there coincident data provided? Looking at the track crossings, it appears that September 24 was cloudy to the point of signal loss, but this isn't explained...signal from September 12th is stronger and appears to cover and cross over S5, so why wasn't this data used? What is the justification and the trade space between small spatial mismatches vs large temporal mismatches? Why March? Having data coincident in both time and space for the validation is ideal, and a strong case with reasoning and justification needs to be made as to why September data was not used and/or was not tractable for use.

While having structure for motion, in situ, and ICESat-2 data all be coincident is ideal, the second best approach is paired validation: coincident UAV and in situ data to validate the method, followed by coincident ICESat-2 and in situ data to validate the scaling to the 1D profile. This is especially appealing since the current work already has separate pairing that is discussed with ICESat-2 and UAV data in boxes A and B in addition to pairing of UAV and in situ data at S5; the only pairing not present is between ICESat-2 and in situ tower measurements. Even if data doesn't simultaneously overlap for all three data sets, finding an overlap between ICESat-2 and the S5 station provides the needed coverage for a compelling validation strategy. I'm unclear on if this is possible, or perhaps why it isn't possible since my expertise is more with ICESat-2 than with tower measurements. My impression is that the most of the weather stations such as S5 collect data in dense time series that are continuous save for maintenance or power outages. Is there a reason why there's not coincidence between S5 and ICESat-2, such as lack of co-occurrence that matches the prevailing wind direction? Some of this is addressed explicitly around line 320, but I'm still skeptical; if wind measurements are occurring in dense (i.e., multi-hertz) time series, brief changes from the prevailing wind direction should still occur, even if they are not sustained on the time scale of hours or days.

Given how hydrologically active the area is, I was surprised by the lack of discussion or mention of water such as lakes and the impact on the retrieval process. Figure 3 shows a profile that appears to have multiple surfaces between 100 and 150 meters that may be water ponding. Around line 195 or 200 would be an appropriate place to discuss this, given the discarding of photons below the median which will help with water surfaces.

The algorithm only uses a single profile; probably fine for this paper, but difficulties in determining the width parameter (or whether a given obstacle meets the width threshold) can likely be improved by examining both of the pairs to assess obstacle persistence in the across track. Similarly, I expect that cross track estimation of surface roughness is feasible at track cross over points given the double beam crossing of the pairs, which would help with the katabatic prevailing wind alignment issues...

The primary roughness retrieval algorithm (i.e., thresholding photons according to confidence class, median filtering, then interpolating with k-nearest neighbor and kriging in constructing profile obstacles) seems reasonably considered, and robust. The alternative formulation which uses the standard deviation of photon spread from the de-trended ATL03 product is less compelling; there is no accounting that I can tell for difference in signal strength or atmospheric conditions; photons for the standard deviation calculation are weighted equally regardless of the per photon quality/confidence flags. While this residual measure is designed primarily to provide an upper bound of the estimated surface roughness, rather than a 'best estimate' of surface roughness, additional corrections and filtering of what photons to consider would improve the metric.

Around lines 375 to 380 there is a discussion of how the ATL03 surface roughness retrieval breaks down at higher elevations...however it is unclear if this is due to sensor tuning for the specific algorithm, or theoretical limits for conceptual mental model that relies on obstacle formalize instead of skin friction parameterization. High resolution DEMs at the higher elevation bands would likely indicate if the formalism adopted in section 2.2 can be scaled to sastrugi in principle, or if the conceptual framework itself is no longer appropriate given the dominance of skin friction related to inherent snow and firn properties. In other words, lack of high elevation UAV DEM coverage such as exisits at the lower S5 or sites 'A' and 'B' does not allow the reader to infer if the Bulk Drag Partitioning method itself is not suited to retrievals at these heights, or if algorithmic implementation as presently tuned is not suited. (Note, this is issue is also raised by the other two referees). Determining this does not require coincident data; simulation of the method on a generic surface with ice hummocks at ~0.6m scale would provide enough context to discuss the issue in the text. 

Technical Corrections:

----------------------

Line 1: Curious if the authors mean 'latent heat' explicitly when they reference moisture in this context

Line 25: 'form drag' is more formally defined later in equation 3; I would include the parameter name (tau_r) here as well to aid readers.

Line 35: This is vague-- are there no physically based drag models that are capable of simulating surface roughness length from an elevation profile period? Or just no models that are used for simulating the exchange between cryospheric surfaces and the atmosphere?

Figure 1: The 'large black box' referenced isn't clear, and is easily mistaken as a graticule; use a different bright color (orange, yellow, red) with higher saturation to highlight the area better.

Line 95: While the fetch footprint is variable, discussion of the range or a small table of the normal values as a function of boundary-layer height / friction velocity would be helpful.

Figure 3: Standard convention is that 'noise' photons are labeled as grey, and 'signal' photons are labeled as black. I realize that color choice here is carried forward with consistency for the figures that follow, so that yellow and grey lines reference the same process/data in figures 4, 7, and 8, but I think that these figures need to have the colors switched as well. The data/noise convention for photon signal/noise is similar in strength to mapping conventions that expect water labeled as blue, or data orientation to point North. If there are concerns for black data dots being too dark in Figure 4 and obscuring the profile, using blue data points is a possible work around (dark blue for signal, cyan or light blue for noise)...but in general, convention and expectation is that the lighter saturation or value assigned in point plots is for noise, and darker points are signal.

Line 120: Can the cut off wavelength be variable? This question is probably related to my comment on lines 375-380 that I discussed at the end of 'specific comments'

Line 145: "...140 km transect AWS..." --> "...140 km transect of AWS..."

Line 200: Some discussion/mention of wet surfaces and standing water is warranted here

Line 225: This should modified to account for signal strength; since the ratio of noise to surface photons varies with signal strength, the standard deviation will be biased between high and low signal strength acquisitions over the same surface. This is true for the background count rate as well, which varies seasonally and between night and day conditions.

Line 240: I'm unclear on exactly what is meant be residual photons here, and if they are weighted or binned by the confidence flags assigned in ATL03.

Line 255: First use of 'L69' I think...a sentence somewhere defining the acronym convention for the various methods would help

Line 270-275: I don't know if it's true to say that there's no relationship between C_d and either H or lambda...especially when the next sentence links increased C_d values with increases in H. I'd change this to say that the there is a weak relationship.

Line 288: I think figure 5 is meant here, not figure 6

Line 313: Eddy covariance measurements are available outside of September, I assume? Even if they aren't available in March specifically, having a date range of the measurement record would be helpful, instead of just the date that the data was pulled for this study.

Figure 8: This would be appropriate to split into a.) and b.) panels. I'm skeptical of the pink 'perfect fit' line; the eddy covariance measurements that the C_d values are inferred from have some spread or standard deviation, so I expect that modeling those uncertainties would produce a flatter line or a bounding envelope.

Line 385: This claim should be tempered a bit. Sure, there is lower contribution to runoff at higher elevations, but the increased surface area relative to the margin means that modeling the high elevation roughness is crucial for understanding and modeling the overall energy exchange between the cryosphere and atmosphere. Also, under changing climate scenarios, run off contribution will increase for high elevations.

Line 430: Capitalization is inconsistent between equation A1 and the following line where the parameters are defined.