Seafloor Density Contrast Derived From Gravity and Shipborne Depth Observations: A Case Study in a Local Area of Atlantic Ocean

Marine gravity data from altimetry satellites are often used to derive bathymetry; however, the seafloor density contrast must be known. Therefore, if the ocean water depths are known, the density contrast can be derived. This study experimented the total least squares algorithm to derive seafloor density contrast using satellite derived gravity and shipborne depth observations. Numerical tests are conducted in a local area of the Atlantic Ocean, i.e., 34°∼32°W, 3.5°∼4.5°N, and the derived results are compared with CRUST1.0 values. The results show that large differences exist if the gravity and shipborne depth data are used directly, with mean difference exceeding 0.4 g/cm3. However, with a band-pass filtering applied to the gravity and shipborne depths to ensure a high correlation between the two data sets, the differences between the derived results and those of CRUST1.0 are reduced largely and the mean difference is smaller than 0.12 g/cm3. Since the spatial resolution of CRUST1.0 is not high and in many ocean areas the shipborne depths and gravity anomalies are much denser, the method of this study can be an alternative method for providing seafloor density variation information.


INTRODUCTION
Shallow seabed density information plays an important role in understanding submarine structure, mineral exploration, military activity, and scientific research (Nagendra et al., 1996;Felix et al., 2002;Hsiao et al., 2011;Sandwell et al., 2014a). For a local ocean area, the shipborne equipment can be used for the detection. However, for the global ocean area, it is very difficult to detect and retrieve the shallow seabed density distribution with high resolution and high precision using shipborne observations.
Currently, there are mainly two methods for the study of the Earth's internal structure, including the density contrast of the seabed, i.e., the density difference between the seawater and the upper crust. One is based on seismic wave data and the other is based on gravity data. For the former, scholars have carried out a lot of investigations. For example, Soller et al. (1982) derived a global crustal thickness model with resolution of 2 • × 2 • based on seismic data; Nataf and Ricard (1996) combined geothermal and seismic data to invert the density structure of the global and upper mantle; Mooney et al. (1998) published a global crustal model crush with a spatial resolution of 5 • × 5 • . Based on seismic, ice sheet and sedimentary datasets, Bassin (2000) derived the crustal model CRUST2.0 with a resolution of 2 • × 2 • and Laske et al. (2013) released the crustal model CRUST1.0 with a spatial resolution of 1 • × 1 • based on CRUST5.1 and CRUST2.0. The above models provide not only the thickness information of each layer, including water layer, ice layer, soft sedimentary layer, hard sedimentary layer, upper crust, middle and lower crust, but also the density information. These models play a significant role in the study of the Earth's internal structure.
However, the above models cannot fully meet all the application requirements. For example, the spatial resolution of the global model published only reaches 100 km, which is obviously insufficient for studying many local problems. In addition, the models mentioned above have great uncertainties in many regions, mainly because of the uneven distribution of seismic observation stations, and thus there are no or few seismic observations in some regions. Therefore, if one wants to retrieve the density distribution of the Earth's interior with higher spatial resolution using information from seismic data, a densification of current network of seismic stations is needed. It would involve a lot of financial and material resources.
Another method is to use gravity field data for density inversion. Many algorithms for inversion of the Earth's internal structure or density distribution using gravity data have been proposed and adopted. Parker (1972) gave the Fourier relationship between the depth of density interface and gravity anomaly in frequency domain. According to the relationship, the depth of density interface can be deduced from gravity anomaly. The related formula has been widely used in inversion of Moho surface depth and seafloor topography, such as Jiang et al. (2015); Hu et al. (2015), and Sandwell et al. (2014a), etc. Indeed, the formula is based on the premise that there is a density anomaly and the density contrast is known. Conversely, if the depth of the density interface is known, the formula can also be used to derive the density contrast. In recent years, with the emergence of new altimetry technologies and the release of new altimetry satellite products, Sandwell et al. (2014a) studied the submarine structure, especially using the gravity gradient data to find and confirm the position of the plate boundary on the seabed, which can also be regarded as a density interface.
Definitely, gravity data can be used to detect the density distribution of the Earth's interior. However, some problems exist in the study of the Earth's interior structure (including the density contrast of the seabed) using gravity data. For example, the inversion results are usually not unique. The fundamental reason is that many sources of gravity anomaly information exist, including not only the seabed density contrast, but also the depth of the density interface (that is, the ocean depth). Therefore, if we can combine the ocean gravity field and the ocean depth information well, it is possible to retrieve the seabed density contrast information with high accuracy.
Indeed, in the inversion of seafloor topography, many methods have been proposed for density contrast determination, such as the iterative search method Ouyang et al., 2014;Xiang et al., 2017). However, the density contrast derived in the bathymetry inversion usually have no physical meaning (Annan and Wan, 2020), since the purpose of these researches is to inverse bathymetry. The idea of this study is to take the density contrasts of the seabed as the inversion product while the water depths are known.
With the advances in space technology, highly accurate marine gravity anomaly products (Hwang et al., 2002;Andersen et al., 2010;Hwang and Chang, 2014;Sandwell et al., 2014b;Zhu et al., 2020) with resolution of 1 × 1 have been derived using satellite altimetry observations. Besides gravity field data, a large amount of shipborne depths have been obtained by many years of ship sounding observations. It is roughly estimated that NOAA (National Oceanic and Atmospheric Administration) has collected more than 50 million shipborne depth observations. Therefore, it is possible to achieve good results by making full use of the shipborne sounding data and gravity data to jointly invert the seabed density contrast. This study presents a case study. Section "Materials and Methods" introduces the inversion algorithm. Study area and data are described in section "Study Area and Data. Results and analysis are presented in section "Results and Analysis" and conclusions are derived finally in section "Conclusion."

Function Model
According to principle of the gravity geological method (GGM) (Ibrahim and Hinze, 1972), we have: whereĝ i denotes the gravity anomaly observation at the i th point; g l denotes the long wavelength part of the gravity anomaly; G means the gravitational constant; ρ represents the density contrast; H is the depth datum andĥ i denotes the water depths at the i th points. Definitely, the gravity anomaly observations have a linear relationship with (H −ĥ). In a local small area, the long wavelength part of gravity anomalyĝ l and density contrast ρ can be seen as constants. Hence, Equation 1 can be rewritten as: Because both gravity anomaly and water depth datasets contain noises, two types of data error should be considered in the determination of the density contrast ρ and long wavelength of gravity anomalyĝ l . Ordinary least squares method only considers the error of dependent variable, i.e., gravity anomaly errors and ignores the error of bathymetric data measured by ships. Hence, the inversion result by traditional least squares method based on Equation 2 would be not accurate enough. Instead, the total least squares (TLS) (Golub and van Loan, 1980;Xu et al., 2012) method was adopted in this study. Considering the two types of observation errors, the following equation can be constructed: In which, v gi and v hi are the noises of gravity anomaly and shipborne depths respectively. Changing Equation 3 into the form of EIV (error in variables) model results in: where A is the coefficient matrix; E A ,E L are the errors of matrices of A and observation vector L; x is the parameter to be estimated, i.e., The final solution can be derived by the following rule.

Solution of TLS
Solution of the TLS can be derived by Singular Value Decomposition (SVD) (Golub and van Loan, 1980). Firstly, augmented matrix C is defined as: And then Orthogonal Trigonometric (QR) decomposition is conducted on matrix C, in which Q is the orthogonal matrix as follows: Based on the singular value decomposition of R-part of matrix C R , the solution of TLS is obtained as Equation 10: where = diag(δ 1 , δ 2 ) with singular values δ 1 ≥ δ 2 ≥ 0. And thus, the estimated parameters can be derived by the following equation.
Finally, we have: The accuracy of TLS solution is evaluated by standard deviation of unit weight (SDUW) as follows:

Filtering Method
In order to test the sensitive band of gravity anomaly to the seafloor topography, which is the basis for deriving the density contrast, band filtering with different pass bands was experimented and the correlation analysis was conducted correspondingly. The band-pass filters proposed by Smith and Sandwell (1994), which is indeed a filter by combining a Gaussian high-pass filter and a Wiener low-pass filter was used. According to Smith and Sandwell (1994), the high-pass filter is defined as Equation 14: where k = 1 λ ,λ denotes the wavelength; s is a parameter of the filter which can be derived by w l (k) = 0.5 at the cutoff wavelength (denoted as λ h ): The low-pass filter is defined as Equation 16: where d is the mean depth of the study area; A is parameter of the filter, which can be derived by w l (k) = 0.5 at the cutoff wavelength (denoted as λ l ) as: Finally, the bandpass filter is obtained as: STUDY AREA AND DATA

Study Area and Observation Data
The study area is a local region in the Atlantic Ocean, i.e., 326 • E-328 • E, 3.5 • N-4.5 • N. The bathymetric data from ship survey were provided by First Institute of Oceanography, Ministry of Nature Resources, and 498,501 shipborne depths are obtained by multibeam echo-sounding with high accuracy. The gravity data used is DTU13 version of gravity anomaly model released by Technical University of Denmark 1 . DTU13 is derived based on DTU10 by adding more satellite observations, such as Cryosat-2 (Andersen et al., 2014). The resolution of DTU13 is 1 × 1 and its precision is better than four mGal (Andersen, 2010;Andersen et al., 2014). The whole study area is divided into eight small areas, denoted as A, B, C, D, E, F, G, and H with at a size of 30 × 30 , as shown in Figure 1, in which water depths are shown by different colors. The three-dimensional distributions of gravity anomaly and shipborne depths are shown in Figure 2. According to this figure, the two maps show very similar characteristics which denote the high correlation between the two types of datasets. . Sea water depth, crustal thickness, and density are derived by the average data in each unit.
Considering that the spatial resolution of CRUST1.0 is 1 • × 1 • and the size of the study area is 1 , only four density contrasts were derived using CRUST1.0 in the study area, for four regions defined as a, b, c, and d corresponding to A and B, C and D, E and F, and G and H shown in Figure 1 respectively. Parameters of CRUST1.0 in each region are presented in Table 1, in which the density contrast ρ 0 is the difference between sea water and the upper crust. It is found that the four density contrasts are almost same except in area CD. According to Table 1, the thickness of the upper crust in region CD is largest compared to other regions. This may cause the differences in the density contrasts with other region, which would be investigated in future research.
It needs to be noted that the effects of sediments are not considered in the results of ρ 0 , but the results derived by the TLS contain the whole signals, including sediments. Hence, it is not accurate to ignore the sediments in the calculation of density contrasts using CRUST1.0 information, especially when the thickness of the sediments is large. In addition, the thickness of the sediments and the upper crust would also lead to the differences with TLS. In order to solve these issues, new density contrasts are derived as follows, in which, ρ sed , ρ c1 are respectively, the densities of the sediment layer and the upper crust; H sed , H c1 are the thicknesses of the sediment layer and the upper crust respectively; ρ w is the density of the sea water. The newly derived results for each region are also given in Table 1.

RESULTS AND ANALYSIS
Initial Results From TLS Figure 3 shows the fitting results of TLS for each region, in which the x-and y-axis denote (H −ĥ i ) and ĝ i respectively. Please note, because the resolution of shipborne depths is much higher than DTU13, the shipborne depths are interpolated to the gird points of DTU13 before the TLS processing in order to make the two datasets have consistent resolution. It is found that the gravity anomalies show linear positive correlations with water depths (i.e.,H −ĥ i ) in each subregion. For example, the larger the water depth, the larger the gravity anomaly. However, the slopes and  intercepts of the fitting lines in each region are not the same. It can be inferred that this is caused by the complexity of the seabed geological structures. The derived density results are shown in Table 2. It is found that the density contrasts of area A-H are in the range of 0.599-1.173 g/cm 3 , and the density contrasts are different due to the unique geological structure in each region. According to the SDUW, the fitting error in area A is the best, and that in area H is the worst. The accuracy of the whole area is between 6.124 and 14.493 mGal, among which the accuracy of area D and G is the closest, which indicates that the submarine geological structures in these two areas are similar. The SDUW of areas A and B have the largest difference in the adjacent areas, which may be caused by the large differences in the complicated seafloor geological structures between the two areas.
Compared to Table 1, the results from TLS are quite different from ρ 0 and ρ 1 . Mean values of the whole study area are given in Table 3. It can be seen from the above table that the density contrasts obtained by CRUST1.0 model are much larger than that obtained by total least squares algorithm. The density contrast obtained by the weighted average method is closer to the result of TLS than that obtained by CRUST1.0 initial results, i.e., ρ 0 . Even so, the difference between the result of weighted average method and that obtained by total least squares is −0.374 g/cm 3 . If the values from CRUST1.0 are true values, the relative error exceed 25%. This indicates that the derived results are not very    reliable. This conclusion is consistent with the fact that the density contrast derived in the GGM for bathymetry usually has no physical meaning (Annan and Wan, 2020). However, it should be noticed that the marine gravity anomalies are correlated with the information of the seafloor topography, including its depth and density contrast. In theory, if the density is known, the depths can be derived; if depths are known, the density can be derived. The possible reason why the large differences exist between the derived results with CRUST1.0 is that gravity anomalies have high correlations with seafloor topography only in a limited wave bandwidth (Marks and Smith, 2012;Wan et al., 2019). Although this issue has been considered in the derivation of bathymetry using GGM by modeling and removing the long wavelength gravity anomaly, the errors of the modeling would certainly cause some errors in the derivation of density contrast. In order to solve this issue, a filtering processing was added both to the gravity anomalies and shipborne depths and the new results are obtained in the next section.

Filtering Processing
The pass bands are designed by referring to studies about bathymetry, such as Smith and Sandwell (1994); Marks and Smith (2012). The sensitive band for bathymetry inversion is usually in the range of 20 ∼ 200 km (Hwang, 1999;Marks and Smith, 2012). This is also true for the density contrast inversion, because the mathematical function relationship is same as that used in bathymetry inversion, i.e., Equation 1. In order to reduce the effect from the non-sensitive bands, we decrease the band width which is usually used for bathymetry inversion and two pass bands are experimented, i.e., (30-100) km and (50-100) km.
As to why 100 km is selected as the maximum cutoff wavelength of the filtering pass-band, it is because the size of the study region is 1 • × 2 • ; and thus, it is difficult to present signals with wavelength larger than 100 km due to the limited size, at least in the latitude direction. It is also because the spatial resolution of CRUST1.0 is 1 • × 1 • , corresponding to approximately 100 km × 100 km in spatial domain. As for the minimum cutoff wavelength of the pass-band, 30 and 50 km are defined as examples. Although, in bathymetry inversion, the minimum cutoff wavelengths of the pass-band are usually lower than 30 km, such as 20 or 15 km (Abulaitijiang et al., 2019), this study enlarges the minimum cut-off wavelength of the passband in order to ensure a higher correlation between the filtered shipborne depths and gravity anomalies. It is difficult to obtain the most sensitive band from gravity anomaly to bathymetry (Wan et al., 2019) like in gravity gradient to bathymetry. It is also true between the gravity anomaly and density contrast, since the mathematical relationship are same in bathymetry and density contrast inversions.
As an example, the filter for pass-band (30-100) km as well as the filtered data are shown in Figures 4, 5 respectively. In Figure 4, the red dotted lines denote the cutoff wavelength, i.e., 30 and 100 km. Compared to Figure 2, the similarity between gravity anomaly and ship-depths shown in Figure 5 is definitely much stronger. This indicates the correlations between the two data sets are higher than that of the data shown in Figure 2. This point is proved further by Table 4, which shows the correlations between the ship-depths and gravity anomalies in each sub-region. Obviously, the correlations are improved largely by the filtering processing. Especially for the region H, the improvements in the correlations arrive at 35% compared to the initial correlations.

New Results
By band-pass filtering gravity anomalies and shipborne depths data, density contrasts are again derived using TLS and shown in Table 5. And the comparisons with weighted mean values of CRUST1.0 are given in Table 6. We named the pass-band (30-100) km as pass-band1 and the pass-band (50-100) km as pass-band2 in this study.  Definitely, the new results are much closer to CRUST1.0 than the initial results and the differences are reduced by more than 60% in the total area. The magnitude of the differences between the new results and CRUST1.0 are smaller than 0.13 g/cm 3 and 0.03 g/cm 3 for the two pass-bands in terms of mean value of the whole area. If the CRSUT1.0's value is the true value, i.e., 1.442 g/cm 3 , the relative errors of the two pass-bands are 8.46% and 1.53%, respectively. In general, in the whole area, the pass-band2, i.e., (50-100) km, yields closer results to CRUST1.0. However, in Region EF, the difference becomes larger. This may   be due to the fact that the spatial resolution of CRUST1.0 is not enough to present the density contrast variation in local regions with size smaller than 1 • × 1 • . This would be investigated further in the future.
Please note, no information of CRUST1.0 is used in the inversion. Since the results are now close to CRUST1.0, it proves that marine gravity anomalies and shipborne depths can be used to derive seafloor density contrasts with a high accuracy. It also proves that the proposed method is effective. It should be emphasized that the spatial resolution of CRUST1.0 is only at 1.0 degree but it is not the case for gravity anomaly and even the shipborne depths. Hence, it is fully possible to derive the density contrast information all over the ocean with a much higher spatial resolution than CRUST1.0, because highly accurate gravity anomalies can be provided by several altimetry missions and a very large amount of shipborne depth observations have been collected by NOAA. The global inversion would be conducted in the future study.

DISCUSSION
In order to show the influence of the bathymetry undulation on density contrast inversion, the mean, deepest, and standard deviation (Std) values of the depths in each subregion are presented in Table 7. Comparing Tables 6, 7, it is found that the deeper the water depths, the lower the inversion accuracy if the region EF is not considered. For example, in the initial shipborne depths, the absolute values of the mean and the deepest depths in GH region are all larger than regions AB and CD. Correspondingly, the inversion accuracy in region GH is much lower than those of regions AB and CD. Even after filtering, the maximum depth in region GH is still larger than those of regions AB and CD; and thus, the accuracy in region GH is still the lowest among the three regions. This phenomenon may be due to the fact: the deeper the water depth, the lower the gravity anomaly signals at the sea surface generated by density body of the seafloor. Since the gravity anomaly accuracies are almost same in the study subregions, the signal to noise ratio (SNR) must be different with different water depths. In general, the SNR of GH is lower than that of other regions. Hence, the density contrast inversion accuracy in the deep ocean area would be poorer than other areas. Because of this, in order to obtain the same inversion accuracy, the gravity anomaly should have a higher accuracy in the deep ocean regions than areas with shallower water depths.
It should be noted that the inversion results of region EF are not consistent with the above phenomenon. This may be caused by the limited spatial resolution of CRUST1.0, which is not enough to present the regional information. In order to check accuracy of CRUST1.0, we compare the sea water depths given by CRUST1.0 with shipborne depths used in this study. The statistics is given in Table 8, in which the differences between shipborne depths and CRUST1.0 are given and the relative error represents the ratio of the difference to shipborne depth value. According to Table 8, mean water depths provided by CRUST1.0 are very close to the shipborne depths and the mean differences are smaller than

CONCLUSION
This study presented a study for seafloor density contrast inversion using gravity anomalies and shipborne depth datasets. The inversion was achieved by a total least squares algorithm, which considers both noises from gravity anomalies and shipborne depths. The results show that if gravity anomaly and shipborne depths are used directly, the derived density contrasts contain large errors, which may be caused by the modeling error of the long wavelength part of gravity anomalies. Hence, a band pass filtering technique was proposed to resolve this issue and the final results show an obvious accuracy improvement of the derived density contrast, and thus verified the effectiveness of the method. As for the sensitive band, in general it is the band of 20 ∼ 200 km. However, in density contrast inversion, it would be better to only use signals in a part of the sensitive band if it is defined as where the correlation is larger than 0.5, because the high correlation between the gravity data and shipborne depths would help improve the inversion accuracy. Therefore, we suggest that the determination of the pass-band of the filter is done by a correlation analysis between the gravity anomaly and shipborne depths, and the sensitive band can be defined as the bands with correlations larger than 0.85 or even higher between the two data sets to ensure high correlation in the pass-band.
Since various altimetry satellite missions (Hwang et al., 2002;Sandwell et al., 2014a;Wan et al., 2019;Zhu et al., 2020) have provided enough observations to derive highly accurate and dense marine gravity products and a large amount of shipborne depths data have been collected by NOAA, high accuracy and resolution of density contrasts can be derived for major areas of the global ocean by the proposed method in this study. In addition, gravity gradient products derived by some institutes (Hwang and Chang, 2014;Sandwell et al., 2014a) may also contribute a lot in the inversion.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: The gravity data can be downloaded from the DTU space center website: https://www.space.dtu.dk/english/ Research/Scientific_data_and_models/Global_Marine_Gravity_ Field. CRUST1.0 is available from the website: https://igppweb. ucsd.edu/~{}gabi/crust1.html. The shipborne depths used in this study are available on the request to the corresponding author.