| 10 Feb 2026
Gravity Topography Modeling of the Denman Glacier Region Using a Geostatistical Approach
Abstract. The Denman Glacier is one of East Antarctica’s most dynamic outlet systems and is modeled to host the deepest continental marine trough, with the potential to contribute up to ∼ 1.5 m of global sea level rise. Yet its bed geometry remains poorly constrained because airborne radar surveys struggle to image steep, narrow, and deeply incised troughs, and existing continental-scale compilations rely heavily on interpolation or mass-conservation assumptions. During the Australian Denman Terrestrial Campaign 2023/24, high-resolution ground-based gravity measurements were collected across the deepest part of the trough, providing short-wavelength constraints that complement ICECAP airborne gravity and radar data.
We use a two-scale, ensemble-based gravity inversion to reconstruct Denman’s bed topography. A geostatistical separation of terrain and non-terrain gravity effects generates multiple plausible regional background fields, and each is explored with a random-walk Metropolis–Hastings Markov Chain Monte Carlo (MCMC) inversion in which small, spatially correlated Gaussian perturbations modify the bed. Within each realization, candidate geometries are jointly evaluated against the gravity signal and radar picks, with gravity providing the dominant constraint in regions lacking radar coverage.
The resulting ensemble reveals a more rugged, spatially variable, and internally segmented subglacial landscape than represented in current bed products. Along cross-trough profiles, the best-fit gravity-derived bed generally falls between BedMachine and Bedmap3 yet exhibits steeper trough walls and greater lateral relief. Depth-to-magnetic-source estimates further support contrasting lithologies across the trough, with crystalline basement to the west and weaker, sedimentary signatures to the east.
The inferred geometry, characterized by steep flanks and an inland-sloping basin, reinforces the susceptibility of Denman Glacier to Marine Ice Sheet Instability. These results highlight the importance of incorporating geophysical inversion methods into future Antarctic bed-mapping efforts and provide an ensemble of bed realizations suitable for ice flow modeling and assessments of grounding line stability.
Lösing et al. used a novel gravity inversion scheme to model to bed topography beneath Antarctica’s Denman Glacier, revealing a steep-flanked and inland-sloping narrow trough. Their approach utilized existing airborne gravity data, and a recently collected ground-based gravity profile. The geostatistical modelling approach, compared to traditional deterministic gravity inversions, provides an avenue for the much-needed estimates of spatial variability in model uncertainties. While I think there is some work to be done to clarify some of the methods and approaches, I think this paper makes a substantial and important contribution to the field.
General comments:
- 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.
- 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.
- 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
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.
- 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?).
Specific comments:
Line 44: You might add (Fig. 1c) after "These gaps".
Line 98: Maybe use ground survey instead of ground line to avoid confusion with grounding line.
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?
Line 131: Did you block-reduce the line data before fitting the equivalent sources? If so, what cell size?
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?
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?
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?
Line 188: "geometry" is a little confusing, maybe "bed topography"?
Line 188: since it's the first mention of it, maybe say "non-terrain gravity effects" for clarity.
Line 188: add "and remove" after isolate
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.
Line 200: It might be helpful to add a reference to Fig. 2A for this section, as that is what it is describing.
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.
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 ...."
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.”
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"
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.
Line 232: short-wavelength residual -> short-wavelength Bouguer disturbance residual
Line 233: Add reference to Fig. 2b center.
Line 238: Add reference to Fig. 2b right side.
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.
Line 243: observations -> isostatically-corrected gravity disturbance
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?
Line 323: the DOI link for Melo and Barbosa is incorrect
Line 336: I think you mean to refer to fig 5c.
Line 414: grounding line zone -> grounding zone, same on line 424 and caption of Fig 10.
Fig 2.
Fig 3.
Fig 5.
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.
Fig 10. The black line is showing the grounding line not the coast line right?
Fig A1. I think adding to the caption that the right side of the x-axis corresponds to ~ 58 days is useful.