Recovering Marine Gravity Over the Gulf of Guinea From Multi-Satellite Sea Surface Heights

A regional gravity field product, comprising vertical deflections and gravity anomalies, of the Gulf of Guinea (15°W to 5°E, 4°S to 4°N) has been developed from sea surface heights (SSH) of five altimetry missions. Though the remove-restore technique was adopted, the deflections of the vertical were computed directly from the SSH without the influence of a global geopotential model. The north-component of vertical deflections was more accurate than the east-component by almost three times. Analysis of results showed each satellite can contribute almost equally in resolving the north-component. This is attributable to the nearly northern inclinations of the various satellites. However, Cryosat-2, Jason-1/GM, and SARAL/AltiKa contributed the most in resolving the east-component. We attribute this to the superior spatial resolution of Cryosat-2, the lower inclination of Jason-1/GM, and the high range accuracy of the Ka-band of SARAL/AltiKa. Weights of 0.687 and 0.313 were, respectively, assigned to the north and east components in order to minimize their non-uniform accuracy effect on the resultant gravity anomaly model. Histogram of computed gravity anomalies compared well with those from renowned models: DTU13, SIOv28, and EGM2008. It averagely deviates from the reference models by −0.33 mGal. Further assessment was done by comparing it with a quadratically adjusted shipborne free-air gravity anomalies. After some data cleaning, observations in shallow waters, as well as some ship tracks were still unreliable. By excluding the observations in shallow waters, the derived gravity field model compares well in ocean depths deeper than 2,000 m.

Marine gravity anomaly is an important geophysical quantity with applications in marine resources exploration and exploitation (Becker et al., 2009), geoid modeling (Olgiati et al., 1995;Hwang, 1998;Soltanpour et al., 2007), delimitation of continent-ocean boundaries (Sandwell et al., 2013), revelation of submarine tectonic structures (Hwang and Chang, 2014;Sandwell et al., 2014;Wan et al., 2020a), and bathymetry augmentation (Hwang, 1999). Marine gravity anomaly can be inverted from geoid heights after SSH is reduced to geoid by the removal of dynamic sea surface topography (Olgiati et al., 1995;Andersen et al., 2010;Nguyen et al., 2020). It can also be derived from the deflection of the vertical, whereby the geoid heights are differenced along satellite orbits to obtain geoid gradients. The along-track geoid gradients are then converted to vertical deflections. This along-track differentiation has the advantage of reducing the impact of long wavelength errors (Hwang et al., 2002;Zhu et al., 2019). Also, the vertical deflections approach does not require crossover adjustment (Hwang et al., 2002).
Since marine surface phenomena mostly have submarine origins, with these submarine activities best revealed through gravity anomalies, it is imperative to derive accurate marine gravity models to aid in understanding marine phenomena. Since the past decade, there have been an increase in the quantity and quality of altimetry data, as well as superior computing hardware. This has led to various researchers to develop gravity field models for their respective countries and economic blocs. For instance: over the South China Sea, Zhang et al. (2017) and Zhu et al. (2020) separately developed 1′ × 1′ marine gravity models by amalgamating multi-satellite SSH data; Nguyen et al. (2020) merged gravity anomaly derived from Cryosat-2 and Saral/AltiKa SSH datasets with shipborne gravity data over the Gulf of Tonkin area of Vietnam. A recent study by Sandwell et al. (2013) has shown that gravity field models derived from modern altimeter satellites are more accurate than those obtained from some governmental shipborne gravity data. This shows the usefulness of satellite-derived gravity models in filling the gaps between tracks of shipborne gravimetry. It must be appreciated that though gravimetry via ships is still reliable, it is expensive, slow and mostly limited to the maritime waters of developed countries. Therefore, SSH datasets from altimetry satellites are the best and cheapest option for developing regions like West Africa.
It is expected that the development of the study region synchronizes with scientific investigations of its largest natural resource (the ocean); however, this is yet to be realised. Even though the region is oil-rich (Osaretin, 2011;Bazilian et al., 2013), the scantiness of reviewed literatures about the maritime waters of West Africa shows that the region's maritime waters is one of the least studied. The study area encompasses 15°W to 5°E, and 4°S to 4°N. Countries of interest include Guinea and Guinea Bissau at the western end, through to the central-south African country of The Congo.
The present study aims to study the gravity field of the region by using vertical deflections derived from the geodetic missions (GM) of Jason-1, Saral/AltiKa, Haiyang-2A; Cryosat-2, and Envisat. We achieve this objective through the remove-restore procedure using EGM2008 as the reference gravity filed model. The use of EGM2008 is due to the fact that it is the most frequently used global geopotential model for marine geodetic and geophysical research (Andersen et al., 2010;Sandwell et al., 2014;Zhang et al., 2017;Andersen and Knudsen, 2019;Zhu et al., 2019;Zhu et al., 2020).

Altimetry Data
The altimetry datasets used in this research comprises GDR (geophysical data record) of eight cycles of Jason-1/GM, seven cycles of Saral/AltiKa, five cycles of Envisat, five cycles of Cryosat-2, and IGDR (interim geophysical data record) of seven cycles of HY-2A/GM. The Jason-1 GDRs were obtained from the Physical Oceanography Distributed Active Archive Center, Jet Propulsion Laboratory, NASA (https://podaac.jpl.nasa.gov/JASON1? sections data). Saral/AltiKa GDRs were obtained from AVISO (https://www.aviso.altimetry.fr/), Envisat and Cryosat-2 GDRs were provided by the European Space Agency (http://earth.esa. int); while HY-2A GDRs were provided by the National Satellite Ocean Administration Service of China (ftp2.nsoas.org.cn). Along-track SSH is computed from each cycle as given by Eq. 1. Table 1 is a summary of the various selection criteria applied to each satellite's dataset during SSH computation. Figure 1 shows the spatial distribution of ground tracks of the different satellites.
where corrections wet troposphere correction + dry troposphere correction+ ionosphere correction + sea state bias correction+ pole tide correction + ocean tide correction+ solid earth tide correction + inverse barometric correction

Reference Gravity Field
In implementing the remove-restore procedure, a global geopotential model (otherwise known as reference gravity field) is required. The reference gravity field of one quantity (which in this case is deflections of the vertical) is removed from its corresponding observed deflections of the vertical to obtain residual deflections of the vertical. The residual deflections of the vertical are then used to construct the residual gravity anomalies. Full gravity anomalies are later obtained by restoring the removed reference gravity field (now in the form of gravity anomalies). The removal of the reference gravity field ensures a more statistically homogeneous and smoother model. Also, it has the effect of ensuring that the gravity field outside the data extent is accounted for (Andersen, 2013). EGM2008 is expressed as coefficients of spherical harmonics and has a maximum degree of 2,190 (Pavlis et al., 2012). It was obtained from the International Centre for Global Earth Models (http://icgem.gfz-potsdam.de/).

Dynamic Topography
In the computation of gravity anomalies, the SSH data must initially be reduced to geoid heights. To achieve this, the timedependent topography of the surface of the sea must be removed. It is well known that the geoid is linearly related to the mean sea surface through the mean dynamic topography (MDT) of the ocean surface. Therefore, to compute the geoid heights, this research used the mean dynamic topography DTU15MDT, obtained from the Technical University of Denmark. It is a product derived from EIGEN-6C4 and DTU15MSS mean sea surface model . The time-independent sea surface topography was removed through filtering similar to the method of Hwang et al. (2002).

Computation of Gravity Anomalies
The SSH i at a point, i, computed from Eq. 1 is transformed into geoid height, N i , as given by Eq. 2: Geoid gradients, β ij , between two consecutive points i and j along the satellite's ground track are then calculated by: where D ij is the spherical distance between the two points i and j. For each geoid gradient, β ij , its corresponding azimuth, α ij , is also computed. To obtain the vertical deflections, the alongtrack geoid gradients and azimuths are gridded using biharmonic spline with a tension of 0.25. The biharmonic spline in tension was executed using the surface module of GMT (Generic Mapping Tools). The north-component, ξ, and  (Hwang et al., 2002), as given by Eq. 4: where v k is the model error.
With reference to the remove-restore technique described in Figure 2, the computed components of deflections of the vertical are reduced to their residual version by subtracting EGM2008simulated deflection components from them. The reference vertical deflections were simulated at the maximum degree/ order of 2,190.
where Δξ and Δη are, respectively, north-component and eastcomponents of residual vertical deflections; ξ EGM2008 and η EGM2008 are the EGM2008-simulated north-and eastcomponents of vertical deflections, respectively. The residual vertical deflection components are then used to invert residual gravity anomaly, δg, through Fourier transform as expressed in Eq. 6: where F −1 is the inverse Fourier transform, c is mean value of gravity, u and v are the spatial frequencies. ΔX and ΔE are, respectively, the Fourier transforms of Δξ and Δη. Finally, using the remove-restore procedure, long wavelength gravity anomaly, Δg EGM2008 , from the reference gravity field at degree/order 2,190 is restored to the residual gravity FIGURE 2 | Remove-restore technique used.
Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 700873 anomaly, δg, to recover the gravity anomaly, Δg. This is simply given as: Reviewed literatures have shown that for most altimetryderived deflections of the vertical, the north component is often times more accurate than the east component. This is attributable to the inclinations of the satellites which usually are towards the north. The effect of this accuracy imbalance propagates to the resultant gravity anomalies derived from the deflections of the vertical (Small and Sandwell, 1992;Zhang, 2017;Wan et al., 2020b). In order to nullify this effect on the derived gravity anomaly model, weights are assigned to the vertical deflection components based on the covariance between them, and their error variances relative to the reference gravity field model. The residual gravity anomaly is derived from each component, δg ξ and δg η as (Small and Sandwell, 1992;Hwang and Parsons, 1996): Using the constraints described in Eq. 9, the weights, w ξ and w η are computed as given by Eq. 10 (Wan et al., 2020b): where δ 2 Δξ and δ 2 Δη are the corresponding error variances of Δξ and Δη, respectively; and δ 2 ΔξΔη is the covariance between Δξ and Δη.

RESULTS
The computed vertical deflection components are shown as Figure 3. They were validated by comparing them with simulated components from EGM2008 at maximum degree/order 2,190. The result of this comparison is summarized in Table 2. With regard to standard deviation of differences, it can be seen from Table 2 that the accuracy of the north component is slightly better than the east component by almost three times. A visual inspection shows that, in general, the positive deflections of the north component ( Figure 3A) have a SW/NE orientation, whereas there is no distinct orientation in the east component ( Figure 3B). With reference to Eqs. 8-10, Table 3 is a compilation of the vertical deflection components' error variances, covariance and their respective weights used for the construction of the residual gravity anomaly. The constructed gravity anomaly model is shown in Figure 4. It compares well with gravity anomalies from EGM2008, DTU13 (obtained from the Technical University of Denmark) and SIOv28 (obtained from Scripps Institution of Oceanography) used as reference models. A summary of this comparison is presented in Table 4. Figure 5 shows similar histograms of gravity anomalies for the inverted model and the reference models; this is an attestation of high accuracy of the derived gravity anomalies. Figure 6 is a visual analysis of differences between the inverted gravity anomaly model and the reference models; this is a further attestation of the results presented in Table 4. On average, the model deviates from the reference models by −0.33 mGal. Further analysis of the histogram indicated that more than approximately 50% of gravity anomalies in the study region falls within the range of −80-80 mGal.

Comparing Derived Gravity Anomaly Model With Shipborne Gravimetry
The accuracy of the constructed gravity anomaly model is further analyzed by comparing it with gravity anomalies observed along tracks of ship cruises obtained from National Centers for Environmental Information of the National Oceanic and Atmospheric Administration (https://maps.ngdc.noaa.gov/ viewers/geophysics/). However, the shipborne gravity data has the disadvantage of being very heterogenous with some erroneous observations made at some points. Some sources of these errors are: the absence of a uniform reference field at the time of measurement, drifting of the gravimeters, absence of base stations for references, etc. (Wessel and Watts, 1988). Therefore, in degree to amalgamate the shipborne gravity anomalies to a common reference, gravity anomalies simulated from EGM2008 were used as reference. With the exception of track RC1312, a second-degree polynomial model is then used to adjust the shipborne gravity anomalies. Track RC1312 was well adjusted by a first-degree polynomial. The adjustment model, given as Eq. 11, is applied separately on each ship track of observations.
where Δg is the difference between shipborne, Δg ship , and EGM2008-simulated, Δg EGM2008 , gravity anomalies; Δt is the difference between observation time, Δt observe , and time of ship's departure, Δt depart . The unknown parameters to be solved are θ 0 , θ 1 , and θ 2 . The adjusted shipborne gravity anomalies are presented in Figure 7. From Figure 7, it can be seen that the northern half of the study region generally has lower gravity anomalies. Higher values can be observed at the southern part, especially in the south-west. The middle section has anomalies ranging approximately from −60 to 45 mGal. These observations generally agree with observations made from the constructed gravity model presented in Figure 4.
For each ship track, statistical summary analysis of the differences between the shipborne anomalies and modelled anomalies is computed; the results are given in Table 5. A total of 288 points were discarded due to their difference Without considering the mean errors, which denote the systematic errors, for major tracks of shipborne gravimetry, the standard deviations of the differences between the shipborne observations and the modelled gravity anomalies are smaller than 5.5 mGal. The ratios of the tracks with standard deviations smaller than 5.5 mGal and 6.0 mGal are about 71.4 and 78.6%, respectively. This index is close to the results of Table 4. Furthermore, from Table 5, it must be appreciated that the adjustment was done using EGM2008 as true values; therefore, the standard deviations of differences between the adjusted shipborne and EGM2008 gravity anomalies are virtually negligible. However, the magnitudes of the mean differences range within 1-20 mGal, which indicates the contribution from altimetry observations. In order to analyze the systematic errors depicted by the large mean errors in Table 5, the accuracy of gravity anomalies is assessed with respect to varying ocean depth. Table 6 is a summary of this analysis. The gravity anomalies are more accurate at depths deeper than 2,000 m. The mean deviation is approximately 5.09 mGal, which signifies the existence of some systematic errors. Therefore, in order to nullify the effects of these errors, observations whose deviation exceeded the mean deviation of the entire dataset were removed. This resulted in removal of values in shallow depths (depth range of 1,001-2,000 m), with the new results compiled in Table 7. It is obvious from Table 7 that the mean and standard deviation of differences escalate in shallow waters (i.e., in depth range of 1,001 ∼ 2,000 m). This shows that observations in shallow waters (depths > 2,000 m) were the main source of systematic differences. Analyzing each satellite's contribution in constructing the gravity anomaly model In order to investigate the contribution of each satellite in the gravity field construction, the north, ξ i , and east, η j , components of deflections of the vertical were computed for each satellite. Summary of each satellite's contribution in establishing the integrated model, (ξ, η), is shown in Table 8. This was achieved by uniformly generating two sets of random numbers (each set contains five elements). The elements of each set are normalized so that their sum equates to one. These normalized values are now used as weights, w, applied to the satellite's vertical deflections as given in Eq. 12. The randomization process was repeated 100 times for each of the two sets; and for each iteration, normalization was performed correspondingly. Using the deflections of integrated model as true values, the set of weights that minimizes the standard deviation of errors was then chosen.
w i ξ i + ... + w 5 ξ 5 w i ξ w j η j + ... + w 5 η 5 w j η w i + ... + w 5 w i 1 i 1, 2, . . . , 5 w j + ... + w 5 w j 1 j 1, 2, . . . , 5 It can be seen from Table 8 that SARAL/AltiKa contributed ∼25% in resolving the north component. This can be attributed to the range accuracy of its Ka-band signal. It is also possible that its good drifting phase spatial resolution (Verron et al., 2021) could be a factor, though this is sometimes lower than that of Cryosat-2. Approximately, 19% proportion of north component was realized from each of HY-2A/GM, Cryosat-2, and Envisat. This can be attributed to the geodetic spatial resolutions of these satellites, except for Cryosat-2 and Envisat which are not geodetic. Even though the datasets from Jason-1/GM are also geodetic, it contributed the lowest north component due to its lower spatial resolution. Jason-1/ FIGURE 6 | Histograms of gravity anomaly deviations from reference models. GM's lower contribution in resolving the north component could also be attributed to its inclination of only 66°. On the other hand, due to this low inclination of 66°, Jason-1/GM contributed ∼28% of the east components, which is higher than that of SARAL/AltiKa. It was only surpassed by Cryosat-2 with ∼34%, which has one of the best spatial resolutions amongst altimetry satellite missions. For Cryosat-2, since it resolved more deflections in east component than north component, it seems that the effect of spatial resolution is larger than the inclination.

CONCLUSION
A regional gravity anomaly model of the Gulf of Guinea has been constructed through the remove-restore method. It was constructed from deflections of the vertical derived from five different altimetry missions. The north-component of vertical deflections was more accurate than the east-component by almost three times. Analysis of the vertical deflections (Table 8) revealed that, on average, SSH data from each satellite can contribute almost equally in resolving the north component. However, Cryosat-2, Jason-1/GM, and SARAL/ AltiKa were the dominant contributors of the east component. The gravity anomaly model was constructed by assigning weights of 0.687 and 0.313 to the north and east components, respectively. A comparison of histograms (refer to Figure 5) of the computed gravity anomaly model with renowned models from EGM2008, SIOv28, and DTU13 showed similar results. The computed gravity anomaly averagely deviated from these reference models by −0.33 mGal. Further assessment was done by comparing it with quadratically adjusted shipborne free-air gravity anomalies. By excluding observations in shallow waters, the derived gravity field model compares well in depths deeper than 2,000 m.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
Conceptualization, RA; methodology, both authors; data curation, both authors; funding acquisition, XW; formal analysis, RA; investigation, both authors; writing-original draft preparation, RA; writing-review and editing, both authors.