EVALUATION OF GGMs BASED ON THE TERRESTRIAL GRAVITY DISTURBANCE AND MOHO DEPTH IN AFAR, ETHIOPIA

. To estimate Moho depth, geoid, gravity anomaly, and other geopotential functionals, gravity data is needed. But, gravity survey was not collected in equal distribution in Ethiopia, as the data forming part of the survey were mainly collected on accessible roads. To determine accurate Moho depth using Global Gravity Models (GGMs) for the study area, evaluation of GGMs is needed based on the available terrestrial gravity data. Moho depth lies between 28 km and 32 km in Afar. Gravity disturbances (GDs) were calculated for the terrestrial gravity data and the recent GGMs for the study area. The model-based GDs were compared with the corresponding GD obtained from the terrestrial gravity data and their differences in terms of statistical comparison parameters for determining the best fit GGM at a local scale in Afar. The largest standard deviation (SD) (36.10 mGal) and root mean square error (RMSE) (39.00 mGal) for residual GD and the lowest correlation with the terrestrial gravity (0.61 mGal) were obtained by the satellite-only model (GO_CONS_GCF_2_DIR_R6). The next largest SD (21.27 mGal) and RMSE (25.65 mGal) for residual GD were obtained by the combined gravity model (XGM2019e_2159), which indicates that it is not the best fit model for the study area as compared with the other two GGMs. In general, the result showed that the combined models are more useful tools for modeling the gravity field in Afar than the satellite-only GGMs. But, the study clearly revealed that for the study area, the best model in comparison with the others is the EGM2008, while the second best model is the EIGEN6C4.


INTRODUCTION
Afar depression is found at the junction of the Gulf of Aden (GOA), Red Sea Rift (RSR), and Main Ethiopian Rift (MER) where earthquakes, plate tectonics, and volcanism happen. The study area is at a critical stage, i.e., the place where the continental rift is changing into oceanic ridge and also it is an ideal locale to study the role of extension and magmatism as rifting progresses to seafloor spreading. Furthermore, the crustal thickness of the Afar rift is attenuated due to the upwelling magma and also the rift is tectonically, volcanically, and 79 seismically active. To estimate the crustal thickness, geoid, gravity anomalies, and other geopotential functionals, the gravity data is needed. Gravimetry is done on, near above (airborne), and far (space-borne) from the Earth's physical surface with different geographical resolutions. The lack of short wavelength (near surface effect) and downward continuation problem occurs on the air-and space-borne data due to weakness of gravitational field. Terrestrial gravity data have full gravity field with varied density of the earth. As Bolkas et al. (2016) and Novák (2010) stated, combined global gravity models (GGMs) are derived from the combination of terrestrial, air-borne gravity, and space-borne gravitation data. But, satellite-only GGMs are derived from satellite observations only. These models characterize the gravitational field mathematically through spherical harmonic functions.
GGM is a global gravity field model of the Earth that can be approximated by mathematical function explicitly describing Earth's gravity field in 3D space, and used for different geodetic and geophysical applications (Barthelmes, 2014). Advancement of space-based technologies and their launches like Gravity field and steady-state Ocean Circulation Explorer (GOCE) (Floberghagen et al., 2011), Gravity Recovery And Climate Experiment (GRACE) (Tapley et al., 2004), and CHAllenging Minisatellite Payload (CHAMP) (Reigber et al., 2002) plays a remarkable role in determining mathematical model of Earth's gravitational field, GGM. These satellite missions improve our understanding of global gravity field with static model and temporal variation by many GGMs (Godah et al., 2017).
The main objective of this study is the evaluation of the accuracy of the following recent satellite-only models: the Sixth release of the GOCE gravity field model through the direct approach (GO_CONS_GCF_2_DIR_R6) , and combined high-degree GGMs, by employing terrestrial gravity data-including eXperimental Gravity field Model 2019 (XGM2019e_2159) (Zingerle et al., 2019), European Improved Gravity model of the Earth by New techniques (EIGEN6C4) , and Earth Gravitational Model of 2008 (EGM2008) (Pavlis et al., 2008)-for approximating the Earth's gravity field. The terrestrial gravity data and the location data (XYZ) were collected by LaCoste-Romberg Gravimeter and Trimble Differential Global Positioning System (DGPS), respectively, in 2017by Geological Survey of Ethiopia (GSE). The ground-based gravity data in the study area were used to assess the performance of GGMs, which best coincides with the study area for gravity field modeling at a local scale, and the comparison results are presented with regard to the range (maximum and minimum), mean, correlation, standard deviation (SD), and root mean square error (RMSE) over the study area. There had never been any previous research on evaluating GGMs utilizing terrestrial gravity data based on gravity disturbance (GD) over Ethiopia, particularly in the Afar region. However, there were three MSc thesis works: Evaluation of Accuracy of Earth Gravity Model 2008 (EGM2008) using Global Positioning System (GPS) and Levelling at Debre Birhan city (Ermias, 2015); central and western Ethiopia (Birbiraw, 2015); and Evaluation of Gravity Field Models: EIGEN-6C4 and GOCO03S combined with EGM08 using GNSS-Levelling at Legedadi and its surrounding areas (Zerihun, 2017).
This study is particularly useful for countries like Ethiopia, where terrestrial gravity data is rare and the region is tectonically, volcanically, and seismically active; earthquakes, plate tectonics, and volcanism all occur in the studied area. This indicates that by identifying the best GGMs from this study or using the similar methods mentioned here, the GGMs can be utilized to produce more precise and reliable results of Moho depth, lithospheric thickness, and geoid. Other researchers can utilize the findings of this study for various geodetic and/or geophysical applications for the entire country or specific places. Furthermore, researchers can utilize this technique to evaluate a large number of GGMs across Ethiopia and recommend the best GGM/s for various applications. 80

DESCRIPTION OF THE STUDY AREA
The study area is found at the junction of the GOA, the RSR and the MER Valley. This rifting has continued up to the present time to form the GOA, the RSR, and the MER; these intersect in the Afar Depression, which covers an area of ~200,000 km 2 . The study area's geographic location points of the corners are 39. 22°E, 14.44°N; 42.35°E, 14.43°N; 42.28°E, 8.63°N; and 39.22°E, 8.66°N. The evaluation procedure of GD refers to a ground measured gravity data set over the study area that is comprised of 1933 data points (Black Points in Figure 1). The map was plotted using a Python interface for the Generic Mapping Tools (PyGMT).

THEORETICAL BACKGROUND
Earth's gravity field is affected by various sources of tidal, radial distance, subsurface density variation and topographic relief, and rotation of the earth field. Some of the sources will be removed using practical Earth's model from the measured gravity and retain only the subsurface effect for geophysical applications. The GD is defined as the magnitude of the gradient of the potential (gravity) at a given point minus the magnitude of the gradient of the normal potential (normal gravity) at the same point: gravity (λ, φ, h) minus normal-gravity (λ, φ, h). In Heiskanen and Moritz (1967), GD, = − is the difference between the magnitudes of the gravity vector, = ‖g ‖ and the normal gravity vector, = ‖γ ‖ at the same point P.
81 Figure 2. The gravity vector , normal gravity vector P, GD vector , and unit vector uP at a point P on the surface of the Earth (Oliveira et al., 2018) An Earth ellipsoid, also known as an Earth spheroid, is a mathematical figure used as a reference frame that approximates the shape of the Earth. It is a sphere with a regular shape. Ellipsoid height (h), also known as geodetic height, is the distance between an ellipsoid and a point on the Earth's surface. Unlike the geoid, the ellipsoid presupposes that the Earth's surface is smooth and that the Earth is homogeneous. Furthermore, ellipsoid assumes that there are no mountains or trenches, and that the ellipsoid surface and mean sea level (msl) are identical.
The geoid is an irregular-shaped "ball" that is used to estimate accurate surface elevations and is a model of global msl. It is the shape where the ocean surface would take under the influence of Earth's gravity and rotation alone. If the effect of winds tides and other factors were ignored, a geoid's surface depicts msl or a supposition of the ocean's surface. The gravitational field of the earth is the single component that influences the shape of the msl. The discrepancy between the ellipsoidal height and the orthometric height is the geoid height (N) above the ellipsoid. Unlike ellipsoidal models, geoid is based locally-or at least more local than the entire surface of the Earth-and it changes and forms an irregular shape. Because the geoid changes in different parts of the planet, the geoid height for the same object or location varies across the globe. The accuracy of the geoid height contributes to the orthometric height's overall accuracy.
The vertical distance along the perpendicular (plumb) line from a point on the Earth's surface to the geoid, the vertical datum that approximates msl, is known as the orthometric height (H). The distance between a place on the Earth's surface and the geoid is the orthometric height of that point. The distance between the point being measured and the Earth's geoid is expressed by the orthometric height, often known as elevation. Because the geoid varies around the world, the orthometric height can be used to identify heights above and below sea level.
The discrepancy between actual gravity and normal gravity obtained at the same point (P) is the GD at any point in space, whereas the gravity anomaly is the discrepancy between observed gravity at P and the normal gravity on the geoid (Q), which is the point where the normal to the ellipsoid at P intersects the geoid.
Gravitational potential at a point P is the potential energy, due to the Earth's gravitational attraction, of a unit mass situated at P. Or it is equal to the work done if a unit mass is brought from infinity to the point P under the influence of the Earth's gravitational field. Normal 82 gravity potential is derived as the potential of a level ellipsoid plus the rotational potential of the ellipsoid.
The gravitational field of a reference ellipsoid of rotation that closely approximates the real Earth is called the normal gravity field. The anomalous field is the discrepancy between the real and the normal fields. As the Earth is approximated by perfect ellipsoid, the gravity field of ellipsoid of revolution is a good approximation to the Earth's gravity field and can thus be considered as a normal terrestrial gravity field. An ellipsoid of revolution-the reference ellipsoid-that is close to the geoid is defined to be an equipotential surface of this normal field, where normal gravity potential is constant.
On and above the level ellipsoid, normal gravity, or the gradient of the normal gravity potential, is shown. Normal gravity is the gravitational force acting on such a "level ellipsoid." The gravity vector represents the force of gravity on a unit mass and is the gradient vector of the gravity potential. The intensity of gravity, or gravity g, is the magnitude of the gravity vector. The force per unit mass, or acceleration, is a dimension of gravity.
Plumb line is a line that, at any given place in space, coincides with gravity's direction. A plumb line's direction corresponds to the direction of a string to which a freely suspended weight, or plumb bob, is attached. In terms of gravitational force, this vector is perpendicular to equipotential surfaces everywhere. The Plumb line is defined by two angles: geographical latitude and longitude, which are determined using astronomical methods to an accuracy of at least 0.1 s of arc. The gravitational vector can be completely determined if g, and are known. Because of the earth's spheroidicity and irregular interior structure, the equipotential surfaces are not parallel to each other.
The value of the normal gravity, γP, is the magnitude of the normal gravity vector. It is the magnitude of the gradient of the normal potential (contains centrifugal potential). It can be calculated at a given point on the ellipsoid (λ, φ, h = 0) and Earth's physical surface (λ, φ, h).
In geodesy, there are two types of anomaly variations: gravity anomalies and GDs, based on gravitational reference model. Gravity anomaly, in geodesy, is the difference between geoidal gravity and normal gravity on the mathematical model ellipsoid (Heiskanen and Moritz, 1967) -but in geophysics, gravity anomaly is understood as the difference between observed gravity field and the field of a reference model. However the GD may be more practical and conceptually more reasonable for geophysical investigations. In general, GD considers the same point and the centrifugal effect can be ignored but the anomaly is not. This means that spherical harmonic function is considered by GD but not by gravity anomaly (Olivera et al., 2018). Both parameters can be either vector or scalar quantity. The harmonic function is the foundation for all processing gravity field. As a result, GD is a valid selection for characterizing gravitational effects produced by subsurface density variation. As Hackney and Featherstone (2003) put it, the GDs are more suitable for geophysical applications.

Global gravity models (GGMs)
GGM is a mathematical expression of the gravity of the Earth based on the combination of terrestrial, air-borne gravity, and space-based gravitational data. Different GGMs (listed in the introduction part) were used for determining the GD at the Earth's topography. The GGMs used in this study described as follows: • The Sixth release of the GOCE gravity field model by means of the direct approach  (v5) and the short wave is from Danish Technical University-13 (DTU13) (Here a smoothed and gridded EGM08 is used).
• Earth Gravity Model of 2008 (EGM2008): It is developed from altimetry, satellite (GRACE), and ground data. It is static global combined gravity field model with d/o of 2190. It has the long wave from early GRACE and the short wave from land gravity and marine and altimetric gravity (Danish National Space Centre (DNSC08)).
• European Improved Gravity model of the Earth by New techniques (EIGEN6C4): It is developed from altimetry, satellite (LAGEOS, GRACE Release03/RL03/ Groupe de RecherchesenGeodesieSpatiale/GRGS/, and GOCE-SGG data), and ground data. The bandlimited combination of normal equations (up to max 370°) is done from the combination of these satellites and terrestrial data, which are generated from observation equations for the spherical harmonic coefficients (Yilmaz et al., 2016). It is static global combined gravity field model with d/o of 2190. It has the long wave from GRACE+GOCE (v4) and the short wave from DTU13, also on land, where a gridded EGM is used.
• eXperimental Gravity field Model 2019 (XGM2019e_2159): It is represented through spheroidal harmonics, corresponding to a spatial resolution of 2' (~4 km). It is developed from altimetry, satellite model (Gravity Observation Combination06s/GOCO06s/), topography, and ground data. It is static global combined gravity field model with d/o of 5399. It has the long wave from GRACE+GOCE v5 and the short wave from land gravity averaged at 1/15° and marine and altimetric gravity (DTU15).
Earth's GGM is determined by the combination of terrestrial, air, marine, ship, and spacebased gravity anomalies derived using spherical harmonic analysis (Rummel et al., 2002). GGMs help to understand the Earth and provide knowledge regarding the Earth, its shape, its interior, and fluid envelope where all gravity fields functional are calculated from these models. The GD ( , ) can be expressed mathematically using spherical harmonic expansion as follows (Barthelmes, 2013): where: , -GD in spherical approximation at point P; -Earth's gravitational constant; radius (geocentric distance); -spherical longitude; spherical latitude; mean radius of the earth; -degree; -order; associated Legendre functions of first kind; and spherical harmonic coefficients; maximum degree of expansion; coefficients of expansion.
The launches of GOCE, GRACE, and CHAMP have led significant achievements in the determination of the Earth's gravitational field. Thus, the technological and scientific developments in artificial satellite techniques and calculation algorithms resulted in releasing high-degree combined GGMs (Yilmaz et al., 2017). In this paper recent satellite-only and high-degree combined models like GO_CONS_GCF_2_DIR_R6, EIGEN6C4, EGM2008, and XGM2019 are studied. The GD of these GGMs was compared with the ground gravity data to check the accuracy of the GGMs.

METHODOLOGY
The ground measured gravity data were collected in the Geodetic Reference System-1980, which is connected to International Gravity Standardization Network-1971 (IGSN71) system found in Addis Ababa. The gravity data and position of the stations were collected by LaCoste-Romberg gravimeter and Trimble DGPS. The collected data point is shown in Figure  1,and the data produced in this figure was collected along the road found in the study area with 2 km stations interval. The GGMs were downloaded from the International Center for Global Earth Models (Ince et al., 2019) (https://www.icgem.gfz-potsdam.de/) data center (GO_CONS_GCF_2_DIR_R6, EGM2008, EIGEN6C4, and XGM2019e_2159) to evaluate with the ground gravity data for the study area.
The ground measured and GGMs were processed to obtain the GD separately at the observation point located anywhere on the topo-surface. The ground data's GD was calculated by determining the normal gravity, which is depending on the latitude and average elevation. The procedure first requires determination of the normal gravity on the ellipsoid (factor of latitude) for each field station, followed by computation of the height dependence normal gravity on the Earth's physical surface (factor of latitude and elevation). The spherical harmonic coefficients of GGMs were processed and then the gridded GDs for the study area were generated; these were calculated on the Earth's topography. Then, the GD is extracted for each ground station from the gridded data, which helps to perform comparative analysis.
The statistical evaluation of GDs was analyzed based on the simple statistics on the GDs and their respective differences. This difference is the difference between the ground GD (δgP,g) and the GD calculated from GGMs (δgP,GGM), which is defined as the residual GD (ΔδgP). It is given by: where: refers to the i th number of data points (i=1−1933).
The statistical analysis of GDs (δgP) and residual GDs (ΔδgP) was carried-out with the range (R) [minimum and maximum], mean (μ), SD (σ), RMSE, and correlation (r) values as the common criteria for the accuracy, relation, and variation. These statistical parameters are given by: 85

Statistical Analysis of GD (δgP)
where: number of data points (field stations); δg -mean of gravity disturbance.

Statistical Analysis of Residual GD (ΔδgP)
where: δg -mean of residual gravity disturbance.
The ground measured data and GGMs were gridded with grid size 0.25 km to enable making a 2D map for each of them. The grid map was used for qualitative interpretation of the comparative evaluation of the GGMs. From the grid maps, the long-wavelength profile curves along and across the rift escarpment (shown in Figure 3a.) were plotted for ground data and GGMs. The plotted profile curve was used for the qualitative interpretation of the comparative analysis of the GDs. The profile curve was analyzed based on the above statistical parameters to check the accuracy of the long-wavelength signal of the GGMs based 86 on the ground measured data. Again, the same procedure (generating 2D maps) was done for the residual GD for qualitative and quantitative interpretation. Here, the evaluation of GGMs focused on the statistical analysis of GD and residual GDs.

GD of the ground data and GGMs
In comparative analysis, the accuracy of GGM was obtained from the GD and its residual generated from the measured gravity discrete data; and the GGM that best fitted with the terrestrial data were selected. The comparative technique for the selection of GGM was based on the statistical analysis of GD and its residual. The statistical values of the GD from each data are given in Table 1. The high correlation coefficient between the terrestrial and GGM GDs in Table 1 does not indicate a high level of accuracy; rather, it indicates a high level of signal structure, which implies whether the GDs from the terrestrial and GGM are in phase or out of phase at each of the field stations (black points in Figure 1). The closer the value of correlation coefficient to ±1, the more accurate the structure of the signal of GGM-derived GD to the terrestrial GD, while the farther the value of correlation coefficient from ±1, the less accurate is the structure of the signal of GGM-derived GD to terrestrial GD at each of the terrestrial data points.
From Table 1, one can conclude that the Sixth release of the GOCE gravity field model by means of the direct approach (GO_CONS_GCF_2_DIR_R6) is not the best model for the study area but the EGM2008 can be selected as the best of all and the next best model is EIGEN6C4. The final best model among the combined model is XGM2019e_2159.
To specify the occurrence of the spatial variation and magnitude of GD, the graphical representations were used for the qualitative interpretations by producing their respective 2D grid map of GD obtained from the ground measured data and the GGM, using Oasis Montaj 8.4 (shown in Figure 3). From Figure 3, one can conclude that the first four maps are almost identical in spatial variation, minimum, and maximum values of gravity distance except for the GO_CONS_GCF_2_DIR_R6 model. It differs from the rest in maximum and minimum values of gravity, as shown in the color bar, and does not contain short wave length as the others,which indicates that it underestimates the ground data.

Profile curve of GD of the ground data and GGMs
From Figure 3, the long-wavelength (for >167 km) profile curve along and across the main rift escarpment in the study area was generated for each GD map. All curves for each profile shown in the above map were plotted in a single map for qualitative interpretation. The statistical values of the GD obtained from the ground measured data and GGMs and its residuals along and across the escarpment are given in Tables 2 and 3.  The EGM2008 and EIGEN6C4 models have the best signal structure relationships with the terrestrial GD compared to the other GGMs, although the XGM2019e 2019 model has a better signal structure compared with the terrestrial GD than the satellite-only model (GO CONS GCF 2 DIR-R6) in Figure 4; this is also applicable for Figure 5.
From Table 2 and Figure 4, it can be ascertained that EGM2008 and XGM2019e_2159 is the best model for the study. From Figure 4, the pattern of the three curves (EGM08, EIGNE6C4, and XGM2019_e2159) is identical with the ground data except that it slightly varies by magnitude in different places, i.e., the signal contains the short and long wavelength effect of the subsurface. This means that the Sixth release of the GOCE gravity field model is not the best model for the study area.   Table 3 and Figure5, it can be ascertained that the GO_CONS_GCF_2_DIR_R6 is not the best model for the study area. As shown in Figures 4 and 5, the majority of the differences are only at the shorter wavelength between the terrestrial data and GGMs, but the long wavelength pattern for the terrestrial data and GGMs are more or less identical, and vary only in magnitude, except for GO_CONS_GCF_2_DIR_R6.
The relatively large values obtained from the RMSEs of the GGM computed/terrestrial GDs could be partly due to the omission and commission errors inherent in the GGMs, possible systematic errors in the observed terrestrial data, and the topographic bias arising from truncation of the topography model and the terrestrially acquired elevation data. This is so because there is always a bias between DEM data and terrestrial data as a result of the great deviations in gradient of the undulating terrain when point values are compared to mean values of DEM. The RMSE of Profile AA' (Table 3) is higher than that of Profile BB' (Table  2), which could be due to topographic and lithological variations in Profile BB'.

Residual GD between the ground data and GGMs i. 2D map
The best approach for comparative analysis is to use the residual GD, which is obtained by subtracting the GGMs gravity disturbance from the ground GD. Here, the evaluation of GGMs is focused on the residual GDs. The statistical values of the residual GDs, calculated from the difference between measured data and GGMs, are given in Table 4. From Table 4, one can conclude that the Sixth release of the GOCE gravity field model by means of the direct approach (GO_CONS_GCF_2_DIR_R6) is not the best model, when compared with the other models for the study area, but the SD and RMSE are small for EGM2008 and EIGEN6C4, respectively; this implies that both the models can best fit the study area.
The closer the value of the RMSE to zero, the more accurate is the GGM-derived GD in relation to the terrestrial GD, while the farther the value of the RMSE from zero, the less accurate is the GGM-derived GD in relation to the terrestrial GD.
The values of the RMSE show the level of accuracy of the GGM-derived GDs in relation to the terrestrial GDs.
The graphical representations were used for the qualitative interpretation of GGMs by producing their respective residual GD 2D maps for each GGM, using Oasis Montaj 8.4 (shown in Figure 6).  Figure 6, the residual GD graphical representations seem identical but the difference between the GD obtained from measured data and the GGM (GO_CONS_GCF_2_DIR_R6) have the largest maximum and smallest minimum of the other differences, as shown in the color bar. One can see its 2D map variation with other maps by focus visualization and it lacks the effect of shallow magma, sedimentary basins, and surface topography. The last map shows that relatively high GD prevails in the western Ethiopian plateau as compared with other regions. The maximum and minimum values of EGM2008 and EIGEN6C4 are almost near each other and they are the best fit model for the study area.

92
ii. Profile curve along or across the escarpment   Figure 7, it can be ascertained that the Sixth release of the GOCE gravity field model by means of the direct approach (GO_CONS_GCF_2_DIR_R6) is not the best model for the study area.  The RMSE values indicate how close the GGM-derived GDs are to the terrestrial GDs. As a result, the EGM2008, EIGEN6C4, and XGM2019e 2159 in Profile BB' (Table 6) are more accurate than the GO CONS GCF 2 DIR R6. Because the terrestrial data has a more dense line across the escarpment than the along one, the three combined model is a better model in Profile BB' than Profile AA'. Because the model lacks ground gravity data, the satellite only model in both profiles has the same level of accuracy.

Radial spectral analysis for depth estimation
Several authors like Bhattacharyya(1978), Gerard and Debeglia (1975), and Spector and Grant(1970) used spectral analysis for depths of gravity and magnetic anomalies. This analysis provides a technique for qualitative studies of large and complex gravity and air-94 borne magnetic data. The logarithm of the radial average of the energy spectrum is plotted vs the radial frequency (Philippe et al., 2006). From this plotted graph, the slopes of the linear are used to separate depth ensembles and provide parameters that are used for the design of numerous filters. As Kivior and Boyd(1998) said the slope of each segment gives information about the depth to the top of an ensemble of gravity bodies.
The value of Lambda (wavelength of the anomalies) which is the inverse of the wave number: the deeper sources produce the largest wavelength. The depth (h) to a statistical ensemble of sources is determined by the ratio between the slopes (s) of the log energy spectrum divided by two. The average spectrum density [ln(E)] for the ground data, EGM08, EIGEN6C4, XGM2019e_2159, and GO_CONS_GCF_2_DIR_R6, is 21.2140.4425, respectively. where: ℎ depth to the source; -depth; power spectrum; -wave number.
The power spectrum estimation was used for periodogram estimate. The periodogram was calculated using the Fast Fourier Transform (FFT) as follows: Truncate the gravity data, and weight the truncated data as the function to be evaluated goes to zero at both ends. The weighing function was normally applied to the first and last 10% of the data. Then, calculate the Discrete Fourier Transform (DFT) and lastly calculate the power spectrum. Finally From Figure 9, the source depth between 28 km and 32 km can be interpreted as the Moho depth beneath the study area; the Moho depth is detected from the long wave length signal of the ground data and GGMs. The estimated Moho depth from the long wave length of the GGMs like EGM08, EIGEN6C4, and XGM2019e_2159 is close to the estimated Moho depth from the ground data, as compared with the GO_CONS_GCF_2_DIR_R6 estimated Moho depth.
The other source of anomaly is found between 7.86 km and 10.25 km, and can be interpreted as the near surface tectonics like uplift of magma (sill and/or dike) as the study area is a seismically and tectonically active region. This layer is not estimated properly using the GO_CONS_GCF_2_DIR_R6 model as it does not have the short wavelength signal, which is the effect of the very near subsurface as compared with the other three GGMs. The last thin source of layer (seen at the very short wave length almost <10 km, which is the effect of the top subsurface density variation) is interpreted as the cover rocks (lava flows or sediments); Afar depression is at the junction of the GOA, the RSR, and the MER Valley, which is largely covered by later volcanics and sediments (Hammond et al., 2011).

CONCLUSIONS
The analysis of the descriptive statistics (minimum, maximum, range, mean, correlation, SD, and RMSE) of the GD and their residual given in Table 1 and 2, respectively, reveals that EGM2008, EIGEN6C4, and XGM2019e_2159 solutions are very close to each other, except for the GO_CONS_GCF_2_DIR_R6. The SD and RMSE of the GD and its residual of the GO_CONS_GCF_2_DIR_R6 are far away from the other three GGMs. This means that the GO_CONS_GCF_2_DIR_R6 is not the best model among the other GGMs. As we compared the four GGMs, the EGM2008 and EIGEN6C4 is the best fit model for the study area as compared with the XGM2019e_2159 and GO_CONS_GCF_2_DIR_R6 model.
The visual interpretation (Figures 3 and 4) of the GD and its residual shows that EGM2008 and EIGEN6C4 have a similar pattern of the GD and its residual, as well as small variation in range (minimum and maximum values of both GDs), as compared with the XGM2019e_2159 and GO_CONS_GCF_2_DIR_R6 model used in this study. 97 The SD of the residual GD is within a range of 20.67-21.27 mGal for the study area but the SD for the GO_CONS_GCF_2_DIR_R6 is the outlier (36.10 mGal). The RMSE of the residual GD is within a range of 25.01-25.65 mGal over the study area but its value for GO_CONS_GCF_2_DIR_R6 is an outlier (39.00 mGal) as compared with others. However, the SD and RMSE for XGM2019e_2159 are larger than the EGM2008 and EIGEN6C4 and smallest from the GO_CONS_GCF_2_DIR_R6. One can conclude from theseobservations that the EGM2008 and EIGEN6C4 model is the best model of the study area.
From the minimum and maximum values in Tables 1 and 2, it is apparent that XGM2019e_2159 and GO_CONS_GCF_2_DIR_R6, strongly, underestimate the GD and its residual. From the visual analysis of the GD and its residual maps (Figures 3 and 4), it can be observed that except for the XGM2019e_2159 and GO_CONS_GCF_2_DIR_R6 (strongly), the values exhibit identical geographical characteristics for the remaining GGMs and their differences from the ground data, but the magnitudes are slightly different.
The comparative results in terms of descriptive statistics of the evaluation of GGM based residual GDresult in the following conclusions in a local scale: the approximation of the residual GD shows that EGM2008 and EIGEN6C4 are almost identical over the study area, and the residuals from the ground data and GD modeling of EGM2008 and EIGEN6C4 are similar and constitute the best model for the study area. XGM2019e_2159 (lower) and GO_CONS_GCF_2_DIR_R6 provide the lowest accuracy in modeling the GD over the study area.
The study shows that the satellite-only model-base is not the best model for the local gravity field model and not good to use for local gravity investigation as compared with the combined gravity model developed from the combination of altimetry, satellite, and topography and ground data, like EGM2008, EIGEN6C4, and XGM2019e_2159, which were also revealed during the literature review conducted by the author. Finally, the qualitative and quantitative analysis results of this study suggest that: due to its better statistics (in terms of correlation, SD, and RMSE), the use of EGM2008 (best) and EIGEN6C4 (better) can be recommended as a feasible GGM for GD modeling tool in geodetic applications in local scales in the study area and also for all over the country. However, the result shows that the three models can be also used for local-scale GD modeling as compared with the XGM2019e_2159 and GO_CONS_GCF_2_DIR_R6; and since the ground data was collected on accessible road by densifying the ground gravity measured data with an equal spatial distribution, the GD can be modeled by GGMs with more accuracy.
The long wavelength of the ground data and GGMs are used for Moho depth estimation on average of 28 km (to the north-eastern part of the study area from the author's literature review) to 32 km (to the western part of the study area from the author's literature review), which is identical with the conclusions of Hammond et al. (2011) and Lavayssière et al. (2018).
The accuracy of GDs derived from four GGMs employing terrestrially gravity data at 1933 points in Afar Region, Ethiopia was statistically evaluated in this study. As previously stated, the accuracy and the resolving power of the data used in the development of a Global Gravity Field Model determine its accuracy and resolution and from this study, we have discovered that EGM2008 and EIGEN-6C4 showed enormous potential to be used as supplements to the terrestrial gravity data. However, these GGMs must be improved in order to improve their accuracy in Afar region, Ethiopia.
The technique and analysis utilized in this work have an impact on the scientific contribution of global geopotential models for geodetic and geophysical applications in the country, such as geoid, quasi-geoid, Moho depth, lithospheric thickness, and other associated parameters.
Because it is used in a variety of geoscience applications, the researcher will continue to apply this methodology and analysis in the future by collecting evenly distributed terrestrial gravity data and conducting a comparison study with GD. Comparison studies can also be made by estimating Moho depth using spectral analysis. A comparison of gravity anomaly analysis and determining geoid using DGPS and Levelling data could be done in the future.