Enhanced gravity-geologic method to predict bathymetry by considering non-linear effects of surrounding seafloor topography

SUMMARY

The gravity-geological method (GGM) is an approach that utilizes marine gravity anomalies (GAs) and shipborne bathymetric data to invert seafloor topography by resolving short-wavelength GAs through the Bouguer Plate approximation. Such an approximation ignores the non-linear effects caused by surrounding seafloor topographical undulations that actually exist in short-wavelength GAs, and thus leaving the space for further modification of GGM. This study thoroughly derives the relationship between seafloor topography and GA, as well as the formula of GGM. Then, we propose a self-adaptive method to improve the accuracy of the inversion significantly: the enhanced GGM (EGGM). The method uses the equivalent mass line method to approximate the non-linear gravitational effects of the surrounding seafloor topography to correct the short-wavelength GAs. By introducing two optimal density contrast parameters, EGGM has been designed to effectively integrate the combined effects of various non-linear factors to a certain extent. The accuracy of the seafloor topography models, produced with a spatial resolution of 1′ × 1′, was evaluated over the study area (132°E–136°E, 36°N–40°N) located in the Sea of Japan. The results indicate that the accuracy of EGGM has a relative improvement of 13.73 per cent compared to that of GGM in the overall study area, while the accuracy of both models is higher than that of the SIO_unadjusted model. The study further investigated the feasibility and stability of EGGM by examining the accuracy of both GGM and EGGM in various water depth ranges and areas with diverse terrain characteristics.

1. INTRODUCTION

Approximately 71 per cent of the surface of the Earth is covered by the water. Seafloor topography data serve as a foundation for understanding the Earth's shape, seafloor tectonic movements and seafloor evolution. It provides essential information for various fields, including maritime navigation, solid earth geophysics, resource development and marine environmental monitoring (Smith & Sandwell 1997). The techniques used to detect seafloor topography primarily include the ship-based sonar sounding, LiDAR sounding, remote sensing image inversion and inversion with satellite altimetry-derived gravity. LiDAR sounding is limited to shallow waters due to its dependence on laser beam penetration. Remote sensing image inversion is affected by radar and spectrometer settings, as well as ocean conditions. Ship-based sonar sounding is known for its high accuracy, but it is less efficient and time-consuming. Therefore, the main technique for obtaining precise global seafloor topography in nowadays and future is likely to be the inversion with satellite altimetry-derived gravity (McMillan et al. 2009; Sandwell & Smith 2009; Sun et al. 2022).

Currently, with satellite altimetry data, marine gravity anomalies (GAs) at the sea surface can be computed with an accuracy of 2–3 mGal (Guo et al. 2022; Zhu et al. 2022). The approach of utilizing altimetry-derived GA to invert seafloor topography has gained widespread applications (Dixon et al. 1983; Vrdoljak et al. 2021; Wei et al. 2021). Observed GA at the sea surface are caused by various factors, such as the seawater mass deficit, uneven distribution of crustal density, and isostatic compensation. Among these factors, the dominance is attributed to the seawater mass deficit (Luo et al. 2002). Seawater mass deficit is caused by filling the space between the sea level and seafloor topography with the seawater while considering the fact that the density of seawater is less than that of the Earth's crust. Therefore, it is possible to obtain the inversion of seafloor topography by clarifying the physical relationship between observed GA at the sea surface and seafloor topography, and developing a proper mathematical model for practical data processing.

Various methods for seafloor topography inversion have been proposed in the past decades, mainly including the admittance function technique (McNutt 1979), the linear regression analysis (S&S, Smith & Sandwell 1997), and the gravity-geological method (GGM) (Kim et al. 2010, 2011; Kim & Yun 2018), etc. Due to the limitations of current geophysical parameter accuracy and theoretical challenges regarding the crustal equilibrium hypothesis, the overall accuracy of the admittance function technique for seafloor topography inversion is not high and requires further development (Jiang et al. 2023). Influenced by the distribution of shipborne bathymetric data and seafloor geological changes, the intermediate product scale factor grid of the linear regression method is prone to anomalies, which can affect the accuracy of the inversion. Therefore, Smith & Sandwell (1997) used shipborne bathymetric data to control the inversion model grid, correcting and adjusting the initial seafloor topographic model of the S&S method. Among these methods, GGM stands out for its ability to provide a simple and efficient prediction of seafloor topography using shipborne bathymetric data and observed GA. It does not require the input of complex geophysical parameters and is applicable to large-scale sea areas, making it highly feasible (Kim et al. 2011; Wei et al. 2021).

GGM was initially proposed for detecting terrestrial bedrock thickness (Ibrahim & Hinze 1972), but it has also been found to be suitable for the inversion of seafloor topography utilizing the satellite altimetry-derived GA (Hsiao 2016; Sun et al. 2017). It inverts the seafloor topography by separating the so-called short-wavelength GA from the observed GA at the sea surface with the help of shipborne bathymetric data at the control points. The short-wavelength GA refers to the portion of the observed GA that has a tight connection to the seafloor topography, specifically arising from the seawater mass deficit. Since the GGM is equivalent to the interpolation of shipborne bathymetric control points on the inverse model, it can incorporate short-wavelength information from shipborne bathymetric control data to a greater extent, resulting in higher accuracy. Kim et al. (2010) compared the power spectral density (PSD) of the GGM model and the S&S model in the study area. The results indicated that the GGM can better represent seafloor topographic information within the wavelength band of less than 12 km.

GGM research primarily focuses on the determination of optimal density contrast parameter and the construction of appropriate short-wavelength GA. Kim et al. (2011) used an iterative approach to determine empirical values for the density contrast parameter. Hsiao et al. (2011) proposed a downward continuation method to estimate the density contrast parameter. Xiang et al. (2017) used a finite element approach based on an adaptive triangular mesh made up of unequal control points to mimic the long-wavelength GA. Xing et al. (2020) simulated the short-wavelength GA using a rectangular prism model, which the density difference parameter was found to have minimal impact on the inversion process and had lost its original physical meaning. By computing the short-wavelength gravity correction at the control points based on the weight parameters added to the Bouguer Plate formulation, An et al. (2022; 2024) refined the long-wavelength gravity field.

In addition to shipborne bathymetric data and satellite altimetry data along with its marine gravity field model, further enhancement of the accuracy in seafloor topography inversion requires the development of non-linear seafloor topography inversion theory and methodology. In 2005, the failure of the U.S. Navy's navigation products to depict a seamount south of Guam resulted in the collision of the U.S. nuclear attack submarine USS San Francisco with the seamount. However, the seamount was detected on the SIO V8.2 undersea terrain model. Compared to the previous version, the V8.2 model incorporates Parker's higher-order non-linear term iteration method proposed by Oldenburg (1974), exhibiting improvements in spectrum, accuracy, and performance across various topographic sea areas (Marks et al. 2010). It shows that one of the avenues for advancing seafloor topography inversion utilizing satellite altimetry and gravity data involves incorporating non-linear terms into the inversion process.

The presence of non-linear effects makes the inversion process non-unique, making it challenging to find an appropriate algorithm to address these non-linearities. With the advancement of computer performance, the use of the artificial intelligence (AI) algorithms that can optimize non-linear effects independently of physical models has become one of major research focus (Annan & Wan 2022; Sun et al. 2023; Zhou et al. 2023; Zhou et al. 2024). However, there is a lack of research grounded in physical model-based methods. Recent studies primarily include the simulated annealing method in the spatial domain (Yang et al. 2018) and higher-order Parker expansion and modified Bott's iteration combined method in the frequency domain (Xu et al. 2023).

The GGM approximates the short-wavelength GA extracted from the observed GA by the Bouguer Plate, which acts as a linear relationship with the seafloor topography for inversion purposes. In fact, the short-wavelength GA at a specific point is influenced not only by the local seafloor topography (i.e. the linear part) but also by the surrounding seafloor topography (i.e. the non-linear part). In this study, the surrounding seafloor topography was represented by rectangular prisms based on the grid size. By utilizing the equivalent mass line approach, the non-linear gravitational effects of the surrounding seafloor topography were derived. This led to the proposal of the enhanced GGM (EGGM), which takes into account the non-linear effects of the surrounding seafloor topography for seafloor topography inversion. Such comprehensive findings aims to enhance researchers’ understanding of the underlying principles of GGM. Combining shipborne bathymetric data and the latest V32.1 satellite altimetry-derived GA published by Scripps Institution of Oceanography (SIO), seafloor topography models with the resolution of 1′ × 1′ were constructed in the Sea of Japan. Experiments regarding different water depth ranges and regions with different terrain features were conducted separately to comprehensively evaluate the inversion results with the aim of verifying the feasibility and stability of EGGM.

2. METHODOLOGY

In GGM, the observed GA at the sea surface can be simply divided into two components: the short-wavelength component and the long-wavelength component. The short-wavelength component is theoretically attributed to the seawater mass deficit resulting from the density contrast between the seafloor crust and seawater. The long-wavelength component is influenced by density anomalies and interface fluctuations within the deep Earth (Ibrahim & Hinze 1972). Therefore, the observed GA at the sea surface $\mathrm{\Delta }{g}_{\mathrm{o}\mathrm{b}\mathrm{s}}$ can be expressed as

where $\mathrm{\Delta }{g}_{\mathrm{s}\mathrm{h}\mathrm{o}\mathrm{r}\mathrm{t}}$ is the short-wavelength GA and $\mathrm{\Delta }{g}_{\mathrm{l}\mathrm{o}\mathrm{n}\mathrm{g}}$ is the long-wavelength GA.

Assuming that the surrounding seafloor topography with the calculation point O as the centre and the radius of R is taken into account, the GA caused by seawater mass deficit at the point O is given by

where G is the universal gravitational constant, $\mathrm{\Delta }\rho$ is the density contrast parameter between the seawater and crust, $\alpha$ is the angle of integration with the calculation point as the centre of the circle, r is the radius of integration taking into account the surrounding seafloor topography, $\xi$ is the integration height of water depth, $\mathrm{\Delta }{g}_O^1$ is the short-wavelength GA commonly used in GGM, which can be regarded as the linear effects caused by a Bouguer Plate with a thickness $h_0$ and a density $\mathrm{\Delta }\rho$⁠, and $\mathrm{\Delta }{g}_O^2$ is the non-linear term resulting from the undulations of the surrounding seafloor topography. For brevity, we refer to it as the correction term for the short-wavelength GA in this study.

2.1. Derivation of the short-wavelength GA

Solving for the integral of $\mathrm{\Delta }{g}_O^1$ yields

Expanding the term $(\sqrt{R^2+h_0^2})$ by the Taylor series expansion and substituting it into eq. (3) yields

As shown in Fig. 1, a reference plane at a certain depth D beneath the oceanic crust is selected. A horizontal plane with a water depth value of X is established, passing through the seafloor point $O^{\prime }$ corresponding to the calculation point O. This configuration forms a Bouguer Plate from the plane X to the plane D. At this point, the Bouguer Plate-induced GA can be calculated as

Assuming that $R\to \mathrm{\infty }$⁠, the linear effects of the seafloor topography can be approximated as the effects of an infinite horizontal Bouguer Plate. Consequently, in such a case, the short-wavelength GA expression can be approximated as

This is the calculation equation of GGM. To validate the assumption $R\to \mathrm{\infty }$ in eq. (5), we tested the accuracy of GGM for different depths D, ranging from 0 to 7000 m in 100 m increments. Theoretically, selecting the D should eliminate some effects of higher-order terms. However, when the assumption holds, the choice of reference plane D should not impact inversion accuracy. The result shows that the choice of D has virtually no effect on inversion accuracy, with differences being less than a submillimetre in magnitude.

It is easy to see that GGM ignores the calculation of the correction terms, and the water depth is directly inverted by the extracted short-wavelength GA, while the correction terms are implicitly absorbed in the long-wavelength GA. In practical calculations, the extraction of the short-wavelength GA proves to be more advantageous when dealing with small topographic relief and a dense distribution of shipborne bathymetric data. However, in the area with large topographic relief or a lack of shipborne bathymetric data, the absorption of correction terms by the long-wavelength GA is uneven and insufficient. At this point, the short-wavelength GA contains a large amount of additional non-linear gravity signals, which leads to reduced accuracy in the inversion process. The calculation of the correction terms of the short-wavelength GA should not be ignored. Therefore, the observed GA should be expressed as

where $\mathrm{\Delta }{g}_{\mathrm{s}\mathrm{h}\mathrm{o}\mathrm{r}\mathrm{t}\mathrm{\_}\mathrm{T}\mathrm{C}}$ is the corrected short-wavelength GA, and $\mathrm{\Delta }{g}_{\mathrm{T}\mathrm{C}}$ is the correction term to the short-wavelength GA.

2.2. Derivation of the correction term to the short-wavelength GA

The correction term to the short-wavelength GA on the calculation point O is the effect of the terrain undulations relative to the plane X. The topography is divided into rectangular prisms based on the grid resolution, with each prism's mass represented by a vertical mass line along its centre. The calculation is conducted using the definite integration.

Given that the gravitational potential of a point mass m with respect to the calculation point is $V=\frac{Gm}{l}$ and the distance from the mass point to the calculation point is l. At the calculation point, the gravitational attraction caused by the mass point is given by

where $V_Z$ is the vertical derivative component of the gravitational potential V in the direction normal to the gravity equipotential surface, and z is the component of the distance from the mass point to the calculation point in normal direction.

If the mass of a rectangular prism is M and the height is $\mathrm{\Delta }H$ (positive above and negative below the plane X), then we have $\mathrm{d}m=\frac{M}{\mathrm{\Delta }H}\mathrm{d}z$⁠. Then the GA caused by this rectangular prism at point O can be derived as

where R is the projected distance from the centre of the prism to the point O on the sea surface, $V=\mathrm{\Delta }x\mathrm{\Delta }y\mathrm{\Delta }H$ is the volume of the prism, $\mathrm{\Delta }x$ and $\mathrm{\Delta }y$ are the grid spacings along the latitude and longitude directions, respectively. The calculation formula for spherical coordinate distance is given by Lin & Denker (2019)

where $\ell$ is the distance from the point A to point B, $(r_A,{\varphi }_A,{\lambda }_A)$and $(r_B,{\varphi }_B,{\lambda }_B)$ are their spherical coordinates, respectively.

The GA caused by all n surrounding rectangular prisms at point O is finally given by

2.3. Calculation formula for practical applications

In general, the actual density contrast between the seawater density of 1.03 g cm⁻³ and the average crustal density of 2.70 g cm⁻³ is 1.67 g cm⁻³. However, in practical applications, the extraction and calculation of the short-wavelength GA may not ideal in theory. The non-linear terms in the calculation of the short-wavelength GA can be affected by factors other than the seafloor topography. For instance, in situations where shipborne bathymetric data is sparsely distributed, the non-linear terms may be influenced by the imperfect construction of the long-wavelength GA and the presence of crustal isostatic compensation. Therefore, the density contrast parameters in GGM are generally regarded as empirical values, and determined using the iterative method with shipborne bathymetric data.

The Airy isostatic hypothesis posits that the Earth's crust is floating over the denser mantle in a state of hydrostatic equilibrium. This equilibrium leads to the Moho surface reflecting the surface topography of the Earth (Calmant & Baudry 1996; Hofmann-Wellenhof & Moritz 2006). According to the hypothesis applied to the seafloor topography, the GA caused by the seawater mass deficit of seamounts is actually smaller compared to that of the seafloor plain. Therefore, in the practical calculations, it is important to account for the effects of various non-linear terms, including crustal isostatic compensation. To enhance the self-adaptability of the non-linear term model and reduce the number of parameters used, this study presents the calculation equation of EGGM for corrected short-wavelength GA, derived from eqs (6), (7) and (11), denoted as

where $\mathrm{\Delta }{\rho }_1$ and $\mathrm{\Delta }{\rho }_2$ are the density contrast parameters for the linear and non-linear terms, respectively.

To account for the combined effects of various factors on the observed GA, the values of $\mathrm{\Delta }{\rho }_1$ and $\mathrm{\Delta }{\rho }_2$ are all determined using the iterative approach. This iterative process ensures that the GA accurately reflects the seafloor topography. The principle of selecting the optimal number of iterations is to observe the trend of the standard deviation (STD) and the correlation coefficient between the inverted bathymetry and the shipborne bathymetry at different iterations and select the optimal density contrast parameter corresponding to the minimum STD and the maximum correlation coefficient (Ouyang et al. 2014). By comparing these statistical measures, the optimal density contrast parameter can be determined, resulting in a better agreement between the inverted bathymetry and the shipborne bathymetry. The algorithm theoretically guarantees that the accuracy is at least equivalent to that of the GGM model, even in extreme cases where the prior model of GGM is not well constructed. This provides a robust solution and is achieved by setting the iteration parameter $\mathrm{\Delta }{\rho }_2$ of the non-linear term to zero.

2.4. Algorithm flow

According to eq. (12), it is known that a prior model needs to be established first to calculate the non-linear terms of the surrounding seafloor topography. The existing seafloor topography models are typically constructed by incorporating shipborne bathymetric data, which implies that the selected check points are also incorporated into these models. Therefore, in order to reduce the impact of the redundant information in the prior model on the comparison between the GGM model and the EGGM model, a prior seafloor topography model based GGM is initially constructed here using the shipborne bathymetic data and the GA at the control points. According to eq. (11), the study uses a radius of approximately 100 km to adequately consider the influence of the surrounding seafloor topography. As each computation point within the study area requires the consideration of a certain radius, the boundary of the GGM-based prior model is extended by 1°, relative to the boundary of the study area, to prevent edge effects. The EGGM algorithm is described as follows, and its flow is depicted in Fig. 2.

1. Use the shipborne bathymetric data and GA to construct a GGM-based prior model.

2. According to eq. (12), set the iterative step size of $\mathrm{\Delta }{\rho }_1$ and $\mathrm{\Delta }{\rho }_2$ to 0.1, respectively, and forward modelling the corrected short-wavelength GA $\mathrm{\Delta }{g}_{\mathrm{s}\mathrm{h}\mathrm{o}\mathrm{r}\mathrm{t}\mathrm{\_}\mathrm{T}\mathrm{C}}^j$ at the control point j with the prior model and different iterative parameters $\mathrm{\Delta }{\rho }_1$ and $\mathrm{\Delta }{\rho }_2$⁠. At this time, X is the water depth at the control point j.

3. Interpolate the observed GA grid $\mathrm{\Delta }{g}_{\mathrm{o}\mathrm{b}\mathrm{s}}^i$ to the control point j to obtain the discrete observed GA $\mathrm{\Delta }{g}_{\mathrm{o}\mathrm{b}\mathrm{s}}^j$ and subtract $\mathrm{\Delta }{g}_{\mathrm{s}\mathrm{h}\mathrm{o}\mathrm{r}\mathrm{t}\mathrm{\_}\mathrm{T}\mathrm{C}}^j$ from $\mathrm{\Delta }{g}_{\mathrm{o}\mathrm{b}\mathrm{s}}^j$ to obtain the discrete long-wavelength GA $\mathrm{\Delta }{g}_{\mathrm{l}\mathrm{o}\mathrm{n}\mathrm{g}}^j$ at the control point j, as given by eq. (7).

4. Grid $\mathrm{\Delta }{g}_{\mathrm{l}\mathrm{o}\mathrm{n}\mathrm{g}}^j$ to form a long-wavelength GA grid $\mathrm{\Delta }{g}_{\mathrm{l}\mathrm{o}\mathrm{n}\mathrm{g}}^i$⁠, and subtract $\mathrm{\Delta }{g}_{\mathrm{l}\mathrm{o}\mathrm{n}\mathrm{g}}^i$ from $\mathrm{\Delta }{g}_{\mathrm{o}\mathrm{b}\mathrm{s}}^i$ to obtain the corrected short-wavelength GA grid $\mathrm{\Delta }{g}_{\mathrm{s}\mathrm{h}\mathrm{o}\mathrm{r}\mathrm{t}\mathrm{\_}\mathrm{T}\mathrm{C}}^i$⁠, as given by eq. (7).

5. Use the prior model and different iterative parameters $\mathrm{\Delta }{\rho }_2$ to calculate the short-wavelength GA correction grid $\mathrm{\Delta }{g}_{\mathrm{T}\mathrm{C}}^i$⁠. At this time, X is the water depth of the prior seafloor topography model.

6. Subtract $\mathrm{\Delta }{g}_{\mathrm{T}\mathrm{C}}^i$ from $\mathrm{\Delta }{g}_{\mathrm{s}\mathrm{h}\mathrm{o}\mathrm{r}\mathrm{t}\mathrm{\_}\mathrm{T}\mathrm{C}}^i$ to obtain the short-wavelength GA grid $\mathrm{\Delta }{g}_{\mathrm{s}\mathrm{h}\mathrm{o}\mathrm{r}\mathrm{t}}^i$ and inverse the bathymetry model by eq. (6). At this time, X is the unknown water depth to be inverted.

7. Select the optimal $\mathrm{\Delta }{\rho }_1$ and $\mathrm{\Delta }{\rho }_2$ according to the criteria of reaching to the minimum STD and the maximum correlation coefficient to inverse the optimal bathymetry model.

3. STUDY AREA AND DATA SOURCES

3.1. Study area

The Sea of Japan, located in the western Pacific Ocean, is the largest marginal sea in the region. It features a relatively narrow continental shelf and is characterized by a predominantly deep-water basin. The northern part of the Sea of Japan, above 40°N latitude, is known as the Japan Sea Basin, where the seafloor is mostly over 3000 m deep and relatively flat. In contrast, the southern part, below 40°N latitude, exhibits a more complex seafloor topography, including basins, ridges, and troughs. Examples of these features include the Yamato Basin in the east and the Tsushima Basin in the southwest.

To assess the performance of EGGM in areas with complex topography, a specific study area within the Sea of Japan has been selected. This study area spans from 132°E to 136°E in longitude and from 36°N to 40°N in latitude, as depicted by the white border in Fig. 3. To mitigate edge effects, the study area's boundary has been extended outward by 1°, resulting in a wider range of 131°E to 137°E in longitude and 35°N to 41°N in latitude. This expanded area ensures that the analysis encompasses the desired region without too many potential edge effects.

3.2. Data sources

3.2.1. Gravity anomaly

In this study, the utilized GA model is the grav_32.1, which was obtained from the Scripps Institution of Oceanography (SIO) with a grid resolution of 1′ × 1′. This model, which is released in August 2022, incorporates updated data from Altika, Cryosat LRM, Cryosat-SAR and Sentinel-3A/B satellites, covering an additional 12 months of observations. To mitigate potential edge effects between land and sea areas, a precise land mask was used to enhance the accuracy of the analysis. The inclusion of the additional data and the utilization of an accurate land mask contribute to a more comprehensive and refined GA model. The data set can be accessed at https://topex.ucsd.edu/pub/global_grav_1 min/.

3.2.2. Global topographic models

In this study, the GEBCO_2022 global terrain model provided by the Nippon Foundation-GEBCO Seabed 2030 Project was utilized for the pre-processing of shipborne bathymetric data (GEBCO Bathymetric Compilation Group 2022 2022). The GEBCO_2022 model, which is released in June 2022, offers a spatial resolution of 15″ and an elevation accuracy of meters (Becker et al. 2009; Mayer et al. 2018). It serves as a comprehensive and reliable resource for analysing and incorporating bathymetric data in research. Further information about the GEBCO_2022 model, including its data and products, can be accessed at https://www.gebco.net/data_and_products/.

Due to the relatively dense distribution of shipborne bathymetric data in the Sea of Japan region, many of final release models use shipborne bathymetry as a proxy for predicting grid bathymetry when constructing models for this region. This hinders the effective assessment of the predictive accuracy of release models. The majority of globally available seafloor topography models incorporate the SIO global topography model, which uses linear regression analysis (S&S) for predictions. Smith & Sandwell (1997) utilized shipborne bathymetric data in the construction process but still adjusted their topography predictions to match the bathymetry-implied heights with the projected topography for grid cells restricted by shipborne bathymetry. In this paper, their unadjusted version of the prediction denoted by SIO_unadjusted model in the following, will be used for comparison with the prediction accuracy results of GGM and EGGM. The corresponding release model for this unadjusted model is topo_18.1.

3.2.3. Shipborne bathmetric data

A total number of 36 371 single-beam bathymetric data were collected from the National Centers for Environmental Information (NCEI) of the National Oceanic and Atmospheric Administration (NOAA) for the study area (https://www.ncei.noaa.gov/maps/bathymetry/). It is worth to note that the quality of shipborne bathymetric data can vary due to different observation conditions. In some cases, data collected earlier may have larger errors, and therefore it is necessary to exclude those data (Smith 1993; Tozer et al. 2019). To identify and remove potentially erroneous data, shipborne bathymetric data that exceeded twice the difference between their depths and the corresponding bathymetry of GEBCO_2022 model were removed, resulting in a data set of 35 578 points. The rejection rate based on this coarse difference rejection process is approximately 2.18 per cent. To reduce the correlation between control points and check points and keep the independence of check data, two independent cruise surveys, namely ‘ODP127JR’ and ‘84 003 111’ were selected as check points. The distribution of various types of points, including control points and check points, is illustrated in Fig. 4.

4. CASE STUDY AND ANALYSIS

4.1. Accuracy analysis in overall study area

The gridding process used the Generic Mapping Tools (GMT) function surface that interpolates or extrapolates data using adjustable tension continuous curvature splines (Smith & Wessel 1990). The used tensor factor is 0.4 in the study. Fig. 5(a) and (b) display the STD and correlation coefficient between the shipborne bathymetry and the inversion of GGM and EGGM, respectively. These plots illustrate the variation of these metrics with different density contrast parameters. The iteration step size for the GGM is set as 0.01, while for the EGGM it is set as 0.1. The reason of doing so is to reduce the computational load in the case of using EGGM. The optimal density contrast parameter $\mathrm{\Delta }\rho$ for the GGM, which corresponds to the minimum STD and the maximum correlation coefficient, is chosen as 1.14 g cm⁻³. The optimal density contrast parameters $\mathrm{\Delta }{\rho }_1$ and $\mathrm{\Delta }{\rho }_2$ of EGGM are chosen as 1.8 and 2.7 g cm⁻³, respectively. As shown in Fig. 5(b), the fluctuation of the curve gradually diminishes with the convergence of $\mathrm{\Delta }{\rho }_1$⁠, where the fluctuation of the curve reflects an iterative process of $\mathrm{\Delta }{\rho }_2$⁠. The influence of $\mathrm{\Delta }{\rho }_2$ on STD and correlation coefficient initially decreases and eventually stabilizes as the iterations progress.

The accuracy of the two inverted models was evaluated using the check points, and the results are shown in Table 1. The evaluation indexes include Max (maximum difference), Min (minimum difference), Mean (mean difference), STD, RMS (root mean square difference) and MAPE (mean average percentage error). As can be seen from Table 1, the RMS of the GGM model is 56.37 m and the RMS of the EGGM model is 48.63 m. Compared with GGM model, the improvement of EGGM model is obvious, with an improvement of 7.74 m in terms of RMS and a relative improvement of 13.73 per cent. The RMS of the SIO_unadjusted model is 225.47 m, much higher than the first two. This illustrates the superiority of the GGM over the S&S in constructing seafloor topography models in areas of dense shipborne bathymetric data. In other indicators, SIO_unadjusted model's performance is significantly inferior to that of GGM and EGGM.

To visually demonstrate the difference distribution of the GGM model, the EGGM model and the SIO_unadjusted model, Fig. 6 shows the corresponding difference histograms and normal distribution curves for the three models. The results show that the majority of the absolute difference between the inverted bathymetry and check points depths for GGM and EGGM models fall within 200 m. Specifically, within the EGGM model, the absolute difference values within 50 and 100 m account for 80.81 and 94.96 per cent, respectively, which are superior to the corresponding values within the GGM model (71.86 and 92.93 per cent). Moreover, upon observing the difference normal distribution curves, it can be observed that the difference distribution for the EGGM model is more concentrated and the difference is smaller. The error distribution of the SIO_unadjusted model is relatively dispersed, with a large portion of the absolute difference remaining above 200 m. The SIO_unadjusted model is not primarily analysed in the following because of its poor performance.

In Fig. 7, the difference between the GGM model and the EGGM model is presented. The difference between the two models is primarily concentrated in areas with strong topographic variations, including the edges of islands and seamounts. The difference distribution allows for a rough distinction of the boundaries of seamounts and islands. Table 2 displays the statistical analysis of the difference between the GGM model and the EGGM model.

Fig. 8 illustrates the seafloor topography for the GGM model and the EGGM model. Both models exhibit a similar overall trend with a complex terrain. The topography of the study area is mainly characterized by the presence of seamounts, which have water depth ranges around 0–1000 m. Additionally, the submarine plains are primarily found in the eastern and northwestern regions, exhibiting water depth ranges around 0–3000 m.

Fig. 9 displays the short-wavelength GA correction calculated using the optimal density contrast parameters $\mathrm{\Delta }{\rho }_1$ and $\mathrm{\Delta }{\rho }_2$ at the discrete control points and gridpoints, respectively. The calculation process follows the steps (2) and (5) of the Algorithm Flow. To visualize the short-wavelength GA correction terms calculated by both types of data, the discrete values calculated at the control points j are first gridded and then plotted. It is evident that the non-linear short-wavelength GA correction exhibits a certain correlation with the seafloor topography, as observed in both the findings of this study and the results presented in the study of Claessens & Filmer (2020).

4.2. Accuracy analysis in different ranges of water depth

To verify the feasibility and stability of EGGM in different water depth ranges, it is important to analyse the correlation between the model accuracy and water depth. The percentage of check points within various water depth ranges was tallied, in which 0–1000 m accounted for 40 per cent, 1000–2000 m accounted for 16 per cent and 2000–3000 m accounted for 44 per cent.

As shown in Table 3, the accuracy indexes of each model in different water depth ranges are compared. It is evident that the EGGM model outperforms the GGM model in all water depth ranges, demonstrating a significant improvement. The RMS is reduced by 5.98, 8.23 and 9.51 m in the water depth ranges of 0–1000 m, 1000–2000 m and 2000–3000 m, respectively. The relative improvements are 15.56, 10.04 and 15.22 per cent, respectively. The SIO_unadjusted model continues to exhibit poor performance in all water depth ranges.

The results indeed verify the feasibility and stability of EGGM in different water depth ranges. The consistent improvements across different water depth ranges imply that EGGM successfully accounts for the non-linear effects of the surrounding seafloor topography. Moreover, the magnitude of improvements does not appear to be influenced by the water depth in the study area. This finding may be attributed to the initial accuracy of the prior model, where a smaller RMS indicates a more accurate non-linear calculation and a higher improvement magnitude.

4.3. Accuracy analysis in regions with different terrain features

To further discuss the feasibility and stability of EGGM in regions with different terrain features, it is important to conduct additional experimental validations (An et al. 2022). Three subregions were selected based on terrain characteristics and the distribution of check points. Specially, subregion A (134°E–135°E, 38°N–39°N), located at the boundary between seamounts and submarine plains, subregion B (134°E–135°E, 36°N–37°N) located at seamounts, and subregion C (135°E–136°E, 38°N–39°N) located at submarine plains, were chosen for the experimental analysis. Fig. 10 shows the locations of the chosen three subregions.

Fig. 11 displays the distribution of control points and check points within the subregion A, B and C. White markers indicate control points, and blue markers indicate check points. In subregion C, an independent cruise survey named ‘HNRO07RR’ was selected for validation. In subregion A, the majority of check points are located on the submarine plains, with some points positioned along the boundary between seamounts and submarine plains. In subregion B, the check points traverse through the seamounts. Check point route in subregion C is located on the submarine plains, and there are scattered submarine hills around. Table 4 provides information on the number of shipborne bathymetric data points available in each area, namely subregion A, B and C.

It has been known that the optimal density contrast parameter obtained from the application of GGM varies depending on the geological and topographic conditions of the seafloor. When comparing the accuracy within the same area, it is often observed that the accuracy of the inversion in the subregion is higher than that that in the entire study area. Therefore, it is necessary to perform separate inversions for each subregion thorough iterative iterations (Wei et al. 2021). The optimal density contrast parameters $\mathrm{\Delta }\rho$ based on GGM in subregion A, B and C are chosen as 1.63, 1.07 and 5.19 g cm⁻³, respectively. The optimal density contrast parameters in subregion A and B based on EGGM are both chosen as $\mathrm{\Delta }{\rho }_1$ = 1.6 g cm⁻³ and $\mathrm{\Delta }{\rho }_2$ = 3.0 g cm⁻³. In subregion C, the optimal density contrast parameters are selected as $\mathrm{\Delta }{\rho }_1$ = 5.2 g cm⁻³ and $\mathrm{\Delta }{\rho }_2$ = 2.8 g cm⁻³. It can be seen that $\mathrm{\Delta }\rho$ of GGM is closer to $\mathrm{\Delta }{\rho }_1$ of EGGM in all three subregions. Therefore, when applying EGGM, the optimal density contrast parameter based GGM can be used as a reference to determine and adjust the iteration interval of $\mathrm{\Delta }{\rho }_1$ and $\mathrm{\Delta }{\rho }_2$ to improve the operation efficiency. To show a more complete iteration curve, the iteration intervals of $\mathrm{\Delta }{\rho }_1$ and $\mathrm{\Delta }{\rho }_2$ are set to (0,3] for subregions A and B, while the iteration intervals of $\mathrm{\Delta }{\rho }_1$ and $\mathrm{\Delta }{\rho }_2$ are set to (0,10] for the subregion C. These iteration intervals can be further reduced in practice.

Fig. 12 shows the iterative curves of the optimal density contrast parameters versus the STD and the correlation coefficient for GGM and EGGM in subregions A, B and C, respectively. In the case of EGGM, the horizontal axis represents the interval of $\mathrm{\Delta }{\rho }_1$⁠, and the fluctuations of curves originate from the iteration of $\mathrm{\Delta }{\rho }_2$⁠. By observing the iterative curves of EGGM, it can be found that the curves tend to converge as $\mathrm{\Delta }{\rho }_1$ increases and the magnitude of the influence of $\mathrm{\Delta }{\rho }_2$ on STD stabilizes.

As indicated in Table 5, EGGM models exhibit enhanced accuracy across different subregions. The RMS of the EGGM model in subregion A is 73.78 m, surpassing that of the GGM model by 3.40 m, resulting in a relative improvement of 4.41 per cent. The RMS of the EGGM model in subregion B is 47.95 m, outperforming that of the GGM model by 6.63 m, with a relative improvement of 12.13 per cent. The RMS of the EGGM model in subregion C is 36.33 m, which is 0.99 m better than that of the GGM model, corresponding to a relative improvement of 2.65 per cent. Overall, all accuracy indices have shown the improvement. The results show that EGGM, which accounts for the non-linear effects of the surrounding seafloor topography in areas with different terrain features, is feasible and stable. Notably, the relative improvement highlights that the subregion B, situated on a seamount, exhibits the greatest improvement, followed by the subregion A, located at the boundary between seamounts and submarine plains, and finally the subregion C, situated on submarine plains.These outcomes agree with theoretical expectations, as regions with greater topographic relief experience stronger non-linear effects from the surrounding seafloor topography, resulting in a more pronounced improvement in the EGGM model. The results confirm that it is reasonable to consider the non-linear effects of the surrounding seafloor topography for bathymetry inversion. Therefore, GGM, which only considers the effects of linear terms, might lead to an unsatisfactory inversion in areas with significant topographic relief. However, EGGM addresses this issue by incorporating the non-linear effects of the surrounding seafloor topography, thus improving the inversion results.

5. CONCLUSIONS

In addressing the limitation of GGM, which considers only the linear terms in the inversion of seafloor topography and neglects the non-linear terms, this study focuses on deriving the non-linear effects of the surrounding seafloor topography on GAs at the sea surface. To address this issue, the self-adaptive algorithm EGGM is introduced, which effectively integrates the combined impact of various factors to a certain extent. The algorithm ensures, theoretically, comparable accuracy to that of GGM even when the GGM-based prior model is not ideal. Moreover, EGGM demonstrates a significant improvement in accuracy while maintaining good model stability. This enhancement in accuracy enables the use of larger iteration step sizes, resulting in reduced running time. Further analysis of both models under different conditions indicates that EGGM performs well across a range of water depths, with the extent of improvement being dependent on the accuracy of the prior model. Greater topographic relief leads to stronger non-linear effects from the surrounding seafloor topography, resulting in greater improvements with EGGM.

In the past, studies on seafloor topography inversion have predominantly focus on modelling the linear relationship between gravity data and the seafloor topography, with less attention given to the influence of non-linear effects. This study demonstrates the importance of non-linear effects in seafloor topography inversion and provides a theoretical foundation to research other methods. In our future works, incorporating more non-linear factors, such as geological conditions, into seafloor topography inversion is planned, with the aim of enhancing the theoretical framework of seafloor topography inversion and integrating it with practical applications.

ACKNOWLEDGEMENTS

The authors would like to thank the reviewers and editors for their valuable comments and suggestions and be grateful to the National Centers for Environmental Information (NCEI) for shipborne bathymetric data, the Scripps Institution of Oceanography for gravity data, the GEBCO Bathymetric Compilation Group 2022 for bathymetric data, David T. Sandwell for providing the unadjusted seafloor topography prediction. The study is supported by the National Natural Science Foundation of China (grant nos. 42192535, 42274006, 42104084, 42004009) and the Scientific Research Fund of Hunan Provincial Education Department (grant no. 23B0173).

DATA AVAILABILITY

The GA data are obtained from the Scripps Institution of Oceanography at https://topex.ucsd.edu/pub/global_grav_1min/. The GEBCO_2022 global terrain model provided by the Nippon Foundation-GEBCO Seabed 2030 Project was utilized for the pre-processing of shipborne bathymetric data at https://www.gebco.net/data_and_products/. The SIO_unadjusted model are obtained from David T. Sandwell. The single-beam bathymetric data were collected from the National Centers for Environmental Information (NCEI) of the National Oceanic and Atmospheric Administration (NOAA) at https://www.ncei.noaa.gov/maps/bathymetry/.

REFERENCES