General Comment:

This manuscript by Seehaus and co-authors presents an application of Gaussian Process Regression (GPR) to interpolate ICESat-2 derived glacier surface elevation change (dh/dt) across two glacierized regions characterized by complex terrain and heterogeneous glacier dynamics: the Larsen-B embayment on the Antarctic Peninsula and central southern Svalbard.

The topic is timely and relevant for the cryosphere community. Satellite altimetry provides highly accurate but spatially sparse measurements, and developing robust approaches to reconstruct spatially continuous elevation change fields remains an important challenge. The use of a probabilistic interpolation framework, such as GPR, is well motivated, and the manuscript demonstrates that the approach can reproduce spatial patterns of glacier thinning and thickening, including signals associated with glacier surge activity.

Overall, the manuscript is clearly written and well structured. The use of two contrasting study regions strengthens the evaluation of the approach, and the sensitivity experiments provide useful insight into the effects of data gaps and terrain complexity. The figures effectively illustrate the spatial patterns of elevation change and interpolation uncertainty. However, several methodological aspects would benefit from clarification or further discussion.

Major Comments:

1. Justification of GRP Kernel

The manuscript adopts a Matérn 5/2 kernel with correlation lengths derived from semivariogram analysis. While this choice appears reasonable, the justification for the specific kernel and parameter values remains somewhat qualitative. It would strengthen the study's methodological rigor if the authors clarified whether kernel hyperparameters were optimized using marginal likelihood within the GP framework or fixed based on the semivariogram analysis. In particular, expanding the discussion of alternative kernel testing and kernel/correlation-length combinations would strengthen the justification for choosing this kernel.

2. Use of Glacier ID as a Predictor

Including the glacier ID as a predictor variable is an interesting approach to reduce spatial leakage between adjacent glaciers. While the examples in the manuscript demonstrate that this approach reduces leakage in practice, the inclusion of this variable also raises some questions. Glacier ID is a categorical variable, rather than a continuous physical variable, like the other predictors. Treating the ID as a numerical feature in the covariance kernel may introduce artificial relationships between glaciers. While this may not be a major concern, it would be helpful for the authors to discuss the statistical implications of this choice and whether alternative approaches were considered.

3. Spatial Correlation Length for Smaller Glaciers

The spatial correlation length used in the model (~11 km) appears relatively large compared to the characteristic scale of many glaciers in the study regions, particularly in Svalbard, where surge-related signals and terminus dynamics may occur over much shorter distances. A brief discussion of how this correlation length affects the model's ability to capture sharp spatial gradients in dh/dt would be useful. In particular, the authors may wish to comment on whether the chosen correlation length could lead to smoothing of localized elevation-change signals.

4. Artificial Gap Experiments

The validation strategy combines comparisons with TanDEM-X elevation change fields and artificial data-gap experiments. These tests are helpful and provide useful insight into the interpolation performance. The artificial gaps used in the sensitivity analysis are circular.  whereas real ICESat-2 sampling gaps are typically elongated and aligned with satellite ground tracks. The authors may wish to briefly discuss whether this difference could influence the interpretation of the gap-filling experiments.

Minor Comments:

1. The detailed discussion of the pre-processing of the ICESat-2 data is well-structured and informative. However, it is difficult to determine exactly how many ICESat-2 observations were used in each region. The addition of a summary sentence with this information (or a table in section 4.1) would aid the readability of the manuscript.
2. When comparing GPR results with the ATL15 product, the authors should note that ATL15 is a coarser-resolution product that applies spatial smoothing for regional elevation change estimates. 
3. Several typographic errors are present throughout the manuscript. Some are noted here, but the manuscript would benefit from detailed proofreading.
   • Figure 1 (inconsistent capitalization in the legends)
   • L. 233 & 390 (study sides rather than study sites)
   • 4.2 Heading (Gaussing rather than Gaussian)
   • L. 396 (acounts rather than accounts)

In this study, the authors have presented work investigating the use of Gaussian Process Regression (GPR) to derive spatially continuous elevation changes over glaciers using ICESat-2 observations. This is a highly relevant and interesting topic, and overall the methods seem to be sound. The manuscript is generally well-written and understandable. I have a few general comments, and a number of specific comments/questions for the authors.

General comments

I think the presentation of the pre-processing steps (§4.1) could be re-arranged or re-worded, to make it clear that you are in effect using three different sets of elevation changes for three different "experiments":

1. Long-term trends using observation pairs separated by at least 20±1 (12±1 for Larsen-B) cycles, so long as the number of cycles separating the two observations is a multiple of 4 - this is the main portion of the results presented and analyzed.
2. Annual height changes, using observation pairs separated by 4±1 cycles.
3. Single cycle (short-term) height changes.

As written, it was not immediately clear that you weren't taking the results of each of these different pairings and using them all in the same GPR regression. Hopefully I have understood this correctly, that you have done the regression using three different ways of slicing the data.

My second general comment regards a different sort of interpolation. In this, you compare your GPR regression to two other forms of interpolation, Inverse Distance Weighting (IDW) and nearest neighbor (NN). Did you also consider comparing to the hypsometric interpolation often used for altimetry surveys on glaciers (e.g., Arendt et al., 2006; Johnson et al., 2013; Nilsson et al., 2015)? I think that might be a useful exercise for demonstrating the utility/improvement of your method against a common approach in glaciology.

I would recommend checking the manuscript carefully for consistency in spelling. For example, Mannerfelt vs Mannerfeld; Moršnevbreen (the version used on NPI maps, for example) vs. Morsjnevbreen, Morsnevbreen, or Morsinevbreen; Matérn vs Matern; Doktorbreen vs Dotorbreen; Liestølbreen vs Liestolbreen; Vallåkrabreen vs. Vallakrabreen, etc. Similarly, double-check the in-text citations and the references list, as there are some inconsistencies.

Specific Comments

l. 11: are the ± values here uncertainty or standard deviation?

l. 30: say here what the repetition rate (or "cycle", later on) is for ICESat-2.

l. 31-35: could combine these statements into a single statement

l. 44: did you explore using slope and aspect as predictors?

l. 61: "Spatially distributed uncertainty estimates are provided" - could also mention that this is a feature of kriging

Fig. 1: label for Morsjnevbreen (Moršnevbreen) could be moved to avoid overlap with Paulabreen

l. 120: what is the vertical datum for the Svalbard DEM provided by NPI?

l. 125: what year(s) do the outlines from Silva et al. correspond to?

l. 126-129: was this step not needed for the Svalbard outlines?

l. 130: previously, you stated December 2024 as your cut-off date - I assume this just means that the data are available until Feb. 2025?

l. 133-134: This threshold (2500 m) would remove obvious clouds, but potentially still include lower-lying clouds. Is there a further filtering step applied? Additionally, what is the "fit quality flag"? And, why only use points within glaciers?

l. 140-145: why the difference in preferred cycle separation for each site?

l. 150-167: maybe state the number of points/observations for each filtering step?

l. 160: is this the absolute deviation from the median of the bin?

l. 171-172: did you use any others? You previously mention slope @ line 153.

l. 177-182: I understand the point here is that you would like to identify areas impacted by surge by comparing the velocity anomaly to "stable" conditions, but couldn't you also use the contemporary velocity over the time separation of your elevation observations as a predictor for the GPR? I think that Hurkmans et al. (2014) used a mean velocity over several years because that was what was available at the time, but I believe that suitable velocity datasets exist to look at this over a smaller time period (e.g., citations at l. 329-330). Hurkmans et al. 2014 is not included in your references.

l. 185: how was this interpolation done?

l. 190: suggest using space rather than period to indicate thousands (e.g., 11 000)

l. 197-200: what implication does this have for the smaller glaciers included in the study area? 

l. 223: wouldn't it be better to do this comparison/correction with the original point measurements, rather than the GPR outputs?

l. 236: I assume this is meant to be Nilsson et al. 2015, which is included in the references list, rather than 2016?

l. 245: why nearest-neighbor and not bilinear?

Fig. 4c), elsewhere: I don't think the units of variance should be m/a

l. 279: I think this suggests that your extreme change rate filter (described at l. 148) should be higher, if observed thinning reaches -30 m/a?

l. 308-311: indicate which glaciers match which pattern/style of propagation

l. 381-382: this also seems to be the case for Hektoria Glacier, but at much lower elevations - why might that be?

Fig. 7: could also include a panel showing a histogram of the offsets/differences, with mean/median, nmad, etc.

l. 395: "offsets peak" rather than "offset peak"

l. 399: remove "revealed" (repeated)

l. 402: volume rather than volumen

l. 415-417: how does this compare to the difference between the GPR output and the observed dh/dt values from the "main" result? That is, if you look at the offset between the GPR output and the original gridded ICESat-2 observations, is the difference observed here similar or greater than the difference between the GPR output and the original observations?

l. 432: I don't see these plots in the Supplement - only the plots for the annual and single cycle GPR outputs.

Fig. 9: suggest using different symbols as well as color here

l. 449: remove figure caption from cross-reference

References