Thank you very much for your helpful comments. We agree that the terminology used in the original manuscript was not always consistent and may have caused confusion. Following your suggestion, we have revised the text to consistently distinguish between gravity disturbances, isostatically corrected gravity disturbances, terrain effects, and non-terrain gravity disturbances throughout the manuscript.

To improve clarity, we have also expanded the description of the gravity-field decomposition workflow and added a dedicated subsection describing the isostatic correction. This section now explains the Airy-compensated Moho model, the density assumptions, and the calculation of the resulting isostatically corrected gravity disturbance. We further revised Figures 2 and 3, their captions, and the corresponding equations to ensure that the terminology is used consistently and that the relationship between the different gravity-field components is explicit.

We believe these changes substantially improve the readability of the manuscript and reduce potential ambiguity regarding the gravity quantities used in the inversion workflow.



Line-by-line responses:

- For the isostatic anomaly calculation, you cite Liebsch 2020, but I could only find some code associated with this. Can you add 1 or 2 sentences in the text or appendix describing this method? Including which bed topography model you used, and which density values. 

We added a separate subsection (Section 2.4) to better explain the isostatic correction, including the underlying method, topographic model, and density assumptions.

- On line 174 you state you will use gravity disturbances instead of FAA. I would change all further references of FAA to gravity disturbance, such as in Fig 2 and its caption, and the equations in sections 3.1. Also I would initially state that the terrain effects are the same thing as the Bouguer correction, and that non-terrain effects are the same thing as the  Bouguer disturbance, but once you state this, stick with one naming convention. 

We added a clarification on the naming of terrain effects and non-terrain disturbances and removed "FAA" from figures, captions, and equations.

- I’ve added lots of suggestions for the text of Steps 1-6. I think this is part is especially important for readers to follow, and I struggled to understand exactly what you did. Hopefully my suggestions can help clarify the workflow you used and help compare/contrast it to the methods of other bathymetry inversion studies to put your methodological advance in context.

- Clarification of gravity anomalies: 

From the code you provided, this seems to be your workflow:

for clarity, bouguer correction = terrain effects, and bouguer disturbance == non-terrain effects

1) subtract isostatic gravity effect from FAA to get an isostatically-corrected FAA

2) calculate the bouguer correction of ice, water, and rock from BedMachine and subtract this from the isostatically-corrected FAA to get the bouguer disturbance

3) calculate a smoothed (regional) bouguer disturbance with a RBF of the data only near constraints

4) Subtract the regional bouguer disturbance from the full bouguer disturbance to get a residual bouguer disturbance.

- up to this point, this should be equivalent to An et al. 2015 DC shift or Tankersley et al. 2025 constraint point regionalization

We added these references in the Methods section to make this conceptual similarity clearer.

5) cut out the residual bouguer disturbance away from radar data, and replace it with SGS to get an ensemble of residual bouguer disturbances.

6) add the regional bouguer disturbance back to each of the residual bouguer disturbance ensemble to get an ensemble of full bouguer disturbances

7) subtract this ensemble of full bouguer disturbances from the isostatically-corrected FAA to get an ensemble of bouguer corrections

I think this differs slightly from the presentation in Figure 2, which omits the important step of calculating the regional bouguer disturbance, and makes it seem like the SGS is applied directly to the full bouguer disturbance, instead of to the residual component. 

I my interpretation of the code / your methodology is correct, I would suggest somewhere you show all the following maps of gravity data, either in figure 3 or a supplementary figure: FAA (Fig 3a right?), moho effect, isostatically-corrected FAA, the mean terrain effect,  the mean Bouguer disturbance, and the regional and residual components of the Bouguer disturbance. I think Fig 3. would be more informative if it showed 1) the isostatically-corrected FAA (instead of the FAA), 2) the mean target terrain effect (current Fig 3b), and the mean gravity misfit (forward effect of Bedmachine - mean target terrain effect). I would then make sure to cross-reference which subplots of Fig 2 correspond to which subplots of Fig 3 / the new supplementary figure.

Thank you for this suggestion. We have substantially revised both Figures 2 and 3 following your recommendation. We have also updated Section 3.1 and added explicit cross-references in the text to the relevant panels in Figures 2 and 3 to clarify the relationship between the gravity-field decomposition and the inversion inputs. We believe these changes make the workflow and interpretation considerably clearer.

- Corrections for ice and water: I read through your code trying to understand how you correct for the gravity effect of ice and water, specifically your python class "PrismGen2". It seems you make 6 sets of prisms, a positive and negative set each for ice, water, and rock, but it was pretty hard to follow. Could you clarify a bit in the text, or supplement, about how you discretized and what density values you used for each of these 6 prism sets, and whether all 6 prism sets were computed at each iteration (or just the bed prisms?). 

In the updated Section, we now explicitly describe how the model domain is discretized and clarify the densities used.

Specific comments:

Line 44: You might add (Fig. 1c) after "These gaps".

Added.

Line 98: Maybe use ground survey instead of ground line to avoid confusion with grounding line.

Changed.

Line 117: "above the ice surface" suggests these altitudes are with respect to the ice surface. Are they? Or is it 1.5-5 km ellipsoidal heights?



We have clarified this in the manuscript as follows:

“Surveys were conducted at varying aircraft heights, typically ranging from approximately 1.5 km to ~4 km ellipsoidal heights.”

Line 131: Did you block-reduce the line data before fitting the equivalent sources? If so, what cell size?

The equivalent sources were fitted directly to the original line data without prior block reduction. The upward-continued field was then predicted onto the BedMachine reference grid in 1 and 2 km resolution for local and regional inversion, respectively.

Line 137: Was the block median here used to go from the 2km resolution gravity grid to the 1km? Why not just predict the equivalent sources at both the 1 and 2km grids? Or does block-reduction essentially give the same result?

Thank you for pointing this out. The wording was unclear. The block-median reduction was not applied to the upward-continued gravity field. Instead, it was used to resample the BedMachine topography onto the 1 km and 2 km reference grids used for the local and regional inversions. The equivalent-source model was then evaluated directly at the nodes of these predefined grids. We have revised the text to clarify this workflow.

Line 138: The fact that you mask the gravity data to 2km is quite novel, as most inversions use full grids without any holes. If you retain the inverted bed solution beyond 2km from gravity data, why limit the gravity to this distance? How would you expect the results to change if you used a fully-interpolated grid of gravity data? 

We have revised the manuscript to better clarify the rationale for the 2 km masking. The mask is chosen to ensure that only gravity values directly supported by observations are used, avoiding artifacts from equivalent-source extrapolation. This reflects a deliberately conservative approach, where we want to avoid forcing the inversion to fit values in regions without observational constraint. Using a fully interpolated gravity grid would likely produce a smoother and more spatially continuous result.

Line 181: Is the ground-gravity data already at the surface elevation, since that is where they were collected? So the AntGG data was only upward continued for the ICECAP comparison, since the AntGG data was already at the surface? Or was the ground-gravity data upward continued to 3.5 km as well?

The ground-based gravity data were collected at (near-)surface elevation and were not upward continued. Upward continuation was applied only to the AntGG2021 field for comparison with the ICECAP data. We have clarified this in the revised manuscript.

Line 188: "geometry" is a little confusing, maybe "bed topography"?

Changed.

Line 188: since it's the first mention of it, maybe say "non-terrain gravity effects" for clarity.

Changed.

Line 188: add "and remove" after isolate

Changed.

Line 197: I think this could be re-worded to be clearer, for example: We therefore (i) forward-model and subtract the ‘known’ terrain effect from the isostatically-corrected gravity disturbance, giving a Bouguer disturbance, (ii) estimate a regional component of the Bouguer disturbance by fitting a smooth trend to the gravity data at carefully selected observations points with topography constraints, (iii) remove this regional component to obtain the residual Bouguer disturbance, and (iv) at locations far from topography constraints, replace these residual values with an ensemble of simulate values. 

Thanks for this comprehensive rewording. We have significantly updated this paragraph, and hope this will improve clarity of the workflow.

Line 200: It might be helpful to add a reference to Fig. 2A for this section, as that is what it is describing.

Added.

Line 213: By “remaining broad signals not explained by mapped topography”, are you referring to the gravity effect of the bed resulting from differences between the true bed and Bedmachine? If so, I would replace the word broad with residual, as these signals are not broad, as there are only short wavelength deviations between the true topography and bedmachine topography. I think adding the below sentence will help clarify: This residual Bouguer disturbance is our desired input to the inversion, so we first need to estimate and remove the regional component.

We removed the word ‘broad’.

Line 216: I suggest changing the sentence to as follows for clarity. "It consists of all bouguer disturbance grid points outside the inversion domain and, within the domain, only bouguer disturbance grid points co-located with radar picks ...." 

Changed.

Line 218: This section is a little confusing, as it's discussing a conditional gravity subset, but ends with saying "Interior points ... treated as unknown topography", which seems to refer to a topography data subset, not a gravity subset. Maybe just remove “and treated as unknown topography.”

Changed.

Line 220: As above, I think a little clarification could help the reader. "broad-scale trend" -> "the regional (or can say long-wavelength) component of the gravity disturbance"

Changed.

Line 225: Up to this point, this technique is the same as the DC shift technique used in An et al. 2015, or the constraint-point-minimization used in Tankersley et al 2025, just with a difference of the interpolation algorithm (radial-basis function, vs min.curvature / biharmonic splines), correct?. This might be good to explicitly state to help put your method in context of other literature.

As mentioned above, we added these references to the Methods section to better place our approach in the context of previous gravity inversion studies.

Line 232: short-wavelength residual -> short-wavelength Bouguer disturbance residual

Changed.

Line 233: Add reference to Fig. 2b center.

Added.

Line 238: Add reference to Fig. 2b right side.

Added.

Line 240: This states B is the sum of T (the regional bouguer disturbance) and r (the residual bouguer disturbance), so I would remove the word regional before “non-terrain effects”, since it’s not the regional component, but the entire bouguer disturbance.

Changed.

Line 243: observations -> isostatically-corrected gravity disturbance

Changed.

Line 308: the gridcell the radar data falls in is replaced with the radar data elevation. Does this not lead to pedestals or holes? How do you achieve a smooth transition? Do you do this for the actually bed inversion depths, or just for computing the density model?

This replacement is applied only within the forward modelling step used to refine the density field, rather than modifying the inversion result itself. Sharp discontinuities are mitigated by the regularisation of the inversion and the spatial sensitivity of the gravity response. We updated this paragraph accordingly.

Line 323: the DOI link for Melo and Barbosa is incorrect

Changed.

Line 336: I think you mean to refer to fig 5c.

Changed.

Line 414: grounding line zone -> grounding zone, same on line 424 and caption of Fig 10.

 Changed.

Fig 2. "Cutout gravity away from known topography" might be clearer than "Cutout unknown Topography", since you are cutting out the gravity effect, not the topography, right?

I think this figure should more clearly relate to your steps 1-6 in the text. For example, it is unclear where, or if, step 3 (regional separation) is shown.

We have significantly revised Figure 2 to improve the clarity of the workflow and its correspondence with the methodological steps described in the text.

Fig 3. a) Is the FAA in ‘a’ adjusted for the Moho gravity effect already? If so, I would change the title to “Isostatically-corrected gravity disturbance”

           b) change the title to target terrain effect, and optionally leave Mean g_target in parenthesis

           c) change title to bouguer disturbance (can leave mean B(x) in parentheses)

 We have changed this Figure significantly.

Fig 5. Can you replace “major geographic boundaries” with the grounding line, or is it different? Same in Fig. 6 and 7.

Thank you for this suggestion. The boundaries are derived from the BedMachine v3 mask and represent the extent of floating ice (mask = 3) and ice-free land (mask = 1). These contours therefore include the ice-shelf front and the land boundary, which corresponds to the line where land is exposed, but also to ice/land–ice-shelf interfaces.

As a result, these boundaries do not strictly correspond to the grounding line. We have therefore retained the more general wording but clarified this explicitly in the figure captions to avoid confusion.

Fig.7. Could you add a small histogram of the residuals for each of these figures?  Then can be stacked vertically to the right side of the figures. 

We added histograms to Figures 7.

Fig 10. The black line is showing the grounding line not the coast line right?

As explained above, we have clarified this in the captions.

Fig A1. I think adding to the caption that the right side of the x-axis corresponds to ~ 58 days is useful.

Changed.

Thank you for your positive assessment of the manuscript and for the thoughtful and constructive comments. We agree that small-scale bed features in regions lacking direct gravity or radar constraints should not be interpreted as uniquely resolved topography. In these areas, the short-wavelength component is generated through conditional stochastic simulations and therefore represents plausible morphological variability rather than deterministic bed structure.

Line-by-line responses:

First – how is the ground gravity data incorporated into the raster of gravity data used in the inversion? L132 states “we upward continued the data to 3.5 km”, and this data was interpolated onto 1 km and 2 km rasters for inversion. However, this section seems to just be referring to the airborne data.  Fig 3a shows distinct negative values (if you really zoom in) which suggest the ground data was incorporated, but how this was done should be explained. It may be better to describe the ground gravity data 1st then the airborne, including how the land data was added. Also – the ground and airborne data directly overlie one another. If a grid was created by a simple block-median reduction the oversampled airborne gravity data would significantly reduce the utility of the nine ground based measurements which do however theoretically have significantly higher resolution due to their lower altitude and without along track filtering.

Thank you for pointing this out. We agree that the previous description was unclear. We have added a paragraph in Section 2.2 to clarify that we added these data at their original heights explicitly after upward continuation of the airborne data.

Second – the method appears to generate a certain amount of noise with wavelengths of 5 to 10 km away from any observational data (gravity or radar). It is not clear what the origin of this noise is and the amplitude is hard to determine. However, it appears to be of a similar form to the Gaussian random perturbation shown in Fig. 4 b and d. My concern is that the method may be introducing noise (bumps and depressions) where there is no real data. It then becomes harder to be confident of the bumps in the bed of the main channel are real or inversion artifacts. One way to convince readers the recovered bed is really due to signal in the gravity data is to plot the gravity data above the topography (like in Fig. 9). Alternatively or in addition it would be good to quantify the amplitude of the noise imposed away from areas with direct constraint, as this could be factored in for example when considering the reliability of water routing across a region where teh bed was derived using this technique.

Thank you for this comment. We agree that in areas far from constraints, the short-wavelength bed features are not directly constrained by gravity. Instead, it is generated by Gaussian random fields. We have therefore clarified that these areas should be interpreted cautiously. At the same time, the standard deviation of the model ensemble (Fig 5a) shows a general trend with an increase in model uncertainty further away from data constraints. However, areas with gravity constraints without radar constraints have the highest uncertainties (e.g. under the main trough and under ice-shelf). It suggests the nature of non-uniqueness of the gravity inversion, and also the interpolated method might underestimate the bed uncertainties. In areas without gravity constraints, our method and the smooth interpolation method are all conditioned by radar measurements. The ensemble-based approach allows for more geologically realistic spatial variability and morphological information than would be obtained from purely smooth interpolation products.

Line by line points:

L132 Is it only the airborne data which is upward continued, or the ground based data as well?

It is only the airborne data. We have clarified this more.

L175 – it is unclear how ground and airborne data are combined.

As in the comment above, we agree that this step was missing from the description, and we added a paragraph for clarification.

L184-186 It is not clear why a different correction was applied along the airborne profile corresponding to the ground data compared to the wider airborne dataset.

We thank the reviewer for pointing this out and agree that the previous wording was unclear. We clarified that a significant local mismatch between ICECAP and both the AntGG2021 reference field, and the ground-based observation persisted after application of this regional correction. Therefore, the additional local adjustment was applied.

L205 What observation altitude was used in the terrain model? I presume it is the same as the previously noted upward continuation altitude, but this should be stated.

The terrain effect was calculated at the respective gravity observation elevations of each dataset. We now clarify this explicitly throughout the manuscript.

L294 I understand you repeat the MCMC approach for each ensemble member to estimate the error, but how do you choose the ‘best’ model, or is it taken as the ensemble average of all the recovered topographies? I think this is reported on L335, but not sure.

That is an important observation. We have clarified this and added that the representative bed model is selected as the individual realisation with the best RMS fit to the gravity data, rather than the ensemble average.

L372-374 Suggests the gravity inversion shows topography “characterised by localised depressions and subdued ridges, which are muted or absent in current continent-scale compilations”. This is true, but are these features seen in the gravity data or artifacts of the inversion? As noted above these undulations are seen away from any data (gravity or radar). It may be that the inversion is generating realistic looking texture. This could be good as the undulations provide realistic texture, but it must be acknowledged that away from gravity and radar observations this is still a form of artificial noise.

Thank you, again, we agree that this distinction should be made more clearly in the manuscript. We also added a paragraph to the Results section.

L376 MAE – this is the first occurrence of this abbreviation so should be defined.

We have added ‘mean absolute error’.

L453-456 These lines are duplicated immediately below (L457-L460).

Thanks for this observation. We removed the redundant lines.