Review of Winstrup et al., “PRODEM: Annual summer DEMs (2019-present) of the marginal areas of the Greenland Ice Sheet”, submitted to ESSD

by Romain Hugonnet, University of Washington

Statement of expertise

I have expertise in the aspects of this paper associated with:

• Remote sensing of surface elevation data for glaciology,
• Laser altimetry with ICESat-2,
• Geostatistical methods used for data fusion, in particular variography and kriging,
• Uncertainty analysis in a spatiotemporal context for the predicted elevations. 

Would be nice to complement with another reviewer being specifically expert in:

• Radar altimetry with CryoSat-2 (though I do have some).

General comment

In their very well-written manuscript, the authors present a rigorously derived dataset of annual (2019-ongoing) gridded surface elevation at 500 m resolution for the Greenland Ice Sheet, predicted from spatiotemporally sparse ICESat-2 and CryoSat-2 elevations.

The study is a truly a great piece of work, and I commend the authors on their attention to detail at all steps of their analysis to produce a robust product for the community. My core expertise is mathematics and DEMs, in particular spatial statistics and uncertainty analysis, and this study is a rare sight in our field in that aspect, and I think the way forward so that other scientists and users can know how reliable a data product really is.

I still have several comments for the authors, which I think could improve some parts of their analysis substantially. Some might be a bit of work, but are important and should be OK given the author’s knowledge in geostatistics, and others are optional.

Sorry if I cite my own work quite a bit in this review, there are not that many other studies looking at these specific problems in detail (geostats + uncertainty quantification + DEMs), and I hope it will be helpful to the authors!

Main comment 1: Refine the modelling of the heteroscedasticity of predicted elevation

Or “variability in error”. This is something the authors have addressed quite a bit in the text and their estimates, in particular in Section 4.1 by constraining the variability in elevation error (CS2 - CS2) with surface roughness (for CS2; Figure 2b), and using the variability of 20 m IS2 points in the 250 m area.

However, Figure 8 shows that this modelling of the variability in error (including subsequent propagation with kriging) is not completely satisfactory:

• Too large predicted errors for the center of distribution,
• Yet does not capture the full range of variability (tails are quite significantly larger than a normal distribution).

Here, I recommend the authors several pists that should improve their reported error.

First, some small but important changes to avoid over-estimating the center distribution:

1. Consider using a Q-Q plot in Figure 8 (or add one), which will represent the tail differences in a clearer fashion,

2. Avoid over-estimating the CS2 measurement error by dividing by square-root of 2 as the error source is combined twice in the cross-validation exercise (e.g., Equation 8, Hugonnet et al. (2022)).

3. Avoid over-estimating the IS2 measurement by taking a simple spread: separate segments will have largely uncorrelated errors, so you need to compute a standard error, i.e. divide by sqrt(N) with N the non-overlapping ICESat-2 segments.

Then, to improve their modelled variability in error:

4. Have a consistent spatial variability dependency with terrain roughness for both IS2 and CS2 as it should be universal (e.g., Hugonnet et al. (2022) for slope + curv), you could apply the same scheme to both CS2 and IS2 (it sounds you did something similar for IS2 on L260-261? These statements are a bit unclear, consider explaining in more details, or adding a Figure if it tells something about the error!).

5. Better constrain the variability of CS2/IS2 error (which might explain Fig. 8 entirely).

I was happy seeing Equation 1, and then disappointed that the sources were not actually individually refined. I do understand that separating the sources is not feasible with the available data, but it is still possible to separate the variability in error mixed within your single error proxy (here, your cross-over differences).

The authors have 40,000 CS2 cross-over and 60,000 IS2 ones at disposal, a good sample to understand different variabilities. For the temporal variability: the authors could bin these differences with the temporal lag to your reconciled summer date represented by your predicted elevation (middle of your June 1st to September 30^th period? This “middle reported date” absolutely needs to be reported in the paper). This should allow to constrain an average variability with the time lag (e.g., Hugonnet et al. (2021), Extended Data Fig. 5ab, which also includies spatial correlations). However, I suspect this variability will likely depends on thinning at the surface at each location, not only on the time lag. You would have to use an estimate of long-term dhdt at the surface (doesn’t need to be super resolved, e.g., Smith et al., 2020, or your own product for multiple years) to combine this into a 2-D (time lag + dhdt) binning. This might explain most of your temporal (and general) error variability.

Main comment 2: Report spatial correlation of errors for your data product

As pointed out by the authors in the short discussion (L617-620), even though your estimates have a reasonable error, it might not be random at all because of the temporal differences in acquisition.

While it would be quite a task to correct these temporal errors in detail (feasible by extending your kriging framework with a dimension of time or seasonality, for instance, or with a more simple models of your multi-annual changes during the summer season), the authors can still model the spatial structure of these errors. And this information would be quite important for downstream applications.

At the local scale, those are systematic errors due to common temporal sampling. But, as long as there is no bias in the time of data acquisition, those will be centered on 0 on average (which the authors validate with their comparison to IS2 in a LOOCV, Table 1). Which means that, at the larger scale of Greenland, they can be interpreted as structured errors with spatial correlations due to the temporal sampling discrepancies (Hugonnet et al. (2021)).

Here, the authors should at least report an “average” variogram that represents the spatial correlation in errors their products due to temporal sampling (this has be done on the LOOCV differences you compute in Section 7.1).

Even better would be to have a varying variogram depending on a map of “time lag to closest elevation used” or similar, that would reliably represent the errors in both space and time!

You could then validate your variogram representing your errors using Block LOOCV, and studying that the spread matches that predicted by your variogram in space (following classical error propagation with a spatial correlation term, e.g. Hugonnet et al. (2022), Equation 18).

However, the “swath” nature of your input data, and hence the temporal sampling, is less than ideal to capture errors in temporal sampling with omnidirectional variograms. It's still an important added value, and I don’t see an easy solution to this directional problem, except trying to come up with a better scheme to correct these temporal errors seasonally to your middle summer date, using your own multi-annual data (as stated at first).

Main comment 3: (Optional) Regionally-varying kriging: could also decompose and standardize your elevation field variance sources to explain the variability at the source

I must have been quite the effort to setup this regionally-variable scheme with GSTools!

First, while the Cressie estimator is reasonably robust to outliers, it is still fairly sensitive for what we get in elevation data (ArcticDEM outliers). You could consider using Dowd’s estimator (Dowd, 1984) that actually matches your NMAD in terms of how it is derived. This estimator is implemented in SciKit-GStat (Mälicke et al., 2022) which has functions to export models to GSTools.

Second, while the regionally-varying kriging is well-implement and definitely valid, it might be over-fitting your data locally by using that regionalized variogram subsample. This would result in an over/under-smoothing of the final interpolated field compared to reality, because of the way the scheme is implemented. Another way to avoid the non-stationarity would be to describe the nature of your spatial covariance and its variability (summed variance components with a certain correlation range, and their individual variability). The variability in error with long correlation range would for instance likely linked to time lag and dhdt (Main comments 1 and 2). Then, using standardization + decomposition (as in done a bit in Hugonnet et al. (2022) with instrument resolution and noise, but was less important there), you could transform your data to reach second-order stationarity, and apply kriging.

I am mostly presenting this option for the authors to realize that there are ways to avoid a regionally-varying variogram implementation. But I think that, given the complexity of this scenario (varying time sampling that is uncorrected), the current approach is completely valid and it would be hard to implement a decomposed + standardized approach.

Main comment 4: Justify the choice of 500 m for the final product

I don’t think it is explained anywhere in the manuscript why the authors chose 500 m for their study, it would be good to add a paragraph on this!

Additionally: ArcticDEM mosaic v4 is released since last month, with significant improvements. If it is not too much work, the authors could update their pipeline with this latest version.

Main comment 5: Drop the nugget for the kriging

Using a nugget in a variogram represents a very specific type of discontinuous spatial pattern at the boundary of the measured location. This is really only relevant to point observations that actually belong to a small spatial ensemble (such as gold nuggets), and is not relevant to inherently continuous elevation data on a grid.



You can actually see on Fig 4. the misfit of the black and red variograms at short spatial lags: they should not have a nugget at all, and it would show if more short lags were sampled.

To improve your spatial lag sampling in the variograms, consider using a non-uniform binning. Those can be long to compute and, for pairwise sampling efficiency, in SciKit-GStat, there are different types of metric spaces just for this, one based on the sampling scheme presented in Hugonnet et al. (2022), Fig. S13. Those are also wrapped more conveniently in xDEM for raster data (xDEM contributors, 2022).

Line-by-line comments:

• Title: Having “present” in the title for the period is probably not a good idea long-term.
• L40-42: This statement unlikely holds true. Typical random errors in modeled ice thickness are +/- 100 meters near flux gates with important spatial correlations as well as systematic errors (Kochtitzky et al., 2022; 2023). While predicted surface elevation errors are typically less than 10 meters (Hugonnet et al., 2021). Even measured (not just modelled) ice thicknesses are usually less precise and with larger correlations in errors than measured surface elevation. The authors could simply make the case that we need more annual data to represent the surface elevation of the ice sheet.
• L48-52: Same remark.
• L81 and onwards: Need to add spaces between numerical values and unit.
• L190: Could directly introduce the Normalized MAD or “NMAD” directly, which has been around for while in elevation data analysis (e.g., Höhle and Höhle, 2009).
• L319: Need to specify: "second-order" stationary!

References:

Mirko Mälicke, Romain Hugonnet, Helge David Schneider, Sebastian Müller, Egil Möller, & Johan Van de Wauw. (2022). mmaelicke/scikit-gstat: Version 1.0 (v1.0.0). Zenodo. https://doi.org/10.5281/zenodo.5970098

Hugonnet, R., Brun, F., Berthier, E., Dehecq, A., Mannerfelt, E. S., Eckert, N., & Farinotti, D. (2022). Uncertainty Analysis of Digital Elevation Models by Spatial Inference From Stable Terrain.

IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 15, 6456–6472. Hugonnet, R., McNabb, R., Berthier, E., Menounos, B., Nuth, C., Girod, L., Farinotti, D., Huss, M., Dussaillant, I., Brun, F., & Kääb, A. (2021). Accelerated global glacier mass loss in the early twenty-first century. Nature, 592(7856), 726–731.

Kochtitzky, W., Copland, L., King, M., Hugonnet, R., Jiskoot, H., Morlighem, M., Millan, R., Khan, S. A., & Noël, B. (2023). Closing Greenland’s mass balance: Frontal ablation of every Greenlandic glacier from 2000 to 2020. Geophysical Research Letters, 50(17). https://doi.org/10.1029/2023gl104095

Kochtitzky, W., Copland, L., Van Wychen, W., Hugonnet, R., Hock, R., Dowdeswell, J. A., Benham, T., Strozzi, T., Glazovsky, A., Lavrentiev, I., Rounce, D. R., Millan, R., Cook, A., Dalton, A., Jiskoot, H., Cooley, J., Jania, J., & Navarro, F. (2022). The unquantified mass loss of Northern Hemisphere marine-terminating glaciers from 2000–2020. Nature Communications, 13(1), 1–10.

Höhle, J., & Höhle, M. (2009). Accuracy assessment of digital elevation models by means of robust statistical methods. ISPRS Journal of Photogrammetry and Remote Sensing: Official Publication of the International Society for Photogrammetry and Remote Sensing , 64(4), 398–406.

xdem contributors. (2021). xdem (v0.0.2). Zenodo. https://doi.org/10.5281/zenodo.4809698

This paper introduces a set of four 500-meter resolution annual Digital Elevation Models (DEMs) representing the Greenland ice sheet marginal zone during the summers of 2019 to 2022, referred to as PRODEMs. Covering a 50km wide band from the ice edge, these DEMs meticulously encompass all outlet glaciers of the Greenland ice sheet. The integration of DEMs derived from ICESat-2 and CryoSat-2 constitutes a commendable initiative, providing an additional valuable source of topographic information for Greenland. While the paper is well-crafted overall, there are certain statements that require precision to ensure accuracy and clarity. It is crucial to address these points before considering the paper for publication.

General comments:

• The authors' decision to generate the DEM specifically for the summer and focus on the marginal zone of the ice sheet raises pertinent questions. What rationale underlies the choice of season, and what significance does the selected ice sheet margin hold in this context? Furthermore, a clarification regarding the specific function of the DEM under these conditions would enhance the paper's comprehensibility.
• The authors employ two types of altimeter data for DEM generation, a method requiring careful consideration of their consistency. While the authors assert that the disparity between the two datasets is negligible, it is crucial to acknowledge that applying these datasets without consistency correction may influence elevation estimations. Particularly noteworthy is the comparison between IS2 and CS2; the latter exhibits uneven spatial distribution and coarse resolution, introducing potential biases in the generated DEM. My concern centers on the relatively lower measurement accuracy and disparate spatial distribution of CS2, aspects that warrant careful attention.
• Regarding data coverage, the authors applied distinct interpolation radii in different regions. To enhance reader understanding, it would be beneficial to include information on data coverage, such as footprint numbers for CS2 and IS2 in each grid. This insight can shed light on the grids effectively covered by the altimeters, offering rationale for the choice of a 500 m resolution. It's noteworthy that the ICESat-2 DEM by Fan et al. (2022), constrained to 10 footprints at 500m resolution, covers only approximately 30% of the entire ice sheet, and a discussion on this limitation could be insightful for the readers.
• The authors should elucidate the rationale behind selecting the Kriging interpolation method. Given the steeper terrain at the ice sheet's edge, it is essential to address whether Kriging interpolation can effectively fulfill the task of estimating elevation in such rugged terrains. Clarification on the suitability of the chosen interpolation method, especially in challenging topographical conditions, would contribute to a more comprehensive understanding for the readers.
• There are reservations regarding the chosen 500m resolution. Satellite altimeters typically yield fewer effective observations at the ice sheet's edge (especially for CS2), and opting for higher resolution may result in a reduced number of usable observation grids. This reduction can potentially impact the interpolation performance. The question arises: Can a DEM heavily reliant on interpolation accurately portray the actual elevations, particularly in the steeper terrains at the ice sheet margins? Additionally, in regions with lower latitudes, where the gaps between satellite orbits are more substantial, concerns linger about ensuring effective spatial interpolation. Addressing these considerations would fortify the paper's discussion on resolution choice and its implications.
• Notably, the entire DEM is subjected to Kriging interpolation, implying a potentially limited actual coverage area by the altimeter. Therefore, it is imperative for the authors to furnish the proportion of space covered by observational data. Elevation derived from Kriging interpolation might deviate from the altimeter-observed elevation, necessitating consideration of this difference, as the potential errors induced by interpolation are a critical aspect. Particularly, given the uneven spatial distribution of CS2, which may introduce additional uncertainty to Kriging interpolation, it becomes crucial to address how the authors accounted for the impact of observation gaps on the results. Discussing these considerations would enhance the paper's transparency regarding the intricacies of the interpolation process and associated uncertainties.
• Addressing DEM uncertainty requires a more detailed exposition. While the spatial and temporal uncertainty is elaborated upon, clarification is needed on how the authors define the instrument and geolocation uncertainties mentioned in equation (1). It would be beneficial to delineate the specific contributions of each factor to the overall uncertainty and provide a proportionate breakdown. Moreover, in discussing the uncertainty of CS2/IS2 elevations, the term 'crossovers' might be prone to misunderstanding. Consider utilizing track-based crossover analysis for data validation, ensuring a clearer understanding of whether the differences approach zero. This adjustment would enhance precision in evaluating the uncertainty associated with CS2/IS2 elevation comparisons.
• The evaluation of the proposed DEM solely through cross-validation may be deemed insufficient. It is recommended that the authors endeavor to gather actual measurements, such as airborne datasets, for Greenland elevation. This would facilitate an independent evaluation of the generated DEM, offering accuracy metrics under various terrain conditions. Such an approach would provide more robust and precise evaluation information, enhancing the credibility of the study's findings.
• In comparing with other DEMs, presenting a distribution map alone may lack depth. Additional comprehensive comparisons, such as showcasing the spatial distribution or numerical values of elevation differences, and establishing relationships between these differences and slope and roughness, are warranted. To enhance the display of the DEM, the authors might consider incorporating a slope map or shaded relief map. These additions would provide a more nuanced evaluation of the DEM data and contribute to a thorough understanding of its effectiveness.
• Consideration should be given to the potential benefits of updating the DEM with the new ArcticDEM mosaic v4, as it has the potential to significantly enhance DEM accuracy and reduce differences in dh (as observed in Figure 10). Such an update could serve as a viable solution to the issues identified and contribute to an overall improvement in the quality of the DEM.

Specific comments:

(p: page, l: line)

Section Introduction: The current paragraph lacks substantial content. To enhance its depth, it is recommended to incorporate the significance of ice sheet elevation monitoring. Elaborating on the importance and implications of monitoring ice sheet elevation would provide readers with a clearer understanding of the broader context and relevance of the study.

p2l44, Supplementing the existing DEM datasets and providing commentary on them is advised. This additional information will not only enhance the comprehensiveness of the study but also offer valuable insights and context for readers evaluating the significance and limitations of the proposed DEM datasets.

p2l53, CS2 is equipped with a Ku-band radar altimeter, which is highly sensitive to moisture changes. Notably, surface melting can introduce significant interference to the echo signal, posing challenges in obtaining accurate elevations during such periods. While the author's emphasis is on capturing surface elevations during summer, it is worth considering the possibility of melting snow, even in areas designated as snow-free. Clarification on the impact of any residual melting snow on elevation data would enhance the understanding of potential influencing factors.

P2l69, what distinguishes the Baseline-D and E datasets, and is there potential for this difference to impact the consistency of the CS2 data? Clarifying the distinctions between these datasets and addressing their potential influence on CS2 data consistency would contribute to a more comprehensive understanding of the data sources employed in the study.

P2l74, in contrast to IS2, CS2 exhibits a relatively limited dataset with irregular coverage, further reduced through quality control measures. It is imperative to elucidate the specific spatial distribution of CS2 data after quality control. Any potential uneven distribution in the data may raise concerns about its impact on subsequent analyses and differences. Addressing these considerations would provide a clearer understanding of the reliability and potential biases associated with the quality-controlled CS2 dataset.

P2l75, Significant disparities exist between LRM data and SARln data. The authors should expound on their considerations regarding the consistency between these two datasets. Providing insights into the strategies employed to address or account for these differences would contribute to a more thorough understanding of the data integration process and potential impacts on the study's outcomes.

P3l83, A higher sampling density often correlates with improved elevation simulation results. In the context mentioned, it would be beneficial to clarify the purpose of avoiding oversampling and elaborate on the potential consequences associated with it. Understanding the reasoning behind this choice and its potential impacts would enhance transparency regarding the sampling strategy employed in the study.

P3l88, However, such downscaling might diminish the inherent advantages of IS2. Additionally, in the case of CS2 data, it is essential to elucidate how the authors addressed the potential impact of observation resolution. Clarifying the strategies employed to mitigate the effects of observation resolution on CS2 data would contribute to a more comprehensive understanding of the data processing methodology.

P3l93, it is evident that IS2 data plays a more central role in DEM generation. Consequently, it raises the question: What specific role does CS2 serve in this context? Considering the non-uniform distribution of CS2, there arises a concern about the potential introduction of uncertainty into elevation estimations. Addressing the distinct contributions and potential uncertainties associated with CS2 data would enhance clarity regarding its significance in the overall DEM generation process.

P4l118, the impact of CS2 signal penetration depth in snow is not explicitly addressed in this section. Providing insights into how the authors considered and accounted for this factor would enhance the completeness of the discussion.

P4l124, the reason for choosing a 500m resolution is not explicitly stated. Including the rationale behind this resolution choice would provide clarity and context for readers seeking to understand the decision-making process.

P5l146, the uncertainty associated with Kriging interpolation is not specified. Offering details on the uncertainties inherent in the Kriging interpolation method would contribute to a more thorough understanding of the potential limitations in the DEM generation process.

P7l236, the stated measurement precision applies specifically to the interior areas of Antarctica. It's crucial to note that for steep terrains, the accuracy is expected to be lower. Emphasizing this distinction would provide a more nuanced understanding of the precision associated with different topographical features.

Section 4.1, The calculation method for DEM uncertainty, particularly how the authors integrated data from the two altimeter datasets, is not clearly presented in this section. Including details on the calculation method would enhance transparency and facilitate a more comprehensive assessment of the uncertainty analysis.

P8l279, Notably, there appears to be a relatively large elevation difference between CS2 and IS2. It is pertinent to inquire whether the authors have considered any correction methods to address or minimize this discrepancy. A discussion on potential strategies or considerations in correcting the observed elevation differences would contribute to a more thorough evaluation of the data integration process.

P9l292, it would be beneficial to expound on the specific steps involved in DETREND. Additionally, clarification on the basis for employing linear interpolation and whether it introduces any uncertainties would be useful.

P9l298, the term "ELEVATION ANOMALIES" requires clarification, and providing details on the proportion of outliers would offer a more precise evaluation of the data characteristics.

Section 5.3, Explicitly mentioning the software used for Kriging interpolation and listing relevant parameter settings would provide essential information for readers seeking to replicate or understand the interpolation methodology.

Section7.1, I strongly recommend incorporating airborne elevation data into the DEM evaluation process. While cross-validation results are valuable, utilizing airborne data for evaluation would provide an additional and crucial perspective on accuracy. This approach ensures a more comprehensive assessment and strengthens the overall validity of the DEM.