On the Feasibility of Seaﬂoor Topography Estimation from Airborne Gravity Gradients: Performance Analysis Using Real Data

: Compared with airborne gravimetry, a technique frequently used to infer the seaﬂoor topography at places inaccessible to ship soundings due to the presence of ice shelf or ice m é lange, airborne gravity gradiometry inherently could achieve higher spatial resolution, thus it is promising for improved inference of seaﬂoor topography. However, its estimation capability has not been demonstrated by real projects. Theoretical analysis through admittance shows that compared with gravity disturbance, gravity gradient is more sensitive to the short-wavelength seaﬂoor topography but diminishes faster with the increase of the distance between the seaﬂoor and airplane, indicating its superiority is recovering short-wavelength topographic features over shallow waters. We present the ﬁrst numerical experiment that estimates seaﬂoor topography from a 0.4-km resolution, real airborne gravity gradients. It is shown that airborne gravity gradiometry can recover smaller topographic features than typical airborne gravimetry, but the estimation accuracy is only ± 17 m due to the presence of subsurface density variations. The long-wavelength e ﬀ ect of the subsurface density variations can be removed with the aid of constraining bathymetry inside the study area, whereas the short wavelengths cannot. This study expands the applications of airborne gravity gradiometry, and helps glaciologists understand its performance in seaﬂoor topography estimation. gradient signal over that caused by the geometric conﬁguration of the seaﬂoor. It is shown that airborne gravity gradiometry can recover small topographic features that are outside of the resolution capability of typical airborne gravimetry, but the estimation accuracy is only ± 17 m due to the presence of large non-terrain e ﬀ ects. The long-wavelength non-terrain e ﬀ ect can be removed if sparse bathymetric data are available inside the study area to provide constraints, whereas the short wavelengths cannot. This suggests that future airborne gravity gradiometry survey should overlap some known bathymetry or land surface elevations, and that more advanced non-terrain e ﬀ ect removing techniques should be investigated in the future to fully realize the potential of the airborne gravity gradients at short wavelengths. We remain optimistic about seaﬂoor topography estimation from airborne gravity gradients, since at some places, the subsurface density variation is signiﬁcant while at some other places the density of terrain is almost constant, where the non-terrain e ﬀ ect is small. Future studies at places where the density of terrain is almost uniform are needed to fully demonstrate the estimation capability of airborne gravity gradients. This study is an important starting point for future work that constantly improves glaciologists’ understanding of the advantage and performance of airborne gravity gradients in seaﬂoor topography estimation. This study could help glaciologists make better decisions when selecting gravity measurement devices for future airborne geophysical missions in polar areas.


Introduction
The airborne gravimetry technique became mature in 1990s, when the appearance of GPS enabled easy access to accurate kinematic acceleration [1]. During the past three decades, airborne gravimetry has been used extensively by the glaciological communities to infer the shape of subglacial lakes [2], the seafloor topography beneath ice shelves and floating ice tongues [3][4][5][6][7][8][9], and the thickness of mountain glaciers [10], etc. At places inaccessible to ship soundings or ice-penetrating radars, airborne gravimetry has become an indispensable tool to infer seafloor topographies that are critical to study ice-ocean interactions, stability of marine-terminating glaciers, etc.

Resolution
The acceleration sensed by an airborne gravimeter is specific force. In the inertial frame, it is the combination of the gravitational acceleration and the kinematic acceleration of the aircraft [15]. Therefore, in addition to the airborne gravimeter, the accuracy of the measured gravity is also affected by the kinematic acceleration derived from aircraft tracking (typically using differentially corrected GPS). Limited by the errors in modeling the numerical derivative and the positioning capabilities of current global navigation satellite systems, the kinematic acceleration can achieve an accuracy of sub-mGal level only with smoothing over intervals of 60s or longer [16]. On the contrary, airborne gravity gradiometry does not need to determine the kinematic acceleration of the aircraft. As a result, the resolution of airborne gravity gradients (typically 0.5 km) is higher than the resolution of airborne gravity measured by a fixed-wing aircraft (typically 2 km). This implies that airborne gravity gradiometry is superior in detecting small topographic features, which is particularly desirable in studies at the individual glacier scale.

Admittances
Parker's infinite series [17] is a formula to calculate the gravitational acceleration caused by a non-uniform volume of material. The regional isostatic compensation mechanism [18] treats the lithosphere as a thin elastic plate that flexes under the load (i.e., topography). Assuming that the measured gravity is only caused by the topography and the corresponding flexure of the lithospheric plate at the Mohorovicic discontinuity (the density interface between the crust and mantle) [19], we combine the first-term approximation of Parker's series and the regional isostatic compensation mechanism to establish the relationship between gravity and topography in the spectral domain.
Remote Sens. 2020, 12, 4092 3 of 14 Assuming the density of the crust is uniform, the gravitational acceleration at sea surface, δg, caused by b, the deviation of the seafloor from the mean sea depth d, is given by (detailed derivations can be found in pp. 23-26 of [20]): where F denotes 2-D Fourier transform, G is the gravitational constant, ∆ρ is the density difference between the crust and sea water, f is the radial frequency, t is the mean thickness of the crust, and Φ e f is the spectral domain function that modifies the Airy response so as to produce the flexure of the lithosphere [19]: where g is the gravitational acceleration, ρ m and ρ c are the density of the mantle and crust, respectively, and D is the flexural rigidity of the plate given by: where E is Young's modulus, υ is Poisson's ratio, and T e is the elastic thickness of the plate. The vertical-vertical gravity gradient (VGG) at sea surface, Γ zz , caused by b, is given by: In Equation (1), the spectral domain parameter that modifies the seafloor topography so as to produce the gravity disturbance is called the gravitational admittance, Z δg f . The theoretical gravitational admittance is: Similarly, based on Equation (4) the gravity gradient admittance Z Γ zz f is defined as: In summary, under certain assumptions (constant density, constant-height measurement surface above all masses, linear approximation), the gravity disturbance and gravity gradients can be approximately related to the seafloor topography in the spectral domain through the admittances: Next, we use the theoretical gravitational admittance and gravity gradient admittance to analyze the advantages and limitations of using the gravity gradients to estimate the seafloor topography. The theoretical admittances (5) and (6) are computed using the values in Table 1 and are shown in Figures 1 and 2. wavelengths. It is one of the reasons why predicting seafloor topography from gravity disturbances based on their theoretical relationship is feasible only at a certain short-to-intermediate, uncompensated wavelength band [22]. As a result, a long-wavelength regional gravity should be removed prior to the inversion when predicting seafloor topography from airborne gravity disturbances. In previous studies, the regional gravity is empirically computed through various methods, such as spherical harmonic models [5,23,24], upward continuation [25,26], matched filters [27], etc. However, a large discrepancy exists between the regional gravity computed by different methods, so does the corresponding gravity-disturbance-estimated seafloor topography. Along with the uniform density assumption to the bedrock, the uncertainty in the regional gravity is one of the largest error sources in gravity-based seafloor topography estimation [28]. The bottom panel of Figure 1 shows that the isostatic compensation has smaller effects on the gravity gradient admittance than on the gravitational admittance. Besides, due to the −2 ̅ factor in the gravity gradient admittance, i.e., Equation (6), the VGG is less sensitive to the long-wavelength topography than the gravity disturbance. This indicates that when predicting seafloor topography from airborne gravity gradients, removing the long-wavelength regional gravity gradient is not necessary, so the uncertainties associated with the computation of regional gravity is avoided. On the other hand, the gravity gradient admittance reaches its maximum at wavelength 75 km, whereas the peak of the gravitational admittance is at the longer wavelengths. This difference is even larger if we use a smaller mean sea depth (see Figure 2). This indicates that gravity gradient has higher relative magnitude at short wavelengths than gravity disturbance. From another perspective, compared to the gravity disturbance, the gravity gradient is more sensitive to short-wavelength The top panel of Figure 1 shows that at long wavelengths, the admittance function is sensitive to the compensation scheme (i.e., the choice of the effective elastic thickness of the plate T e ). As can be found by comparing the uncompensated (solid curve) and compensated (dashed curves) admittances in the top panel of Figure 1, at long wavelengths, the gravity disturbance generated by seafloor topography is partially cancelled out by isostatic compensation, which in fact is uncertain to a large extent. Besides, the density variations within the mantle, if they exist, affect the gravity disturbance at the sea surface mainly in the long wavelengths because the short-wavelength gravity disturbance is attenuated in the upward continuation [21]. In summary, the theoretical relationship between gravity disturbance and seafloor topography contains large uncertainties at long wavelengths. It is one of the reasons why predicting seafloor topography from gravity disturbances based on their theoretical relationship is feasible only at a certain short-to-intermediate, uncompensated wavelength band [22]. As a result, a long-wavelength regional gravity should be removed prior to the inversion when predicting seafloor topography from airborne gravity disturbances. In previous studies, the regional gravity is empirically computed through various methods, such as spherical harmonic models [5,23,24], upward continuation [25,26], matched filters [27], etc. However, a large discrepancy exists between the regional gravity computed by different methods, so does the corresponding gravity-disturbance-estimated seafloor topography. Along with the uniform density assumption to the bedrock, the uncertainty in the regional gravity is one of the largest error sources in gravity-based seafloor topography estimation [28]. element and the measurement point, while gravity disturbance attenuates with the square of distance. So, for deep ocean areas, the magnitude of the gravity gradient at the sea surface may be small, which then also implies a small signal-to-noise ratio. For these areas, the superiority of the gravity gradient at short wavelengths may be cancelled by the low signal-to-noise ratio because the attenuation rate of gravity disturbance is slower. The variation of admittance with respect to mean sea depth is shown in Figure 2. As the mean sea depth increases, the gravity gradient admittance decreases fast, and the peak of the admittance gradually moves to longer wavelengths. Although the gravitational admittance also decreases when mean sea depth increases, the change rate is comparatively slower. Figure 2. The theoretical gravitational admittance (top) and gravity gradient admittance (bottom) with respect to mean sea depth .

Numerical Experiments
As analyzed in Section 2, compared with airborne gravity disturbance, airborne gravity gradient avoids uncertainties associated with the regional gravity, is more sensitive to the short-wavelength seafloor topography, and thus is considered superior when predicting seafloor topography over shallow waters. This section uses high-resolution (0.4 km) airborne gravity gradient data in a small region (12.6 km × 18.4 km) to infer the seafloor topography, for the purpose of investigating its ability to detect small topographic features in offshore areas.

Airborne Gravity Gradients and Multibeam Bathymetry
In December 2012, Bell Geospace Inc. (BGI), Houston, TX, USA conducted a gravity gradiometry survey over St. George's Bay [29]; see Figure 3a. It is located off the southwestern shore of Newfoundland, Canada. The FTG full tensor gradiometer is installed on a Basler Turbo BT-67 aircraft, which was flown to gently drape the land topography in order to minimize the distance between measurement altitude and ground. The reported accuracy of the FTG is about 5 to 7 Eötvös [29]. The The bottom panel of Figure 1 shows that the isostatic compensation has smaller effects on the gravity gradient admittance than on the gravitational admittance. Besides, due to the −2π f factor in the gravity gradient admittance, i.e., Equation (6), the VGG is less sensitive to the long-wavelength topography than the gravity disturbance. This indicates that when predicting seafloor topography from airborne gravity gradients, removing the long-wavelength regional gravity gradient is not necessary, so the uncertainties associated with the computation of regional gravity is avoided.
On the other hand, the gravity gradient admittance reaches its maximum at wavelength 75 km, whereas the peak of the gravitational admittance is at the longer wavelengths. This difference is even larger if we use a smaller mean sea depth (see Figure 2). This indicates that gravity gradient has higher relative magnitude at short wavelengths than gravity disturbance. From another perspective, compared to the gravity disturbance, the gravity gradient is more sensitive to short-wavelength topography. Because the prediction band is in the short-to-intermediate wavelengths, we may conclude that the gravity gradient is superior to the gravity disturbance in seafloor topography estimation.
One should note that the gravity gradient diminishes with the cube of distance between the mass element and the measurement point, while gravity disturbance attenuates with the square of distance. So, for deep ocean areas, the magnitude of the gravity gradient at the sea surface may be small, which then also implies a small signal-to-noise ratio. For these areas, the superiority of the gravity gradient at short wavelengths may be cancelled by the low signal-to-noise ratio because the attenuation rate of gravity disturbance is slower. The variation of admittance with respect to mean sea depth is shown in Figure 2. As the mean sea depth increases, the gravity gradient admittance decreases fast, and the peak of the admittance gradually moves to longer wavelengths. Although the gravitational admittance also decreases when mean sea depth increases, the change rate is comparatively slower.

Numerical Experiments
As analyzed in Section 2, compared with airborne gravity disturbance, airborne gravity gradient avoids uncertainties associated with the regional gravity, is more sensitive to the short-wavelength seafloor topography, and thus is considered superior when predicting seafloor topography over shallow waters. This section uses high-resolution (0.4 km) airborne gravity gradient data in a small region (12.6 km × 18.4 km) to infer the seafloor topography, for the purpose of investigating its ability to detect small topographic features in offshore areas.

Airborne Gravity Gradients and Multibeam Bathymetry
In December 2012, Bell Geospace Inc. (BGI), Houston, TX, USA conducted a gravity gradiometry survey over St. George's Bay [29]; see Figure 3a. It is located off the southwestern shore of Newfoundland, Canada. The FTG full tensor gradiometer is installed on a Basler Turbo BT-67 aircraft, which was flown to gently drape the land topography in order to minimize the distance between measurement altitude and ground. The reported accuracy of the FTG is about 5 to 7 Eötvös [29]. The designed (nominal) flight altitude during the entire survey is 80 m above the water. The flight lines are 400 m apart, with data sampled along-track at a one second interval, which translates to a 59.7-m spatial sampling interval using the 215 km/h speed. As part of the data processing, BGI projected the gravity gradients onto the Universal Transverse Mercator (UTM) zone 21 N.
Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 14 designed (nominal) flight altitude during the entire survey is 80 m above the water. The flight lines are 400 m apart, with data sampled along-track at a one second interval, which translates to a 59.7-m spatial sampling interval using the 215 km/h speed. As part of the data processing, BGI projected the gravity gradients onto the Universal Transverse Mercator (UTM) zone 21 N. High horizontal resolution (10 m × 10 m) multibeam bathymetry data for this area were obtained from the Geological Survey of Canada (Atlantic) [30]; see Figure 3b. The reported accuracy of the multibeam echo sounding is usually in the level of 10 cm [31]. Besides, the accuracy of the multibeam sounding is generally evaluated according to the International Hydrographic Organization Standards for Hydrographic Surveys [32], which requires the uncertainty of the reduced depths to be smaller than 0.25 + 0.0075 × meter, where b is the seafloor depth (the mean depth of the study area is 38.56 m). It is much smaller than the accuracy of the gravity estimated seafloor topography. Therefore, the bathymetric data are treated as true depths to evaluate the accuracy of the gravitygradient-inferred seafloor topography.

Inversion Methods
The seafloor topography is estimated from airborne gravity gradients using very fast simulated annealing (SA) [33], a global optimization method that can process high-dimensional nonlinear inversion problems. The SA algorithm iteratively searches over a set of candidate seafloor topographies. During each iteration, the gravity gradients are forward computed from seafloor topography and compared with the measured gravity gradients. The SA algorithm can asymptotically find an optimal seafloor topography that minimizes the discrepancy between the measured and forward-computed gravity gradients. Interested readers are referred to [8,34] for a more detailed description of the SA method. In the remainder of this subsection, we only describe the key steps of SA, namely the forward model and the cost function (defined in the next paragraph).
We discretize the study region (inside the dashed lines shown in Figure 3b) into adjacent right rectangular prisms with a horizontal dimension of 400 m × 400 m. Each prism is assigned a uniform density representing either sea water (1030 kg/m 3 ) or bedrock (2200 kg/m 3 [29]). The gravity gradient at flight altitude, treated as the effect of mass deficiency caused by the sea water between the sea surface and the seafloor, is forward computed from the right rectangular prisms model using Nagy et al.'s formula [35]. For each computation point, only the prisms within 23 km of horizontal distance from the point are used in the forward computation. According to formulas given by Jekeli [36], the errors introduced by neglecting the masses outside the 23-km inner zone is 0.5 Eötvös, which is one order smaller than the accuracy of the gravity gradient measurements. Note that the seafloor High horizontal resolution (10 m × 10 m) multibeam bathymetry data for this area were obtained from the Geological Survey of Canada (Atlantic) [30]; see Figure 3b. The reported accuracy of the multibeam echo sounding is usually in the level of 10 cm [31]. Besides, the accuracy of the multibeam sounding is generally evaluated according to the International Hydrographic Organization Standards for Hydrographic Surveys [32], which requires the uncertainty of the reduced depths to be smaller than 0.25 + 0.0075 × b meter, where b is the seafloor depth (the mean depth of the study area is 38.56 m). It is much smaller than the accuracy of the gravity estimated seafloor topography. Therefore, the bathymetric data are treated as true depths to evaluate the accuracy of the gravity-gradient-inferred seafloor topography.

Inversion Methods
The seafloor topography is estimated from airborne gravity gradients using very fast simulated annealing (SA) [33], a global optimization method that can process high-dimensional nonlinear inversion problems. The SA algorithm iteratively searches over a set of candidate seafloor topographies. During each iteration, the gravity gradients are forward computed from seafloor topography and Remote Sens. 2020, 12, 4092 7 of 14 compared with the measured gravity gradients. The SA algorithm can asymptotically find an optimal seafloor topography that minimizes the discrepancy between the measured and forward-computed gravity gradients. Interested readers are referred to [8,34] for a more detailed description of the SA method. In the remainder of this subsection, we only describe the key steps of SA, namely the forward model and the cost function (defined in the next paragraph).
We discretize the study region (inside the dashed lines shown in Figure 3b) into adjacent right rectangular prisms with a horizontal dimension of 400 m × 400 m. Each prism is assigned a uniform density representing either sea water (1030 kg/m 3 ) or bedrock (2200 kg/m 3 [29]). The gravity gradient at flight altitude, treated as the effect of mass deficiency caused by the sea water between the sea surface and the seafloor, is forward computed from the right rectangular prisms model using Nagy et al.'s formula [35]. For each computation point, only the prisms within 23 km of horizontal distance from the point are used in the forward computation. According to formulas given by Jekeli [36], the errors introduced by neglecting the masses outside the 23-km inner zone is 0.5 Eötvös, which is one order smaller than the accuracy of the gravity gradient measurements. Note that the seafloor topography model was extended by 23 km from each of the four sides of the study region to enable the forward computation at the study region margin. The resolution of the extended seafloor topography is still 400 m by 400 m. The extended seafloor topography in regions not covered by multibeam bathymetry is interpolated from the global topography model version 18.1 released by the Scripps Institution of Oceanography [22], whose spatial resolution is about 1230 m. The cost function is an indicator that measures the discrepancy between the forward computed and measured gravity gradients. In this study, it is defined as the L 2 norm of the difference between the forward computed and measured gravity gradients.

Analysis and Removal of the Non-Terrain Effects
The high-resolution high-accuracy multibeam bathymetric data shown in Figure 3b enable us to analyze the gravity gradients caused by non-terrain factors (called the non-terrain effects). Firstly, the gravity gradients at flight altitude were forward computed from the multibeam bathymetry represented by adjacent right rectangular prisms. They are called the terrain effects and are shown in Figure 4a. Note that for each computation point, only prisms within its inner zone of radius, 23 km, were used in the forward computation. The gravity gradient due to one prism was computed using formulas in [35]: where ρ is the density of the right rectangular prism, (x, y, z) are the coordinates of the computation point, (ξ, η, ζ) are the variables of integration over the right rectangular prism, [x 1 , x 2 ], [y 1 , y 2 ] and [z 1 , z 2 ] are the x axis, y axis, and z axis boundaries of the right rectangular prism, respectively. Next, the difference between the measured gravity gradients and the terrain effects is defined as the non-terrain effects and shown in Figure 4b. The non-terrain effects are caused by the sub-topography density anomalies, regional isostatic compensation, and model errors in the forward computation due to the neglected masses outside the 23-km inner zone. Since the model errors are in the level of 0.5 Eötvös [36], and the isostatic compensation is a regional effect that is not significant over small regions [18], the non-terrain effects shown in Figure 4b are mainly attributed to the sub-topography density anomalies. density anomalies, regional isostatic compensation, and model errors in the forward computation due to the neglected masses outside the 23-km inner zone. Since the model errors are in the level of 0.5 Eötvös [36], and the isostatic compensation is a regional effect that is not significant over small regions [18], the non-terrain effects shown in Figure 4b are mainly attributed to the sub-topography density anomalies. Comparison between Figure 4a and Figure 3a shows that the terrain effect contributes less than 50% of the gravity gradient measurements (the RMS of the terrain effect shown in Figure 4a is about 49.13% of the RMS of the directly measured gravity gradients shown in the black dashed rectangle in Figure 3a). Topography estimation from gravity gradients in this area is feasible yet inaccurate since Comparison between Figures 4a and 3a shows that the terrain effect contributes less than 50% of the gravity gradient measurements (the RMS of the terrain effect shown in Figure 4a is about 49.13% of the RMS of the directly measured gravity gradients shown in the black dashed rectangle in Figure 3a). Topography estimation from gravity gradients in this area is feasible yet inaccurate since sub-topography density variation introduces large errors. So, the non-terrain effects shown in Figure 4b should be removed prior to the inversion. Recent studies [4,7,8,37,38] tend to remove the non-terrain effects by (1) calculating the difference between the measured and forward-computed VGGs over areas constrained by bathymetric data or land surface elevation, (2) interpolating this difference into areas without constraining data using a minimum curvature algorithm [39], (3) and then subtracting the interpolated difference from the measured VGGs. This method is practical in offshore areas, where land surface elevation (e.g., from the GIMP model [40]) or sparse ship soundings are usually available inside or in the surrounding of the study area to provide the constraints [3][4][5][6][7][8]. We examined the performance of this method by designing 5 experiments with different constraints. Experiment No. 1 directly estimates the seafloor topography from the measured VGGs without removal of the non-terrain effects. Experiment No. 2 assumes that the bathymetric depths in the surrounding of the study area (see Figure 5a) are available to the researchers, estimates the non-terrain effects (see Figure 5d) using constraints from known bathymetry, and removes them from the measured VGGs prior to the inversion. Experiment No. 3 assumes that sparse ship soundings inside the study area (see Figure 5b) are also available to the researchers. Experiment No. 4 increases the amount of available ship soundings inside the study area (see Figure 5c). Besides, we also added experiment No. 0 to examine the reliability of the simulated annealing method. Experiment No. 0 estimates seafloor topography from the synthesized VGGs shown in Figure 4a that are forward computed from multibeam bathymetric data. This means the non-terrain effects were completely removed, so ideally the estimated seafloor topography should be consistent with the multibeam bathymetry.
Comparison between the actual non-terrain effects (shown in Figure 4b) and the estimated non-terrain effects (shown in Figure 5d-f) indicates that experiment No. 2 to No. 4 only removes the long-wavelength component of the non-terrain effects from the measured VGGs. The difference between them is the residual error in the input to the inversion model estimating the seafloor topography, and is shown in Table 2. Experiment No. 1 directly uses the measured VGGs, so its residual error is the actual non-terrain effects. Experiment No. 0 uses the terrain effect as input to the inversion model. This is equivalent to a complete removal of the non-terrain effects, so its residual error is 0. Comparison between the actual non-terrain effects (shown in Figure 4b) and the estimated nonterrain effects (shown in Figure 5d-f) indicates that experiment No. 2 to No. 4 only removes the longwavelength component of the non-terrain effects from the measured VGGs. The difference between them is the residual error in the input to the inversion model estimating the seafloor topography, and is shown in Table 2. Experiment No. 1 directly uses the measured VGGs, so its residual error is the actual non-terrain effects. Experiment No. 0 uses the terrain effect as input to the inversion model. This is equivalent to a complete removal of the non-terrain effects, so its residual error is 0. (d-f) Non-terrain effects estimated using constraints from bathymetric data shown in (a-c), respectively; they are removed from the measured vertical-vertical gravity gradients prior to the inversion. Table 2. Statistics of the residual error, namely the difference between the actual non-terrain effects shown in Figure 4b and

Results, Evaluation, and Discussion
The SA programs of the five experiments were run on the supercomputer Oakley at the Ohio Supercomputer Center. Each program exited after 1 × 10 7 iterations. The estimated seafloor topographies are shown in Figure 6. The multibeam bathymetry is treated as the true depth to evaluate the accuracy of the seafloor topography estimated from airborne gravity gradients using SA. The statistics of the difference between them are shown in Table 3. The statistics of the difference between the measured VGGs and the one forward computed from the seafloor topography estimation are shown in Table 3

Results, Evaluation, and Discussion
The SA programs of the five experiments were run on the supercomputer Oakley at the Ohio Supercomputer Center. Each program exited after 1 × 10 iterations. The estimated seafloor topographies are shown in Figure 6. The multibeam bathymetry is treated as the true depth to evaluate the accuracy of the seafloor topography estimated from airborne gravity gradients using SA. The statistics of the difference between them are shown in Table 3. The statistics of the difference between the measured VGGs and the one forward computed from the seafloor topography estimation are shown in Table 3 too.  Figure 6a shows that the seafloor topography was successfully recovered from the synthesized VGGs shown in Figure 4a using SA, with estimation error of ±3 m. The error is mainly caused by the large dimensionality of the inversion, and the comparatively small number of iterations used. This confirms the reliability of simulated annealing in recovering the seafloor topography, so that the estimation error in the other four experiments is mainly attributed to the sub-seafloor density variations. Table 3. (a) Statistics of the difference between the multibeam bathymetry and the seafloor topography estimated from gravity gradients using simulated annealing; (b) Statistics of the difference between the measured gravity gradients and the one forward computed from the seafloor topography estimation.  Figure 6b shows that experiment No. 1 recovers the general shape of the seafloor topography, including the sharp drop at 375 km Easting, the 2-km-wide Y-shaped trough system at the southeast corner and its 0.8 km wide branch at km [378. 5,5365]. Considering the typical resolution of airborne gravimetry is 2 km (capable of recovering topographic features with wavelength longer than 4 km), it is shown that the airborne gravity gradient is capable of recovering shorter-wavelength topographic features than gravity disturbance. However, the estimation error is as large as ±24 m. This is because the airborne gravity gradient measurements used in this study were conducted for geologic purposes in areas with complex sub-topography geologic structures, where the assumption of uniform density does not hold. The large non-terrain effects (22 Eötvös RMS, mainly due to sub-seafloor density variations) prevents experiment No. 1 from achieving high estimation accuracy. This experiment highlights the need to remove the non-terrain effects. Note that this is not the drawback of the airborne gravity gradiometry. In seafloor topography estimations using gravity disturbance, the non-terrain effect is one of the largest error sources too [22,28].
Experiments No. 2 to No. 4 remove the non-terrain effects using constraints from known bathymetry. In experiment No. 2, the non-terrain effects are estimated from constraining bathymetry in the surrounding of the study area. Figure 6c and Table 3 indicate that compared with experiment No. 1, experiment No. 2 does not improve the estimation accuracy. This can be explained through the residual errors shown in Table 2. The STD of the residual error in experiment No. 1 is in the same level with that in experiment No. 2. Experiment No. 3 has constraining bathymetry both inside and in the surrounding of the study area. Figure 6d and Table 3 indicate that the estimation accuracy is improved to ±17 m. Experiment No. 4 increases the amount of the bathymetric data inside the study area, but surprisingly, Figure 6e and Table 3 do not show improvements in the estimation accuracy. Probably it is because the randomly selected constraining bathymetry does not capture the proper feature points of the seafloor topography. Combining these results and the comparison between Figure 4b with Figure 5d-f, it is shown that the sparse bathymetric constraints could be used to remove the long-wavelength component of the non-terrain effects, which highlights the need to make future gravity gradiometry survey overlap the known bathymetry or land surface elevations.

Conclusions
Airborne gravimetry has been used extensively by the glaciological community to estimate the seafloor topography at places covered by ice shelves, iceberg, or ice mélange, which are inaccessible to traditional ship soundings. Airborne gravity gradiometry is an alternative source of gravity data. It inherently has the potential to achieve higher spatial resolution for improved inference of seafloor topography, and thus is in need of demonstration. By analyzing the theoretical admittances, we demonstrated the superiority of airborne gravity gradient in recovering small topographic features over shallow waters. We also present the first numerical experiments that use real airborne gravity gradient data to estimate the seafloor topography. Understandably, we achieved only moderate success since in our study area, the unknown variations in subsurface mass density dominate the gravity gradient signal over that caused by the geometric configuration of the seafloor. It is shown that airborne gravity gradiometry can recover small topographic features that are outside of the resolution capability of typical airborne gravimetry, but the estimation accuracy is only ±17 m due to the presence of large non-terrain effects. The long-wavelength non-terrain effect can be removed if sparse bathymetric data are available inside the study area to provide constraints, whereas the short wavelengths cannot. This suggests that future airborne gravity gradiometry survey should overlap some known bathymetry or land surface elevations, and that more advanced non-terrain effect removing techniques should be investigated in the future to fully realize the potential of the airborne gravity gradients at short wavelengths. We remain optimistic about seafloor topography estimation from airborne gravity gradients, since at some places, the subsurface density variation is significant while at some other places the density of terrain is almost constant, where the non-terrain effect is small. Future studies at places where the density of terrain is almost uniform are needed to fully demonstrate the estimation capability of airborne gravity gradients. This study is an important starting point for future work that constantly improves glaciologists' understanding of the advantage and performance of airborne gravity gradients in seafloor topography estimation. This study could help glaciologists make better decisions when selecting gravity measurement devices for future airborne geophysical missions in polar areas.