Gravity inversion for sub-ice shelf bathymetry: strengths, limitations, and insights from synthetic modeling

Sub-ice-shelf bathymetry strongly influences ice shelf stability by guiding melt-inducing water masses and through pinning points that resist the flow of the overriding ice. Collecting sub-ice-shelf bathymetry data using active source seismic surveying or direct observations is accurate but time-consuming and often impractical. Gravity methods provide a pragmatic, but more uncertain, alternative, by which observed variations in Earth's gravitational field are used to estimate the underlying bathymetry. We utilize a new open-source gravity inversion algorithm developed specifically for modeling sub-ice-shelf bathymetry and estimating the spatially variable uncertainty in the results. The inversion is tested on a suite of models created with high-resolution multibeam bathymetry data. These tests enable (1) determination of the best practices for conducting bathymetric inversions, (2) recognition of the limitations of the inversion and uncertainty quantification, and (3) identification of where community efforts should be focused for the future determination of Antarctica's sub-ice-shelf bathymetry. With an airborne gravity survey with 10 km spacing, 3 mGal of errors, a distribution of known bathymetry measurements, and a regional gravity field representative of the average Antarctic ice shelf, we achieve a root mean squared error of the inverted bathymetry of 23 m. We find that estimating and removing the regional component of gravity before the inversion is the largest source of error in the resulting bathymetry model, but this error can be greatly reduced with additional known bathymetry points. We analyzed Antarctic ice shelves and predict that, if high-resolution gravity data were available, gravity inversion could potentially improve bathymetry models for 95 % of them compared to interpolated products like Bedmap2.

1 Introduction

Over the past few decades (1992–2020), Antarctica has been losing a significant amount of ice, averaging 92 billion tons per year, adding over 7 mm to global sea-level rise (Otosaka et al., 2023). Most of this ice loss is happening in places where relatively warm ocean waters come into contact with the floating ice shelves that fringe the continent (Rignot et al., 2019). These ice shelves exert a resistive force, buttressing, slowing the flow of inland ice into the ocean. This buttressing comes from lateral drag and resistive stresses where the ice touches the sea floor at pinning points (Dupont and Alley, 2005; Matsuoka et al., 2015). But in the past 20 years, about 37 % of these pinning points have shrunk due to ice shelf thinning (Miles and Bingham, 2024), reducing this buttressing effect and increasing ice flux into the ocean (e.g. Pritchard et al., 2012; Scambos et al., 2004). Because of Antarctica's cold climate and limited surface melting, most ice shelf thinning occurs from below, through basal melting. The largest ice shelves, like Ross, Ronne-Filchner, and Amery, are considered “cold-water” shelves, based on the ocean conditions beneath them (Adusumilli et al., 2020).

In these cold-water shelves, basal melt tends to happen at the ice front and along the deep grounding zones, the region where ice transitions from grounded to floating. This grounding zone melt occurs because the melting point of ice drops under pressure, so even cold, salty waters can melt ice at depth (Alley et al., 2015). This cold and salty water, such as High Salinity Shelf Water (formed during sea ice formation), travels along seafloor troughs towards the grounding zones (e.g. Tinto et al., 2019). This highlights sub-ice shelf bathymetry as a strong control on the future of the Antarctic Ice Sheet through its routing of melt-inducing waters and the control on pinning point geometry.

Techniques for collecting sub-ice-shelf bathymetry data are limited. Radar (ground or airborne) cannot image through the ocean cavity due to signal attenuation caused by water. Conventional over-ice seismic surveys provide accurate measurements, but can only cover small spatial scales (10's of kilometers), unless significant funding and time are available (e.g., the RIGGS survey; Bentley, 1984). Similarly, direct observations through drill holes or Autonomous Underwater Vehicles (AUVs) are impractical for adequate coverage of large ice shelves. However, because the seafloor is much denser than water and ice, it creates a detectable signal in Earth's gravity field. This signal can be observed from airborne or ground-based gravity surveys and then used in a method called gravity inversion to model the shape of the sub-ice shelf seafloor. This technique is commonly used throughout Greenland (e.g. Millan et al., 2018; An et al., 2019a) and Antarctica (e.g. Vaňková et al., 2023; Locke et al., 2025; Constantino and Tinto, 2023; Charrassin et al., 2025) for modeling bathymetry where it is otherwise impractical to measure.

There are two main types of gravity inversions. Density inversions estimate how density changes underground, often used in mineral exploration. Geometry inversions (also called structural inversions) estimate the shape of a surface which represents a change in density, like the seafloor beneath an ice shelf. They are commonly used to model the boundary between Earth's crust and mantle (e.g Uieda and Barbosa, 2017; Borghi, 2022), or the contact between sediment and basement rock (e.g. Santos et al., 2015; Barbosa et al., 2007). Density inversions are common, well-supported by open-source software, and widely used in industry (Oldenburg and Pratt, 2007; Cockett et al., 2015; Rücker et al., 2017). Geometry inversions, on the other hand, don't have reliable, freely available tools. Researchers often rely on expensive, proprietary software that is hard to customize. These tools also struggle to handle uncertainty or use known measurements to improve accuracy.

We use an open-source Python-based geometry inversion tool called Invert4Geom. It is designed specifically to model topography from gravity data and to estimate spatially variable uncertainty in the results. We run a set of synthetic tests to see how well it performs under different conditions, such as noisy gravity data, different spacing between measurements, and the presence of a long-wavelength geologic component of the gravity field. We analyze all Antarctic ice shelves to predict inversion performance based on the results of these synthetic tests. We discuss the general assets and limitations of using gravity inversions to recover bathymetry and provide guidance for field data collection. Finally, we highlight ice shelves for which bathymetry would be improved using a gravity inversion, and ice shelves where inversion would degrade the bathymetry reliability.

2 Methods

2.1 Geophysical inversion

Geophysical inversion involves determining a set of physical parameters that, when forward modeled, produce predicted data that closely match observed measurements. Forward modeling is the calculation of what data would be observed given a model of physical parameters. In this study, we perform a least-squares inversion to solve the matrix equation

$$\begin{array}{}\text{(1)}& \mathbf{J}\mathit{x}=\mathit{m},\end{array}$$

where J is the Jacobian (or sensitivity) matrix, x represents the perturbations applied to the model's physical parameters, and m are the data misfit values. Each element of the Jacobian quantifies the sensitivity of a predicted data value to a small change in a physical parameter, i.e., ${\displaystyle \frac{\partial f_j}{\partial p_i}}$, which is the partial derivative of the jth predicted data with respect to the ith physical parameter. In our application, the physical parameter of interest is seafloor topography, and the observed data are variations in Earth's gravity field, specifically the residual component of the gravity misfit, described below. Solving Eq. (1) for x using an iterative least-squares approach yields corrections to the initial bathymetry model that reduce the mismatch between observed and modeled gravity. A detailed mathematical derivation is provided in Appendix A1.

2.2 Gravity reduction

Prior to inversion, gravity data must be reduced to isolate the signal relevant to the density contrast of interest. The reduction and preprocessing workflow (Fig. 2a) comprises the following steps:

1. Field corrections: Measured gravity is corrected for instrument drift, flight line leveling inconsistencies, aircraft maneuvers, and the effects of measuring gravity on a moving platform (Eötvös correction) (Hinze et al., 2005, 2013). These corrections yield observed gravity (g_obs), which reflects the combined effect of Earth's gravitational field and its rotation.

2. Normal gravity correction: We then subtract the gravity of a reference Earth model (denoted γ), which accounts for variations due to latitude and elevation. This yields the gravity disturbance (δg), sometimes referred to as the free-air anomaly, which contains only gravity effects resulting from deviations between the actual Earth and the reference model (Fig. 1).

3. Terrain correction: To isolate internal density structure, the effects of known topographic masses are removed from the disturbance. These terrain mass effects reflect all topographic deviations between the normal Earth and the real Earth, where the normal Earth is defined by air above the ellipsoid and crust below it (Fig. 1a). Once calculated and removed (see Appendix A2 for details), the remaining gravity signal is the topography-corrected disturbance (δg_TC), sometimes referred to as the complete Bouguer anomaly. This contains only gravity effects from density variations within the Earth.

Lastly, Invert4Geom uses gravity data on a regular grid, as opposed to scattered observation points. We use the equivalent source technique to interpolate the scattered gravity data onto a regular grid. Equivalent sources find a set of point sources at depth which, when forward-modeled, reproduce the measured gravity data (Soler and Uieda, 2021; Dampney, 1969). These sources can then be used to predict the gravity anomaly at any desired point, such as each location on an even grid.

Figure 1 Defining the terrain mass effect (a–d) and discretizing the components (e–h). (a) The normal Earth model, (b) the simplified real Earth model, with air, ice, water, and earth, (c) terrain masses separated into components above and below the ellipsoid, (d) real Earth topographic masses which are anomalous relative to the normal Earth and their corresponding assigned density contrasts, which, combined, make the Terrain Mass Effect, (e) prisms between the ellipsoid and the ice surface, (f) prisms between the ellipsoid and the water surface (ice base), (g) prisms between the ellipsoid and the rock surface (topography onshore and bathymetry offshore). (h) the combination of prisms of (e)–(g), where overlapping sets of prisms create the different shades and corresponding densities.

For conventional gravity inversions targeting internal interfaces (e.g., the Moho or sediment-basement contact), δg_TC is typically used. However, for bathymetry inversions, we retain the gravitational signal from the topography itself, and therefore apply only a partial terrain correction. The resulting gravity anomaly, which we refer to as the partial topography-corrected disturbance (δg_{TC_partial}), is given by:

$$\begin{array}{}\text{(2)}& \mathit{\delta }{\mathit{g}}_{\text{TC\_partial}}={\mathit{g}}_{\mathrm{obs}}-\mathit{\gamma }-{\mathit{g}}_{\mathrm{ice}}-{\mathit{g}}_{\mathrm{water}}.\end{array}$$

This removes the gravitational effects of ice and water, while retaining the effects of the rock surface (onshore topography and offshore bathymetry). A rigorous discussion of this formulation and its implementation is provided in Appendix A2.

Figure 2 Schematic workflow diagram for (a) gravity processing steps and (b) running the inversion. Blue boxes show the input datasets, green boxes show the key processes, and red lines signify iterative components of the inversion.

2.2.1 Forward modeling

Both gravity reduction and inversion require forward modeling of the gravity field associated with a given topographic surface. The surface is discretized into an array of adjacent vertical rectangular prisms, and the total gravitational field at any given observation point is computed as the sum of the individual contributions of each prism. For prisms of constant density, gravity is calculated using the analytical solutions of Nagy et al. (2000). We use a density-contrast approach (Fig. 3e) rather than using absolute densities (Fig. 3b). This approach reduces computational costs (using one rather than two prism layers per topographic interface) and minimizes gravitational edge effects due to thinner prisms and lower absolute density values. Further details are provided in Appendix A4.

Figure 3 Equivalent methods of discretizing a density contrast (a) across a surface. (b) Two prism layers, each assigned absolute densities. One prism layer either above (c) or below (d) the surface, using density contrast values. (e) One prism layer around a reference surface, with prisms above having a positive density contrast and prisms below a negative contrast.

2.2.2 Gravity misfit

The gravity misfit used in the inversion is defined as

$$\begin{array}{}\text{(3)}& \mathit{m}=\mathit{\delta }{\mathit{g}}_{\text{TC\_partial}}-f\left({\mathit{p}}_{\mathrm{0}}\right)\end{array}$$

where f(p₀) is the gravity predicted from the initial bathymetry model. Note that as the initial model approaches the truth, the misfit is equivalent to the topography-corrected disturbance (δg_TC), or the complete Bouguer anomaly. This misfit comprises two components: (1) a residual component caused by differences between the actual and modeled bathymetry, and (2) a regional component due to deeper subsurface density variations not accounted for in the model. Because the inversion targets corrections to bathymetry alone, it is essential to isolate the residual component by estimating and removing the regional component before inversion. Further details are provided in Appendix A3.

2.2.3 Regional separation

The regional component of the gravity misfit can be estimated using several methods, most of which assume that the regional field consists of longer-wavelength anomalies, while the residual field dominates the short-wavelength components. This assumption underpins techniques such as upwards continuation or the use of high-altitude gravity data to estimate the long-wavelength signals (e.g. Tinto et al., 2015; Muto et al., 2016), as well as low-pass filtering of the gravity data (e.g. Eisermann et al., 2020). Alternatively, the regional field can be forward-modeled from a geologic model, developed using other geophysical data (Hodgson et al., 2019), prior geologic knowledge (Tinto and Bell, 2011), or crustal density models derived from gravity inversion (Tinto et al., 2019; Cochran et al., 2020; An et al., 2019b).

The last method, which we utilize in this paper, leverages locations with known bathymetry depths, referred to as constraint points. The inversion defines the residual misfit as the gravity effect due to differences between the true and starting bathymetry. At constraint points, where the two bathymetries are nearly identical (within measurement error), the residual misfit should be zero. Consequently, the entire misfit at these points is attributed to the regional component. By sampling the misfit at each constraint point and interpolating these values, we estimate the regional field.

This technique has key advantages. It imposes a form of smallness regularization, preventing significant changes to the depth at constraint points. Aside from the choice of interpolation method, it avoids the subjective parameter selection required by the other techniques, such as low-pass filter widths or upward continuation heights. However, it assumes the wavelength of the regional field is larger than the spacing between constraint points; shorter-wavelength regional anomalies will not be captured and may be absorbed into the residual (i.e., the inverted bathymetry). This approach, often referred to as a “DC-shift” (Millan et al., 2020; Yang et al., 2021; An et al., 2019b), is here termed Constraint Point Minimization to emphasize that it estimates a spatially variable regional field. Estimating the regional field is among the most challenging steps of a bathymetric inversion, particularly in areas of sparse constraint coverage (Brisbourne et al., 2014). Over or underestimating the regional field directly biases the inversion. Once the regional field is estimated, it is removed from the total misfit to yield the residual misfit:

$$\begin{array}{}\text{(4)}& {\mathit{m}}_{\mathrm{residual}}=\mathit{\delta }{\mathit{g}}_{\text{TC\_partial}}-f\left({\mathit{p}}_{\mathrm{0}}\right)-{\mathit{m}}_{\mathrm{regional}}.\end{array}$$

Combining with Eq. (2):

$$\begin{array}{}\text{(5)}& {\mathit{m}}_{\mathrm{residual}}={\mathit{g}}_{\mathrm{obs}}-\mathit{\gamma }-{\mathit{g}}_{\mathrm{ice}}-{\mathit{g}}_{\mathrm{water}}-f\left({\mathit{p}}_{\mathrm{0}}\right)-{\mathit{m}}_{\mathrm{regional}},\end{array}$$

where g_obs is the observed gravity, γ, g_ice, g_water are forward-calculated, f(p₀) is the forward-calculate effect of the starting model, and m_regional is the estimated regional field. The remaining residual ideally reflects only the discrepancy between the initial and true bathymetries.

2.3 Running the inversion

2.3.1 Jacobian matrix

The first step of the inversion (Fig. 2b) is computing the sensitivity matrix, or Jacobian (Eq. 1). Each entry in the Jacobian is the partial derivative of gravity produced by a specific prism at a given observation point, with respect to the prism's thickness. Although an analytical solution exists for this derivative (Nagy et al., 2000), we use an approximation to reduce the computational cost. Each vertical prism is approximated as a sector of an annulus (a cylindrical shell), whose vertical derivative can be calculated more efficiently (McCubbine, 2016). See Appendix A5 for a derivation of this approximation.

2.3.2 Least-squares solver

With the Jacobian and the gravity residuals, Eq. (1) is solved to find the set of surface corrections (Δp) that minimizes the objective function (Fig. 2b). We use a damped iterative least squares solver (LSQR; Paige and Saunders, 1982), which finds the minimum-norm solution, that is, the one with the smallest $\left|\right|\mathbf{\Delta }\mathit{p}|{|}^{\mathrm{2}}$ among all solutions that fit the data equally well.

2.3.3 Regularization

Inversion problems involving potential fields are inherently ill-posed: many solutions can satisfy the data. Regularization helps stabilize the solution by incorporating prior knowledge or desirable properties (Aster et al., 2018). We use two forms: (1) smallness, which constrains the inversion to remain close to known bathymetry measurements, and (2) smoothness, which suppresses unrealistic roughness in the inverted surface.

Smallness

Invert4Geom provides two smallness regularization techniques. The first comes from the Constraint Point Minimization method used for estimating the regional field. This method ensures that near constraint points, most of the gravity misfit is accounted for in the regional field, resulting in minimal bathymetric correction near those locations. An optional second method uses a spatial weighting grid to taper surface corrections towards zero at constraint points. The grid assigns a weight of 0 at constraint points, increasing smoothly to 1 with distance. This weight is multiplied element-wise with each iteration's surface correction grid before it is applied (Fig. 2b), ensuring a smooth transition from constrained to unconstrained areas. While effective, this increases both the number of interactions and total runtime. Both methods require constraint points, but in most applications, at least some are available from radar, outcrops, or shipboard surveys.

Smoothness

Smoothness regularization is enforced through damping in the LSQR solver. The damping parameter restricts the magnitude of surface corrections in each iteration, reducing noise and promoting realistic bathymetry. While this gives smoother results, it requires extra iterations. Section 2.5 outlines how to choose an optimal damping value.

After each iteration, the surface corrections (possibly weighted) are applied to update the bathymetry. The new model is re-discretized into prisms, forward-modeled, and a new residual is computed. This process repeats until one of the termination conditions is met.

2.4 Stopping criteria

The inversion terminates when the first of the following conditions is met:

1. The maximum number of iterations (N) is reached;

2. The ℓ²-norm of the residual ($\left|\right|\mathit{m}|{|}_{\mathrm{2}}$, Eq. 4) falls below a defined threshold;

3. The changed in ℓ²-norm between iterations (Δℓ²-norm) is below a defined tolerance.

2.5 Hyperparameter estimation

Hyperparameters are parameters whose values need to be chosen (not calculated) and that also have a significant effect. Two key hyperparameters affect inversion quality: the regularization damping parameter and the uniform density contrast value (if it is not previously known). Below, we describe how to determine their optimal values.

2.5.1 Damping value cross-validation

The damping parameter balances solution stability and data misfit. Too high a value, the solution under-fits the data; too low a value, the data are overfit and noisy results are produced. To find an optimal value, we use cross-validation of the input gravity data (Fig. 4a; Uieda and Barbosa, 2017). The gravity data are split into training and testing sets. The inversion is run using the training set for a candidate damping value. The root mean square (RMS) misfit between the unused testing gravity points (the partial topography-corrected disturbance minus the estimated regional gravity) and the forward-modeled gravity effect of the resulting bathymetry gives a score. This is repeated for multiple damping values, and the one with the lowest score is selected.

To construct the training and testing sets, the data are interpolated onto a grid at half the target resolution using equivalent sources (Soler and Uieda, 2021). The training set is those points that fall on the desired-resolution grid, and the remaining interleaved points become the testing set.

Figure 4 Workflows for selecting the optimal values for the two inversion parameters: (a) the regularization damping value and (b) the density contrast value. Both workflows follow the general methods of Uieda and Barbosa (2017). Blue boxes show input datasets and green boxes show the key processes.

2.5.2 Density contrast estimation

If no seafloor density contrast is known, it can be estimated via the inversion (Fig. 4b). For a candidate density contrast value, the full inversion is run, including the creation of the starting bathymetry and the regional estimation, and the output bathymetry is compared to known constraint point depths. The RMS of these differences is the score. The optimal density contrast minimizes this score. If the Constraint Point Minimization method is used for regional separation, these scores become biased, since Constraint Point Minimization enforces a zero residual at constraint points. Two alternatives are:

• Cross-validation: Split constraints into K-folds of training and testing subsets (Hastie et al., 2009), using spatial blocks to minimize autocorrelation (Roberts et al., 2017). For each fold, the training set is used in the inversion, and the testing set is used for scoring.

• Alternative regional estimation: Temporarily use a different regional field method during optimization, then revert to Constraint Point Minimization once the best density contrast is found.

2.5.3 Optimization routines

Both hyperparameter estimation methods rely on computing a score for each candidate value. Invert4Geom supports two search strategies: (1) a grid search, where values are tested uniformly across a defined range, and (2) an adaptive optimizing approach, which uses prior scores to refine the search space, converging more efficiently on the optimal value (Akiba et al., 2019).

2.6 Uncertainty quantification

We use a Bayesian framework, specifically Monte Carlo simulation, to assess the spatial uncertainty of the inverted bathymetry (Jansen et al., 1994). For a generalized problem of inputs x and a function y(x), Monte Carlo simulation estimates the probability distribution of y(x) when you know the probability distribution of x (Helton et al., 2006). For a set of parameters, we draw N samples from their input distributions (Fig. 5b) and run the full inversion workflow for each set (Fig. 5c). The resulting N bathymetry models are used to compute a cell-wise weighted standard deviation (Fig. 5d). Each model is weighted by the inverse square of its RMSE with constraint points (Appendix A6). This standard deviation grid serves as our spatial uncertainty estimate. To evaluate each parameter’s individual effect, separate Monte Carlo runs can isolate their influence by sampling one parameter at a time.

Figure 5 Schematic workflow diagram for the Monte Carlo uncertainty analysis. (a) The Monte Carlo simulation, consisting of (b) sampling the inputs from their respective distributions, and (c) implementing the various components of the workflow, depending on which inputs have been sampled. This process is repeated N times, yielding N inverted bathymetry grids. (d) Weighted cell-specific statistics are computed. Each inverted bathymetry grid (large purple boxes) is weighted by the respective inversion's constraint depth RMSE. These weights are used to compute the weighted mean and weighted standard deviation of each grid cell (bold purple squares).

2.6.1 Parameter distributions

Defining input parameter distributions is critical. For point data (gravity, constraints), we assume normal distributions based on user-provided or estimated uncertainties. For gravity, uncertainties often come from survey crossover errors; for constraints, uncertainty may be based on data source or be proportional to depth. For hyperparameters (density contrast and damping), the results of the optimizations can be used to define each parameter's distribution. For a well-defined best score, a normal distribution centered on that value with a small standard deviation can be used. For an even distribution of similar scores, a uniform (Boxcar) distribution can be used with the width defined by the range of the similar scores. To assess uncertainty in the regional field, we also assign distributions to parameters used in its estimation (e.g., equivalent source depth or damping in Constraint Point Minimization). Probability distributions for these parameters can be derived from the respective optimization procedures.

2.6.2 Sampling strategy

Since full inversion workflows are computationally intensive, we use Latin Hypercube Sampling to efficiently cover the input parameter space with a small number of runs (Jansen et al., 1994; Helton et al., 2006). In Latin Hypercube Sampling, each input distribution is divided into N equal-probability intervals, and one value is sampled from each. The N samples for each parameter are randomly combined into N total input sets. For point-wise-sampled inputs (e.g., gravity data), regular random sampling suffices.

2.7 Metrics

In this study, we use a few unique metrics to describe the distribution of data and assess model performance. There are many ways to describe how data points are spaced out, whether it is airborne gravity data or over-ice seismic survey points, but not all metrics work equally well for all survey types. For example, spacing (usually the average distance between points) only makes sense if your data are evenly spread across a grid. If data follow tracks, like flight lines or ship surveys, the mean distance between points is skewed toward unrealistically small values. Spatial density (how many points fall within a certain area) has a similar issue; it is strongly affected by the sampling frequency along those tracks. Line spacing, the average distance between survey lines, only works well for grid-based surveys and does not handle gaps, tie lines, or irregular coverage very well. To enable robust comparisons of data distribution with varying collection types and spatial distributions, we use a metric we call median proximity, which describes the median distance of each grid point in a region to the nearest data point. To calculate this, we create a grid for the region of interest with a small cell size (∼ 100 m) and for each cell, calculate the distance to the nearest data point. We then calculate the median value of all the grid cells. This approach avoids the pitfalls of sampling frequency, uneven line spacing, and coverage gaps, and gives a more consistent picture of how well an area is constrained by data. We use median proximity to describe how bathymetry constraints (Fig. 6d) and gravity data (Fig. 11e) are distributed, which helps us compare the spatial distributions of data between various Antarctic ice shelves and synthetic datasets created here.

To evaluate the performance of various techniques in our synthetic tests, we use root mean square error (RMSE). RMSE is useful because it ignores the sign of the errors and is less sensitive to outliers than the mean. For example, the RMSE between a “true” bathymetry grid and an inverted grid tells us how far off our model is, on average. Finally, we introduce a topographic improvement score to measure whether the inversion gives better results than just interpolating the sparse bathymetry points. This score is the difference between the RMSE of the interpolated starting model and the RMSE of the inverted model: RMSE_{starting_bathymetry}−RMSE_{inverted_bathymetry}. A score of 0 means no improvement, while a positive value shows that the inversion helped reduce error and improve the bathymetry.

Lastly, throughout this study, we discuss the relative amplitudes of regional gravity fields. Mathematically, we quantify this as the standard deviation of the regional field, but in text we use the term variability. This refers to the amplitude of variation in the signal, regardless of the mean level, allowing a comparison between different regional gravity fields.

2.8 Antarctic ice shelf analysis

To set up our synthetic tests realistically, we needed to choose values for several key parts of the bathymetry inversion process, such as how accurate the data are, how the data are spread out, and how variable the regional gravity field is. To choose which values to test, we analyzed all 164 Antarctic ice shelves to figure out the typical ranges for these variables. For each shelf, we collected nearby bed measurements from IBCSO v2 (Dorschel et al., 2022b) and BedMachine v3 (Morlighem, 2022). These included data from multibeam and single-beam sonar, over-ice seismic surveys, airborne and surface-based radar, and exposed rock outcrops. We left out any points derived from indirect methods like gravity inversion, mass conservation, or streamline diffusion. From the compiled points, we calculated the constraint proximity as described earlier.

For each ice shelf we used AntGG-2021 gravity data (Scheinert et al., 2024) and Bedmap2 (Fretwell et al., 2013) ice surface, ice base, and bed topography to calculated the topography-corrected gravity disturbance, δg_TC, using densities of 917, 1025, and 2670 kg m⁻³ for ice, water, and earth, respectively. We chose not to use Bedmap3 (Pritchard et al., 2025) as it included gravity-inversion derived bathymetry models for some shelves. Using our compiled bed measurements, we then used the Constraint Point Minimization method to estimate the regional gravity field. For each shelf, we created grids of constraint proximity, topography-corrected gravity disturbance, and regional field, and then masked those to the ice shelf outlines based on MEaSUREs (Mouginot et al., 2017). To use realistic values of the error and spatial distribution of gravity data, we analyzed a range of published sub-ice shelf bathymetry studies (Tinto et al., 2019; Yang et al., 2021; Eisermann et al., 2021; Constantino et al., 2020; Jordan et al., 2020; Eisermann et al., 2020, 2024; Hodgson et al., 2019; Wei et al., 2020; Cochran et al., 2020; Millan et al., 2020; Bird et al., 2025).

2.9 Synthetic testing

To assess the performance of the inversion software, showcase its use, and explore the sensitivity of the inversion to certain parameters, we ran a series of tests using a semi-synthetic model. This model is based on bathymetry data from Antarctica's Ross Sea (Fig. 6a), collected using shipborne sonar and compiled by IBCSO (Dorschel et al., 2022a). In this region, the median distance to the nearest bathymetry data point is ∼ 400 m (Fig. 6d). IBCSO originally gridded the data at 500 m spacing, but for our purposes we filtered and regridded it at 2 km spacing. We forward modeled its gravity effect using a constant density contrast of 1476 kg m⁻³, based on typical densities of seawater ∼ 1024 kg m⁻³ and seafloor rock (∼ 2500 kg m⁻³) (Bamber and Bentley, 1994). This gave us a synthetic “observed” gravity dataset, calculated at 2 km spacing and 1 km observation height.

Figure 6 Ross Sea topographic data used for synthetic inversions. (a) 500 m IBCSO bathymetry (Dorschel et al., 2022a), re-gridded at a 2 km spacing, (b) difference between a and c, (c) starting bathymetry model from the interpolation of constraint points, shown by black dots, with smaller dots denoting the constraint outside the ice shelf boundary, (d) distance to the nearest IBCSO data point (black dots), giving a median proximity of ∼ 400 m. (e) Crystalline basement topography from the ANTOSTRAT shipborne seismic data compilation (Brancolini et al., 1995). Inset map shows the study site location in the Ross Sea.

Across different tests, we used the full-resolution gravity grid, added synthetic noise, or sampled the data along simulated flight lines to mimic an airborne survey. In some cases, we also added a regional gravity component by forward modeling the effect of basement topography using ANTOSTRAT seismic data (Fig. 6e; Brancolini et al., 1995). Because we know the “true” bathymetry and the “true” regional gravity field in these synthetic setups, we can directly assess how well the inversion performs under different conditions, as well as how reliable our uncertainty estimates are.

We ran five core inversion tests:

• Test 1: the ideal case, using full-resolution, noise-free gravity data with no regional component.

• Test 2: adds noise to the gravity data.

• Test 3: uses simulated airborne gravity data instead of the full grid.

• Test 4: includes a regional gravity component.

• Test 5: combines all the above factors for a more realistic test case.

We also ran four inversion ensembles, each with 100 runs:

• Ensemble 1: explores the relationship between the number of constraint points and the variability of the regional field. We test 10 different amounts of the number of constraint points, from 4 points to 650 points, and 10 different levels of the regional gravity field variability, from very low variability (∼ 4 mGal standard deviation) to a very high variability (∼ 40 mGal standard deviation).

• Ensembles 2–4: look at how gravity noise and flight line spacing affect results. For these ensembles, we test 10 values of gravity noise, from none to 5 mGal of random noise, and 10 numbers of gravity flight lines, from 6 to 60 lines across the region. These values are tested with no regional gravity field, and medium variability field, and a highly variable field (ensembles 2, 3, and 4, respectively).

To make sure our tests reflect real-world conditions, we set all these parameters (e.g., constraint proximity, gravity noise, and flight line spacing) based on the ranges we found across Antarctica's ice shelves and in published gravity surveys. Each test and ensemble has a corresponding Jupyter notebook, available in the supplementary material.

2.10 Software implementation

The inversion algorithm we use in this study is implemented in Python as the open-source package Invert4Geom (Invert4Geom Community and Tankersley, 2025), released under the MIT license. It builds on many scientific Python libraries. For gravity modeling, both with prisms and equivalent sources, we use Harmonica (Fatiando a Terra Project et al., 2024b; Soler and Uieda, 2021). For spatial operations like gridding and interpolation, we rely on Verde (Uieda, 2018). To generate figures and download geospatial data, we use PyGMT and PolarToolkit (Uieda et al., 2021; Tankersley, 2024; PolarToolkit Community and Tankersley, 2025). We handle least squares regression with SciPy (Virtanen et al., 2020), optimize hyperparameters with Optuna (Akiba et al., 2019), and perform uncertainty sampling using UQpy (Olivier et al., 2020).

All of the experiments were run in Jupyter notebooks (Pérez and Granger, 2007), which are documented to explain the details of using the software. These notebooks were used to create a website; https://mdtanker.github.io/synthetic_bathymetry_inversion/ (last access: 18 November 2025​​​​​​​). The full source code, including all results and figures used in this paper, is publicly available in a GitHub repository, https://github.com/mdtanker/synthetic_bathymetry_inversion (last access: 18 November 2025​​​​​​​), and the state of the code used in this analysis is preserved at Zenodo (Tankersley et al., 2025).

3 Results

3.1 Antarctic ice shelf analysis

To assess gravity inversion suitability across Antarctica, we compiled all available bathymetry measurements and created grids of three key variables for each ice shelf: constraint proximity, topography-corrected gravity disturbance, and regional gravity misfit. To avoid bias from the large number of small shelves, we restricted our analysis to the 50 largest ice shelves (those exceeding 1000 km²) out of the 164 named shelves in Antarctica. Among these, the median constraint proximity ranged from 2.5 km (Thwaites) to 32 km (Conger-Glenzer), with an overall mean of 10 km and median of 8 km (Fig. 7c). Most shelves had median proximities between 5 and 10 km. The maximum distance to the nearest constraint ranged from 14 km (Drygalski) to 84 km (Shackleton). The skewness of constraint proximity, which reflects how evenly distributed the constraint points are, ranged from near zero (0.008) for Conger-Glenzer to strongly positive (1.5) for Ronne-Filchner (Fig. 7f – x-axis).

We also evaluated the regional gravity field variability (amplitude of variation) using two measures: the standard deviation of the topography-corrected gravity disturbance and the standard deviation of the regional gravity misfit. We found the former to be a more reliable indicator, as it is less influenced by the number and distribution of constraint points. Although it can be affected by the resolution of the topographic model, our focus is on long-wavelength signals, which are well captured by the BedMap2 topography grids. Across the 50 shelves, regional field variability (based on topography-corrected gravity disturbance) ranged from 7 mGal (Holmes) to 66 mGal (Ekström), with a mean of 23 mGal and median of 18 mGal (Fig. 7d). Most ice shelves exhibited values between 10 and 20 mGal.

To better understand typical survey parameters, we reviewed published studies on airborne gravity data collection over Antarctic shelves. Reported flight line spacings typically range from 5 to 20 km, with 10 km being the most common. Estimating gravity uncertainty, however, is more complicated. Most studies use crossover analysis, but methods vary widely: some report the mean, others the standard deviation or RMS of crossover differences; some apply leveling (of varying orders), while others do not, making comparisons difficult. In general, reported uncertainties range between 1 and 5 mGal. To provide a consistent metric, we used the AntGG-2021 uncertainty grid to calculate the RMS uncertainty across each ice shelf. The lowest shelf-wide RMS uncertainty we found was 3.8 mGal (Pine Island), highlighting the inconsistency in reporting practices.

Figure 7 Antarctic ice shelf analysis; (a) proximity to nearest bed measurements and (b) topography-corrected gravity disturbance for all Antarctic ice shelves using AntGG-2021 gravity disturbances (Scheinert et al., 2024) and Bedmap2 grids of ice surface, ice base, and bedrock topography (Fretwell et al., 2013). Black outlines in (a) and (b) show ice shelf boundaries from MEaSUREs (Mouginot et al., 2017). (c–d) histograms of median constraint proximity and standard deviation of topography-corrected gravity disturbance for ice shelves larger than 1000 km², respectively; (e) scatter plot of median constraint proximity and standard deviation of topography-corrected gravity disturbance, with the range of values tested in Ensemble 1 shown as a blue box and values used in Test 5/Ensemble 3, and Ensemble 4 shown as blue and red triangles, respectively. (f) scatter plot of skewness of constraint proximity and standard deviation of topography-corrected gravity disturbance. Ice shelves are labeled numerically (see legend for names), and the 24 previously individually-inverted shelves are in red.

Figure 8 Results of the inversions and uncertainty analyses for each of the five scenarios (each column is a scenario). Thefirst row (1a–6a) shows the inverted bathymetry results, all with the same colorscale, second row (1b–6b) shows the inverted bathymetry error (row 1 minus true bathymetry (Fig. 6a), third row (1c–5c) shows the uncertainty assessment for each inversion. All plots or rows 2 and 3 share the same colorscale to enable direct comparisons between error and uncertainty, and between results of each Test. Last column shows the starting bathymetry model from the interpolation of constraint points, and its error. White dots depict locations of points of known bathymetry depth (constraints). Profile A-A', location in (1a), shows the inverted bathymetry for each scenario. Scalebar in (1a).

Figure 9 Comparison of the five inversion scenario results. Each scenario has an associated color, shown in the legend. (a) The RMSE of the inverted bathymetries (triangles) and the RMS of the uncertainty estimates (squares). (b) Radially averaged power spectral density (PSD) of bathymetry errors showing the dominant wavelengths of the inversions' errors. (c) Results of the damping cross-validations, with optimal found values shown by squares and +/- the estimated standard deviations shown by vertical lines. (d) Results of the density contrast optimizations, with optimal found values shown by squares, +/- the estimated standard deviations shown by vertical lines, and the true density contrast shown by the vertical black line. Dashed lines show the results of the alternative optimization approach tested in Tests 4 & 5.

3.2 Starting model

We construct a low-resolution bathymetry model that serves as both the initial model and the basis for generating synthetic gravity data in Tests 1–5. This model is derived by interpolating individual bathymetry point measurements (constraint points). To simulate realistic situations, we include 42 evenly spaced interior constraints representing an over-ice seismic survey, along with a surrounding border of constraints mimicking denser measurements near the ice shelf periphery, such as airborne radar over grounded ice or multibeam bathymetry offshore the calving front (Fig. 6c – small black dots). For this study area, the 42 interior points yield an average spacing of 23 km, and a median constraint proximity of 8.5 km within the interior region, closely matching the Antarctic-wide median for ice shelves larger than 1000 km² (Fig. 7c). The resulting starting bathymetry model has an RMSE of 26 m and a maximum absolute error of 164 m (Fig. 6b). Spectral analysis (Fig. 9b – black line) indicates that the interpolation error is dominated by wavelengths near 40 km, consistent with twice the point spacing.

3.3 Test 1: Best case scenario

3.3.1 Setup

This test represents an idealized, best-case scenario, establishing a baseline for inversion performance. The resulting errors are expected to be minimal and primarily reflect limitations in estimating the two key inversion hyperparameters: the regularization damping parameter and the density contrast. We generate a synthetic gravity disturbance dataset by forward-modeling the full-resolution gravity response of the real bathymetry. No noise or regional gravity component is added. To calculate the gravity misfit used in the inversion, we apply the starting bathymetry model using a prism-layer method with an incorrect density contrast of 1350 kg m⁻³, while the synthetic data is generated using the true value of 1476 kg m⁻³.

3.3.2 Test 1 Results

We identify the optimal damping value by performing cross-validation over eight candidate values (0.001 to 0.1), finding an optimal value of approximately 0.01 (Fig. 9c – blue line). Using this damping value, we then find the optimal density contrast by testing 10 values between 1000 and 2400 kg m⁻³, yielding an optimal value of 1480 kg m⁻³, just 4 kg m⁻³ above the true value (Fig. 9d – blue line). This optimized inversion (Fig. 8-1a) achieved an RMSE of 1.8 m and a maximum absolute error of 25 m relative to the true bathymetry (Fig. 9a – blue triangle). The inversion accurately recovered the true bathymetric features with only minor deviations (Fig. 8-1a and 1b). The radially averaged power spectral density (Fig. 9b – blue line) shows uniformly low errors across all wavelengths, with slight peaks near 5 and 80 km.

To assess uncertainty, we performed a Monte Carlo analysis by sampling damping and density contrast values from probability distributions (see Appendix A7). Using Latin Hypercube Sampling, we ran 20 inversions. The resulting uncertainty estimate (Fig. 8-1c) had an RMS of 1.4 m (Fig. 9a – blue square), close to the true inversion RMSE of 1.8 m. Spatial patterns of uncertainty also align well with the actual error distribution (Fig. 8-1b vs 1c), indicating a well-calibrated uncertainty estimate under ideal conditions.

3.4 Test 2: Errors in gravity data

3.4.1 Setup

In this scenario, we expand upon Test 1 by incorporating observational noise into the synthetic gravity disturbance data to simulate realistic airborne gravity survey conditions. As in Test 1, the baseline gravity disturbance is generated by forward-modeling the true bathymetry. To simulate data collection errors, we add random Gaussian noise with a standard deviation of 3 mGal (Fig. 10b), resulting in peak noise amplitudes of ∼ 12 mGal. As standard practice in airborne gravity processing, to mitigate this noise while preserving the underlying gravity signal, we applied Gaussian low-pass filters with widths ranging from 0 to 32 km (in 2 km increments).

For each filter, we computed the RMSE relative to the original noise-free gravity data (Fig. 10a). A filter width of 24 km achieved the best balance between noise reduction and signal retention, lowering the RMSE from 3 to 0.61 mGal and reducing the maximum error from ∼ 12 to ∼ 4 mGal with minimal loss of gravity signal (Fig. 10c). Following noise reduction, we proceeded as in Test 1: using an incorrect density contrast value, calculating the gravity misfit, and performing both regularization damping cross-validation and density contrast optimization.

Figure 10 Errors in the gravity data for Test 2. (a) Noise-free gravity disturbance data from forward modeling of the true bathymetry, (b) pseudo-random noise added to (a), (c) difference between (a) and (d), representing the data loss from noise and filtering, with the same colorscale as (b), (d) noisy gravity data after low-pass filtering, (e) determining the optimal low-pass filter width to reduce noise while minimizing data loss. noise-contaminated gravity disturbance data,

3.4.2 Test 2 Results

The density contrast optimization identified an optimal value of 1471 kg m⁻³ (Fig. 9d – orange line), 5 kg m⁻³ below the true value. Using this and the optimal damping value, the resulting inversion (Fig. 8-2a) produced an RMSE of 16 m and a maximum absolute error of 133 m when compared to the true bathymetry (Fig. 9a – orange triangle). This represents a notable increase in error compared to the ∼ 2 m RMSE of the idealized scenario in Test 1. Despite this, the inversion still successfully recovered all major bathymetric features, although it introduced visible short-wavelength noise into the result and increased error content in the 0–15 km wavelength range, relative to the starting model error (Fig. 9b – orange line).

For uncertainty quantification, we again used a Monte Carlo simulation. This time, in addition to sampling the damping and density contrast values, we also added a new realization of Gaussian noise (3 mGal standard deviation) to the gravity data in each iteration. After noise addition, we applied the 24 km Gaussian filter as in the main inversion. The resulting uncertainty estimate (Fig. 8-2c) had an RMS of 11 m (Fig. 9a – orange square).

3.5 Test 3: Airborne gravity survey

3.5.1 Setup

This scenario builds upon Test 1 by simulating a realistic airborne gravity survey, in contrast to the full-resolution synthetic gravity disturbance data used previously. The synthetic airborne survey follows a typical Antarctic design, found from the analysis of many published Antarctic airborne gravity surveys. From this analysis, we use flight lines spaced 10 km apart, orthogonal tie lines every 50 km, and data points sampled every 500 m along-track at a constant altitude of 1 km. The resulting gravity data have a median proximity of 2.3 km (Fig. 11e). To generate the input for inversion, we first sample the forward-modeled gravity effect of the true bathymetry (Fig. 11a) along the synthetic flight lines. We then interpolate these sampled values to a regular grid using the equivalent source method (see Appendix A8). The interpolated grid (Fig. 11c) shows an RMSE of 0.3 mGal and a maximum error of 4.3 mGal relative to the full-resolution gravity (Fig. 11a). The difference with the full-resolution gravity (Fig. 11b) highlights localized interpolation errors, primarily near the gravity high and domain edges.

To estimate the gravity data uncertainty introduced by this interpolation, we performed a Monte Carlo simulation of the equivalent source method, sampling the interpolation parameters (e.g., source depth and damping; see Appendix A7). The resulting uncertainty estimate (Fig. 11d) has a mean of 0.09 mGal, lower than the true error but capturing its general spatial distribution. Using the interpolated gravity data, we followed the same inversion workflow as in Test 1: starting with an incorrect density contrast, computing the gravity misfit, and optimizing the regularization damping and density contrast values.

Figure 11 Airborne gravity survey for Test 3. (a) Full-resolution gravity disturbance data from forward modeling of the true bathymetry (b) absolute difference between (a) and (c), (c) the equivalent source interpolated airborne gravity disturbance data, sampled along flight paths (black lines in a), (d) the estimated uncertainty of the equivalent source interpolation, and (e) the distance to the nearest gravity data point.

3.5.2 Test 3 Results

The density contrast optimization yielded an optimal value of 1464 kg m⁻³ (Fig. 9d – green line), which is 12 kg m⁻³ below the true value. With this and the optimal damping value, the inversion result (Fig. 8-3a) achieved an RMSE of 9.1 m and a maximum absolute error of 148 m relative to the true bathymetry (Fig. 9a – green triangle). The inversion successfully recovered all major bathymetric features, though it exhibited small localized deviations and greater error at wavelengths below 10 km (Fig. 9b – green line).

To assess uncertainty, we repeated the Monte Carlo analysis as in Test 1, this time including variability in three components: the damping parameter, density contrast, and equivalent source interpolation parameters. The resulting uncertainty estimate (Fig. 8-3c) had an RMS of 4.6 m (Fig. 9a – green square), and effectively captured the spatial distribution of errors observed in the final inversion (Fig. 8-3b vs 3c).

3.6 Test 4: Regional gravity field

3.6.1 Setup

This fourth scenario extends Test 1 by incorporating a regional component into the gravity disturbance data. The regional component of the gravity field is generally the long-wavelength component created from subsurface geologic density variations. To simulate this, we use crystalline basement topography from the ANTOSTRAT shipborne seismic compilation (Fig. 6e; Brancolini et al., 1995). We forward modeled its gravity effect and normalized the resulting field between 0 and 90 mGal (Fig. 12f), yielding a standard deviation of approximately 18 mGal. This value was chosen as it matches the median regional field variability for Antarctic ice shelves larger than 1000 km² (Fig. 7d).

The regional gravity was then added to the forward-modeled gravity effect of the true bathymetry to form the synthetic gravity disturbance data (Fig. 12a). As in Test 1, we calculated the starting gravity using an incorrect density contrast of 1350 kg m⁻³, and then computed the gravity misfit (Fig. 12c). We estimated and removed the regional component from this misfit using the Constraint Point Minimization technique. To interpolate the regional field from the constraint points, we tested multiple methods and found the equivalent source technique to be the most accurate, despite being more computationally demanding than minimum curvature interpolation. During the equivalent source fitting, we applied no damping and optimized the source depth using cross-validation (Fig. 12h; see Appendix A8).

Following regional field removal, the residual gravity misfit (Fig. 12d) was used to find the optimal damping parameter value. To find the optimal density contrast, we used both methods described in Sect. 2.5.2: K-Folds cross-validation with Constraint Point Minimization and a constant regional correction approach. For the cross-validation approach, the regional estimation was repeated for each fold of constraint points. The inversion was then rerun using the optimal density contrast. To isolate the performance of the regional gravity estimation from density contrast estimation errors, we also performed a separate test using the true density contrast value of 1476 kg m⁻³. This allowed us to directly assess the accuracy of the regional component removal by comparing the estimated field to the known regional input.

Figure 12 Gravity anomalies for Test 4: regional gravity component. (a) Gravity disturbance, the combination of the forward modeled gravity effects of the bathymetry and basement surfaces, (b) the forward gravity of the starting bathymetry model, (c) the gravity misfit, calculated as the difference between (a) and (b), (d) the residual component of the gravity misfit, calculated as the difference between (c) and (e), and used as the input to the inversion, (e) the regional component of the gravity misfit, estimated with the Constraint Point Minimization technique, (f) the true regional component of the misfit, which is attempted to be reproduced by (e), (g) the error in the estimated regional component, calculated as the difference between (e) and (f), (h) the determination of the optimal source depth during the fitting of equivalent sources during Constraint Point Minimization. Profile A–A', location in (a), shows the values of (c)–(f) along the profile. Distance to the nearest constraint (white dots) is shown by gray-scale color above the profile.

3.6.2 Test 4 Results

Using the true density contrast, the estimated regional component (Fig. 12e) had an RMSE of 0.47 mGal relative to the true regional field (Fig. 12f). The K-Folds cross-validated density contrast optimization yielded an optimal value of 2099 kg m⁻³ (Fig. 9d – red line), which is 623 kg m⁻³ higher than the true value. Despite this overestimation, the final inversion (Fig. 8-4a) achieved an RMSE of 7.7 m (Fig. 9a – red triangle) and a maximum absolute error of 32 m. All major bathymetric features were recovered, although errors of up to ∼ 20 m appeared in various locations. These errors predominantly occurred at wavelengths greater than 40 km (Fig. 9b – red line).

To evaluate uncertainty, we ran the same Monte Carlo simulation as in previous tests. The resulting uncertainty estimate (Fig. 8-4c) had an RMS of 5.1 m (Fig. 9a – red square), lower than the true error, but it successfully captured most regions of elevated error.

3.7 Test 5: Realistic scenario

3.7.1 Setup

This final scenario combines all components from Tests 1–4 to simulate a realistic inversion scenario, incorporating all expected sources of complexity: (1) unknown true density contrast, which must be estimated, (2) errors in gravity disturbance data, (3) gravity data acquired along airborne flight paths, and (4) a regional gravity component that must be estimated and removed. To generate the gravity disturbance dataset, we combined the full-resolution forward-modeled gravity effect of the true bathymetry with the regional gravity field produced by the basement topography. We then added 1 mGal of pseudo-random noise, sampled the data along simulated flight lines, applied a low-pass filter, and re-interpolated it onto a grid using equivalent sources (Appendix A8). The resulting gravity disturbance data are shown in Fig. 13c. All other parameters, such as the number of constraint points, gravity error levels, flight line spacing, and regional field variability, were consistent with those used in previous tests.

Figure 13 Errors in the gravity data for Test 5. (a) Noise-free full resolution gravity disturbance data from forward modeling of the true bathymetry and basement surfaces, (b) observed gravity after noise-contamination, sampling along flight lines (black lines in a), low-pass filtering, and interpolation. (c) the difference between (a) and (b), represents the data loss from noise, airborne surveying, filtering, and interpolation, (d) estimated uncertainty from equivalent source interpolation, and (e) determining the optimal low-pass filter width to reduce noise while minimizing data loss.

As in Test 2, we applied low-pass filtering to reduce short-wavelength gravity noise. Here, the noise was first added to gridded data, which were then sampled along flight paths. We applied a range of filter widths to the flight-line data and calculated the RMSE between the filtered and noise-free signals. The optimal filter width was 24 km (Fig. 13e), representing the best compromise between noise reduction and signal preservation. After interpolation, the resulting dataset had an RMSE of 1.1 mGal relative to the noise-free full-resolution model, and a maximum absolute error of 8.2 mGal. The uncertainty introduced by the interpolation was estimated via a Monte Carlo simulation of the equivalent source fitting parameters (Appendix A7). The RMS of this estimated uncertainty was 0.11 mGal (Fig. 13d), lower than the actual error but broadly capturing its spatial distribution.

With the interpolated gravity data, we followed the same steps as in Tests 3 and 4: we computed the gravity misfit using an incorrect density contrast, estimated and removed the regional component using the Constraint Point Minimization method, then performed damping cross-validation and density contrast optimization. As before, we also repeated the regional separation using the true density contrast to assess the performance of regional field estimation independently from density estimation.

Figure 14 Regional gravity estimation for Test 5; (a) the true regional gravity from the forward modeling of acoustic basement topography, (Fig. 6e), (b) errors in the regional field estimation, as the difference between (a) and (c), (c) estimated regional gravity field from Constraint Point Minimization, (d) the determination of the optimal source depth during the fitting of equivalent sources during Constraint Point Minimization. Profile A–A', location in (a), shows the gravity misfit, regional and residual components, and the true regional field along the profile. Distance to the nearest constraint (white dots) is shown by gray-scale color above profile figure.

3.7.2 Test 5 Results

Using the true density contrast, the estimated regional gravity field (Fig. 14c) had an RMSE of 1.2 mGal compared to the true regional component (Fig. 14a). Using the K-fold cross-validation approach, the density contrast optimization returned an optimal value of 2017 kg m⁻³ (Fig. 9d – purple line), which is 541 kg m⁻³ higher than the true value. Despite this large discrepancy, the inversion (Fig. 8-5a) resulted in an RMSE of 23 m with respect to the true bathymetry (Fig. 8-5b). While the major bathymetric features were well recovered, localized errors up to 100 m were introduced. Overall, the inversion improved upon the starting model, with lower errors at wavelengths greater than 20 km, but slightly higher errors at shorter wavelengths (Fig. 9b – purple line).

We assessed the uncertainty of this inversion using a Monte Carlo simulation that varied four components: damping, density contrast, gravity noise, and interpolation effects (Fig. 15). In total, 20 inversions were run, each sampling from the distributions of these parameters. The resulting bathymetry uncertainty (Fig. 8-5c) had an RMS of 10 m (Fig. 9a – purple square), compared to a true error RMS of 23 m. While the uncertainty values were lower, they correctly identified most regions of elevated error.

To understand the relative contribution of each uncertainty source, we performed additional sensitivity tests, each including only a single parameter (10 runs per test):

• Damping parameter: RMS 1 m, spatially uniform;

• Interpolation of flight-line data: RMS 4.8 m, with localized high values near topographic ridges and model edges;

• Gravity data: RMS 8.5 m, generally spatially uniform;

• Density contrast: RMS 6.9 m, especially concentrated along topographic ridges.

Among all components, uncertainty from density contrast had the largest effect and dominated the spatial distribution of the total uncertainty (Fig. 15b).

Figure 15 Sensitivity analysis for Test 5; realistic scenario. (a) Deterministic inverted bathymetry error (same as Fig. 8-5b), (b) inversion uncertainty, same as Fig. 8-5c, (c) uncertainty from the damping parameter value, (d) uncertainty from the density contrast value, (e) uncertainty from the gravity data, and (f) uncertainty from the interpolation of flight line data. All plots share the same colorscale.

3.8 Ensemble 1 – Comparing the number of constraints and the regional gravity variability

As shown in Test 4, accurately estimating the regional component of the gravity misfit is a vital part of the inversion process. The reliability of regional field estimation depends on both the variability of the regional field itself and the number and spatial distribution of constraint points within the inversion domain. To test this, we ran an ensemble of 100 synthetic inversions, covering all combinations of 10 levels of regional field variability and 10 varying amounts of constraint points. Regional field variability, quantified as the standard deviation of the regional misfit, varied linearly from ∼ 4 to ∼ 40 mGal. The number of constraint points (excluding the bordering constraint points) ranged from 4 to 650, logarithmically spaced. This range corresponds to median constraint proximities of 2.5 to 18 km, capturing the conditions encountered on most Antarctic ice shelves larger than 1000 km² (Fig. 7e – blue box).

These tests were performed in idealized conditions: full-resolution synthetic gravity data with no added noise or flight line interpolation. To reduce computational cost, we omitted the damping parameter optimization and instead used a constant damping value of 0.025 for all runs. For each inversion, we first estimate the optimal density contrast using 8 trial values ranging from 1400 to 3300 kg m⁻³, assuming an initial regional field of 0 mGal (Fig. 4b). After selecting the best-fitting density contrast, we re-estimate the regional field via Constraint Point Minimization and repeat the inversion. We evaluated the performance of each inversion using two metrics: (1) inverted bathymetry RMSE, between the inverted bathymetry and the true synthetic model, and (2) topographic improvement, the difference between the RMSE of the initial interpolated bathymetry and the RMSE of the inverted bathymetry. This metric indicates how much improvement the gravity inversion offers beyond a simple interpolation of the constraints. A score of 0 implied no net improvement, while a negative score indicates the inversion made the bathymetry worse.

Figure 16 Ensemble 1; number/proximity of constraint points vs regional gravity field variability. (a) Ensemble of inversions colored by each inverted bathymetry's RMSE, and (b) by each inversion's topographic improvement, black line is the 0 m contour. Points plotted on (a) and (b) represent the predicted values for large (>1000 km²) Antarctic ice shelves whose values were within the bounds of Ensemble 1, with values and names in the legend. Red stars and labels indicate shelves with published individual inversions. Red and blue triangles denote the constraint proximity and regional variability used in Ensemble 4 and Test 5/Ensemble 3, respectively. x-axis for (a) and (b) is on a linear scale for constraint proximity, but a log scale for the number of interior points. (c) Bathymetry RMSE vs number/proximity of constraint points grouped by regional gravity variability, each colored line corresponds to a row of (a). Solid and dashed lines are plotted on lower and upper x-axes, respectively. (d) Bathymetry RMSE vs regional variability grouped by constraint proximity, each line corresponds to a column of (a). Red lines show trend lines for the min and max slopes. (e) Histogram of predicted inverted bathymetry RMSE and topographic improvement scores for all 110 ice shelves included in the analysis.

The starting bathymetry RMSE values across Ensemble 1 ranged from 6 to 46 m, while the inverted bathymetry RMSEs ranged from 2 to 121 m. Corresponding topographic improvements ranged from −74 to +32 m. As shown in Fig. 16a, the best inversion performance was achieved in cases with many constraints (low constraint proximity) and low-variability regional gravity fields, consistent with expectations. Figure 16c further illustrates that inversion RMSE decreases exponentially with the number of constraints (solid lines). When plotted against constraint proximity (dashed lines, RMSE also decreases, though more gradually. In both cases, the effect is stronger when the regional field variability is high, indicating that dense constraints are especially important under conditions with highly variable regional fields.

Figure 16d reveals a roughly linear increase in inversion RMSE with regional field variability, for all levels of constraint proximity. The slope of this relationship ranges from 0 to 2.2 m RMSE per 1 mGal increase in regional field variability. With few constraints (i.e., high constraint proximity), inversion accuracy is highly sensitive to the variability of the regional field. In contrast, when constraint density is high, the inversion becomes less sensitive to the variability of the regional field.

Figure 16b shows that, in most scenarios, the inversion improves topography relative to simple interpolation of constraint points. However, in the upper-right corner of the parameter space, where constraint proximity is high and regional fields are highly variable, the inversion performs worse than interpolation, resulting in negative topographic improvement. This highlights that errors introduced by poor regional field estimation can outweigh the benefits of inversion in data-sparse geologically complex settings. On the other hand, with many constraints, the initial bathymetry is already relatively accurate, leaving limited room for improvement. In such cases, while the inversions still add value, the gains are modest.

We used these ensemble results to make best-case predictions of inversion RMSE and topographic improvement for 110 Antarctic ice shelves. Remaining ice shelves were not included because their median constraint proximities and/or estimated regional field variability (calculated as the standard deviation of the topography-corrected gravity disturbance) were outside the limits of the values tested in Ensemble 1. To avoid clutter, only ice shelves larger than 1000 km² are shown in Fig. 16a and b. These predictions assume perfect gravity data, i.e., full-resolution and noise-free, and should be interpreted as lower bounds of inversion error. Under these ideal conditions, we estimate that 90 % of ice shelves would yield inverted bathymetry RMSEs below ∼ 25 m (Fig. 16e). Among the largest ice shelves (>1000 km²), Thwaites and Pine Island show the lowest predicted errors (2 and 3 m, respectively), while Riiser-Larsen and Amery show the highest (121 and 83 m, respectively). In terms of topographic-improvement, we predict that 95 % of ice shelves would benefit from gravity inversion over simple interpolation (Fig. 16e). The greatest predicted improvements are for Lazarev and Getz (23 m each), while Riiser-Larsen (−75 m) and Amery (−37 m) are predicted to perform worse than interpolation, reinforcing the challenge of resolving bathymetry in regions with sparse constraints and highly-variable regional fields.

3.9 Comparing gravity noise and flight line spacing

To evaluate the relative importance of the gravity data errors (e.g., noise, leveling artifacts) and gravity data resolution (i.e., flight line spacing), we ran three ensembles of 100 inversions. Each ensemble tested 10 gravity noise levels (0–5 mGal) and 10 numbers of gravity flight lines (3–30 lines in both E–W and N–S directions), corresponding to average line spacings of 6–60 km, or median proximities to gravity data points of 0.8–8 km. All ensembles used a constant along-line observation spacing of 500 m. The three ensembles differ in regional gravity field variability; Ensemble 2 (Fig. 17 – column 1) had no regional field, Ensemble 3 (Fig. A1 – column 2), had a medium-variability regional field (18 mGal standard deviation), and Ensemble 4 (Fig. 17 – column 2) had a highly-variable regional field (41 mGal standard deviation). These variabilities span the range estimated across most Antarctic ice shelves (Fig. 7d). All ensembles used the same constraint locations and starting bathymetry as in Tests 1–5. Constraint proximities and regional field variability for these ensembles are shown as triangles in Fig. 7e.

Observed gravity was low-pass filtered with an optimal filter width estimated by the same method as Test 5. Ensemble 2 used the true density contrast, while Ensembles 3 and 4 used the density optimization procedure from Ensemble 1. To reduce computational cost, we applied a fixed damping value of 0.025.

Figure 17 Gravity noise vs quantity of gravity data. 1st column; Ensemble 2 with a no regional field and 2nd column; Ensemble 4 with a high-variability regional field (41 mGal standard deviation). (a–b) Grid cell color indicates each inversion's RMSE with the true bathymetry. Cyan line is the 26 m contour corresponding to the RMSE of the starting bathymetry model. x-axis is on a linear scale for line spacing and data proximity, but a logarithmic scale for total flight distance. (c–d) ensemble results grouped by average gravity line spacing, each colored line corresponds to a column from the top plot (a-b), (e–h) ensemble results grouped by gravity noise level, each colored line corresponds to a column from the top plot (a–b), (e–f) x-axis units are flight km on a linear scale, and (g–h) x-axis units are median gravity data proximity on a linear scale. For (e)–(f) the gray dashed lines show the starting bathymetry RMSE and for (c), (d), (g), and (f), the red lines and slope values show the line of best fit and its slope for each subplot. All rows share the same colorscales and y-axes.

3.9.1 Ensemble 2 (No regional field)

Inverted bathymetry RMSEs ranged from 3 to 49 m, with topographic improvements (relative to the 26 m RMSE starting model) from −23 to 23 m. As expected, inversions with more gravity data and lower noise performed better (Fig. 17a – lower left). All inversions with median gravity proximity below ∼ 5–6 km (line spacing <35 km) outperformed the starting model. Figure 17c shows a strong linear relationship between inversion RMSE and gravity noise for all line spacings. Slopes ranged from 5.3 m RMSE per mGal noise (many flight lines) to 1.7 m RMSE per mGal noise (few flight lines), indicating that reducing noise consistently improves inversion performance. Grouping results by noise level, we observed an exponential decay in RMSE with increasing amount of gravity data (Fig. 17e). Due to the exponential relation between total flight distance and line spacing or gravity data proximity, this equates to a roughly direct linear relation between inversion RMSE and median proximity or flight spacing (Fig. 17g). At low noise, RMSE increased ∼ 5 m per km increase in median gravity data proximity. At higher noise, this relation became more exponential: initially, added gravity data greatly improved RMSE, but additional data gave diminishing returns.

3.9.2 Ensemble 4 (High variability regional field)

Results for Ensemble 3 were similar to Ensemble 4, and are described in the Appendix A8.4. Inverted bathymetry RMSEs ranged from 12 to 33 m, with topographic improvements from −7 to 14 m. Only inversions with median gravity proximity <7 km (line spacings <50 km) showed improvement over the starting model. The effect of gravity noise was still linear (Fig. 17d), but weaker than in Ensemble 2: slopes ranged from 0.2 to 3.5 m RMSE per mGal noise. Similarly, the increase in RMSE with data proximity (Fig. 17h) was ∼ 2.5 m km⁻¹ proximity, reduced from Ensemble 2. These results indicate that a regional field dampens the benefits of additional gravity data or lower noise.

This analysis demonstrates that both reducing gravity noise and increasing flight-line density improve inversion performance. However, their relative benefits diminish in the presence of highly-variable regional fields. In low-regional-field settings, improvements from denser or cleaner data are more pronounced.

4 Discussion

Our synthetic inversions and Antarctic ice shelf analysis highlight key factors governing the accuracy of performing gravity inversions for bathymetry.

4.1 Quality of gravity data

Introducing 1 mGal (3 mGal before filtering) of random noise (Test 2) degraded bathymetry accuracy by ∼ 14 m RMSE. Though this sensitivity of gravity noise to inversion performance weakens when gravity coverage is sparse or the regional field is highly-variable (Fig. 17c–d), the relationship remains linear in all cases. In practice, this means every reduction in gravity noise directly improves inversion results, making gravity data quality a top priority, especially for dense surveys over regions with regional gravity fields with low variability, as in these situations, the dominant error source is likely gravity data errors.

4.2 Quantity of gravity data

Simulating an airborne survey (Test 3) added ∼ 7 m RMSE to the inverted bathymetry, notably lower than the effect of typical noise levels. Across Ensembles 2–4, inversion error falls exponentially with increasing total flight distance (Fig. 17e–f), particularly for regional fields with low variability. Practically, this means in regions of limited gravity data, inversion error can be greatly reduced by collecting additional gravity data. However, for regions where sufficient gravity data exists, the benefits of collecting additional data quickly diminish. The point at which there are diminishing benefits is specific to the geometry of an ice shelf, and thus, there is no specific value for optimal flight line spacing. Median gravity data proximity has a robust, nearly linear relation to RMSE in low-noise scenarios (Fig. 17g–h). This means any reduction in median gravity data proximity will result in a direct improvement to the inversion. This improvement is strongest with low noise and regional fields with low variability. Practically, this means collecting additional gravity data will offer significant improvement to inversion performance if it notably reduces overall median gravity proximity (i.e, filling a large data gap). For two shelves with similar geometries, the shelf with the less-variable regional field will benefit more from additional gravity data collection.

4.3 Importance of constraint points

Point measurements of bathymetry beneath and surrounding an ice shelf (constraint points) underpin every stage of our inversion: initial model creation, regional-field estimation, regularization, and density contrast estimation. Even with sparse, noisy gravity data, well‐distributed constraints keep inversion errors low; conversely, with few constraints, inversions only perform well if there is a regional field with exceptionally low variability. Ensemble 1 shows inversion RMSE decays exponentially with the number of constraints (Fig. 16c), and transitions from linear to exponential dependence on median constraint proximity as regional field variability increases. Practically, this means while any additional constraint point will reduce inversion error, there are only certain situations where additional constraints will offer notable improvement; (1) the addition of that constraint point will notably reduce overall median constraint proximity (i.e. filling a large data gap) or (2) a highly-variable regional gravity field is present. As with gravity flight line spacing, since each ice shelf geometry is unique, we cannot offer guidelines on the ideal number or median proximity of constraint points.

To address this, we identified which Antarctic ice shelves would likely benefit most, or least, from additional bathymetric measurements. The skewness of the constraint proximity values offers insight into the spatial distribution of existing constraints. Low or negative skew suggests that large areas of the shelf are far from existing constraints. Adding a constraint in one of these areas would reduce constraint proximity for a large portion of the ice shelf, thus notably improving inversion performance. High positive skew indicates that only a small portion of the shelf is distant from existing constraints. In these cases, adding a new constraint would have only a minor effect on median constraint proximity and inversion performance. We have calculated the skewness of constraint proximity for all large (>1000 km²) Antarctic ice shelves (Fig. 7f – x-axis). As a demonstration, we have added a single constraint to the largest data gap for two ice shelves, Cosgrove (low skew, +0.32) and Ronne-Filchner (high skew, +1.5), both with similar median constraint proximities. For Cosgrove, the additional point reduced the median constraint proximity by ∼ 1 km (∼ 10 % improvement), whereas for Ronne-Filchner, the reduction was just 6 m (∼ 0.1 % reduction) (Fig. 18).

Figure 18 The impact of adding a constraint point (green triangle) on median constraint proximity for two ice shelves; Cosgrove (a–b) and Ronne-Filchner (c–d). (a) and (c) shows the current state of each ice shelf's constraint proximity. (b) and (d) show each shelf's constraint proximity recalculated after adding a single constraint point (green triangle). Background shows MODIS-MOA imagery, red line shows ice shelf extent from MEaSURES, and black dots show existing constraint points.

From this, we infer that ice shelves with a low skewness in constraint proximity would benefit significantly from each new bathymetry measurement (left side of Fig. 7f). Those to benefit the most are likely Conger Glenzer, Rennick, Prince Harald, and Tracy Tremenchus. Conversely, shelves with high skewness will likely only benefit significantly from additional bathymetry measurements if there is a highly-variable regional field (upper right side of Fig. 7f). These include Fimbul, Ekstrom, Brunt-Stancomb, and Riiser-Larsen. Shelves which likely will not benefit significantly from additional constraints are those with high skewness and low regional fields (lower right side of Fig. 7f), such as Ronne-Filchner, Holmes, Nickerson, Sulzberger, or Thwaites.

4.4 Data collection recommendations

Our synthetic tests yield these guidelines:

• Regional gravity fields with low variability:

  • Prioritize filling large bathymetry gaps: add a few well-placed constraint points.

  • If constraints are adequate, focus on acquiring initial gravity lines, even at lower precision. Once a minimal gravity grid is in place, shift to improving gravity data quality through increasing numbers of tie lines, reducing flight speed, or re-flying noisy lines.

• Regional gravity fields with high variability:

  • Inversion errors are dominated by regional‐field misestimation, so additional gravity data or noise reduction have limited impact.

  • Instead, concentrate on bathymetry constraints until median proximity gains plateau, balancing field effort and model improvement.

We recommend using median proximity, rather than point/line spacing or count, to assess and compare spatial coverage of both gravity and bathymetric data. Calculating median constraint proximity for all ice shelves (Fig. 7a and e) has revealed some interesting findings. For example, the Ross Ice Shelf is often viewed as well-constrained concerning bathymetry measurements due to the extensive array of ∼ 223 constraint points measurements evenly distributed across the shelf from the RIGGS survey (Bentley, 1984). Yet, based on the median constraint proximity we calculate (17.3 km), it is in the worst ∼ 5 % of all Antarctic ice shelves (9th out of 164), and also has the 6th highest maximum constraint distance (62 km). This is likely partially due to the roughly circular shape of the Ross Ice Shelf. For shelves which are more elongate, while total area may be similar, all locations are closer to either the grounded ice or the open ocean, where there are more likely to be available bed/bathymetry measurements.

Future work is needed to test the optimal survey designs for constraint points. Measuring constraints on an even grid naturally provides the most uniform spatial coverage, and adding constraints at the locations of maximum data proximity is optimal for reducing median constraint proximity and reducing errors at the locations where they are highest. However, we hypothesize that focusing constraint point placement according to the gravity disturbance data, either at the maxima or minima, or at locations with high gradients, may provide improvement in estimation of the regional gravity field, and thus provide an effective means for reducing inversion errors. In addition to the ice shelves listed in the above section based on their constraint proximity skew, the following shelves have very large median constraint proximities and would benefit greatly from additional constraints: Shackleton, Roi-Baudouin, Wilkins, Ross, and Amery.

4.5 When are inversions worthwhile?

Gravity inversion generally outperforms simple interpolation of constraints, except when combining sparse constraints, high gravity noise, and highly-variable regional gravity fields. By mapping predicted RMSE and topographic improvement for each Antarctic ice shelf, we can identify ideal inversion targets (large positive improvement scores) and shelves where inversion may degrade results (large negative scores). Recently, (Charrassin et al., 2025) have performed the first continent-wide inversion for bathymetry beneath all ice shelves. For shelves which we predict to have a large positive topographic improvement score, the models produced by Charrassin et al. (2025) are likely improvements over simple interpolation of constraints (Bedmap2). In contrast, shelves with large negative scores, such as Riiser-Larsen, Amery, and Ross, should be interpreted with caution, as the inversion may have introduced more errors than it resolved. This is especially true for shelves without dedicated gravity surveys within the AntGG-2021 compilation used by (Charrassin et al., 2025) (e.g., Riiser-Larsen and Ross). Gravity inversions in Antarctica are typically restricted to the extent of ice shelves. However, from our findings of inversions outperforming interpolation in many settings, we propose that gravity inversion could offer significant improvement to our understanding of subglacial topography, or bathymetry of the open ocean.

4.6 Uncertainty and parameter sensitivity

Our Monte Carlo-based uncertainty estimates match true RMSE within ∼ 5 m for Tests 1–4 and correctly highlight most high‐error areas (e.g., steep ridges). In the semi-realistic Test 5, uncertainties are underestimated but still highlight the most uncertain regions. This provides a useful tool for determining where results are robust or guiding where future data collection efforts should be focused. For Test 5, we were able to partition the inversion uncertainty into components resulting from various parts of the inversion. The component with the lowest uncertainty (<1 m RMS) was the regularization damping parameter, owing to the accurate value found through the cross-validation optimization routine. The 3 mGal of gravity noise resulted in a bathymetry uncertainty RMS of 8.5 m. The bathymetry uncertainty caused by the interpolation of the airborne gravity data gave an RMS of 4.8 m. Lastly, the bathymetry uncertainty attributed to the choice in density contrast value resulted in an RMS of 6.9 m. While this high uncertainty is attributed to the density contrast choice, the density contrast optimization procedure is strongly affected by an incorrectly estimated regional gravity field, as was the case in Test 5. For this reason, a portion of these uncertainties should be attributed to the presence of a regional field. These insights direct efforts toward the most impactful improvements, whether refining inversion parameters, acquiring seafloor density knowledge, or densifying gravity/constraint datasets.

4.7 Hyperparameter estimation

Our cross-validation procedures consistently recovered optimal values for both the damping parameter and the density contrast (Fig. 9c–d). Damping parameter selection was robust across all scenarios, showing little sensitivity to gravity noise, sparse flight-line sampling, or the presence of a regional gravity field. Density contrast optimization also performed well in most cases, accurately recovering the true value despite added gravity noise and flight-line sampling. However, when the regional gravity field was misestimated, the optimized density contrast became biased, tending toward values that were too high. This behavior reflects the trade-off in the inversion process: a higher density contrast can achieve the same gravity misfit with a smaller modification to the topography. Consequently, in the presence of regional signals, the optimization favors larger density contrasts to avoid introducing bathymetric artifacts. Importantly, this means that the density contrast optimization does not necessarily identify the true density contrast, but instead selects the value that yields the lowest bathymetric error. This overestimation is evident in Tests 4 and 5 (Fig. 9d – red and purple boxes), where the estimated density contrast was significantly higher than the true value. Despite being biased, the optimized density contrast led to better inversion results. In Test 4, for example, using the optimized (but overestimated) density contrast resulted in a lower bathymetric RMSE of 7.7 m, compared to 10 m when using the true density contrast, the same value used to generate the synthetic gravity data. This outcome illustrates that in scenarios with a misestimated regional field, the optimization process prioritizes minimizing inversion error over recovering the physically correct parameter. A similar pattern is seen in the ensemble tests. Ensemble 2, which lacked a regional gravity field and thus did not include density contrast optimization, sometimes underperformed relative to Ensembles 3 and 4. Although Ensembles 3 and 4 had more variable regional gravity fields, their inclusion of density contrast optimization allowed them to choose the density contrast that best reduced inversion errors.

5 Conclusions

Gravity inversion is a viable technique for understanding the shape of the seafloor beneath floating ice shelves. We present an open-source gravity inversion algorithm specifically designed to model bathymetry and its uncertainty. Using synthetic tests based on real bathymetry data from the Ross Sea, Antarctica, we evaluated the algorithm's performance sensitivity to various factors. Our tests show that the estimation of the regional gravity field is a dominant source of inversion error. However, this error can be significantly reduced with sufficient point measurements of bathymetry. Interpolation errors from standard survey flight-line spacing and typical gravity data uncertainty (∼ 1 mGal after filtering) introduced additional bathymetric errors of comparable wavelength, but the magnitude of errors from gravity data uncertainty was higher. We demonstrate that our Monte Carlo-based uncertainty assessment effectively estimates both the mean error in the inverted bathymetry and its spatial variability. This approach also highlights which components of the inversion contribute most to the uncertainty, enabling a more nuanced understanding of the often large sub-ice shelf bathymetry errors. Overall, gravity inversion improves upon simple interpolation of sparse bathymetry measurements (e.g., as used in Bedmap2) in nearly all scenarios, except where there are very high variablility regional gravity fields. This indicates that the majority of Antarctic ice shelves, and possibly areas of grounded ice or open ocean, are suitable targets for gravity-based bathymetry modeling. To guide future applications, we estimate for each Antarctic ice shelf the median distance to the nearest bathymetric constraint and the local variability of the regional gravity field. This allows us to predict the best-case inversion errors and identify which shelves are likely to benefit most from additional bathymetry measurements or are currently most suitable for gravity inversions. While we provide general guidelines on data collection, the optimal amount of gravity and constraint data is shelf-specific, due to each shelf's unique geometry and the exponential relation between inversion performance and data density. Since bathymetric point measurements, typically from over-ice seismic surveying, are essential at every step in the inversion procedure, we emphasize the need for expanded efforts to collect them. More broadly, we expect these findings to inform the planning and prioritization of future field campaigns aimed at improving sub-ice shelf bathymetry models.

Appendix A

A1 Geophysical inversions

In general, inversions consist of three components;

1. Forward operator (f): This is a mathematical means of linking a physical Earth property to the expected geophysical response.

2. Physical Earth property model (p): This is a representation of the physical Earth property of interest. While real-world properties are continuous (i.e., a smoothly varying topography), computation often necessitates the discretization of the continuous properties into discrete components. The starting model (p₀) typically will include the prior geologic knowledge of the region.

3. Fit function (ϕ): This function describes the similarity between the measured (observed) data and the predicted data (resulting from the forward calculation of the property model (f(p))). The goal of the inversion is to minimize this function so that f(p) is as close to the observed data as possible while maintaining realistic characteristics.

Here, a mathematical derivation of a generalized discrete inversion problem is shown. The inversion problem can be framed as finding the set of physical parameters (p) which when modeled with the forward operator (f) produce predicted data (d^pred) as close to the observed data (d^obs), as possible. This inversion is discrete because the topography is represented with a series of grid cells, instead of a continuous function. The difference between the observed and predicted data gives the misfit m, where

$$\begin{array}{}\text{(A1)}& \mathit{m}={\mathit{d}}^{\mathrm{obs}}-{\mathit{d}}^{\mathrm{pred}}={\mathit{d}}^{\mathrm{obs}}-f\left(\mathit{p}\right).\end{array}$$

Here, the ℓ²-norm (mean squared error) of the gravity misfit is the metric used to define the closeness between the predicted and observed data, where

$$\begin{array}{}\text{(A2)}& \left|\right|\mathit{m}|{|}_{\mathrm{2}}=\sqrt{\sum _{i=\mathrm{1}}^N[{\mathit{d}}_i^{\mathrm{obs}}-f_i(\mathit{p}){]}^{\mathrm{2}}},\end{array}$$

where N is the number of observed data points.

This is defined as the fit function, ϕ(p). The fit function can also be expressed as

$$\begin{array}{}\text{(A3)}& \mathit{\varphi }\left(\mathit{p}\right)=\mathit{m}\cdot \mathit{m}=[{\mathit{d}}^{\mathrm{obs}}-f(\mathit{p}\left)\right]\cdot [{\mathit{d}}^{\mathrm{obs}}-f(\mathit{p}\left)\right].\end{array}$$

In this context, the inversion problem is to determine a vector of parameters $\stackrel{\mathrm{̃}}{\mathit{p}}$ of length M which minimizes ϕ(p). This minimum occurs where the gradient of ϕ(p) at the vector $\stackrel{\mathrm{̃}}{\mathit{p}}$ has zero length. The gradient of ϕ(p) evaluated at any p is an M-dimensional vector is defined as

$$\begin{array}{}\text{(A4)}& \mathrm{\nabla }\mathit{\varphi }\left(\mathit{p}\right)=\left[\begin{array}{c}{\displaystyle \frac{\partial \mathit{\varphi }\left(\mathit{p}\right)}{\partial {\mathit{p}}_{\mathrm{1}}}}\\ {\displaystyle \frac{\partial \mathit{\varphi }\left(\mathit{p}\right)}{\partial {\mathit{p}}_{\mathrm{2}}}}\\ \mathrm{⋮}\\ {\displaystyle \frac{\partial \mathit{\varphi }\left(\mathit{p}\right)}{\partial {\mathit{p}}_M}}\end{array}\right].\end{array}$$

Evaluating the gradient at the ith element, using Eq. (A3) gives

$$\begin{array}{}\text{(A5)}& \begin{array}{rl}{\displaystyle \frac{\partial \mathit{\varphi }\left(\mathit{p}\right)}{\partial {\mathit{p}}_i}}& {\displaystyle }={\displaystyle \frac{\partial }{\partial {\mathit{p}}_i}}\sum _j[{\mathit{d}}_j^{\mathrm{obs}}-f_j(\mathit{p}\left)\right]\cdot [{\mathit{d}}_j^{\mathrm{obs}}-f_j(\mathit{p}\left)\right]\\ {\displaystyle }& {\displaystyle }=\sum _j{\displaystyle \frac{\partial }{\partial {\mathit{p}}_i}}[{\mathit{d}}_j^{\mathrm{obs}}-f_j(\mathit{p}\left)\right]\cdot [{\mathit{d}}_j^{\mathrm{obs}}-f_j(\mathit{p}\left)\right]\\ {\displaystyle }& {\displaystyle }=-\mathrm{2}\sum _j{\displaystyle \frac{\partial f_j\left(\mathit{p}\right)}{\partial {\mathit{p}}_i}}\cdot [{\mathit{d}}_j^{\mathrm{obs}}-f_j(\mathit{p}\left)\right]\end{array}\end{array}$$

Substituting Eq. (A5 into the elements of Eq. (A4) gives

$$\begin{array}{}\text{(A6)}& \mathrm{\nabla }\mathit{\varphi }\left(\mathit{p}\right)=-\mathrm{2}\mathbf{J}\left(\mathit{p}\right)[{\mathit{d}}^{\mathrm{obs}}-f(\mathit{p}\left)\right],\end{array}$$

where p is the parameter vector of dimension M, d^obs is the observed data vector of dimension N and J is the Jacobian matrix of dimension M×N, and is given by

$$\begin{array}{}\text{(A7)}& \mathbf{J}\left(\mathit{p}\right)=\left[\begin{array}{ccc}\frac{\partial f_{\mathrm{1}}\left(\mathit{p}\right)}{\partial {\mathit{p}}_{\mathrm{1}}}& \mathrm{\dots }& \frac{\partial f_N\left(\mathit{p}\right)}{\partial {\mathit{p}}_{\mathrm{1}}}\\ \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ \frac{\partial f_{\mathrm{1}}\left(\mathit{p}\right)}{\partial {\mathit{p}}_M}& \mathrm{\dots }& \frac{\partial f_N\left(\mathit{p}\right)}{\partial {\mathit{p}}_M}\end{array}\right]\end{array}$$

Each element of the Jacobian matrix is given by

$$\begin{array}{}\text{(A8)}& {\mathbf{J}}_{ij}={\displaystyle \frac{\partial f_j\left(\mathit{p}\right)}{\partial {\mathit{p}}_i}},\end{array}$$

which is the partial derivative of the jth predicted data with respect to the ith physical parameter. The Jacobian (J) is a sensitivity matrix that describes how much the predicted data changes for an infinitesimal change in the physical parameter. For example, in the case of a gravity inversion attempting to recover the density of subsurface cubes, if the model consists of 100 cubes, and there are 10 observation points, the Jacobian would be a 100×10 matrix and J_1,1 (the first entry) would be the 1st cube's gravitational derivative with respect to density at the first observation point. The Jacobian, therefore, describes how sensitive each pair of cubes and observations is to a change in density.

Once the Jacobian is created, a matrix equation of the form Ax=b is set up, where A is the Jacobian, x is the unknown variable, and b is the data misfit.

$$\begin{array}{}\text{(A9)}& \mathbf{J}\mathit{x}=\mathit{m}\end{array}$$

A solution to this linear equation can be found directly, by finding the inverse or pseudo-inverse matrix of the Jacobian, or iteratively via a least-squares solver (Jacoby and Smilde, 2009). The solution gives the correction which, when applied to the physical parameters of interest, minimizes the misfit between the observed data and the starting model.

The derivative of the forward operator with respect to the parameter of interest determines the linearity of an inversion. For a gravity inversion attempting to recover rock density, the derivative of gravity with respect to density is constant, and thus the inversion is considered linear (Aster et al., 2018). Conversely, a gravity inversion attempting to recover the topography of a surface, or the thickness of a layer, is non-linear since the vertical derivative of gravity (derivative with respect to depth) is dependent on the depth. Non-linear inversions cannot be solved directly, with matrix inversion (Jacoby and Smilde, 2009). They must be linearized, which is the purpose of the Jacobian matrix. The Jacobian enables a minimum-norm solution to be found with an iterative solver, where at each iteration, the calculated corrections (x) are applied to the starting model, the residuals (m) are recalculated, the Jacobian is updated, and Eq. (A9) is solved again, giving a new set of corrections. This is repeated until ϕ(p) is suitably low.

A2 Gravity reduction

For a geometric gravity inversion, we need to process the gravity data (Fig. 2a) to isolate the gravity signal resulting from the topography of interest. This involves initial processing of the measured gravity data, such as drift corrections, flight line leveling, and accounting for tidal effects, aircraft maneuvers, and the effects of measuring gravity on a moving platform (Eötvös correction) (Hinze et al., 2005, 2013). This results in observed gravity data, which we take to be the signal that is produced by (1) the gravitational attraction of all massive bodies in Earth and (2) the rotation of the Earth. Observed gravity is defined as;

$$\begin{array}{}\text{(A10)}& g_{\mathrm{obs}}\left(q\right)=\mathit{\gamma }\left(q\right)+\mathit{\delta }g\left(q\right),\end{array}$$

where, for an observation point q, g_obs(q) is the observed gravity, γ(q) is the attraction of the simplified theoretical model of the Earth (the normal Earth, Fig. 1a) and δg(q), the gravity disturbance, encompasses the gravity effects of all deviations between the normal Earth and the real Earth.

To calculate the gravity effect of the normal Earth at each observation, γ(q), we use the analytical solution of the gravity effect of the WGS84 ellipsoid (Li and Götze, 2001), implemented within the Python package Boule (Fatiando a Terra Project et al., 2024a). Once removed from the observed gravity, we are left with the gravity disturbance, δg(q). For further discussion on the derivation of the gravity disturbance and the relation to the free-air anomaly, see Blakely (1995). From here on, we omit the (q) for clarity.

The gravity disturbance consists of two components;

$$\begin{array}{}\text{(A11)}& \mathit{\delta }\mathit{g}={\mathit{g}}_{\mathrm{TME}}+{\mathit{g}}_{\mathrm{subsurface}},\end{array}$$

where g_TME is the terrain mass effect and g_subsurface is the gravity effect of anomalous subsurface geologic bodies (gray bodies in Fig. 1). The terrain mass effect reflects all topographic deviations between the normal Earth and the real Earth. The normal Earth model is shown in Fig. 1a, as defined by air above the ellipsoid and crust below. A simplified model of the true Earth is shown in Fig. 1b, and contains air, ice, water, crust, and various geologic bodies within the crust. These masses are separated in Fig. 1c into components above and below the ellipsoid. From this, the masses that are anomalous with respect to the normal Earth can be distinguished, as shown in Fig. 1d. These anomalous masses include water, ice, and crust above the ellipsoid, and air, ice, water, and geologic bodies below the ellipsoid. The gravitational effect of these anomalous masses can be approximated by assuming each mass's density and setting it relative to the density of the component of the normal Earth that the mass is replacing Fig. 1d.

g_TME can be separated into its components, as defined by;

$$\begin{array}{}\text{(A12)}& {\mathit{g}}_{\mathrm{TME}}={\mathit{g}}_{\text{ice\_surface}}+{\mathit{g}}_{\text{water\_surface}}+{\mathit{g}}_{\text{rock\_surface}}.\end{array}$$

Each component of the terrain mass effect can be calculated using a layer of prisms, as shown in Fig. 1e–g.

1. g_{ice_surface}; prisms between the ice surface and the ellipsoid are assigned densities of ρ_ice−ρ_air for prisms above the ellipsoid, and ρ_air−ρ_ice for prisms below the ellipsoid (Fig. 1e).

2. g_{water_surface}; prisms between the water surface (ice base) and the ellipsoid are assigned densities of ρ_water−ρ_ice for prisms above the ellipsoid, and ρ_ice−ρ_water for prisms below the ellipsoid (Fig. 1f).

3. g_{rock_surface}; prisms between the rock surface (topography onshore and bathymetry offshore) and the ellipsoid are assigned densities of ρ_earth−ρ_water for prisms above the ellipsoid, and ρ_water−ρ_earth for prisms below the ellipsoid (Fig. 1g).

This configuration of prisms correctly discretizes the terrain mass effect because the ice surface is equal to the water surface in areas of no ice thickness, and the water surface is equal to the rock surface in areas of no water thickness. Due to this, the overlap between the three sets of prisms, shown in Fig. 1h, gives the appropriate densities, as defined in Fig. 1d. From these three sets of prisms, the terrain mass effect can be calculated at the gravity observation points using the analytical solutions of Nagy et al. (2000).

Subtracting this terrain mass effect from the gravity disturbance yields the topography-corrected disturbance, δg_TC. In theory, this the topography-corrected disturbance is the gravity effect caused by, and only by, anomalous subsurface bodies, which have a density contrast different from the assumed density contrast of the crust (gray bodies in Fig. 1). In reality, there is an additional component resulting from inaccurate topographic data used in the terrain mass effect calculation. In most geophysical applications, this component is sufficiently small to be ignored. However, for the case of a bathymetry inversion, this is the component of interest which is used as input to the inversion.

A3 Gravity misfit

In this gravity inversion, as implemented in Invert4Geom, the gravity misfit is defined as the topography-corrected disturbance (δg_TC) calculated from the initial bathymetry model, p₀.

$$\begin{array}{}\text{(A13)}& \begin{array}{rl}{\displaystyle }\mathit{m}& {\displaystyle }=\mathit{\delta }{\mathit{g}}_{\mathrm{TC}}\left({\mathit{p}}_{\mathrm{0}}\right)\\ {\displaystyle }& {\displaystyle }=\mathit{\delta }\mathit{g}-{\mathit{g}}_{\mathrm{TME}}\left({\mathit{p}}_{\mathrm{0}}\right).\end{array}\end{array}$$

Substituting Eq. (A11) for δg gives;

$$\begin{array}{}\text{(A14)}& \mathit{m}={\mathit{g}}_{\mathrm{TME}}+{\mathit{g}}_{\mathrm{subsurface}}-{\mathit{g}}_{\mathrm{TME}}\left({\mathit{p}}_{\mathrm{0}}\right),\end{array}$$

and substituting in the expanded forms of g_TME gives;

$$\begin{array}{}\text{(A15)}& \begin{array}{rl}\mathit{m}& =({\mathit{g}}_{\text{ice\_surface}}+{\mathit{g}}_{\text{water\_surface}}+{\mathit{g}}_{\text{rock\_surface}})\\ & +{\mathit{g}}_{\mathrm{subsurface}}-({\mathit{g}}_{\text{ice\_surface}}+{\mathit{g}}_{\text{water\_surface}}+f({\mathit{p}}_{\mathrm{0}}\left)\right)\\ & ={\mathit{g}}_{\mathrm{subsurface}}+{\mathit{g}}_{\text{rock\_surface}}-f\left({\mathit{p}}_{\mathrm{0}}\right).\end{array}\end{array}$$

where f(p₀) is the forward gravity effect of the initial bathymetry model. This starting bathymetry model should contain all available knowledge of the density distribution of the seafloor and its depth. If there is no available information on depth, a flat starting model must be created at a reasonable depth. If there is no available information on the seafloor density, a reasonable value must be chosen, or the optimal density value can be inferred from the data, as described in Sect. 2.5. Note that in geometric inversions for subsurface layers (Moho, basement), there is an additional parameter that must be chosen, the reference level for the prism layer. For bathymetry inversions, since the input anomaly is the topography-corrected disturbance, the reference level should be set to 0 meters (the ellipsoid surface).

Before the inversion, we need to separate the gravity misfit, m, into a regional, m_reg, and a residual component m_res. m_reg is equivalent to g_subsurface, and is the gravity effect of subsurface geologic variations (gray bodies in Fig. 1). m_res is equivalent to g_{rock_surface}−f(p₀), defined as the gravity effect resulting from inaccuracies between the true bathymetry and the starting bathymetry model.

A4 Discretization and edge effects

There are several options for discretizing a topographic surface, representing a density contrast between two materials, as a series of adjacent vertical right-rectangular prisms (Fig. 3). The most conceptually simple method is two use two layers of prisms, one to represents the true density and geometry of the material on either side of the surface (Fig. 3b). In this method, the topographic surface is represented by the bottoms of the upper prism layer and the tops of the lower prism layer. Alternatively, the density contrast can be represented by a single layer of prisms either above or below the surface, ending at an arbitrary reference, typically the minimum or maximum value of the surface (Fig. 3c–d). For this, the prism’s density values are relative to the medium on the other side of the surface (i.e., if prisms are below the surface, their densities are ρ_below−ρ_above). The last option is to choose a flat reference level (z_ref; typically the mean value of the surface or sea level) and create one layer of prisms between this and the surface (Fig. 3e). Prisms above the reference level are assigned a positive density contrast (Δρ=ρ_below−ρ_above), and prisms below the surface a negative density contrast (Δρ=ρ_above−ρ_below). All of these methods result in similar forward gravity calculations. Note that the absolute density option (Fig. 3b) has twice as many prisms as the other options, significantly increasing its computational expense, regardless of whether or not a buffer zone is used.

Any forward gravity calculation of a layer of prisms will include gravitational edge effects. These effects are typically a decay in the gravitational field towards the edges of the prism model due to the void space (density of 0) beyond the edge of the model. The magnitude of decay is dependent on both the average thickness and density of prisms near the edge. Thicker and or denser prisms create a larger contrast with the void space, resulting in a greater amount of decay. Therefore, the absolute density method of discretization (Fig. 3b) produces the largest edge effects due to the high-density values and large total thickness of the prisms. The relative density methods (Fig. 3c–d) produce a smaller edge effect, since the prisms are only above or below the surface, reducing their total thickness. The density values are also less than in the absolute density method. The reference surface method (Fig. 3e) produces the smallest edge effects due to both (a) the mean prism thickness being smaller (depending on the chosen reference level value), and (b) positive and negative density prisms having opposite edge effects, partially canceling each other out. Note that with a reference level completely above or below the surface, this reference method becomes identical to methods (c) and (d) of Fig. 3.

A buffer zone can be used to reduce the impact of these edge effects by calculating forward gravity confined to an inner region, while the prism model extends beyond. This limits the majority of the edge effects to outside the region of calculation. Too large of a buffer zone will add an unnecessary amount of computation (more prisms), while too small a zone will introduce unacceptably large edge effects to the calculations.

A5 Annulus approximation

A vertical prism can be approximated as a sector of an annulus (a cylindrical shell). This approximation greatly simplifies the calculation of the vertical derivative of the prisms' gravity. The vertical component of gravity of an entire annulus is defined by

$$\begin{array}{}\text{(A16)}& \mathit{\delta }g=\mathrm{2}\mathit{\pi }G\mathit{\rho }(R-r+\sqrt{r^{\mathrm{2}}+h^{\mathrm{2}}}-\sqrt{R^{\mathrm{2}}+h^{\mathrm{2}})}\end{array}$$

where ρ is the density, G is the gravitational constant, r and R are the inner and outer radii, respectively, and h is the average height of the cylinder (Hammer, 1939). r, R, and h are relative to the observation point, located in the center bottom of the annulus. McCubbine (2016) provides an adaptation of Eq. (A16) which isolates the gravity of just a portion of this annulus. This is accomplished by adding factor f, relating the area of the annulus to the area of a sector of it:

$$\begin{array}{}\text{(A17)}& f={\displaystyle \frac{{\mathit{\alpha }}^{\mathrm{2}}}{\mathit{\pi }(R^{\mathrm{2}}-r^{\mathrm{2}})}}\end{array}$$

where

$$\begin{array}{}\text{(A18)}& {\displaystyle }r=\sqrt{x^{\mathrm{2}}+y^{\mathrm{2}}}-\sqrt{{\displaystyle \frac{{\mathit{\alpha }}^{\mathrm{2}}}{\mathrm{2}}}}\text{(A19)}& {\displaystyle }R=\sqrt{x^{\mathrm{2}}+y^{\mathrm{2}}}+\sqrt{{\displaystyle \frac{{\mathit{\alpha }}^{\mathrm{2}}}{\mathrm{2}}}}\end{array}$$

where α is the prism width, x, y, and z are coordinates of the prism center, relative to the observation point. This factor f is then multiplied by Eq. (A16) to give the vertical component of gravity for a section of an annulus, δg_annulus=fδg.

Taking the derivative of this with respect to h gives

$$\begin{array}{}\text{(A20)}& {\displaystyle \frac{\partial g}{\partial h}}=\mathrm{2}f\mathit{\pi }G\mathit{\rho }h\left({\displaystyle \frac{\mathrm{1}}{\sqrt{r^{\mathrm{2}}+h^{\mathrm{2}}}}}-{\displaystyle \frac{\mathrm{1}}{\sqrt{R^{\mathrm{2}}+h^{\mathrm{2}}}}}\right).\end{array}$$

While the approximation of a prism as a section of an annulus introduces some errors, the calculation is efficient and simple to implement.

A6 Uncertainty quantification

During the calculation of cell-wise statistics of the inversion ensembles resulting from the Monte Carlo analysis, the statistics can be weighted by the inverse square of each inverted bathymetry's root mean square difference with the constraint point depths, if they are available (Fig. 5d). This reduces the bias of inversions in the ensemble, which performed poorly, as defined by not adhering to the constraints (Schnaidt and Heinson, 2015).

A7 Estimating parameter probability distributions

A7.1 Test 1: Best case scenario

In Test 1, two parameters were included in the Monte Carlo analysis: (1) the damping parameter and (2) the density contrast. For the damping parameter, we use a standard deviation (0.24) calculated from half of the standard deviation of the damping values in the cross-validation (Fig. 9c – blue), omitting any outliers ($\left|Z\mathrm{score}\right|>\mathrm{2}$). Similarly, the density contrast optimization resulted in a well-defined best solution, so for its probability distribution we manually chose a standard deviation of 5 kg m⁻³ and a mean of the best-found value, 1480 kg m⁻³ (Fig. 9d – blue).

A7.2 Test 2: Errors in gravity data

In Test 2, three parameters were included in the Monte Carlo analysis: (1) the damping parameter, (2) the density contrast, and (3) the gravity data points. For the damping parameter, we use a standard deviation (0.24) calculated from half of the standard deviation of the damping values in the cross-validation (Fig. 9c – orange, omitting any outliers ($\left|Z\mathrm{score}\right|>\mathrm{2}$)). As in Test 1, we manually chose a standard deviation of 5 kg m⁻³ for the density contrast (Fig. 9d – orange). Each gravity observation point is sampled from a probability distribution with a mean of the point's value (contaminated with noise) and a standard deviation of 3 mGal.

A7.3 Test 3: Airborne gravity survey

In Test 3, three parameters were included in the Monte Carlo analysis: (1) the damping parameter, (2) the density contrast, and (3) the gravity data points. For the damping parameter, we use a standard deviation (0.24) calculated from half of the standard deviation of the damping values in the cross-validation (Fig. 9c – green), omitting any outliers ($\left|Z\mathrm{score}\right|>\mathrm{2}$). We manually chose a standard deviation of 20 kg m⁻³ for the density contrast (Fig. 9d – green). Each gravity observation point is sampled from a probability distribution with a mean of the point's value (the equivalent source gridded value) and a standard deviation. The standard deviation is unique for each point and was determined from a Monte Carlo analysis of the equivalent source gridding technique. For this, the source depth and damping parameters were each sampled from their respective probability distributions and an ensemble of interpolated gravity data was created. The optimal values for source depth and damping found from the equivalent source gridding (Appendix A8) were used as the means of each probability distribution, ∼ 10 km and ∼ 0.002, respectively. A standard deviation of 2.5 km (25 % of 10 km) was used for the source depth, and a standard deviation of ∼ 3 (the absolute value of the optimal damping value's base 10 exponent) was used for the base 10 exponent of the damping value. An ensemble of 10 equivalent source interpolated gravity disturbance grids was created with the source depth and damping values sampled from their uncertainties. This resulted in a gravity uncertainty, shown in Fig. 11d, with a mean of ∼ 0.1 mGal and a maximum value of 0.8 mGal. Comparing this with the true error in the gravity data (Fig. 11b) shows that this uncertainty technique performed well for estimated both the mean value and the spatial distribution of the gravity errors.

A7.4 Test 4: Regional gravity field

In Test 4, two parameters were included in the Monte Carlo analysis: (1) the damping parameter and (2) the density contrast. For the damping parameter, we use a standard deviation (0.24) calculated from half of the standard deviation of the damping values in the cross-validation (Fig. 9c – red), omitting any outliers ($\left|Z\mathrm{score}\right)|>\mathrm{2}$). The density contrast optimization resulted in a poorly defined best solution (Fig. 9d – red), so we use a large standard deviation of 500 kg m⁻³ for sampling the density contrast values.

A7.5 Test 5: Realistic scenario

In Test 5 three parameters were included in the Monte Carlo analysis; (1) the damping parameter, (2) the density contrast, and (3) the gravity data. The damping parameter was assigned an standard deviation of 0.24, calculated from half the standard deviation of the cross-validation results (Fig. 9c – purple), the density contrast was assigned a standard deviation of 500 kg m⁻³, the gravity data were assigned a constant standard deviation of 3 mGal.

A8 Equivalent source interpolation

A8.1 Test 3: Airborne gravity survey

In Test 3, the full-resolution gravity data is sampled along the paths of an airborne survey. We then use equivalent sources to interpolate this data back to a full grid. We block-mean the data at a 2 km spacing to reduce aliasing effects due to oversampling along flight lines (Uieda, 2018). Two parameters are needed to fit the equivalent sources: the source depth and a damping value. The optimal value for each of these can be estimated through a cross-validated optimization. For each set of parameter values, a five-fold cross-validation is performed where the scoring metric is the mean of the five scores, and each score is the RMSE between the gridded gravity data and the fold's omitted gravity data. To find the set of parameter values that gives the optimal score, we perform 12 trials of parameter values with dampings between $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{3}}$ and 10 and source depths (relative to observation heights) between 100 m and 100 km. This resulted in an optimal source depth of ∼ 10 km and a damping of ∼ 0.002. We use these fitted sources to predict the gravity data across the entire grid, with a 2 km spacing (same as the spacing of the forward-modeled gravity effect of the true bathymetry).

A8.2 Test 4: Regional gravity field

In Test 4, we introduced a regional component to the observed gravity data. To estimate and remove this component, we use Constraint Point Minimization, which requires interpolation of the gravity data. For this, we use the equivalent source technique and find the optimal source depth with the same cross-validated optimization as above. For this application, we use no damping. We perform 20 trials of source depth values (relative to observation heights) between 100 m and 400 km. We use these fitted sources to predict the gravity data across the entire grid, with a 2 km spacing (Fig. 12d). Once the optimal density contrast value was found, we repeated this procedure and found the optimal equivalent source depth to be ∼ 214 km.

A8.3 Test 5: Realistic scenario

In Test 5, the same cross-validated optimization as Test 3 was used to determine optimal parameter values for interpolating the flight line data. This resulted in an optimal source depth and damping of 57 km and $\mathrm{1.54}\times {\mathrm{10}}^{-\mathrm{20}}$, respectively. With the optimal density value, the optimal equivalent source depth and damping parameters for the regional estimation, found with the same technique as Test 4, were found to be ∼ 10 km and 0.0002, respectively.

A8.4 Ensemble 3 (Medium variability regional field)

Inverted bathymetry RMSEs ranged from 9 to 31 m, with topographic improvements from −5 to 17 m. Only inversions with median gravity proximity <8 km (line spacings <50 km) showed improvement over the starting model. The effect of gravity noise was still linear (Fig. A1e), but weaker than in Ensemble 2: slopes ranged from 0.3 to 4.2 m RMSE per mGal noise. Similarly, the increase in RMSE with data proximity (Fig. A1k) was ∼ 2.6 m per km proximity, reduced from Ensemble 2.

Figure A1 Ensembles 2–4; gravity noise vs quantity of gravity data, with varying levels of regional field variability. 1st column; Ensemble 2 with a no regional field, 2nd column; Ensemble 3 with a medium-variability regional field (18 mGal standard deviation), and 3rd column; Ensemble 4 with a high-variability regional field (41 mGal standard deviation). 1st row; (a–c) Grid cell color indicates each inversion's RMSE with the true bathymetry. Cyan line is the 26 m contour corresponding to the RMSE of the starting bathymetry model. x-axis is on a linear scale for line spacing and data proximity, but a logarithmic scale for total flight distance. 2nd row; (d–f) ensemble results grouped by average gravity line spacing, each colored line corresponds to a column from the top plot (a–c), 3textsuperscriptrd and 4th rows; (g–l) ensemble results grouped by gravity noise level, each colored line corresponds to a column from the top plot (a–c), (g–i) x-axis units are flight km on a linear scale, and (j–l) x-axis units are median gravity data proximity on a linear scale. For rows 2–4 the gray dashed lines show the starting bathymetry RMSE and for rows 2 and 4, the red lines and slope values show the line of best fit and its slope for each subplot. All rows share the same colorscales and y-axes.