the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Highprecision 1′ × 1′ bathymetric model of Philippine Sea inversed from marine gravity anomalies
Dechao An
Xiaotao Chang
Zhenming Wang
Yongjun Jia
Xin Liu
Valery Bondur
Heping Sun
The Philippine Sea, located at the edge of the northwestern Pacific Ocean, possesses complex seabed topography. Developing a highprecision bathymetric model for this region is of paramount importance, as it provides fundamental geoinformation essential for Earth observation and marine scientific research, including plate motion, ocean circulation, and hydrological characteristics. The gravity–geologic method (GGM), based on marine gravity anomalies, serves as an effective bathymetric prediction technique. To further strengthen the prediction accuracy of conventional GGM, we introduce the improved GGM (IGGM). The IGGM considers the effects of regional seafloor topography by employing weighted averaging to more accurately estimate the shortwavelength gravity component, along with refining the subsequent modeling of longwavelength gravity component. In this paper, we focus on seafloor topography modeling in the Philippine Sea based on the IGGM, combining shipborne bathymetric data with the Scripps Institution of Oceanography (SIO) V32.1 gravity anomaly. To reduce computational complexity, the optimal parameter values required for IGGM are first calculated before the overall regional calculation, and then, based on the terrain characteristics and distribution of sounding data, we selected four representative local sea areas as the research objects to construct the corresponding bathymetric models using GGM and IGGM. The analysis indicates that the precision of the IGGM models in four regions is improved to varying degrees, and the optimal calculation radius is 2^{′}. Based on the above finding, a highprecision ${\mathrm{1}}^{\prime}\times {\mathrm{1}}^{\prime}$ bathymetric model of the Philippine Sea (5–35° N, 120–150° E), known as the BAT_PS model, is constructed using IGGM. Results demonstrate that the BAT_PS model exhibits a higher overall precision compared to the General Bathymetric Chart of the Oceans (GEBCO), topo_25.1, and DTU18 models at singlebeam shipborne bathymetric points.
 Article
(16422 KB)  Fulltext XML
 BibTeX
 EndNote
The Philippine Sea serves as a convergence zone for continental and oceanic plates, resulting in frequent and intense plate tectonic activity (Richter and Ali, 2015; Lallemand, 2016; Holt et al., 2018). The Philippine Sea has complex topography, including trenches, island arcs, ridges, seamounts, basins, and rifts, exhibiting typical characteristics of a trench–arc–basin system. Additionally, it stands among the highly dynamic areas for geospatial exploration globally, attracting significant attention as a focal area for international scientific research. As fundamental geoinformation for marine scientific research, the accurate acquisition and application of highprecision seafloor topography are crucial (Kunze and Smith, 2004; Jena et al., 2012; Hirt and Rexer, 2015; Hu et al., 2020). They provide necessary information to guarantee an indepth comprehension of the marine environment and support marine development and governance (Ryabinin et al., 2019; Wolfl et al., 2019).
The Seabed 2030 project, a collaboration under the General Bathymetric Chart of the Oceans (GEBCO) and the Nippon Foundation, strives to achieve the creation of a comprehensive global mapping of the seafloor topography by 2030 (Mayer et al., 2018). In May 2023, the International Hydrographic Organization (IHO) announced that Seabed 2030 project had collected mapping data covering 24.9 % of the seabed. Additionally, multiple GEBCO models (such as GEBCO_2019, 2020, 2022, etc.) have been released, which integrate singlebeam shipborne bathymetric data, highresolution multibeam shipborne bathymetric data, and predicted bathymetry based on satellite altimetry data. For unexplored areas that have not been covered by bathymetric data, the marine gravity field derived from satellite altimetry data (Zhu et al., 2022; Zhou et al., 2023; Wei et al., 2023) is currently the main source for constructing highprecision and highresolution seafloor bathymetric models, compared to alternative observation methods (Bondur and Grebenyuk, 2001; Smith, 2004; Hilldale and Raff, 2008; Yeu et al., 2018; Tozer et al., 2019; Yu et al., 2022, 2023; Xu et al., 2023).
The prediction of seafloor topography based on satellite altimetry data has undergone a progression from onedimensional to twodimensional approaches, and various methods have been developed, such as the gravity–geology method (GGM), the frequency domain method (Parker, 1973; Watts, 1978; Smith and Sandwell, 1994), and the simulated annealing method (Yang et al., 2018, 2020). The GGM, proposed by Ibrahim and Hinze (1972), was initially used to measure the height of bedrock beneath sediments in land areas and subsequently employed in the inversion studies of seabed topography. Based on gravity anomalies derived from satellite altimetry missions, Smith and Sandwell (1994) constructed the corresponding bathymetric model and observed a significant correlation between seafloor topography and gravity anomaly within the wavelength band of 15 to 160 km. Yang et al. (2018) predicted seafloor topography in the western Pacific Ocean using vertical gravity gradient data through simulated annealing. Compared with the shipborne measured depth, the root mean square of the prediction result was 236 m, representing a 22 % improvement in accuracy over SIO (Scripps Institution of Oceanography) bathymetric model. However, the simulated annealing in the spatial domain includes forward and inverse modeling, requiring significant computational resources for modifying the initial model through continuous iteration. While the GGM and S&S (Smith and Sandwell, 1994) method use the linear relationship between seafloor topography and gravity data to construct empirical functions, the accuracy in regions with complex topography needs improvement. The frequency domain method is based on the spectral relation between seafloor topography and marine gravity field, using the firstorder term of the Parker formula to approximate bathymetry. This method omits the effect of higherorder terms and includes several complex geophysical parameters.
As one of the commonly used methods for bathymetric prediction, GGM can use shipborne bathymetry and gravity anomalies to build a bathymetric model with high accuracy (Nagarajan, 1994; Wei et al., 2021; Kim et al., 2010; Hsiao et al., 2011; Annan and Wan, 2020). The GGM divides the gravity anomaly into a shortwavelength component and longwavelength component. The shortwavelength gravity can be calculated using the Bouguer plate formula, while the estimation of the longwavelength gravity is particularly important. Annan and Wan (2020) adopted an adaptive mesh form to approximate the longwavelength gravity with a prediction accuracy of 180.20 m in the Gulf of Guinea, comparable to the accuracy of the ETOPO1 model and SIO model. In the conventional GGM, the process of calculating shortwavelength gravity from sea depth according to the Bouguer plate formula is a onetoone mapping, which ignores the gravity effect caused by surrounding seafloor topography. To overcome this limitation, An et al. (2022) proposed the improved GGM (IGGM), which considers the effect of regional seafloor topography. The IGGM finely calculates the shortwavelength gravity by weighted averaging and refines the subsequent longwavelength gravity modeling, significantly enhancing the overall accuracy of the bathymetric model.
Due to its prominent strategic position and complex marine environment, constructing an accurate bathymetric model of the Philippine Sea is of utmost importance. Therefore, this paper focuses on constructing a higherprecision bathymetric model for this region, utilizing the IGGM in combination with shipborne bathymetric data and the SIO V32.1 gravity anomaly. To effectively enhance the calculation efficiency of the bathymetric model and be guided by the topographic characteristics and the distribution of shipborne bathymetric data, four representative areas are chosen as the research objects before the overall solution. This preliminary stage involves exploring the optimal values of required parameters and assessing the applicability of IGGM for the Philippine Sea. Subsequently, a ${\mathrm{1}}^{\prime}\times {\mathrm{1}}^{\prime}$ bathymetric model of the Philippine Sea (5–35° N, 120–150° E) is obtained using a priori values of certain parameters. Finally, the accuracy of the bathymetric model is evaluated by comparing with existing models, including the GEBCO_2022, topo_25.1, and DTU18 models.
The Philippine Sea, located between the two island chains in the western Pacific Ocean, consists of the Philippine Basin, Parece Vela Basin, Shikoku Basin, Kyushu–Palau Ridge, Izu–Ogasawara Trench, Mariana Trench, and Mariana Trough. Recognizing the importance of the seafloor topography as essential geoinformation within this region, this paper focuses on the Philippine Sea (5–35° N, 120–150° E) as the study area, which presents significant challenges for shipborne bathymetry and the inversion of seafloor topography based on satellite altimetry data due to its complex terrain and geomorphology.
This paper used singlebeam shipborne bathymetry data obtained from the National Centers for Environmental Information (NCEI). However, due to the large time span over which the bathymetric data were acquired, some data collected before GPS became operational were often poorly localized and contained significant measurement errors (Smith, 1993). Therefore, it is necessary to use a bathymetric model with higher accuracy as a reference (e.g., GEBCO_2022 model) to ensure the quality of the shipboard bathymetric data. This involves interpolating the GEBCO_2022 model to the shipborne bathymetry points, calculating the difference with the shipborne depth, and then excluding points where the difference exceeds 3 times the standard deviation (SD) of the whole difference. The gravity anomalies were obtained from the SIO V32.1 model with a resolution of 1 arcmin, which was released in August 2022. The V32.1 model in the study area is shown in Fig. 1.
The GEBCO_2022 Grid, released in June 2022, is a continuous global terrain model with a grid spacing of 15 arcsec (GEBCO, 2022), briefly referred to as the GEBCO model. Over time, GEBCO has developed a range of bathymetric datasets and products. GEBCO provides a Type Identifier (TID) grid that specifies the type of source data used for each grid cell, including singlebeam shipborne bathymetry, multibeam shipborne bathymetry, lidar bathymetry, prediction depths derived from satellite gravity, and grid interpolation depths. By continually collecting and incorporating the latest and most relevant bathymetric data, the GEBCO series of models represents the cuttingedge level of accuracy for terrain models in the world.
Furthermore, the DTU18 and topo_25.1 bathymetric models were used in this paper, both of which have a grid resolution of 1 arcmin. The DTU18 model was released by the Technical University of Denmark (DTU) in 2019. SIO has been continuously updating its bathymetric models for an extended period, and the topo_25.1 model, released in January 2023, is the latest version of the current series.
3.1 Principle of the improved gravity–geologic method
In fact, there is a nonlinear relationship between seafloor topography and marine gravity anomalies. In many geodetic calculations, the nonlinear issue can be linearized by employing a suitable reference field (Hwang, 1999). The principle is to decompose the gravity anomaly into two main components: the regional gravity field representing the longwavelength gravity, and the residual gravity field representing the shortwavelength gravity. Among them, the longwavelength component is generated by deeper mass variations and the shortwavelength component is derived from variations in the local bedrock under sediment (Kim et al., 2010; Hsiao et al., 2016). The gravity anomaly (Δg) is composed of the longwavelength component (Δg_{reg}) and the shortwavelength component (Δg_{res}):
Using the known depth at control points to obtain the shortwavelength gravity component ($\mathrm{\Delta}{g}_{\mathrm{res}}^{{j}_{n}})$:
where G is the gravitational constant ($\mathrm{6.672}\times {\mathrm{10}}^{\mathrm{11}}$ m^{3} kg s^{−2}); Δρ is the density contrast (kg m^{−3}) between seawater and bedrock; ${E}_{{j}_{n}}$ is the depth at j_{n} point; and D is a reference datum, which is generally the deepest depth of control points.
The GGM uses the Bouguer slab Eq. (2) to compute the shortwavelength component of the sea surface point, which represents a singlepoint calculation form as shown in the left part of Fig. 2 and is not considered rigorous. In response to this problem, IGGM refines the longwavelength gravity model by considering the shortwavelength gravity effect of regional seafloor topography, as shown in the righthand side of Fig. 2. Using the control point (j_{n}) as the center, R is the calculated radius for estimating the effect of seafloor topography, m is the number of sounding points within the calculation range, and ${j}_{n}^{m}$ is the encompassing shipborne sounding point.
Based on Eq. (2), introducing cos ^{k}θ_{m} as the weight parameter, the shortwavelength gravity effect ($\mathrm{\Delta}{g}_{\mathrm{res}}^{{j}_{n}^{m}})$ of the surrounding shipborne points on the control point is
where θ is calculated by the arctangent value, based on the depth of the control point and the horizontal distance between the control point and the surrounding sea surface points. k is an unknown parameter determined by the iterative algorithm and is used in conjunction with cosθ_{m} as a weighting parameter to quantify shortwavelength gravity effects caused by the surrounding points.
Then, the shortwavelength component ($\mathrm{\Delta}{g}_{\mathrm{res}}^{{j}_{n}}$) is corrected using weighted averaging:
Subtracting the refined shortwavelength component calculated by Eq. (4) from the gravity anomaly to obtain the longwavelength component ($\mathrm{\Delta}{g}_{\mathrm{reg}}^{{j}_{n}}$) at the control point j_{n}:
The longwavelength component at control points is gridded using a tension spline function to obtain the longwavelength gravity field. The longwavelength component ($\mathrm{\Delta}{g}_{\mathrm{reg}}^{i}$) at any point i can be calculated through cubic spline interpolation. Subsequently, the shortwavelength component ($\mathrm{\Delta}{g}_{\mathrm{res}}^{i}$) at the prediction point i is
Finally, based on a variation in the Bouguer formula, the predicted depth (E_{i}) at point i is inversely calculated using the shortwavelength component:
3.2 Calculation process of improved gravity–geologic method
The key of IGGM is how to determine and calculate the optimal values of Δρ, k and calculation radius R, so that gravity anomalies can better characterize the basic geoinformation of seafloor topography. In this study, referring to the process of determining Δρ in GGM, the parameters Δρ, R and k are also calculated by iterative method in IGGM. The value ranges of Δρ, R and k are set as 0.1×10^{3} kg m^{−3}–5×10^{3} kg m^{−3}, 0^{′}–30^{′}, and 0.1–15, and their corresponding iteration steps are 0.1×10^{3} kg m^{−3}, 0.1^{′}, and 0.1, respectively. The correlation coefficient and SD between the predicted depth and the shipborne bathymetry at points i are analyzed under the conditions of different parameter values. When the correlation coefficient is the largest and the SD error is the smallest between predicted depths and measured depths (i.e., ${max}_{\mathrm{1}\le i\le n}C{C}_{i}$ and ${min}_{\mathrm{1}\le i\le n}{\mathrm{SD}}_{i}$), the corresponding parameters are the optimal values. In the iterative process of IGGM, the parameters are independent of each other. The optimal values of Δρ, R, and k are determined after all the iterations of the parameters have been completed. Figure 3 presents the flowchart of IGGM, providing a detailed visual representation of the process.
4.1 Parameter determination and applicability verification of IGGM
Due to the large size of the Philippine Sea and the substantial amount of bathymetric data, improving the calculation efficiency of the bathymetric model becomes essential. If certain unknown parameters can be determined before the overall model calculation, it can effectively achieve the purpose. In view of this, according to the different topographic features and bathymetric data distribution, four representative local areas within the Philippine Sea were first selected as the research subregions for analyzing and determining the unknown parameters in the IGGM. Then the optimal parameter values chosen afterwards were used to construct the final bathymetric model of the Philippine Sea. The four selected areas were as follows: A (21–27° N, 130–137° E); B (24–28° N, 138–145° E); C (8–15° N, 128–138° E); and D (10–16° N, 140–148° E), as shown in Fig. 4a. Area A is situated in the Daito Basin and the northern part of the Philippine Basin. Area B includes parts of the Shikoku Basin, the Izu–Ogasawara Islands arc and Trench, along with other adjacent sea areas. Area C is mainly situated in the Philippine Basin and the Parece Vela Basin, with less topographic relief, except for the Palau Ridge at the junction of the two basins. Area D is similar to area B and has complex terrain, including the Parece Vela Basin, and the Mariana Trench with depths exceeding 10 000 m. Table 1 provides the removal statistics of the four areas, and Fig. 4b illustrates the distribution of shipborne bathymetric data after gross error removing.
To verify the applicability of IGGM in the Philippine Sea, corresponding bathymetric models were constructed using GGM and IGGM in four areas. Multiple independent cruises were selected within each area, which were not involved in the model calculations, for independent verification purposes (represented by the black points in Fig. 4b) to ensure the reliability of the accuracy evaluation between GGM and IGGM. In both GGM and IGGM, the control points were divided into two parts. The first part was utilized for model calculation under different parameter values, while the second part was used for iteratively selecting the local optimal solutions of required parameters. Once the optimal values of each parameter were determined through iteration, the final bathymetric model was calculated using all the control points. This twostep process allowed for a comprehensive evaluation of the parameter values and ensured that the final bathymetric model was calculated by the best possible combination of shipborne bathymetric data, gravity anomalies, and optimized parameter values. For the division of control points in the iterative method, the conventional GGM still adopted the original proportional selection method, where one control point for iteration was selected at an interval of three points on each cruise, as shown in Fig. 5. On the other hand, IGGM divided control points using the cruise selection method, based on the distribution of shipborne cruises, as shown in Fig. 6. The proportional selection method tended to consider the entire area for obtaining the optimal values of the parameter. Contrastingly, the cruisedistributionbased selection weakened the effects caused by the interpolation of neighboring points.
This section mainly involves the following steps. First, the bathymetric models corresponding to the initial values of each parameter were calculated using green control points and were interpolated to obtain the predicted depth at white control points. Second, the correlation coefficient and SD of predicted depths and shipborne measured depths at white control points were calculated. The parameter values that yield the largest correlation coefficient and smallest SD (i.e., ${max}_{\mathrm{1}\le i\le n}C{C}_{i}$ and ${min}_{\mathrm{1}\le i\le n}{\mathrm{SD}}_{i}$) were considered optimal. Subsequently, the final bathymetric model was constructed by combining all control points and using the optimal parameter values. The predicted depths at checkpoints were then obtained through cubic spline interpolation, and the prediction accuracies for both GGM and IGGM were evaluated based on the shipborne bathymetric data. The optimal parameter values for GGM and IGGM in the four regions were presented in Table 2. It reveals that the optimal computational radius R of IGGM is 2^{′} in different areas, but there are differences in the optimal values of other Δρ and k. Additionally, the density contrast value in area D shows a significant difference. Figure 7 displays the total number of surrounding points within a 2^{′} radius centered on each shipborne point, with black points indicating no other singlebeam bathymetry point distribution. Statistical analysis shows that the numbers of the black points in the four areas are 0, 454, 131, and 130, respectively. The percentages of surrounding shipborne points exceeding 20 are 72.1 %, 59.1 %, 72.9 %, and 92.8 % for the four areas, respectively. The optimal selection of density contrast in both GGM and IGGM enables the gravity anomalies to better characterize the distribution of seafloor topography (Nagarajan, 1994; Hu et al., 2012; Kim and Yun, 2018). However, it should be noted that the density contrast is used as an empirical parameter only for obtaining the optimal accuracy of the bathymetry model, thus diminishing its original physical significance.
The corresponding bathymetric models constructed by GGM and IGGM were interpolated to obtain the predicted depths at the checkpoints. Table 3 displays the accuracy comparison between the predicted depth and the shipborne bathymetry in each area. Results indicate that the IGGM models in the four areas show varying degrees of accuracy improvement. The most significant enhancement is observed in areas A and D, where the SD is reduced by approximately 16.36 % (25.88 m) and 11.05 % (33.39 m). In areas B and C, the improvements of IGGM are limited to within 10 m, which is only slightly better than GGM models. The shipborne bathymetric data are evenly distributed in areas A and B, but area A is located in the Daito Basin with relatively slow topographic changes and exhibits the most substantial accuracy improvement. In contrast, there are Shikoku Basin and Izu–Ogasawara Island arc in area B, whose eastern side is located at the junction of the Philippine Sea plate and the Pacific Plate, and the topographic drop can reach up to thousands of meters. Terrain fluctuation is the main factor that affects the inversion accuracy, so the complex terrain in area B leads to limited accuracy improvement of IGGM. Area C, positioned in the western Philippine Basin and the Parece Vela Basin, has relatively gentle terrain. However, the scarcity of sounding bathymetric data in this region, results in a large error when predicting checkpoint depths using gridding and interpolation processing. Although area D includes the Mariana Trench with drastic topographic variations, it has a large number of shipborne bathymetric data.
The mean absolute percentage error (MAPE) is regarded as an indicator of relative accuracy and defined as the average value of the ratio (positive) of the prediction error to the measured depth:
where H_{i} is the shipborne measured depth, and E_{i} represents the predicted depth. A smaller MAPE indicates a higher accuracy of the bathymetric model. The result shows that IGGM models in the four regions have higher relative accuracy and can better characterize the geoinformation of seafloor topography. In conclusion, compared to GGM models, IGGM models effectively improve accuracy, with the degree of improvement affected by the distribution of shipborne bathymetry and terrain fluctuation. The significant improvement is observed in areas with flat terrain or sufficient and uniform distribution of bathymetric data. Conversely, less data and larger terrain fluctuation result in reduced accuracy improvement. Additionally, gridding and interpolation also introduce varying degrees of errors.
In areas A, B, and C, the GGM model and IGGM models show relatively consistent bathymetry performance overall, as shown in Fig. 8. However, notable differences are observed in area D, particularly concentrated in the vicinity of the Mariana Trench and Mariana Trough, such as 10–12° N, 143–145° E; 12–13° N, 140–141° E; 10–11° N, 140–141° E; and 15–16° N, 143–144° E. The GEBCO model, which incorporates the latest shipborne bathymetric data with high accuracy, was used as the standard to evaluate the accuracy levels of GGM and IGGM models. Figure 9a shows the distribution of seafloor topography in the GEBCO model. Figure 9b and c illustrate the differences between this model, the GGM model, and IGGM model, respectively. Since the GEBCO model has a grid interval of 15 arcsec, it needs to be interpolated to the corresponding grid points using the cubic spline. It is evident that the IGGM model exhibits smaller errors than the GGM model and is closer to the GEBCO model within four ranges marked by black boxes in Fig. 9a. Especially, in the vicinity of the Mariana Trench marked in the black box of Fig. 9a, the GGM model shows significant errors, further indicating the higherprecision performance of the IGGM model in area D as depicted in Fig. 8. Within the red box range of Fig. 9c, the IGGM model exhibits larger errors than the GGM model. This problem may be attributed to the lack of shipborne bathymetric data in this range, resulting in the optimal parameters selected by IGGM being less applicable to a certain local area. According to statistics, differences between the GEBCO model and the GGM and IGGM models within the range of 0–300 m account for 86.2 % and 88.8 %, while the proportions of the differences above 500 m are 6.1 % and 3.5 %, respectively. The GEBCO model, despite having a grid interval of 15 arcsec, exhibits lower true spatial resolution in marine regions that lack of measured bathymetric data, as illustrated in the basemap of area C in Fig. 6.
Based on the above statistics, it is evident that the IGGM, which considers the effect of regional seafloor topography, demonstrates higher accuracy at checkpoints and achieves a significant improvement in the Philippine Sea. These results also confirm that the IGGM has better applicability. For all parameters required in IGGM, the optimal calculation radius in the four areas remains consistent at 2^{′}, which is not affected by the distribution of seafloor topography and shipborne bathymetric data, but the optimal values of Δρ and the index k in the weight factor differ in each area. Given these findings, a predetermined calculation radius of 2^{′} is adopted for the final bathymetry modeling of the Philippine Sea (5–35° N, 120–150° E).
4.2 Bathymetric model based on IGGM in the Philippine Sea
Considering the huge amount of shipborne bathymetric data in the entire Philippine Sea and the limitations in computer capacity, the region (5–35° N, 120–150° E) is divided into nine subregions of 10°×10° each, as shown in Fig. 10a. To ensure data quality, the shipborne bathymetric data within each subregion are preprocessed by removing gross errors that exceed 3σ, with the GEBCO model used as a reference for comparison. Following the procedure of IGGM, the bathymetric models of each subregion are inverted. Finally, the ${\mathrm{1}}^{\prime}\times {\mathrm{1}}^{\prime}$ bathymetric model of the Philippine Sea (BAT_PS) is finally obtained, as depicted in Fig. 10b. During the construction of the BAT_PS model, the optimal parameters values required in IGGM are given in Table 4.
The accuracy of the BAT_PS model was evaluated using three reference models: GEBCO, topo_25.1, and DTU18. The BAT_PS, GEBCO, topo_25.1, and DTU18 models were interpolated to all shipborne bathymetry points within the study area. As all the singlebeam bathymetric data were involved in the model calculation, they could serve as an evaluation of internal consistency accuracy. The statistical results were presented in Table 5, with two categories of assessment: primary check (PC) and second check (SC). The primary check is the analysis of the difference between the predicted depth and the measured depth at all shipborne measured points. In the second check, the points with poor prediction quality based on the 3σ criterion are removed before conducting the statistical analysis. The results demonstrate that the BAT_PS model exhibits better agreement with the shipborne bathymetric data in the primary check, with SD and MAPE values of 89.56 m and 9.79 %, outperforming the other three models. MAE represents the absolute magnitude of prediction errors, while MAPE reflects the percentage magnitude of prediction errors relative to the shipborne bathymetry. In the second check after removing gross errors, the accuracy levels of the BAT_PS model and the GEBCO_2022 model are approximately equal, surpassing topo_25.1 model and DTU18 model. Overall, the accuracy ranking from high to low is BAT_PS model, GEBCO model, topo_25.1 model, and DTU18 model.
In the statistics of two check accuracies, both before and after gross error elimination, the DTU18 model outperforms the topo_25.1 model, with MAPE values of 20.88 % and 20.26 %, compared to the topo_25.1 model's values of 29.97 % and 27.19 %. In the second checks, the topo_25.1 model has excluded poorquality points, with both the SD and RMS showing better performance compared to the DTU18 model. However, it is noteworthy that the MAPE of topo_25.1 is still worse. The topo_25.1 model may exhibit a large relative error in shallow seas. Figure 11 shows the relationship between the relative errors between four models and different depths. The relative error is described as the percentage of the prediction error and measured depth. MAPE can be regarded as a comprehensive statistical measure of relative error. As the depth increases, the relative errors within each model exhibit a decreasing trend.
Table 6 presents the statistics of MAPE at different depths after the second check. Specifically, the topo_25.1 model has a large relative error within a depth of 2000 m, which explains the phenomenon discussed above. After reaching a depth greater than 2000 m, the relative error in the topo_25.1 model gradually becomes smaller than that of the DTU18 model. The construction of the BAT_PS model only relies on singlebeam data without the inclusion of multibeam bathymetric data. Consequently, the accuracy in some deepsea areas is slightly lower. Additionally, the relative errors in each model increases when the depth exceeds −8000 m. This can be attributed to most of these areas being near trenches, where both the accuracy of shipborne bathymetry and the precision of inversion methods face challenges due to the abrupt terrain.
To further improve the accuracy of the BAT_PS model and scale down its prediction error, we use shipboard bathymetry data for the final constraints in the following steps (SRTM, 2019): (1) interpolate the BAT_PS model to obtain predicted depths at ship measurement points and calculate the difference from the actual measured depths; (2) supplement grid points located 5^{′} away from the ship measurement points with zerodepth differences; (3) use the GMTs (Generic Mapping Tools) module surface to generate a corrected grid (Fig. 12a) by combining the depth differences at ship measurement points and the zero values at grid points; and (4) restore the corrected grid to the PS model to obtain the CBAT_PS (constrained BAT_PS) model, as shown in Fig. 12b. About 96.29 % of the corrected values were counted to be within the range of −100–100 m, indicating that the overall accuracy of the BAT_PS model is still reliable. Finally, after being constrained by shipborne bathymetry, the CBAT_PS model has also been uploaded.
Topography, especially bathymetric topography, has been a significant research focus, with unique applications in the field of geoscience. Although multiple techniques for bathymetric inversion and their corresponding improvement methods are temporarily constructed, a fully rigorous bathymetric inversion theory is yet to be established. Consequently, it is necessary to refine and perfect the original method based on the current theory, serving as the foundation and focus of future bathymetric inversion research. In this paper, the IGGM, which considers the shortwavelength component effect of regional seafloor topography, is used to predict bathymetry in the Philippine Sea. The primary conclusions are summarized as follows:

Four local areas in the Philippine Sea were selected as research subjects, and corresponding bathymetric models were calculated using both GGM and IGGM. Overall, the accuracy of IGGM at checkpoints outperformed that of GGM, with improvements of 16.36 %, 3.31 %, 3.80 %, and 11.05 %. This verifies the applicability of IGGM in the Philippine Sea. As a result, the optimal calculation radius of IGGM was set at 2^{′} for the subsequent construction of the overall bathymetric model of the Philippine Sea.

A highprecision BAT_PS model, with a grid of 1 arcmin, was established for the Philippine Sea (5–35° N, 120–150° E) using IGGM. The SD error between the predicted depths of the BAT_PS model and shipbornemeasured depths reached 89.56 m at the singlebeam shipborne points, and the accuracy of the BAT_PS model was 59.91 m without poor quality points. This accuracy is essentially equivalent to that of the GEBCO_2022 model and is significantly better than the topo_25.1 and DTU18 models. In addition, we use the shipborne bathymetric data to constrain the BAT_PS model and make the final CBAT_PS model.

In this study, only singlebeam bathymetric data were used for the construction of the BAT_PS model, while multibeam bathymetric data are abundant in certain areas, such as the western Philippine Sea basin, the Mariana Trench, and Mariana Trough. Future research will explore integrating multibeam data into the model construction to further enhance the accuracy of the BAT_PS model. Additionally, the accuracy of seafloor terrain inversion heavily relies on the ocean gravity field. The observation data from the Surface Water and Ocean Topography (SWOT) swath altimetry satellite are hoped to further enhance the resolution and accuracy of the marine gravity field, thereby greatly improving the prediction precision of bathymetry.
All data used in this study are publicly available through the NCEI (https://www.ncei.noaa.gov/maps/bathymetry/, NOAA National Centers for Environmental Information, 2023), GEBCO (https://doi.org/10.5285/e0f0bb80ab442739e0536c86abc0289c, GEBCO Bathymetric Compilation Group, 2022), SIO (https://topex.ucsd.edu/pub/global_grav_1min/, Scripps Institution of Oceanography, 2023a; https://topex.ucsd.edu/pub/global_topo_1min/ Scripps Institution of Oceanography, 2023b), and DTU (https://ftp.space.dtu.dk/pub/DTU18/1_MIN/, Technical University of Denmark, 2023). The BAT_PS bathymetric model, CBAT_PS model, and source codes are available from https://doi.org/10.5281/zenodo.10370469 (An, 2023).
Conceptualization: DA and JG. Methodology: DA and JG. Validation: DA, JG, XC, ZW, YJ, XL, VB, and HS. Writing – original draft: DA. Review and editing: DA, JG, XC, ZW, YJ, XL, VB, and HS. All authors contributed to writing and revising the paper.
The contact author has declared that none of the authors has any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
We express our gratitude to the following organizations: NCEI for providing shipborne bathymetric data, GEBCO for offering bathymetric model, SIO for sharing gravity anomaly model and bathymetric model, and the DTU for providing bathymetric model. We also extend our appreciation to the editors and reviewers for their invaluable comments.
This research has been supported by the National Natural Science Foundation of China (grant nos. 42192535, 42274006, and 42242015) and the National Key Research and Development Program of China (grant no. YS2018YFE011371).
This paper was edited by Dan Lu and reviewed by Richard Fiifi Annan and one anonymous referee.
An, D.: Highprecision 1′×1′ bathymetric model of Philippine Sea inversed from marine gravity anomalies, Zenodo [code], https://doi.org/10.5281/zenodo.10370469, 2023.
An, D., Guo, J., Li, Z., Ji, B., Liu, X., and Chang, X.: Improved gravitygeologic method reliably removing the longwavelength gravity effect of regional seafloor topography: A case of bathymetric prediction in the South China Sea, IEEE T. Geosci. Remote, 60, 4211912, https://doi.org/10.1109/TGRS.2022.3223047, 2022.
Annan, R. F. and Wan, X.: Mapping seafloor topography of Gulf of Guinea using an adaptive meshed gravitygeologic method, Arab. J. Geosci., 13, 301, https://doi.org/10.1007/s12517020052978, 2020.
Bondur, V. G. and Grebenyuk, Y. V.: Remote sensing methods for determining the bottom relief of coastal zones of seas and oceans, Mapp. Sci. Remote Sens., 38, 172–190, https://doi.org/10.1080/07493878.2001.10642174, 2001.
GEBCO Bathymetric Compilation Group 2022: The GEBCO_2022 Grid – A continuous terrain model of the global oceans and land, NERC EDS British Oceanographic Data Centre NOC [data set], https://doi.org/10.5285/e0f0bb80ab442739e0536c86abc0289c, 2022.
Hilldale, R. C. and Raff, D.: Assessing the ability of airborne LiDAR to map river bathymetry, Earth Surf. Proc. Land., 33, 773–783, https://doi.org/10.1002/esp.1575, 2008.
Hirt, C. and Rexer, M.: Earth2014: 1 arcmin shape, topography, bedrock and icesheet models – Available as gridded data and degree10,800 spherical harmonics, Int. J. Appl. Earth Obs., 39, 103–112, https://doi.org/10.1016/j.jag.2015.03.001, 2015.
Holt, A. F., Royden, L. H., Becker, T. W., and Faccenna, C.: Slab interactions in 3D subduction settings: The Philippine Sea Plate region, Earth Planet. Sc. Lett., 489, 72–83, https://doi.org/10.1016/j.epsl.2018.02.024, 2018.
Hsiao, Y. S., Kim, J. W., Kim, K. B., Lee, B. Y., and Hwang, C.: Bathymetry estimation using the gravity geologic method: An investigation of density contrast predicted by the downward continuation method, Terr. Atmos. Ocean. Sci., 22, 347–358, https://doi.org/10.3319/TAO.2010.10.13.01(OC), 2011.
Hsiao, Y. S., Hwang, C., Cheng, Y., Chen, L., Hsu, H., Tsai, J., Liu, C., Wang, C., Liu, Y., and Kao, Y.: Highresolution depth and coastline over major atolls of South China Sea from satellite altimetry and imagery, Remote Sens. Environ., 176, 69–83, https://doi.org/10.1016/j.rse.2016.01.016, 2016.
Hu, M., Li, J., and Jin, T.: Bathymetry Inversion with GravityGeologic Method in Emperor Seamount, Geomat. Inf. Sci. Wuhan Univ., 37, 610–612, https://doi.org/10.13203/j.whugis2012.05.008, 2012.
Hu, Q., Huang, X., Zhang, Z., Zhang, X., Xu, X., Sun, H., Zhou, C., Zhao, W., and Tian, J.: Cascade of internal wave energy catalyzed by eddytopography interactions in the deep South China Sea, Geophys. Res. Lett., 47, e2019GL086510, https://doi.org/10.1029/2019GL086510, 2020.
Hwang, C.: A bathymetric model for the South China Sea from satellite altimetry and depth data, Mar. Geod., 22, 37–51, https://doi.org/10.1080/014904199273597, 1999.
Ibrahim, A. and Hinze W. J.: Mapping buried bedrock topography with gravity, Ground Water, 10, 18–23, https://doi.org/10.1111/j.17456584.1972.tb02921.x, 1972.
Jena, B., Kurian, P. J., Swain, D., Tyagi, A., and Ravindra, R.: Prediction of bathymetry from satellite altimeter based gravity in the Arabian Sea: Mapping of two unnamed deep seamounts, Int. J. Appl. Earth Obs., 16, 1–4, https://doi.org/10.1016/j.jag.2011.11.008, 2012.
Kim, K. B. and Yun, H. S.: Satellitederived bathymetry prediction in shallow waters using the gravitygeologic method: A case study in the west sea of Korea, KSCE J. Civil Eng., 22, 2560–2568, https://doi.org/10.1007/s122050170487z, 2018.
Kim, K. B., Hsiao, Y. S., Kim, J. W., Lee, B. Y., Kwon, Y. K., and Kim, C. H.: Bathymetry enhancement by altimetryderived gravity anomalies in the East Sea (Sea of Japan), Mar. Geophys. Res., 31, 285–298, https://doi.org/10.1007/s1100101091100, 2010.
Kunze, E. and Smith, S. G. L.: The role of smallscale topography in turbulent mixing of the global ocean, Oceanography, 17, 55–64, https://doi.org/10.5670/oceanog.2004.67, 2004.
Lallemand, S.: Philippine Sea Plate inception, evolution, and consumption with special emphasis on the early stages of IzuBoninMariana subduction, Prog. Earth Planet. Sci., 3, 15, https://doi.org/10.1186/s4064501600856, 2016.
Mayer, L., Jakobsson, M., Allen, G., Dorschel, B., Falconer, R., Ferrini, V., Lamarche, G., Snaith, H., and Weatherall, P.: The Nippon foundation – GEBCO seabed 2030 project: The quest to see the world's oceans completely mapped by 2030, Geosciences, 8, 63, https://doi.org/10.3390/geosciences8020063, 2018.
Nagarajan, R.: Gravitygeologic investigation of buried bedrock topography in northwestern Ohio, Ohio, The Ohio State University, 1994.
NOAA National Centers for Environmental Information: Seafloor Mapping [data set], https://www.ncei.noaa.gov/maps/bathymetry/, last access: 16 September 2023.
Parker, R. L.: The rapid calculation of potential anomalies, Geophys. J. Int., 31, 447–455, https://doi.org/10.1111/j.1365246X.1973.tb06513.x, 1973.
Richter, C. and Ali, J. R.: Philippine Sea Plate motion history: Eocenerecent record from ODP Site 1201, central West Philippine Basin, Earth Planet. Sc. Lett., 410, 165–173, https://doi.org/10.1016/j.epsl.2014.11.032, 2015.
Ryabinin, V., Barbiere, J., Haugan, P., Kullenberg, G., Smith, N., McLean, C., Troisi, A., Fischer, A., Arico, S., Aarup, T., Pissierssens, P., Visbeck, M., Enevoldsen, H. O., and Rigaud, J.: The UN decade of ocean science for sustainable development, Front. Mar. Sci., 6, 470, https://doi.org/10.3389/fmars.2019.00470, 2019.
Scripps Institution of Oceanography: V32.1 gravity anomaly model [data set], https://topex.ucsd.edu/pub/global_grav_1min/ (last access: 16 September 2023), 2023a.
Scripps Institution of Oceanography: topo_25.1 bathymetric model [data set], https://topex.ucsd.edu/pub/global_topo_1min/ (last access: 16 September 2023), 2023b.
Smith, W. H. F.: On the accuracy of digital bathymetric data, J. Geophys. Res.Sol. Ea., 98, 9591–9603, https://doi.org/10.1029/93JB00716, 1993.
Smith, W. H. F.: Introduction to this special issue on bathymetry from space, Oceanography, 17, 6–7, https://doi.org/10.5670/oceanog.2004.62, 2004.
Smith, W. H. F. and Sandwell, D. T.: Bathymetric prediction from dense satellite altimetry and sparse shipboard bathymetry, J. Geophys. Res.Sol. Ea., 10, 21803–21824, https://doi.org/10.1029/94JB00988, 1994.
Technical University of Denmark: DTU18 bathymetric model [data set], https://ftp.space.dtu.dk/pub/DTU18/1_MIN/, last access: 16 September 2023.
Tozer, B., Sandwell, D. T., Smith, W. H. F., Olson, C., Beale, J. R., and Wessel, P.: Global bathymetry and topography at 15 arc sec: SRTM15+, Earth Space Sci., 6, 1847–1864, https://doi.org/10.1029/2019EA000658, 2019.
Watts, A. B.: An analysis of isostasy in the world's oceans 1. HawaiianEmperor seamount chain, J. Geophys. Res.Sol. Ea., 83, 5989–6004, https://doi.org/10.1029/JB083iB12p05989, 1978.
Wei, X., Liu, X., Li, Z., Chang, X., Luo, H., Zhu, C., and Guo, J.: Gravity anomalies determined from mean sea surface model data over the Gulf of Mexico, Acta Oceanol. Sin., 42, 52–70, https://doi.org/10.1007/s1313102321786, 2023.
Wei, Z., Guo, J., Zhu, C., Yuan, J., Chang, X., and Ji, B.: Evaluating accuracy of HY2A/GMderived gravity data with the gravitygeologic method to predict bathymetry, Front. Earth Sci., 9, 636246, https://doi.org/10.3389/feart.2021.636246, 2021.
Wolfl, A.C., Snaith, H., Amirebrahim, S., Devey, C. W., Dorschel, B., Ferrini, V., Huvenne, V. A. I., Jakobsson, M., Jencks, J., Johnston, G., Lamarche, G., Mayer, L., Millar, D., Pedersen, T. H., Picard, K., Reitz, A., Schmitt, T., Visbeck, M., Weatherall, P., and Wigley, R.: Seafloor mapping – the challenge of a truly global ocean bathymetry, Front. Mar. Sci., 6, 283, https://doi.org/10.3389/fmars.2019.00283, 2019.
Xu, C., Li, J., Jian, G., Wu, Y., and Zhang, Y.: An adaptive nonlinear iterative method for predicting seafloor topography from altimetryderived gravity data, J. Geophys. Res.Sol. Ea., 128, e2022JB025692, https://doi.org/10.1029/2022JB025692, 2023.
Yang, J., Jekeli, C., and Liu, L.: Seafloor topography estimation from gravity gradients using simulated annealing, J. Geophys. Res.Sol. Ea., 123, 6958–6975, https://doi.org/10.1029/2018JB015883, 2018.
Yang, J., Luo, Z., and Tu, L.: Ocean access to Zachariæ Isstrøm glacier, northeast Greenland, revealed by OMG airborne gravity, J. Geophys. Res.Sol. Ea., 125, e2020JB020281, https://doi.org/10.1029/2020JB020281, 2020.
Yeu, Y., Yee, J., Yun, H. S., and Kim K. B.: Evaluation of the accuracy of bathymetry on the nearshore coastlines of western Korea from satellite altimetry, multibeam, and airborne bathymetric LiDAR, Sensors, 18, 2926, https://doi.org/10.3390/s18092926, 2018.
Yu, J., An, B., Xu, H., Sun, Z., Tian, Y., and Wang, Q.: An iterative algorithm for predicting seafloor topography from gravity anomalies, Remote Sens., 15, 1069, https://doi.org/10.3390/rs15041069, 2023.
Yu, L., Du, J., Zhai, R. Wu, F., and Qian, H.: A fast generalization method of multibeam echo soundings for nautical charting, J. Geovisualization Spat. Anal., 6, 2, https://doi.org/10.1007/s41651021000965, 2022.
Zhou, R., Liu, X., Li, Z., Sun, Y., Yuan, J., Guo, J., and Ardalan, A. A.: On performance of vertical gravity gradient determined from CryoSat2 altimeter data over Arabian Sea, Geophys. J. Int., 234, 1519–1529, https://doi.org/10.1093/gji/ggad153, 2023.
Zhu, C., Guo, J., Yuan, J., Li, Z., Liu, X., and Gao, J.: SDUST2021GRA: global marine gravity anomaly model recovered from Kaband and Kuband satellite altimeter data, Earth Syst. Sci. Data, 14, 4589–4606, https://doi.org/10.5194/essd1445892022, 2022.