Moho Geometry Estimation Beneath the Gulf of Guinea From Satellite Gravity Data Based on a Regularized Non‐Linear Gravity Inversion

In this work, we present a new gravimetric Moho model over the Gulf of Guinea computed from the inversion of the global gravity EIGEN‐6C4 model. The computation procedure contains a regularized non‐linear inversion algorithm. Before applying the inversion algorithm, the gravity data from the EIGEN‐6C4 model have been corrected from gravimetric effects of topography, bathymetry, and sediments. To calculate these corrections, a forward modeling has been applied using tesseroids to take into account the earth curvature, and via the Legendre Quadrature integration technique. The inversion procedure of the Bouguer disturbances integrates the use of hyperparameters in particular the regularization parameter, the anomalous Moho density contrast and the reference Moho depth. This regularized non‐linear inversion method has been tested on synthetic data contaminated by a normally distributed noise. Through this test, the inversion method showed its effectiveness in estimating a Moho model. The application of this method on the Gulf of Guinea has produced a high‐resolution gravimetric model, whose depths vary between 7.75 and 44.50 km. The greatest depths (greater than 40 km) are located along the eastern part of the study area. The comparison of our gravimetric model with the isostatic Airy, GEMMA (GOCE Exploitation for Moho Modeling and Applications) and CRUST1.0 Moho models, globally shows good similarities. The comparison of our gravimetric Moho model to the GEMMA gravimetric model exhibits large residuals around the continental margin, thus reflecting the limit of the inversion method for areas with great Moho depth gradients. The differences observed between our gravimetric model and the CRUST1.0 model in certain regions such as southern Nigeria show the impact of the neglect of anomalous crustal and mantle sources in the accuracy of the solutions obtained after inversion. However, our gravimetric model yields an improved and high‐resolution representation of the Moho undulations over the Gulf of Guinea compared to previous gravimetric Moho models.

(Gravity Recovery and Climate Experiment; Tapley et al., 2004) and GOCE (Gravity field and steady-state Ocean Circulation Explorer; Floberghagen et al., 2011), many global gravitational models have been developed, showing the Earth's external gravitational field at a global scale. With homogeneous and (almost) global coverage, the accuracy and resolution of these models have been improved by combining satellite gravity data with airborne, terrestrial and sea-borne gravity data. Increasingly, these gravitational models are showing their performance for Moho modeling with regard to previous studies in the World (Chen et al., 2018;Chisenga & Yan, 2019;Ghomsi et al., 2020;Uieda & Barbosa, 2017). For an optimal estimation of the Moho geometry using gravity data, some stripping gravity corrections are necessary. Indeed, Bouguer gravity disturbances must be devoid from density heterogeneities within the Earth's crust, upper mantle and deep mantle. Then, the resulting refined gravity disturbances are mainly caused by the Moho undulations beneath the Earth's crust.
For more than a decade, several methods for estimating Moho undulations using gravity disturbances have been revolutionized (Abrehdary et al., 2017;Abrehdary & Sjöberg, 2019Ashagrie et al., 2020;Bagherbandi & Eshagh, 2012;Eshagh, 2014;Eshagh & Hussain, 2016). Among the inversion methods most used by the scientific community, we can include the VMM (Vening Meinesz-Moritz) method (Eshagh, 2016), the direct gravity inversion method (Chen & Tenzer, 2017a), and the Generalized Parker-Oldenburg method (Chen & Tenzer, 2017b). Subsequently, Rathnayake et al. (2021) showed the importance of evaluating the performance of these gravity inversion methods with independent Moho models (such as Airy's isostatic and seismic models). In this study, we mainly use a non-linear gravity inversion method to estimate Moho depth undulations over the Gulf of Guinea. This inversion method developed by Uieda and Barbosa (2017), is based on the discretization algorithm of Uieda et al. (2016) for the forward modeling and on a smoothness regularization to ensure the stability of the inversion. The determination of the regularization parameter is carried out by the hold-out cross validation technique of Kim (2009). Several studies around the world have evaluated and validated the relevance of this inversion method to Moho depth recovery (Deng & Shen, 2019;Li et al., 2022;Sobh et al., 2019). So, we first test the inversion method through a synthetic example and then we estimate the Moho undulations over the Gulf of Guinea using a high-resolution gravity model. This work is organized into six sections. Section 2 describes geological and structural settings of the Gulf of Guinea. The study area and the data used are briefly presented in Section 3. Section 4 provides the methodology applied to estimate the Moho relief. A synthetic example is shown in Section 5. Section 6 presents the different results followed by a discussion. This study is concluded in Section 7. It is worth noting that this work mainly focuses on the validation of the inversion method on the Gulf of Guinea and the computation of a high-resolution gravimetric Moho model at the scale of the study area.

Geological and Structural Settings of the Gulf of Guinea
Figure 1 presents a simplified geological map with some major structures of the Gulf of Guinea. The Gulf of Guinea has been formed between the Upper Jurassic and Lower Cretaceous as a result of the separation of the African and South American plates (Binks & Fairhead, 1992;Guiraud & Maurin, 1992). This phase was characterized by many tectonic activities thus creating block faulting and transform faults. Some volcanic rocks dating from the Middle Jurassic period are present in the western part of the study area. So, the tectonics which have affected the Gulf of Guinea is estimated to have occurred no later than the Middle Jurassic (Dumestre, 1985;Kjemperud et al., 1992). The Gulf of Guinea contains three major fracture zones including: The Saint Paul fracture zone located in the north-west, the Romanche fracture zone and the Chain fracture zone located in the east. The Cameroon Volcanic Line (CVL) could be linked to the various islands and seamounts of the Gulf of Guinea.
From a geotectonic point of view, the Gulf of Guinea region can be subdivided into five provinces.
• The North Gabon rifted terrane: it has been formed during the Upper Cretaceous. This geological structure includes a large rift zone on the continental crust (Karner & Driscoll, 1999;Karner et al., 1997). The faults which have affected this zone have a NW-SE trend and the transfer faults in the NE-SW direction are continuous with the oceanic fracture zones (Teisserenc & Villemin, 1989). The northern limit of this rifted zone is marked by a fracture zone. • The Rio Muni margin: it is located in the south-western part of Cameroon. This margin is oriented in the NNE-SSW direction. It contains the Kribi-Campo and Rio-Muni basins. The continental rifting affecting this margin appear to have the same direction as the Rio-Muni margin. Meyers et al. (1998) have shown that the crust of this part of the Gulf of Guinea is made up of dislocated continental blocks and crossed by dykes, mafic magma and thick volcanic accumulations. Moreover, in this same zone, Rosendahl and Groschel-Becker (1999) identified the Kribi fracture zone following the NE-SW direction. The north of the margin is a whole complex made up of a thick layer of sediment dating from the Tertiary. • The Benue Trough: this structure has a great impact on the geological history of the study area. The Benue Trough has been created after the occurrence of a fracture in the east central Africa, which has begun at the beginning of the Cretaceous (145-66 Ma). This major geological structure progresses toward the Gulf of Guinea in a NE-SW direction. The transition zone between the Benue trough and the Gulf of Guinea is marked by the presence of the Niger Delta, where sediments dating from the Tertiary and more recent are found. It is an aborted rift that is part of the West and Central African rift system (WCARS). • The CVL progresses in the Gulf of Guinea with a NE-SW direction. The southern part of the Cameroon line contains many seamounts, probably born from the Cameroon hot spot. It also contains several volcanic centers namely the Mt Cameroon, the volcanic islands of Bioko, Principe, São Tomé and Annobon, and several submarine volcanoes. Several theories have been proposed on the origin of the Cameroon volcanic line. Burke (2001) thinks that the CVL would have formed following the reactivation of a mantle plume. This plume would start from a new short-wavelength, shallow mantle convection pattern. The submarine hills observed at the southwest of the CVL follow the same direction and are said to have the same source with the CVL. • The oceanic domain is made up of several fracture zones which have been highlighted using satellite gravity and seismic data. For example, we have the Kribi fracture zone. It is a transform fault zone of 75 km wide which progresses in a NE-SW direction.

Study Area and Data
In this section, after presenting the topographic and bathymetric data of the study area, the characteristics of the global gravity field model utilized in the computation process are detailed. The different software used in the computation process and design of the figures are also presented.

Topographic/Bathymetric Data
The study area covers a large part of the Gulf of Guinea and extends from 6° to 13° east longitude and −7° to 8° north latitude. Figure 2 also presents the topography and bathymetry of the study area. The accuracy and resolution of elevation data are essential to ensure good corrections of free-air gravity disturbances. Up to now, one of the bathymetric models with a good resolution at global scale is the 15 arc-second Shuttle Radar Topography Mission (USGS, 2017), that is, SRTM15+ (Shuttle Radar Topography Mission). Figure 2 shows high bathymetric values up to about 6,000 m depth. Maximum altitudes are mainly observed on the CVL and its extension into the Gulf of Guinea. The maximum altitude is around 4,000 m. This altitude is very close to the height of Mount Cameroon (4,037 m) estimated by Kamguia et al. (2015) and which happens to be the highest peak in Central Africa. In the southwest of the CVL, the bathymetry reveals several irregularities attributed to hills or seamounts.
To  active seismic studies, receiver function data and gravimetric constraints. Figure 2 shows a set of 54 points with an estimate of the depth of the Moho. The set of points will be considered as seismic data in the computation procedure.

Gravity Data
Gravity data used to estimate the Moho geometry over the Gulf of Guinea was extracted from the EIGEN-6C4 model. The EIGEN-6C4 is a combined global geopotential model available up to the degree and order 2190 (Förste et al., 2014). This model was developed jointly by the GFZ Potsdam and the GRGS Toulouse (Groupe de Recherche de Géodésie Spatiale). This model was deducted from the combination of data from the Lageos, GRACE, GOCE and DTU missions. EIGEN-6C4 contains the complete SGG (Satellite Gravity Gradients) data of the GOCE mission. The satellite-derived EIGEN-6C4 gravity model can be downloaded from the ICGEM (International Center for Global Gravity Field Models) database at GFZ Postdam.
The satellite-derived EIGEN-6C4 model provides high-resolution gravity data worldwide. This model is widely and successfully used for tectonic and crustal studies (Sahoo & Pal, 2021), or to estimate gravimetric Moho model (Avellaneda-Jiménez et al., 2022;Chisenga & Yan, 2019). In addition, the work of Gautier et al. (2022) have shown that the EIGEN-6C4 model offers a good performance in the representation of the Earth's gravity field over the Gulf of Guinea. The raw gravity data was downloaded with a grid spacing of 0.1° and at a height of 50 km above the ellipsoid. The choice of these parameters is a compromise to calculate the data with optimal resolution and precision.

Software Materials
The procedure to estimate a gravimetric Moho model requires some numerical codes. The main programming language used in this work is Python. The various results obtained were computed from the following open-source libraries: Fatiando a Terra (Uieda et al., 2013; http://www.fatiando.org) to implement forward modeling codes using tesseroids as well as the procedure for inversion, SciPy and NumPy (Jones et al., 2001; http:// scipy.org) to get any results from arrays, matplotlib (Hunter, 2007; http://matplotlib.org) and Seaborn (Waskom et al., 2014; http://seaborn.pydata.org/) to provide a high-level interface for drawing attractive and informative statistical graphics. The scipy.sparse package is applied to certain codes to provide functions to deal with sparse data as linear systems and matrices. The numerical codes containing the different libraries have been compiled from Jupyter notebook (Pérez & Granger, 2007; http://jupyter.org/) in order to visualize results as data, graphics or models. In this work, several figures have been designed in PyGMT (Python interface for Generic Mapping Tools; Uieda et al., 2021). The numerical codes applied to this project have been modified from the source codes freely available (Uieda & Barbosa, 2017; or https://github.com/pinga-lab/paper-moho-inversion-tesseroids).

Methodology
The determination of Moho undulations goes through preliminary corrections applied to the observed gravity. The gravity observed on a point P is actually the addition of the potential created at this point by the reference ellipsoid representing the Earth and all kind of gravity disturbances. In fact, the Earth's density anomalies caused by the variation of Moho relief are contained in gravity disturbances. The other effects are caused by the topography, the sedimentary cover, the mass deficiency of the oceans, the density variations within the crust and the upper mantle… In this work, the corrections applied to gravity disturbances ∆g to obtain the full Bouguer disturbances are the effects of topography ∆g topo , bathymetry ∆g bath and sediments ∆g sed . Other effects are assumed to be negligible.
The sediment-free Bouguer disturbance can defined as follows: So compute these gravity corrections, the forward modeling has been applied using tesseroids via the Gauss-Legendre Quadrature integration technique (Asgharzadeh et al., 2007) improved by Uieda et al. (2016).
The methodology of computation the Moho geometry begins with a discretization of the Moho into tesseroids. A numerical function is therefore established between the gravity of each tesseroid and the Moho depth. Then, we pose an inverse problem to estimate the depth parameters of each tesseroid they fit as much as possible to the observed gravity. Assuming that the vector of depths to be estimated is v, and g 0 the gravity disturbance observed, the gradient-based iterative optimization method enables to find the solution which minimizes the non-linear function φ(v): To ensure a better stability of the inverse problem, we integrated a regularization in the inversion procedure. The most widely used regularization technique is the one established by Tikhonov and Arsenin (1977). The general expression of the smooth regularization function, noted f(v) is: where M is a matrix of finite size n × m, n being the number of tesseroids juxtaposed in the defined grid, and m the number of first-order differences between the depths of the juxtaposed tesseroids.
The estimation of the best regularization parameter is important to smooth the obtained solutions with respect to reality. So, to inverse the full Bouguer disturbance, Bott's method (Bott, 1960) has been required. To improve the stability of this method, a regularization has been applied using Tikhonov's formalism. At any iteration i, the difference vector ∆v i in the Tikhonov regularization is estimated by solving the following nonlinear equation: with A i being a Jacobian matrix of size n x m and which terms are defined by: The regularized inversion procedure involves the estimation of some hyperparameters. These hyperparameters have a real impact on the optimization of the solutions obtained after inversion. The first hyperparameter is the regularization parameter (μ) which smoothes and stabilizes the solutions obtained after inversion. In this work, the estimation of this parameter has been done using the cross-validation hold-out method (Kim, 2009). The other two hyperparameters are the reference depth of the Normal Earth's Moho (z ref ), and the density contrast of the anomalous Moho (∆ρ). The estimation of these two parameters requires the information on Moho depth at certain points in the study area, for example, Moho depths derived from seismological data.
For more details on the equations applied in the inversion method and the theory behind the estimation of inversion hyperparameters, the ready is kindly referred to the work of Uieda and Barbosa (2017).

Synthetic Example
In this section, the inversion method is tested to evaluate its performance on the study area. The global CRUST1.0 is chosen to apply the inversion method. In fact, the Moho depths extracted from the CRUST1.0 model are used as basic data. A forward modeling is applied to get the gravity information on the Gulf of Guinea and then, the inversion is applied to estimate the anomalous Moho depth. The optimal regularization parameter (μ) is estimated by cross-validation while the Moho density-contrast (∆ρ) and the mean reference depth (z ref ) are estimated using seismic points. This synthetic example also aims to evaluate the procedure of estimation of hyperparameters and their effectiveness in optimizing the solutions obtained by the inversion method. Figure 3 shows a simple sketch illustrating a section of the lithosphere. The reference depth z ref of the Normal Earth's Moho is materialized in Figure 3.
Thus, a discretization is applied to the global CRUST1.0 model to obtain a grid of tesseroids (n lat × n lon = 25 × 33, i.e., a total of 825 elements). The discretized model is presented in Figure 4a with a spatial sampling step at 0.5° × 0.5°. For the forward modeling of the gravimetric information (synthetic gravity data), we assume that the contrast parameters of density and reference depth of the Moho are respectively 500 kg/m 3 and 20 km. A gravimetric perturbation (noise) of normal distribution with a zero mean and a standard deviation of 5 mGal is introduced in the synthetic model. Figure 4b presents the noise-corrupted synthetic gravity data generated from the discretized model. The hyperparameters are estimated by cross-validation and from seismic data. The cross-validation consists of calculating the MSE (mean square error) statistical value over a range of parameters and choosing the optimal parameter that minimizes the MSE. The determination of the optimal regularization parameter is done over a set of values ranging from 10 −7 to 10 −2 . As shown in Figure 5a, the MSE is calculated for each value of μ and the optimal value is obtained for μ = 2.15 × 10 −4 . From this optimal value and the generated seismic points (Figure 2 Thereby, the hyperparameters estimated have been used in the synthetic data inversion procedure. Figure 6a presents the estimated Moho depths with a spatial sampling of 0.5° × 0.5°. We notice that the computed solutions are quite smooth with a good continuity in the arrangement of the tesseroids, which indicates the efficiency of the optimal regularization parameter determined by cross-validation. Figure 6b shows the differences between the  (Figure 5a) which is about 21.79 km, the uncertainty of Moho depths can be estimated. So, an uncertainty of 1.07 km over 27 km represents a relative uncertainty of about 5%. This value means an approximate accuracy of 5% on the depth of the Moho computed. These statistical results testify of the performance of the inversion method in estimating Moho depths. In fact, the largest differences are observed near high-gradient areas on the Moho depth map (Figure 4a), especially around continental margin areas. Figure 6.c presents the gravity residuals, representing differences between observed and predicted gravity disturbances. The gravity residuals are obtained by applying a forward modeling to the Moho depth residuals. Figure 5c shows a normal distribution of these gravity residuals, with a mean of 0.05 mGal and a SD of 5.29 mGal. These statistical values correlate well with those of the gravimetric perturbation introduced above (0.00 mgal; 5.00 mgal, respectively). Considering that the inversion process acts as a filter noise on the observed gravity data, these results suggest that the disturbing noise introduced in the simulation has been correctly removed by the inversion method. Similarities are observed between signatures of the Moho depth residuals and those of the gravity residuals. Figure 6d shows the differences between seismic Moho depths ( Figure 2) and the estimated Moho depths. The differences range from −1.4 to 1.9 km, with a mean of 0.11 km and a SD of 0.64 km. Despite the contamination of the synthetic data with a normally distributed noise, the inversion method applied in this study proves to be effective in estimating a Moho depth model.

Results and Discussion
In this section, the results of the computation procedure of Moho depths in the Gulf of Guinea are presented. As we pointed out before, the tesseroids are used for the discretization of the models to respect the curvature of the Earth (we assume a radius of 6,378.137 km as we are close to the equator in this study area). Crustal and mantle heterogeneities are neglected and the lateral variation of Moho density contrast is assumed to be zero. Before applying the inversion method, gravity corrections must be properly calculated and applied to obtain full gravity disturbances.

Computation of Full Gravity Disturbances
The gravity data used here are extracted from the EIGEN-6C4 model. The gravity disturbances in the Gulf of Guinea are shown in Figure 7a. The topographic effects (Figure 7b) are calculated using a tesseroid model and the digital terrain model SRTM15+. Conventional densities are used in the computational algorithm, 2670 kg/m 3 for the topography and −1630 kg/m 3 for the bathymetry. Then, the computed topographic corrections are applied to the gravity disturbances to obtain the Bouguer disturbances (Figure 7c). The topographic corrections are greatly significant in the study area (they vary from −324 to 85 mGal). The differences of trends between gravity disturbances ( Figure 7a) and Bouguer disturbances (Figure 7c) are noticeable. Also, the study area is covered by a large sedimentary layer. Figures 7d, 7e, and 7f show the upper, middle and lower sediment thicknesses respectively. The gravimetric effects of the three sedimentary layers have been calculated using a discretized tesseroid model. Non-negligible values ranging from −2 to −102 mGal are obtained. After correcting the gravimetric effect of sediments, the complete Bouguer disturbances are obtained. These Bouguer disturbances will be introduced into the inversion algorithm to estimate the anomalous Moho depth.

Estimation of Hyperparameters and Inversion of Gravity Disturbances
The estimation of Moho depths starts with the determination of the optimum hyperparameters. The regularization parameter μ is estimated using cross-validation. By fixing the Moho reference depth and density-contrast parameters respectively to 20 km and 500 kg/m 3 , the regulation parameter is estimated on a logarithmic scale of values ranging from 10 −10 to 10 −2 . As shown in Figure 8a Figure 9, we observe a great lateral gradient along the continental margin. On the continental plate, the computed gravimetric model estimates the greatest depths located at the west of the study area (the Moho depth values oscillate averagely between 40 and 45 km). On the other hand, in southern Nigeria, an uplift of the Moho can be observed down to an average depth of 30 km. Around the continental margin, the Moho geometry varies brutally (between 25 and 30 km in depth).
In the oceanic domain, the lowest depths of the estimated gravimetric Moho are located on either side of the prolongation of the CVL in the Gulf of Guinea.

Comparison of Moho Models and Discussion
To further analyze the gravimetric Moho undulations computed in this study, a comparison is made with some classic Moho models, in particular a Moho model according to Airy's isostatic hypothesis (Airy, 1855;Heiskanen, 1931), a gravity-based Moho model extracted from GEMMA (Reguzzoni & Sampietro, 2015) and a CRUST1.0 seismic Moho model. According to the law of isostasy, the anomalous Moho undulations are compensated by the variations of the relief on the surface of the Earth. Large topographic masses would be compensated by great Moho depths. The comparison of an isostatic Moho model with a gravimetric Moho model is important to better understand the isostatic compensation phenomena that dominate a region. Figure 10.a shows the Moho depth calculated according to the Airy isostatic theory. To obtain these isostatic depths, the density contrast parameter (ρ = 450 kg/m 3 ) and the average reference Moho depth (z ref = 32.5 km) evaluated above have been used. The SRTM15+ model (Figure 2) provides topographic and bathymetric data to apply Airy isostatic hypothesis on the   study area. Table 1 presents the statistics of the isostatic Moho model obtained. The isostatic Moho depths range from 12.6 to 40.60 km. At first glance, the isostatic depths seem to have the same trends with our gravimetric Moho model (Figure 9). Strong gradient zones are observed during the transition from the continental plate to the oceanic plate. The compensation of the surface relief seems to be respected in view of the great depths observed on the CVL. The differences between the isostatic model and our gravimetric model range from −9.32 to 9.47 km, with a mean of −2.30 km and a standard deviation of 2.67 km (Figure 10c). These differences are actually signatures of an isostatic compensation. Therefore, we can observe three possible situations. When the residual values are positive, there is overcompensation (the relief of the gravimetric Moho is below the relief of the isostatic Moho). When the residual values are negative, there is undercompensation (the relief of the gravimetric Moho is above the relief of the isostatic Moho). The isostatic compensation is considered to be complete for zero residual values. In view of the signatures modeled in Figure 10c, the isostatic overcompensation occurs mainly on the continental plate. The areas most affected by this overcompensation are located along the eastern part of the study area. The undercompensation would affect the continental margin in a particular way. The similarities between the areas most marked by an overcompensation and a dense sedimentary cover (Figure 7e) are quickly made. The oceanic part of the study area is generally undercompensated with respect to the continental part. However, Figure 10c also shows areas where the isostatic compensation is almost total.
In addition, an evaluation of our gravimetric model has been made compared to the global gravimetric model GEMMA. The GEMMA Moho model has a spatial sampling of 0.5° × 0.5° and has been designed from GOCE satellite and seismic data. Figure 10b shows the spatial distribution of Moho depths according to the GEMMA model. The Moho depths range between 7.35 and 37.51 km, with a mean of 17.73 km and a standard deviation of 9.02 km (Table 1). Unlike the isostatic model, the GEMMA model shows shallower Moho depths in some areas. For example, we observe Moho depths less than 10 km below areas marked by the Chain and Romanche fractures, and also below the rifting tectonic plate observed at the west of Gabon. Also, contrary to the isostatic model, the GEMMA model shows an amplification of Moho depths (over 35 km) under some areas, particularly in southern Ghana, and along the eastern part of the study area. The differences between our gravimetric Moho model and the GEMMA model are shown in Figure 10c. These differences range between −5.72 and 11.19 km, with a mean of 2.95 km and a standard deviation of 2.32 km ( Table 1). The greatest differences are mainly observed along the margin of the continental plate, with some negative residual values in the northwest and positive residual values at the west of the study area. However, these differences are less marked in the oceanic part of the study area.
The gravimetric Moho model estimated in this study is also compared to the CRUST1.0 seismic Moho model. The differences between the two models are shown in Figure 10e. Table 1 presents the statistics of these differences. The CRUST1.0 seismic Moho model closely agree with our gravimetric Moho model compared to the isostatic and GEMMA models. This could be due to the use of seismic points generated from the CRUST1.0 model in the inversion algorithm of our gravimetric Moho model. Nevertheless, the differences ranging from −11.60 to 9.20 km are observed in Figure 10e. This could be the consequence of using a constant density contrast in the inversion procedure. Thus, positive residual values could be observed in areas where the density contrast is less than 450 kg/m 3 , and negative residual values for areas where the density contrast is greater than 450 kg/m 3 . However, Chen and Tenzer (2017b) showed that applying a variable density contrast in the inversion algorithm of gravity disturbances could lead to very unrealistic solutions. Furthermore, as also claimed by Chen et al. (2018) in their work, these differences could also be due to the low precision of data extracted from the CRUST1.0 model, in particular the large errors in the uppermost density data. The differences observed between the CRUST1.0 seismic model and our gravimetric model also show us potential areas where crustal and mantle anomalies, assumed to be zero initially, have an impact on the solutions obtained from the inversion algorithm. For example, the CVL and its extension in the Gulf of Guinea is dominated by positive residual values. This would be the consequence of the presence of anomalous structures in this part of the crust or upper mantle. Also, the lowest residual values are observed in southern Nigeria, thus showing an underestimation of the Moho depth in this area by our gravimetric model. This could be explained by the presence of bodies much denser than the crust. The estimated gravimetric model would be even more accurate if all anomalous sources were subtracted from the gravity disturbances, leaving only those caused by the Moho undulations.
Despite these differences caused by the complexity of the Moho geometry in the study area, our gravimetric Moho model generally shows good similarities with the other models. The RMS values are 3.52, 3.75, and 2.98 km in comparison with respectively the isostatic, GEMMA and CRUST1.0 models. In addition, the inversion algorithm used to present the undulations of the Moho seems satisfactory. In addition to having a good resolution, the gravimetric model obtained shows superficial Moho depths (up to 7.35 km) absent in the isostatic and CRUST1.0 models. Also, the greatest depths are developed in our models (up to 44.5 km), while only the GEMMA model shows a maximum depth of 37.51 km in the study area. Nevertheless, the main aim of this work was to compute a high-resolution gravimetric Moho model over the Gulf of Guinea based on a regularized non-linear gravity inversion.

Conclusion and Outlook
In this work, a high-resolution gravimetric Moho model has been computed beneath the Gulf of Guinea based on a regularized nonlinear inversion algorithm. One of the advantages of this inversion method is that the Earth's curvature is taken into account through the use of tesseroids in the discretization of models. Moreover, the method used incorporates a regularization technique in the inversion procedure, to ensure the stability of solutions and to obtain a realistic Moho depth model. The test of this inversion method in a synthetic example has shown its effectiveness in producing a smooth Moho model, with a density contrast assumed constant. Also, the cross-validation technique to determine the optimal regularization parameter proves to be efficient for the algorithm. Despite the injection of noise in the gravity disturbances, the inversion method applied in this study has been able to show its stability. However, the inversion method seems to have a poor performance for areas with high variations in Moho depths (e.g., the continental margin in our study area).
The application of the regularized non-linear inversion method to the Gulf of Guinea yielded a high-resolution gravimetric Moho model. Our computed model shows that the Moho undulations in the Gulf of Guinea vary between 7.75 and 44.50 km in depth. Despite the differences observed with the Airy's isostatic, GEMMA and CRUST1.0 Moho models, our gravimetric Moho model shows good similarities with these models. The resolution and accuracy of our gravimetric model are assessed with regard to its ability to estimate the most superficial Moho depths in the study region, not identified by the Airy's isostatic and CRUST1.0 models. The differences established between the CRUST1.0 seismic model and our gravimetric model showed that the neglect of anomalous crustal and mantle sources in the correction of gravity disturbances could have a great influence on the accuracy of the computed Moho model. For example, the negative residual values observed in southern Nigeria could reflect the presence of high-density structures in the crust. The influence of the use of a constant density contrast could also affect the quality of the solutions obtained after inversion. Nevertheless, some previous works have shown that the use of variable density contrast could lead to very unrealistic solutions.
Up to now, the estimation of a gravimetric Moho model remains an alternative for geodynamic and tectonic studies given the limitations of the developed inversion methods. Therefore, the most effective option remains the densification of the study area with seismic data. Improving inversion methods using gravity data involves to take into account anomalous sources in the crust and in the upper mantle. Also the use of a variable and accurate density contrast is essential for a more realistic estimation of the Moho anomalous relief.