This manuscript describes a new method for characterizing the vertical heterogeneity of Greenland’s firn using VHF radar sounding and coincident laser altimetry measurements. Using full-wave electromagnetic modeling, the authors show that the offset in the apparent surface elevation measured by these two instruments increases significantly when one or more ice layers are present within the range resolution of the radar system. They demonstrate that this is due to constructive and destructive interference amongst layers that shifts the apparent location of the surface echo in the radar sounding data. The authors then map this offset across Greenland using OIB MCoRDS and ATM data collected between 2011-2019 and discuss the spatial and temporal trends in the context of regional climate. In particular, they show that their method captures both the inland expansion of ice layers following extreme melt in 2012 and the later replenishment of pore space in Southwest Greenland during a period of cooler summers from 2014-2018.

This is a technically complete and well-written paper that introduces an interesting new method for studying the heterogeneity of Greenland’s firn. It is exciting to see that that the apparent reduction in shallow ice content since 2012 holds at a larger scale than the initial firn core studies in Southwest Greenland, and this seems like a promising method for tracking trends in percolation zone extent and structure. I have some comments on framing and small details, but overall, I think this is an excellent piece of work.

Major Comments:

[1] The motivation for the study, and particularly why the use of VHF radar sounding together with altimetry, did not come across clearly to me in the introduction.  The authors list a slew of methods that are currently applied to study firn properties, but do not clearly establish what gap their method is then trying to fill. I think this is particularly important in Greenland because UHF Accumulation Radar data that can directly resolve the firn structure is available in all but one of the years that this study presents, so it is not immediately clear why one would want to use the VHF and laser altimeter combination instead. I think there may be some compelling reasons – for example, portability to Antarctica, where almost no Accumulation Radar data was collected – but the paper would be strengthened if the authors laid those motivations out in the introduction.

[2] Might surface roughness also modulate dz? I could image that for the radar, coherent integration of reflections from facets with different heights could set up constructive or destructive interference patterns (similar to layers at different depths) that might shift the surface height from the expected average of the facet heights. The anomalously high dz values in a few spots in the interior make me wonder if features like sastrugi might impact the retrievals at all.  

[3] What is the typical spread of values within the each 10km grid cell in Figure 7? Does the variability in dz provide any additional information, particularly in terms of the quality of the retrieval? To me, it seems like a potentially useful quantification of uncertainty in the boundaries between radar facies, which is somewhat lacking in the current version of the paper. I could also imagine that dz would be highly variable in the percolation zone, but maybe the dz anomalies in the interior would show a more consistent bias if they’re formed by some localized by consistent density anomaly.

Minor Comments:

Line 79: Note that the MCoRDS documentation gives an empirical windowing factor of 1.53 for a 20% Tukey window  on transmit/receive plus a frequency domain hanning window (https://data.cresis.ku.edu/data/rds/rds_readme.pdf). That seems to be more consistent with the current processing parameters for the MCoRDS data – Rodriguez-Morales et al. (2014) reports a processing chain that using a much more aggressive Blackman window on receive, which does not match up with what is in the current processing spreadsheets.

Lines 81 – 82: this seems not quite right for the various horizontal resolutions. The formula given is for the pulse-limited cross-track resolution, but needs to be noted as such, because the MCoRDS data have undergone synthetic aperture focusing and have a much smaller footprint in the along-track direction. The nominal along-track resolution should be ~25 m after focusing and incoherent summation, and typically the trace spacing is about half that (~14 m as correctly noted here). (Again, see the data documentation: https://data.cresis.ku.edu/data/rds/rds_readme.pdf.)

Line 115: Since you upsample by a factor of 10 before picking the surface, I would think that the vertical picking error should be +/- ~3 ns which would translate to height offsets of 0.32-0.27 m.

Line 117: Considering my comments on the along-track resolution, perhaps specify here what the total length of your moving average window actually is.

Figure 2: it’s hard to see the difference between black and blue in the figure, perhaps do something with more contrast like black and red?

Line 134: give the cutoff number for the maximum deviation from the straight line here for reproducibility.

Line 145: Are there significant differences in the ablation zone dz offset between different flights within the same season? If so, is only using the west coast for calibration adequate? For example, I might expect some flights in North Greenland not to cross this calibration area.

Line 148: Is the distribution generally close to gaussian in these 10 km grids cells? Is there much difference if you look at the median vs mean?

Line 161: Is a citation to a paper missing here? Courville & Perry (2021) seems to be a software, so it’s not clear what “Equation 15” is in that context. The software also appears to use csv input of the ideal transmit waveform, rather than an equation.

Line 195: does variability in the dry firn density matter? In the dry snow zone that can lead to a standard deviation of small scale density variations on the order ~35 kg/m^3 near the surface, which maps to plausible density contrasts of up to ~100 kgm^3 just from hoar formation and wind scour (Hörhold et al., 2011) (see the North Greenland Traverse cores high resolution measurements as an example). At least one previous statistical model of typical firn profiles did try to take this into account, but was also looking primarily at ice lenses impacts (Culberg & Schroeder, 2020).

Line 246: Okay, generally answers my question above, but given that later on your median offset in the percolation zone is only 2.2 m, it does make me wonder dry firn variability matters (and might be another explanation for the strange areas of large dz in the dry snow zone).

Section 3.1: not suggesting that you should actually do this for this particular paper, but it is possible to reprocess the OIB MCoRDS data with a higher vertical sampling rate. The posted data has been downsampled by at least a factor of 2 in the vertical and is not the native sampling rate of the radar. It might be worth mentioning in the discussion or future work for others considering using this method.

Figure 9: can the snow radar tell you anything about what is going on here with its higher resolution? For example, rule in or out very thin layers of higher density near the surface due to wind scour or other such processes?

Line 402: might be an appropriate place to also point to (Rennermalm et al., 2021)?

Line 465: I am not sure that this mapping method adds much to this beyond showing the limits of the percolation zone.

Line 473: Summit and NEEM both had evidence of thin ice crust/layer formation after 2012, so if it is indeed being buried slowly, it could maybe explain some of these anomalies (though doesn’t help with what’s already there in spring 2012) (Nghiem et al., 2012).

Line 493: Given that dz often has something to do with multiple peaks in the radar data, I expect this might work really well for capturing the same kinds of properties and would let you get around the need for coincident laser altimetry collection, which would be great.

References:

Culberg, R., & Schroeder, D. M. (2020). Strong Potential for the Detection of Refrozen Ice Layers in Greenland’s Firn By Airborne Radar Sounding. IGARSS 2020, 7033–7036.

Hörhold, M. W., Kipfstuhl, S., Wilhelms, F., Freitag, J., & Frenzel, A. (2011). The densification of layered polar firn. Journal of Geophysical Research: Earth Surface, 116(1), Article 1. https://doi.org/10.1029/2009JF001630

Nghiem, S. V., Hall, D. K., Mote, T. L., Tedesco, M., Albert, M. R., Keegan, K., Shuman, C. A., DiGirolamo, N. E., & Neumann, G. (2012). The extreme melt across the Greenland ice sheet in 2012. Geophysical Research Letters, 39(20), Article 20. https://doi.org/10.1029/2012GL053611

Rennermalm, Å. K., Hock, R., Covi, F., Xiao, J., Corti, G., Kingslake, J., Leidman, S. Z., Miège, C., Macferrin, M., Machguth, H., Osterberg, E., Kameda, T., & McConnell, J. R. (2021). Shallow firn cores 1989–2019 in southwest Greenland’s percolation zone reveal decreasing density and ice layer thickness after 2012. Journal of Glaciology, 1–12. https://doi.org/10.1017/jog.2021.102

Within, “Mapping the vertical heterogeneity of Greenland’s firn from 2011-2019 using airborne radar and laser altimetry” the authors devise a radar surface echo retrieval for shallow firn facies classification. The manuscript is well developed and utilizes modeling and in-situ data to improve radar signal interpretation. Much thought is given to the electromagnetic modeling set-up regarding the presence of snow cover on ablated ice and ice layers within the percolation zone of the firn. Modeling results provide a compelling argument for the divergence of radar surface echo waveforms and the misclassification of surface picks based solely on maximum amplitude. When the surface echo retrieval is applied ice sheet wide, interesting and confirmatory information is discovered, such as the replenishment of homogeneous firn within the percolation zone after the significant 2012 melt event. Results from this radar retrieval approximately confirm facies boundaries established by reanalysis model results, but I am left wanting a bit more analysis on how closely the radar and reanalysis resemble each other in a spatiotemporal context, as well as some discussion as to how this new radar information can inform or validate firn modeling. Those points aside, this manuscript is well-written, complete from an investigative standpoint, and deserving of publication. I recommend minor revisions which I have annotated in the attached .pdf.

Major Comments:

Throughout the analysis, I gained a more intimate understanding of how the VHF radar signal interacts with the near-surface, kudos to the authors for this. The efforts put forth in the modeling suggest that a given ice layer in the percolation zone must be buried greater than 8 m as to not interfere with the wave form surface echo. This is approximately 10 years’ worth of snow accumulation. However, empirically it was determined and the authors state that “between 2013 and 2014, areas of heterogeneous firn in the dry-snow zone reverted back to homogeneous firn.” This disagreement between the modelled information and radar retrieval is not reconciled within the manuscript discussion. Explanation for this phenomenon should be provided.

The correct identification of the “true” radar surface required an algorithm which upsamples the MCoRDS data to finer fast-time resolution. The author’s choice of resampling algorithm has a significant effect on the outcome of their results. If for example a piecewise linear upsampling was applied, the original MCoRDS signal and the upsampled signal would appear identical. Using a Fourier method as you have described, introduces "Gibb's Phenomenon" which can be seen as the oscillations at the tail end of the surface reflection signal. I suspect also that such phenomena are occurring to produce the troughs in the signal seen in Figure 6d & e (i.e., the Fourier series is struggling to represent the flat discontinuity of the signal between samples 3&4 of 6d). I appreciate the approach you have taken, and it is fortuitous that the upsampled signal reasonably recreates the modelled signal. However, acknowledging the origins of this "spurious" waveform bulge, be that an effect of the resampling algorithm or your choice in sample lags, should be considered and described more thoroughly.

Explanations of why dz values near the summit of Greenland are significantly higher than the surrounding dry-snow zone data retrievals remains unsubstantiated. Answers to these concerns are not supported in Scanlan et al. (2023), which pertain to higher frequency radar systems and show both consistently low retrieved density and high density within the same region of Greenland’s summit. Through tighter analysis, or a fleshed-out hypothesis, effort should be considered to reconcile the work of these authors.

Minor Comments: See the attached annotated .pdf