A New × Global Seaﬂoor Topography Model Predicted from Satellite Altimetric Vertical Gravity Gradient Anomaly and Ship Soundings BAT_VGG2021

: In this paper, we construct a new 1 (cid:48) × 1 (cid:48) global seaﬂoor topography model, BAT_VGG2021, using the satellite altimetric vertical gravity gradient anomaly model (VGG), SIO curv_30.1.nc, and ship soundings. Approximately 74.66 million single-beam depths and more than 180 GB of multibeam grids were downloaded and adopted from the National Centers for Environmental Information (NCEI), Japan Agency Marine-Earth Science and Technology (JAMSTEC), and Geosciences Australia (GA). The SIO curv_30.1.nc model was used to predict seaﬂoor topography at 15~160 km wavelengths, and ship soundings were used to calibrate topography to VGG ratios. The accuracy of the new BAT_VGG2021 model was assessed by comparing it with ship soundings and existing models. The results indicate that the standard deviation of differences between the predicted model and ship soundings is about 40~80 m, and ~93% of the differences are within 100 m, similar to that of the SIO topo_20.1.nc model. The new BAT_VGG2021 model shows better accuracy than the DTU18BAT, ETOPO1, and GEBCO_08 models, and has been improved signiﬁcantly from our last model, BAT_VGG2014. ((cid:1863)) . The results show that the transform functions work like band-pass filters. The transform function, (cid:1852) (cid:3047)(cid:3042)(cid:3043)(cid:3042)(cid:2879)(cid:3034)(cid:3045)(cid:3028)(cid:3031) ((cid:1863)) , suppresses the effect of isostasy and enlarge signal at wavelengths shorter than ~100 km. topography model at wavelengths longer than 160 km. The ship soundings were cleaned by comparing with topo_20.1.nc. About 90% of the cleaned ship points were applied to constrain topography to VGG ratios at 15~160 km wavelength bands, and 10% were used to assess the model accuracy. The VGG model, curv_30.1.nc was used to predict seaﬂoor topography model at 15~160 km wavelength bands. The VGG model, curv_30.1.nc was used to predict seafloor topography model at 15~160 km wavelength bands.


Introduction
The seafloor covers~71% of the solid earth and has diverse tectonic features. Knowledge of the seafloor topography (seawater depth) plays a pivotal role in geosciences research, resource exploitation, and environmental protection, etc. However, mapping the seafloor requires a significant investment of labor and money.
Traditionally, seafloor topography can be measured directly by echo sounder systems equipped on a ship. However, it is time-consuming work to measure the global seafloor due to the slow speed of the ship. Scientists have suggested that people have surveyed only~18% of the seafloor, at an effective resolution of~1 km, in the past centuries [1]. In the latest SRTM+V2.1 seafloor topography model at 15 arc-second resolution, only~10.84% of the seafloor is directly constrained by acoustic surveys [2]. It is assessed that hundreds of ship-years and billions in terms of financial cost are required to map the global seafloor using echo sounding systems [1,3]. Thus, it is very difficult to map the oceans at a proper resolution directly by ship currently.
Fortunately, the development of satellite altimetry technology provides an indirect way to recover seafloor features at a moderate resolution and accuracy. Altimeters measure the height of the ocean's surface, which can be used to derive gravity anomaly. Then the gravity anomaly can be used to predict seafloor topography at short-to-middle 2 of 17 wavelengths [4][5][6]. Since the launching of the Skylab in 1973, many satellite altimetric missions, such as Geosat, ERS, and Jason, have been carried out successfully. The altimetric derived gravity models have been updated constantly with the accumulation of satellite altimetry data and improvements in the data processing technology [7][8][9][10][11][12]. At the same time, a series of global seafloor topography models, constructed with altimetric data and ship soundings, were published by the Scripps Institute of Oceanography (SIO), Danmarks Tekniske Universitet (DTU), Wuhan University (WHU), and the International Hydrographic Organization (IHO), etc. [2,[4][5][6][13][14][15][16].
At present, most of the publicly available seafloor topography models, such as SIO topo_20.1.nc and DTU18BAT, have been predicted from altimetric gravity anomalies. While scientists have suggested that the vertical gravity gradient anomaly (VGG) may be used to strengthen seafloor topography at short wavelengths [17], few papers have been published [18][19][20][21]. For more than ten years, we have engaged in seafloor topography model construction using ship soundings and VGG and published a global model, BAT_VGG2014, in 2014 [22]. Presently, new altimetry data from AltiKa, CryoSat, and Sentinel-3A/B have been collected to derive a new version of the VGG model, SIO curv_30.1.nc. In this paper, we construct a new 1 × 1 global seafloor topography model using ship soundings and the latest version of the VGG model. Approximately 74.66 million single-beam depths and 6.6 GB of multibeam grids were downloaded from the National Centers for Environmental Information (NCEI),~120 GB of multibeam grids were downloaded from the Japan Agency for Marine-earth Science and Technology (JAMSTEC), and~54 GB of multibeam grids were downloaded from Geosciences Australia (GA). The accuracy of the new model was assessed by comparing it with ship soundings and existing models such as SIO topo_20.1.nc, DTU18BAT, ETOPO1, GEBCO_08, and BAT_VGG2014 [23,24].

Theory and Methods
According to the lithospheric flexural isostasy theory [25], seamounts loading will introduce flexure of the oceanic crust (Moho discontinuity) where R(k) and H(k) represent the Fourier transform of Moho flexure and seafloor undulations; ρ m , ρ c , ρ w are densities of the mantle, crust, and water, respectively; k = 2π/λ is the wavenumber, and λ is the wavelength; Φ e (k) is the lithospheric flexural response function [26] Φ e (k)= Dk 4 (ρ m − ρ c )g + 1 −1 (2) where D is lithospheric flexural rigidity, D = ET 3 e /[12(1 − υ 2 )], and E is Young's modulus, T e is lithospheric effective elastic thickness; υ is Poisson's ratio; g is the average gravity acceleration.
Seamounts and the corresponding Moho flexure introduce most of the gravity anomalies on the sea surface. Based on Parker's formula [27], these gravity anomalies, ∆G(k), can be given by where G is gravitation constant, d is mean water depth, and t is mean crustal thickness. This series expansion formula converges very quickly. Substituting Equation (1) into Equation (3), and considering n = 1, Equation (3) can be simplified as In the frequency domain, the VGG, ∆G z (k), will be Thus, the transform functions between seafloor and gravity anomaly, These transform functions are composed by coefficient, 2πG(ρ c − ρ w ), exponential decay function, exp(k), isostatic response function, ϕ(k), and wavenumber, k, where Using parameters in Table 1, the exponential decay function, exp(k), works like a lowpass filter (thick blue line in Figure 1), and suggests high-frequency topography signal is suppressed due to upward-continuation of seawater depth. The isostatic response function, ϕ(k), works as a high-pass filter (red lines in Figure 1), and indicates the gravity signal is weakened by isostatic compensation mass in the oceanic crust. The combined effect of these functions results in a band-pass filter of transform function, as shown in Figure 2. Figure 2 indicates that compared to Z topo−grav (k), the transform function between seafloor topography and VGG, Z topo−g rad (k), suppresses the effect of isostasy and enlargesthe signal at shorter wavelengths (<~100 km). where G is gravitation constant, d is mean water depth, and t is mean crustal thickness. This series expansion formula converges very quickly. Substituting Equation (1) into Equation (3), and considering n = 1, Equation (3) can be simplified as In the frequency domain, the VGG, ∆ ( ), will be Thus, the transform functions between seafloor and gravity anomaly, These transform functions are composed by coefficient, 2 ( − ), exponential decay function, exp(k), isostatic response function, ( ), and wavenumber, k, where Using parameters in Table 1, the exponential decay function, exp(k), works like a lowpass filter (thick blue line in Figure 1), and suggests high-frequency topography signal is suppressed due to upward-continuation of seawater depth. The isostatic response function, ( ), works as a high-pass filter (red lines in Figure 1), and indicates the gravity signal is weakened by isostatic compensation mass in the oceanic crust. The combined effect of these functions results in a band-pass filter of transform function, as shown in Figure 2. Figure 2 indicates that compared to ( ), the transform function between seafloor topography and VGG, ( ), suppresses the effect of isostasy and enlargesthe signal at shorter wavelengths (<~100 km).   Isostasy analysis of seafloor features suggested that the seafloor topography shows high coherence with gravity anomalies at certain short-to-middle wavelength bands [28,29]. For example, in the northwestern Pacific (144°~180°E, 0°~36°N), Figure 3 indicates that the seafloor topography and gravity anomaly or VGG show high coherence at certain wavelength bands. Thus, the satellite altimetric gravity anomalies or VGG were used to constrain seafloor topography at 10~160 km or 20~200 km wavelength bands [4,5,19].  0°~36°N). This indicates that the seafloor topography and gravity anomaly or VGG show high coherence at certain wavelength bands. At long wavelengths (>500 km), the topography-VGG coherence is lower than the topography-GA coherence.
A band-pass filter was designed to process gravity anomaly or VGG, and seafloor topography at certain wavelength bands can be predicted by where BF indicates a band-pass filter. Isostasy analysis of seafloor features suggested that the seafloor topography shows high coherence with gravity anomalies at certain short-to-middle wavelength bands [28,29]. For example, in the northwestern Pacific (144 •~1 80 • E, 0 •~3 6 • N), Figure 3 indicates that the seafloor topography and gravity anomaly or VGG show high coherence at certain wavelength bands. Thus, the satellite altimetric gravity anomalies or VGG were used to constrain seafloor topography at 10~160 km or 20~200 km wavelength bands [4,5,19].  Isostasy analysis of seafloor features suggested that the seafloor topography shows high coherence with gravity anomalies at certain short-to-middle wavelength bands [28,29]. For example, in the northwestern Pacific (144°~180°E, 0°~36°N), Figure 3 indicates that the seafloor topography and gravity anomaly or VGG show high coherence at certain wavelength bands. Thus, the satellite altimetric gravity anomalies or VGG were used to constrain seafloor topography at 10~160 km or 20~200 km wavelength bands [4,5,19]. The coherence between seafloor topography and gravity anomaly (black dots) or VGG (red dots) in the northwestern Pacific (144°~180°E, 0°~36°N). This indicates that the seafloor topography and gravity anomaly or VGG show high coherence at certain wavelength bands. At long wavelengths (>500 km), the topography-VGG coherence is lower than the topography-GA coherence.
A band-pass filter was designed to process gravity anomaly or VGG, and seafloor topography at certain wavelength bands can be predicted by where BF indicates a band-pass filter.
. This indicates that the seafloor topography and gravity anomaly or VGG show high coherence at certain wavelength bands. At long wavelengths (>500 km), the topography-VGG coherence is lower than the topography-GA coherence.
A band-pass filter was designed to process gravity anomaly or VGG, and seafloor topography at certain wavelength bands can be predicted by where BF indicates a band-pass filter.

Data Sources
In this paper, the VGG model, existing seafloor topography model, and ship soundings were used to construct a new global seafloor topography model. The latest version of the SIO VGG model, curv_30.1.nc, which includes new altimetric data from AltiKa, CryoSat, and Sentinel-3A/B satellites, was used to constrain seafloor topography at 15~160 km wavelength bands, as shown in Figure 4. The SIO topography model, topo_20.1.nc, was used to control seafloor topography at wavelengths longer than 160 km. The ship soundings were used to calibrate topography-to-VGG ratios at 15~160 km wavelength bands.

Data Sources
In this paper, the VGG model, existing seafloor topography model, and ship soundings were used to construct a new global seafloor topography model. The latest version of the SIO VGG model, curv_30.1.nc, which includes new altimetric data from AltiKa, CryoSat, and Sentinel-3A/B satellites, was used to constrain seafloor topography at 15~160 km wavelength bands, as shown in Figure 4. The SIO topography model, topo_20.1.nc, was used to control seafloor topography at wavelengths longer than 160 km. The ship soundings were used to calibrate topography-to-VGG ratios at 15~160 km wavelength bands. The ship soundings, including single-beam points and multibeam grids, were downloaded from NCEI, JAMSTEC, and GA, as show in Figure 5. The NCEI provided ~74.66 million single-beam points (black dots) in more than 3800 cruise surveys collected since 1960Sand ~6.6GB of multibeam grids (blue dots). Approximately 120GB of multibeam grids were downloaded from the JAMSTEC (yellow dots), with most of the data distributed around Japan and the northwestern Pacific. Multibeam grids, ~54GB, around Australia were obtained from GA (purple dots). These ship depths distribute unevenly over the world, and most of them are in the northern hemisphere. All data sources are listed in Table 2. The ship soundings, including single-beam points and multibeam grids, were downloaded from NCEI, JAMSTEC, and GA, as show in Figure 5. The NCEI provided~74.66 million single-beam points (black dots) in more than 3800 cruise surveys collected since 1960Sand~6.6 GB of multibeam grids (blue dots). Approximately 120 GB of multibeam grids were downloaded from the JAMSTEC (yellow dots), with most of the data distributed around Japan and the northwestern Pacific. Multibeam grids,~54 GB, around Australia were obtained from GA (purple dots). These ship depths distribute unevenly over the world, and most of them are in the northern hemisphere. All data sources are listed in Table 2.
The ship soundings show an uneven distribution of data quality all over the world oceans, mostly due to the use of analog echo sounders and poor positing before the availability of satellite navigation [14,30]. Smith (1993) assessed the accuracy of ship soundings collected between 1955 and 1992 in Lamont-Doherty Earth Observatory and found the least accurate data were in the southern oceans where the median of the crossover errors at intersecting ship tracks were 100-250 m [30]. The modern multibeam grids from the NCEI, JAMSTEC, and GA have a relative accuracy of about 1% [1]. The single-beam depths were cleaned and edited by scientists at the NCEI before being provided to the public, but there are still a few significant errors. However, we focused on constructing a global seafloor topography model, rather than editing and correcting ship soundings. Thus, we cleaned the ship soundings and removed ship data with obvious errors by comparing them with the SIO topo_20.1.nc model, before constructing the new seafloor topography model. We calculated the standard deviation (STD) of differences between the ship soundings and the SIO topo_20.1.nc model. In each 2 • × 2 • segment, ship depths with ship-model differences larger than twice of the STD were deleted. In the global oceans (−180 •~1 80 • E, −75 •~7 0 • N), 5% of the single-beam points were deleted.The multibeam grids were re-sampled to cover 160 million 15 arc-second grids using the GMT blockmedian [31]. Approximately 90% of the cleaned ship points were applied to construct the new model, and 10% were used to assess the model accuracy.
2020. This version includes an additional year of AltiKa, CryoSat, and Sentinel-3A/B data than the previous version.
The ship soundings, including single-beam points and multibeam grids, were downloaded from NCEI, JAMSTEC, and GA, as show in Figure 5. The NCEI provided ~74.66 million single-beam points (black dots) in more than 3800 cruise surveys collected since 1960Sand ~6.6GB of multibeam grids (blue dots). Approximately 120GB of multibeam grids were downloaded from the JAMSTEC (yellow dots), with most of the data distributed around Japan and the northwestern Pacific. Multibeam grids, ~54GB, around Australia were obtained from GA (purple dots). These ship depths distribute unevenly over the world, and most of them are in the northern hemisphere. All data sources are listed in Table 2.

Data Processing Procedure
The following processing steps were applied to manipulate the data and construct a new 1 × 1 global seafloor topography model, as shown in Figure 6.

1.
The SIO topo_20.1.nc model was filtered to construct a reference model at wavelengths longer than 160 km, h long (x). Then, the reference depths at the ship points, h re f _ship (x ), were interpolated from h long (x).

2.
At ship points, the residual depths, h resi_ship (x ), can be calculated by subtracting h re f _ship (x ) from the observed depths, h ship (x ).
3. The SIO curv_30.1.nc model was band-pass filtered, downward continued, and divided by k to construct VGG at 15~160 km wavelength bands, ∆G z_down (x), and then was used to interpolate VGG at the ship points, ∆G z_ship (x ).

4.
The topography-to-VGG ratios at the ship points were calculated by The ratios were then gridded to a 1 × 1 grid, S(x).
The final seafloor topography model becomes  topography model at wavelengths longer than 160 km. The ship soundings were cleaned by comparing with topo_20.1.nc. About 90% of the cleaned ship points were applied to constrain topography to VGG ratios at 15~160 km wavelength bands, and 10% were used to assess the model accuracy.
The VGG model, curv_30.1.nc was used to predict seafloor topography model at 15~160 km wavelength bands.
For example, Figure 7 shows seafloor topography prediction results in the northwestern Pacific (144 •~1 80 • E, 0 •~3 6 • N). The seamounts, such as Marcus-Wake Guyots, were recovered very well. Compared with ship depths, the differences between the predicted model and ship depths have a mean of −0.182 m and a standard deviation of 35.561 m. That means the data processing procedure was properly performed and VGG can be used to construct a seafloor topography model with high accuracy. For example, Figure 7 shows seafloor topography prediction results in the northwestern Pacific (144°~180°E, 0°~36°N). The seamounts, such as Marcus-Wake Guyots, were recovered very well. Compared with ship depths, the differences between the predicted model and ship depths have a mean of −0.182 m and a standard deviation of 35.561 m. That means the data processing procedure was properly performed and VGG can be used to construct a seafloor topography model with high accuracy.  0°~36°N). The seafloor topography at wavelengths longer than 160 km, ℎ ( ) (a), the VGG grid at 15~160 km wavelength bands, ∆ _ ( ) (b), the topography-to-VGG ratios at ship points ' (c), and the predicted 1′×1′ seafloor topography model, ℎ( ) (d).

The New 1′ × 1′ Global Seafloor Topography Model
The new 1′ × 1′ global seafloor topography model, BAT_VGG2021, was predicted from VGG and ship soundings, as shown in Figure 8. The model clearly revealed seafloor features, such as mid-ocean ridges, seafloor plateaus, and seamount chains. Figure 9 shows the global distribution of differences between BAT_VGG2021 and ship soundings. The results indicate that the new predicted model fits ship measurements very well. The standard deviation of model-ship differences is 45.464 m, and ~93% of the differences are within 100 m. Figure 10 shows the differences between BAT_VGG2021 and the SIO topo_20.1.nc model. The standard deviation difference of these two models is 80.732 m, ~84% of the differences are within 100 m, and ~95.8% of the differences are within 200 m. The seafloor topography at wavelengths longer than 160 km, h long (x) (a), the VGG grid at 15~160 km wavelength bands, ∆G z_down (x) (b), the topography-to-VGG ratios at ship points s(x ) (c), and the predicted 1 × 1 seafloor topography model, h(x) (d).

The New 1 × 1 Global Seafloor Topography Model
The new 1 × 1 global seafloor topography model, BAT_VGG2021, was predicted from VGG and ship soundings, as shown in Figure 8. The model clearly revealed seafloor features, such as mid-ocean ridges, seafloor plateaus, and seamount chains. Figure 9 shows the global distribution of differences between BAT_VGG2021 and ship soundings. The results indicate that the new predicted model fits ship measurements very well. The standard deviation of model-ship differences is 45.464 m, and~93% of the differences Remote Sens. 2021, 13, 3515 9 of 17 are within 100 m. Figure 10 shows the differences between BAT_VGG2021 and the SIO topo_20.1.nc model. The standard deviation difference of these two models is 80.732 m, 84% of the differences are within 100 m, and~95.8% of the differences are within 200 m. The frequency distribution histograms of the differences between BAT_VGG2021 and ship soundings, and the differences between BAT_VGG2021 and the SIO topo_20.1.nc model, are shown in Figure 11. These results indicate that VGG can be used to constrain seafloor topography at 15~160 km wavelength bands, and the data processing method proposed in this paper is correct.
Remote Sens. 2021, 13, x FOR PEER REVIEW 9 of 16 The frequency distribution histograms of the differences between BAT_VGG2021 and ship soundings, and the differences between BAT_VGG2021 and the SIO topo_20.1.nc model, are shown in Figure 11. These results indicate that VGG can be used to constrain seafloor topography at 15~160 km wavelength bands, and the data processing method proposed in this paper is correct.   The frequency distribution histograms of the differences between BAT_VGG2021 and ship soundings, and the differences between BAT_VGG2021 and the SIO topo_20.1.nc model, are shown in Figure 11. These results indicate that VGG can be used to constrain seafloor topography at 15~160 km wavelength bands, and the data processing method proposed in this paper is correct.   (a) (b) Figure 11. The frequency distribution histogram of the differences between BAT_VGG2021 and ship soundings (a), and differences between BAT_VGG2021 and SIO topo_20.1.nc model (b).

Accuracy Evaluated by Comparing with Ship Soundings and Existing Models
We can assess the accuracy of the BAT_VGG2021 model by comparing it with ship soundings and existing models, such as SIO topo_20.1.nc, DTU18BAT, ETOPO1, GEBCO_08, and our last model BAT_VGG2014. The topo_20.1.nc model is the latest version of the global seafloor topography model released by SIO and is being predicted from ship soundings and satellite altimetric gravity anomalies. DTU18BAT is the latest version of the seafloor topography model released by DTU and is being constructed by a combination of the GEBCO model and satellite altimetric gravity anomalies. ETOPO1 is a global topography model released by National Geophysical Data Center, and an old version of the SIO seafloor topography model was used to fill depths in world oceans [24]. GEBCO_08 model is gridded from digital contours. BAT_VGG2014 is a model predicted from single-beam soundings and the satellite altimetric VGG model [22].
We compared existing models with ship soundings in the north Pacific (120°~280°E, 0°~70°N), south Pacific (120°~300°E, −75°~0°N), north Atlantic (280°~360°E, 0°~70°N), south Atlantic (−60°~20°E, −75°~0°N), and Indian Ocean (20°~120°E, −75°~26°N), respectively. Table 3 summarizesthe differences between existing global seafloor topography models and ship soundings in these five regions. The results show that the STDs of differences between the new BAT_VGG2021 model and ship soundings are about 40~80 m, with 39.6 m in the north Pacific and 77.6 m in the south Atlantic. The northern hemisphere shows better accuracy than in the southern hemisphere, mainly due to there being more ship depths covered in the north. The results suggest that the BAT_VGG2021 model has accuracy similar to the SIO topo_20.1.nc model and significantly better than the ETOPO1 (a) (b) Figure 11. The frequency distribution histogram of the differences between BAT_VGG2021 and ship soundings (a), and differences between BAT_VGG2021 and SIO topo_20.1.nc model (b).

Accuracy Evaluated by Comparing with Ship Soundings and Existing Models
We can assess the accuracy of the BAT_VGG2021 model by comparing it with ship soundings and existing models, such as SIO topo_20.1.nc, DTU18BAT, ETOPO1, GEBCO_08, and our last model BAT_VGG2014. The topo_20.1.nc model is the latest version of the global seafloor topography model released by SIO and is being predicted from ship soundings and satellite altimetric gravity anomalies. DTU18BAT is the latest version of the seafloor topography model released by DTU and is being constructed by a combination of the GEBCO model and satellite altimetric gravity anomalies. ETOPO1 is a global topography model released by National Geophysical Data Center, and an old version of the SIO seafloor topography model was used to fill depths in world oceans [24]. GEBCO_08 model is gridded from digital contours. BAT_VGG2014 is a model predicted from single-beam soundings and the satellite altimetric VGG model [22].
We compared existing models with ship soundings in the north Pacific (120°~280°E, 0°~70°N), south Pacific (120°~300°E, −75°~0°N), north Atlantic (280°~360°E, 0°~70°N), south Atlantic (−60°~20°E, −75°~0°N), and Indian Ocean (20°~120°E, −75°~26°N), respectively. Table 3 summarizesthe differences between existing global seafloor topography models and ship soundings in these five regions. The results show that the STDs of differences between the new BAT_VGG2021 model and ship soundings are about 40~80 m, with 39.6 m in the north Pacific and 77.6 m in the south Atlantic. The northern hemisphere shows better accuracy than in the southern hemisphere, mainly due to there being more ship depths covered in the north. The results suggest that the BAT_VGG2021 model has accuracy similar to the SIO topo_20.1.nc model and significantly better than the ETOPO1 Figure 11. The frequency distribution histogram of the differences between BAT_VGG2021 and ship soundings (a), and differences between BAT_VGG2021 and SIO topo_20.1.nc model (b).

Accuracy Evaluated by Comparing with Ship Soundings and Existing Models
We can assess the accuracy of the BAT_VGG2021 model by comparing it with ship soundings and existing models, such as SIO topo_20.1.nc, DTU18BAT, ETOPO1, GEBCO_08, and our last model BAT_VGG2014. The topo_20.1.nc model is the latest version of the global seafloor topography model released by SIO and is being predicted from ship soundings and satellite altimetric gravity anomalies. DTU18BAT is the latest version of the seafloor topography model released by DTU and is being constructed by a combination of the GEBCO model and satellite altimetric gravity anomalies. ETOPO1 is a global topography model released by National Geophysical Data Center, and an old version of the SIO seafloor topography model was used to fill depths in world oceans [24]. GEBCO_08 model is gridded from digital contours. BAT_VGG2014 is a model predicted from single-beam soundings and the satellite altimetric VGG model [22].
We compared existing models with ship soundings in the north Pacific (120  Table 3 summarizesthe differences between existing global seafloor topography models and ship soundings in these five regions. The results show that the STDs of differences between the new BAT_VGG2021 model and ship soundings are about 40~80 m, with 39.6 m in the north Pacific and 77.6 m in the south Atlantic. The northern hemisphere shows better accuracy than in the southern hemisphere, mainly due to there being more ship depths covered in the north. The results suggest that the BAT_VGG2021 model has accuracy similar to the SIO topo_20.1.nc model and significantly better than the ETOPO1 and GEBCO_08 models. The accuracy has been improved obviously from our last BAT_VGG2014 model, which has STD differences of about 120~160 m, due to more altimetric data and multibeam depths being used. Figure 12 shows the frequency distribution of the differences between BAT_VGG2021 and ship soundings in the five regions. It indicates that more than 93% of the differences are within 100 m, except for 87.74% in the south Atlantic. Table 3. The statistics of differences between global seafloor topography models and ship soundings.

Region Model Minimums (m) Maximums (m) Mean (m) STD (m)
North Pacific (120 and GEBCO_08 models. The accuracy has been improved obviously from our last BAT_VGG2014 model, which has STD differences of about 120~160 m, due to more altimetric data and multibeam depths being used. Figure 12 shows the frequency distribution of the differences between BAT_VGG2021 and ship soundings in the five regions. It indicates that more than 93% of the differences are within 100 m, except for ~87.74% in the south Atlantic.

Model Evaluated by Independent Multibeam Grids of MH370 Searching Area
The satellite altimetric data play a pivotal role in filling the blanks between ship tracks. We can discuss the important role of satellite data by comparing with ship-alone grids and independent multibeam depths. For example, in the northeastern Indian Ocean (76 •~1 20 • E, −44 •~− 6 • N), the ship soundings cover is sparse except for west coast of Australia, as showed by the black lines in Figure 13. The ship soundings shown by the black lines were used to construct the seafloor topography model with satellite data, noted as BAT_Pre. The colored depths were multibeam grids of the MH370 searching area, and were used to assess the predicted model, BAT_Pre, as independent data.
Using a combination of ship soundings (black lines in Figure 13) and the satellite altimetric VGG model, a 1 × 1 seafloor topography model (BAT_Pre, as shown in Figure 14) in the northeastern Indian Ocean was constructed using the data processing procedure proposed in Section 3.2 of this paper. For comparison, a ship-alone grid (BAT_Ship, as shown in Figure 15) was also constructed using GMT tools [31]. Comparing BAT_Pre ( Figure 14) to BAT_Ship (Figure 15), the VGG and ship predicted model reveals more seafloor details. To the east of the Ninetyeast Ridge, the BAT_Pre model reveals more seamounts around the Raiit Rise, and other ridges and troughs south of the Wharton Basin. Along the southeast Indian Ridge, the BAT_Pre model recovers many troughs relating to transform faults. tracks. We can discuss the important role of satellite data by comparing with ship-alone grids and independent multibeam depths. For example, in the northeastern Indian Ocean (76°~120°E, −44°~−6°N), the ship soundings cover is sparse except for west coast of Australia, as showed by the black lines in Figure 13. The ship soundings shown by the black lines were used to construct the seafloor topography model with satellite data, noted as BAT_Pre. The colored depths were multibeam grids of the MH370 searching area, and were used to assess the predicted model, BAT_Pre, as independent data. Figure 13. Distribution of ship soundings in the northeastern Indian Ocean. The black lines indicate ship soundings used to construct seafloor model with satellite data, noted as BAT_Pre. The colored depths were multibeam grids of MH370 searching area, and were used to assess the predicted model as independent data. Using a combination of ship soundings (black lines in Figure 13) and the satellite altimetric VGG model, a 1′ × 1′ seafloor topography model (BAT_Pre, as shown in Figure 14) in the northeastern Indian Ocean was constructed using the data processing procedure proposed in Section 3.2 of this paper. For comparison, a ship-alone grid (BAT_Ship, as shown in Figure 15) was also constructed using GMT tools [31]. Comparing BAT_Pre ( Figure 14) to BAT_Ship (Figure 15), the VGG and ship predicted model reveals more seafloor details. To the east of the Ninetyeast Ridge, the BAT_Pre model reveals more seamounts around the Raiit Rise, and other ridges and troughs south of the Wharton Basin. Along the southeast Indian Ridge, the BAT_Pre model recovers many troughs relating to transform faults.     The multibeam grids of the MH370 searching area ( Figure 13) were not used in the BAT_Pre model ( Figure 14) and BAT_Ship model ( Figure 15). Thus, we used these grids to assess model accuracy quantitatively. Table 4 gives the differences between the BAT_Pre model, BAT_Ship model, and multibeam grids in the MH370 searching area. The results show that the STD difference between the BAT_Pre and multibeam grids is 106.2 m, better than BAT_Ship which is 148.5 m. Figure 16 shows the differences between the BAT_Pre model and multibeam grids. Figure 17 shows the differences between the BAT_Ship model and multibeam grids. Comparing these two figures, we found that the BAT_Pre model is better, especially around the Diamantina Fracture Zone and East Indiaman Ridge where seafloor topography changes rapidly. The multibeam grids of the MH370 searching area ( Figure 13) were not used in the BAT_Pre model ( Figure 14) and BAT_Ship model ( Figure 15). Thus, we used these grids to assess model accuracy quantitatively. Table 4 gives the differences between the BAT_Pre model, BAT_Ship model, and multibeam grids in the MH370 searching area. The results show that the STD difference between the BAT_Pre and multibeam grids is 106.2 m, better than BAT_Ship which is 148.5 m. Figure 16 shows the differences between the BAT_Pre model and multibeam grids. Figure 17 shows the differences between the BAT_Ship model and multibeam grids. Comparing these two figures, we found that the BAT_Pre model is better, especially around the Diamantina Fracture Zone and East Indiaman Ridge where seafloor topography changes rapidly.

Conclusions
We constructed a new 1 × 1 global seafloor topography model, BAT_VGG2021, from ship soundings and satellite altimetric VGG model (SIO curv_30.1.nc). The SIO topo_20.1.nc model was used to constrain seafloor topography at wavelengths longer than 160 km. The VGG model was used to constrain seafloor topography at 15~160 km wavelength bands. The ship soundings were used to calibrate the topography-to-VGG ratios at 15~160 km wavelength bands. The accuracy of the new model was assessed by comparing it with ship soundings and existing models. The results show that the STD differences between BAT_VGG2021 and ship soundings are about 40~80 m. The north hemisphere (north Pacific, north Atlantic) shows higher accuracy than the south hemisphere (south Pacific, south Atlantic), mostly due to more coverage of ship soundings. The new model has an accuracy similar to the SIO topo_20.1.nc model, better than the ETOPO1, DTU18BAT, and GEBCO_08 models, and improved significantly on our last BAT_VGG2014 model, especially in the northern hemisphere. These improvements may be due to the application of the new VGG model and the collection of more ship soundings. The distinct role of the satellite VGG model in seafloor topography modeling was discussed by comparison with a ship-alone grid and independent multibeam grids. The predicted model can reveal more details of seafloor features.

Data Availability Statement:
The model produced in this article may be provided by the author (M. Hu, mzhhu@whu.edu.cn), but it is prohibited to be used for commercial purposes.