GONGOLA BASIN USING Gongola Basin Geoid Determination using Isostatic Models and Seismic Reflection Data and Geophysical Interpretation

The application of Stokes’ formula to create geoid undulation requires no masses outside the geoid. Usually, a constant density of 2.67 g/cm 3 is used in the determination of the geoid which introduces error in the reduced gravity anomalies (Helmert’s condensation) and consequently the geoid. In this paper, isostatic models of Airy— Heiskanen and Pratt—Hayford were utilized in the determination of the geoid by considering the planar and spherical approximation models the indirect effect of the topographic lateral density variation on the geoid was computed as additive correction for the improvement of the accuracy of the computed geoid. Additional density information deduced from seismic and well log information were considered for the variable density computation. The geopotential geoid undulations were computed from the EGM 2008 model. The residual geoid was obtained by subtracting the local isostatic geoid from the geopential geoid. Geoid and gravity admittance studies were also carried out to complement the results from the residual geoid. The planar and spherical approximation results showed similar characteristics; but a change in magnitude in both models. Our results suggest that the effects of topographic lateral density variations in geoid determination are significant and should be considered in rift basins. The geophysical analysis of the geoid results show that the north-east has positive residual geoid which indicates the presence of high density intrusive igneous rocks, while the south-east has negative residual geoid which indicates the dominant presence of low density sedimentary rocks. The results also show that the radial distribution of the anomalous mass obtained using the geoid/residual geoid anomaly uniquely matched that obtained using the seismic reflection data which inferred the presence of hydrocarbon accumulation in the southeast zone of the project area. The gravity and geoid admittance studies corroborated the residual geoid and seismic reflection results.

Fraser and colleagues have developed a GIS based system to calculate terrain corrections using the real topographical rock density values [Fraser et al., 1998]. The results show that in the Skeena region of British Columbia Canada, the terrain corrections to gravity can change by a few mGals when real topographical density is used. Further, [Pagiatakis et al., 1999] showed that the effect of lateral density variations on the geoid can reach about 10 cm in the Skeena region and several millimeters in New Brunswick, where the terrain is moderate. Also, [Sjöberg, 2004] showed that the total effect on the geoid from lateral density for the deepest lake on Earth (Lake Bajchal) and the highest mountain on Earth (Mt. Everest) can reach up to +1.5 cm and +1.78 cm respectively. These differences are the reason for using improved density model in geoid determination based on Stokes' function [Kuhn, 2003]. Several studies have previously investigated the use of lateral varying topographic density models in gravimetric geoid determination (e. g. [Martinec et al., 1995;Marti, 1997;Kuhtreiber et al., 1998;Pagiatakis et al., 1999;Tziavos, Featherstone, 2000;Haung et al., 2001;Hunegnaw, 2001;Kiamehr, 2006a,b]). In the Nigerian geoid determination [Nwilo et al., 2007] modeled the local geoid of Lagos (Nigeria) based on geometrical interpolation approach. By using the orthometric height and ellipsoidal height, the empirical geoid height were computed. The surface interpolation utilized the kriging approach. Isioye and colleagues utilized a five to eight parameter model to fit the GPS/Leveling to the EGM 2008 model to improve the determination of orthometric height observed from GPS in the study area (Portharcourt, Nigeria) [Isioye et al., 2011]. Ezeigbo and colleagues examined some factors that affect the accuracy of gravimetric geoid determination us-ing mean gravity anomalies over geographically defined grids [Ezeigbo et al., 2007]. Gravity anomaly data obtained from satellite altimetry mission were used in the evaluation of the Stokes' and Vening Meinesz's integral formulae. The result shows that the most significant parameters that affect the accuracy of gravimetric geoid determination are the minimum spherical distance from the computation point, the size and distribution of the observed gravity anomaly data. Okiwelu and colleagues determined geoid undulation for Nigeria using the spherical harmonics expansion employed in the Earth Gravitational Model 2008 (EGM 2008) referenced to the WGS 84 (World Geodetic System 1984) [Okiwelu et al., 2011]. In their study, they found that the highest geoid undulations are centered over the North Central region of Nigeria with relatively lower values (16-20 m) confined to the Nigerian sedimentary basins (Bornu basin and Benue Trough). In this study, the Pratt-Hayford and Airy-Heiskanen isostatic models as well as density data derived from seismic and well log observations (Kolmani River-1 log) will be utilized in the determination of the geoid. The main advantage of using isostatic models in geoid determination is their small indirect effect, together with a smooth field of gravity anomalies [Kuhn, 2003[. Isostatic geoid are also long wavelength in nature and as such are dominated by the signature of deepmantle anomalies (e. g. [Hager, Clayton, 1989]) and/or subducted slabs (e. g. [Ricard et al., 1993]). The introduction of the additive lateral density variation indirect effect and the primary indirect effect terrain effects will improve the accuracy of the computed geoid and hence, the accuracy of the interpreted geophysical structures.
Variations in the height of the geoidal surface are related to density anomalous distributions within the Earth and the geoid undulations help to understand the internal structure of the earth. Fig. 1 [ Lowrie, 2007] shows that positive geoid indicates the presence of high density excess mass, while negative geoid indicate regions of mass deficiency or low density mass deposits. The use of geoid in geophysics has been presented by [Vanicek, Christou, 1994;Featherstone, 1997] in which; the relationship between the geoid and deep Earth mass density anomaly structure, strain and stress fields, tectonic forces, isostatic state of ocean lithosphere, Earth rotation, geophysical prospecting and ocean circulation are discussed. However, the importance of geoid in the determination of the anomalous density distribution has also been recognized by [Kaula, 1967;Chase, 1985;Lambeck, 1988]. Brown in his work [Brown, 1983] recognized the correlation between the geoid and the deep-Earth mass density anomalies; while [Christou et al., 1989] showed its correlation with near surface mass density anomalies. Also correlations between geoid and earth mantle convection is established by [Runcorn, 1967] and with plate tectonic features and seismic tomography by [Silver et al., 1988]. Featherstone used geoid to determine the lateral extent of known geological structures [Featherstone, 1997]. The geoid approach in this research, serves as a complementary way of studying the earth interior and mass density distribution in the Gongola basin. Most previous studies have relied heavily on geophysical approach; for example [Okereke, 1988;Osazuwa et al., 1992;Ugbor, Okeke, 2010;Okiwelu et al., 2010;Okiwelu et al., 2011].

Airy-Heiskanen Model (Planar and spherical approximation)
The model is applied under the following assumptions [Rummel et al., 1998;Kuhn, 2003;Ilk, Witte, 2007] that the isostatic compensation takes place completely and locally, i.e. the compensation mass is directly under the regarded topographic mass which makes the compensation depth variable (Fig. 2).
The root thickness t for planar approximation is obtained by [Heiskanen, Moritz, 1967;Grant, West, 1987] as: where t р -root thickness (planar), ρ cr -crustal density (2.67 g/cm 3 ), Δρ -the variable density contrast between the lower crust and the upper mantle, H -topographic height. The root thickness t sp for spherical approximation is obtained by [Rummel et al., 1998;Ilk, Witte, 2007] where R -radius of the earth (6371 km), T -normal thickness.
Usually it is applied assuming: -a constant compensation depth (D) of 100 km at which the hydrostatic equilibrium is achieved; -that due to constant compensation depth, the condition of equilibrium leads to a laterally variable mass density. Fig. 3 illustrates the Pratt-Hayford model.
The variable density for the planar model is expressed as: ρ is given as [Rummel et al., 1998;Ilk, Witte, 2007]: sp L ρ -density variation on land (spherical), ρ t -weathered tertiary density value, R -radius of the Earth.

Spherical Harmonics Representation of Geoid
Spherical harmonics are often used to approximate the shape of the geoid. The geoid undulations are evaluated from spherical harmonic coefficients at the surface of the ellipsoid, not taking into account the difference between height anomalies and geoid undulations [Heiskanen, Moritz, 1967, p. 325] nor the effect on the geoid of the downward continuation of gravity from the surface [Sjöberg, 1998a;Kaban et al., 2004]. Geoid undulation (N) can be expressed in spherical harmonics as [Rapp et al., 1991]: The corresponding free air gravity anomaly can be obtained from the anomalous potential as [Heiskanen, Moritz, 1967, p. 97].
Δg F -free air anomaly. The Earth Gravitational Model EGM 2008 [Palvis et al., 2008] is complete to spherical harmonic degree and order 2159 and contains additional coefficients extending to degree 2190 and order 2159. The harmonic degree implies that short wavelength anomalous features can be studied which is relevant in the determination of the near surface mineral depth and its mass density distribution [Featherstone, 1997]. The free air anomaly derived from the spherical harmonic model is stated as follows:

Isostatic Geoid Anomalies 2.4.1 Co-Geoid and Geoid undulations
The geoid undulation N is the vertical separation between the geoid and the ellipsoid. Fig. 4 shows the geometric separate of the geoid and the co-geoid. In gravimetric geoid determination, the computed geoid does not match with the actual geoid. The difference between the actual geoid and the computed geoid is called the indirect effect.
The formula for computing co-geoid undulation from a 3D density model Δρ(x, y, z) is expressed as [Turcotte, Schubert, 1982]: N c -co-geoid, γ 0 -normal gravity, G -universal gravitational constant, p m -mantle density, ρ crcrustal density. For Airy-Heiskanen compensation with depth of the normal T, surface topography H and the root t. From Eq. 1, taking downward continuation, the co-geoidal undulation is given for planar approximation [Crovetto et al., 2008] as: In the Airy-Heiskanen formula, we assume a perfect isostatic balance using T=30 km, density 2.67 and density contrast ρ m −ρ cr =0.6 g/cm 3 . For spherical approximation with respect to Eq. 2, it is derived as: с m -upper mantle density (3.27 g/cm 3 ).
In the use of Pratt-Hayford model for planar earth (Eq. 3), with depth compensation D and normal density ρ cr , H -elevation, the co-geoid undulation associated with positive topography is derived from Eq. 7 as: For spherical approximation with respect to Eq. 4, it is derived from Eq. 7 as After computing the co-geoid undulation, the geoid undulation N is computed by adding a number of additive corrections (which sum up as the indirect effects) to the co-geoid N -geoid undulation, δN -total indirect effect, which must use the same density model as the gravity anomalies, δN Δp -lateral isostatic density anomaly indirect effect on the geoid, δN DWC -indirect effect due to downward continuation. The contribution of this effect cancels in their sum.

Primary Indirect Terrain Effect
The primary indirect topographical effect (PITE) is the separation between the geoid and the co-geoid caused by the condensation of the topography and the atmosphere δV -change of potential at the geoid which depends on the reduction method used.

The Effect of Lateral Density Variation on the Geoid
Sjöberg [2004] showed that the total effect of the geoid due to the lateral density anomaly could be represented as a simple correction proportional to the lateral density anomaly and the elevation of the computation point square. The combined topographic effect on the geoid [Sjöberg, 2001;Kiahmehr, 2006b] including the zero and first degree terms is well approximated by where G -gravitational constant, ρ=ρ(θ, λ) is the laterally variable topographic density at the colatitude (θ), and longitude (λ), H -height, γ 0normal gravity. The lateral density variations indirect effect on geoid determination was computed using Eq. 14. If the density of the topography at the computational point is where ρ cr is the standard density (2.67 g/cm 3 ), Δρ(ϕ, λ) is the lateral density anomaly with respect to the standard density. The total effect of Δρ on the geoid undulation becomes For the Airy-Heiskanen model for geoid undulation computation, the lateral density variation indirect effect on the geoid for planar approximation with respect to Eq. 1 is given as Based on the concept above, the indirect effect in spherical approximation with respect to Eq. 2, is derived as From the Pratt-Hayford model for geoid undulation computation, the lateral isostatic density variation indirect effect for planar approximation on the geoid with respect to Eq. 3, is given as [Kiamehr, 2006b] The indirect effect in spherical approximation with respect to (Eq. 4), is derived as Fig. 4. The Relationship between geoid and co-geoid [Kuhn, 2003].

Geoid Undulations computations from isostatic models
From the above formulations for the co-geoid, primary indirect terrain effect and lateral density variation indirect effect, the final expression for geoid undulation obtained from the isostatic models are:

Residual Geoid Undulations
In order to reveal the short wavelength geoidal features which are assumed to reflect crustal and lithospheric structures, the long wavelength component of the geoid assumed to originate in the mantle is removed by the process termed detrending [Featherstone, 1997]. Residual geoid undulation is obtained as the difference between the geoid obtained from EGM 2008 model and the geoid undulations obtained from the isostatic models: where ΔN I -residual undulation for the various isostatic models, N EGM 2008 -geoid undulation for EGM 2008, N I -geoid undulation for the various isostatic models.

Geoid and Gravity Admittance Evaluation
Gravity and geoid anomalies reflect lateral heterogeneities in the earth's density structure. Because such anomalies often correlate with topography, it became a standard approach to use the relationship between bathymetry and the gravity or geoid to gain information about the subsurface density structure and the style of isostatic compensation of the topographic load assuming that compensation occurs on a regional basis. In the following, the geoid admittance (Geoid to Topographic Ratio) as a spectral function is calculated by [Keifer, Hager, 1991]: A geoid (λ) -geoid admittance, N(λ) -geoid undulation, H(λ) -topography.
In relating the geoid undulations to free air gravity anomalies in the spectral domain, the following expression is used [Keifer, Hager, 1991] where A gravity (λ) -gravity admittance, Δg FA (λ)free air anomaly.
In the case of the Airy-Heiskanen model, the topographic height is considered in terms of surface topography and mantle (dynamic topography) to determine the deep earth mass density anomalies [Bowin, 1983] and near-surface mass density distribution [Christou et al., 1989].

Methodology
The following data were utilized in the gravimetric geoid determination.
1. Digital terrain model (DTM) for modeling the shape of topography. 2. Digital Density Model (DDM) for modeling the spatial distribution of density in the topography and deeper masses. The information from seismic and well log observations were used to derive the density model. 3. Isostatic models for analytically modeling the Earth's outer masses. 4. EGM 2008 geoid undulation data. 5. Residual geoid obtained as difference between EGM 2008 and the local isostatic geoid. Fig. 5 shows the flowchart for the geoid modeling and its geophysical analysis.

The Study Area and Gravity Data Acquisition
Gongola basin of Northern Nigeria is one part of a series of Cretaceous and later rift basins in Central and West Africa whose origin is related to the opening of the South Atlantic [Obaje et al., 2006]. Many authors have noted that the Benue rift (which includes Gongola basin) have many features in common with other intra-continental East African rift such as Baikal rift and the Rio Graade rift, for example [Logatchev, 1993;Shemang et al., 2001;Ugbor et al., 2010]. These rift systems are associated with volcanism and regional up-lift. The basin contains thick sediment accumulations (mainly Cretaceous) in excess of 7 km in the North-Eastern part of the study area; deposited under varying environments. Seismic studies in the area have led the insight into the structure of the crust and mineral viability of the basin.
The gravitational data for this investigation is a set of Bouguer gravity and isostatic residual gravity anomalies observed at 1813 shot points with station interval of 500 m within Gongola basin.

Results
The following results shown in Fig. 6-10 il- lustrate the geoid and geophysical findings within the basin.

Geoid and Residual Geoid Anomaly Maps
The spatial behavior of the geoid undulation for P-H and A-H and that of the EGM 2008 (Fig. 6) have similar characteristics with respect to the density variation in west-east direction; but differ in magnitude. However, the results from the P-H spherical approximations are better than the other results in the planar and spherical approximations in terms of producing the minimum sum of squares of the residual geoid values (Fig. 7).
The isostatic geoid undulations are found to have long wavelength features between 4000<λ<8600 m. These features of the geoid are due to density variations in the lower mantle and resulting deformations of the core mantle boundary and other boundaries in the mantle [Richards, Hager, 1984;Okiwelu et al., 2011]. The residual geoid undulations for the P-H isostatic models range between -7 m to +11 m in the spherical approximation model. The negative residual geoid undulation values correspond with the region of less dense intrusive rocks in the SE, while the positive residual geoid values correspond to regions of high density subsurface intrusive rock deposits zones of the project area. The spatial behavior of the P-H isostatic residual gravity anomaly, the P-H geoid undulation and the P-H residual geoid (as shown in Fig. 7, a, b) correlates and this helps in the definition of the basin's subsurface density distribution. The large bias in the residual geoid undulation is based on the use of spherical harmonics coefficients in the EGM 2008 geoid model. Gravity disturbances [Heiskanen, Moritz, 1967], rather than gravity anomalies, are computed from the spherical harmonics so as to account for the masses between the reference ellipsoid and geoid [Kaban et al., 2004]. Geoid undulations are evaluated from the spherical harmonic coefficients at the surface of the ellipsoid, not taking into account the difference between height anomalies and geoid undulations [Heiskanen, Moritz, 1967] nor the effect on the geoid of the downward continuation of gravity from the surface [Sjöberg, 1998a]. Correlation exists in the geoid undulation between basement complex zones and the EGM 2008 geoid undulations. In comparison with the geology, there is a substantial agreement with respect to the sedimentary zones with large negative residual geoid undulations which is due to lateral density anomaly. This is due to the large density contrast which has a direct influence on the computed geoid. The area with high sedimentary thickness between 4 km and 7 km in the North-East has high density igneous intrusive rocks and as such, depicts positive residual geoid undulation values as shown in Fig. 7, c. The de- Fig. 6. Pratt-Hayford, Airy-Heiskanen and EGM 2008 geoid undulations. crease in geoid undulations from west to east within the Gongola sedimentary basin shows the presence of depression in the NE and SE zones which is consistent with the evolution of the trough. The high geoid undulation transiting between the NE and SE (Fig. 6, a, b) shows the presence of a ridge Fig. 7. Correlation of the isostatic residual gravity, geoid and residual geoid maps.  [King, 1950;Cratchley, Jones, 1965;Wright, 1976;Benkhelil, 1982;Whiteman, 1982]. It is a rift basin with plate dilation leading to the opening of the Gulf of Guinea [Benkhelil, 1989;Fairhead, Binks, 1991]. Benkhelil [1989] suggested that the evolution trough could also be as a result of tension resulting in a rift or wrench related fault basin. Mesozoic to Cenezoic magmatism has accompanied the evolution of the tectonic rift as it is scattered all over and throughout in the trough [Coulon et al., 1996;Abubakar et al., 2010]. A magmatic old rift was also suggested for the Gongola basin by [Shemang et al., 2001] while [Abubakar et al., 2010] suggested the evolution as a combination of mantle upwelling or rise of a mantle plume which resulted in crustal stretching and thinning and the emplacement of basic igneous material within the basement and sediment which resulted in rifting. The reflection seismic time/depth structural maps also define the SE zone as a zone with favourable structural geometry for hydrocarbon accumulation. Fig. 8 shows the correlation between the Pratt-Hayford isostatic residual gravity anomaly, the seismic depth map (at prospect horizon) and the Pratt-Hayford residual geoid undulation. The reflection seismic depth map corroborates the Pratt-Hayford geoid and residual geoid structural pattern and hence, the favourable locations for hydrocarbon accumulation (oval shapes in the maps) shown in the residual geoid are justified.

Analysis of Geoid and Gravity Admittance Maps
Short wavelength geoid provides information about the near surface features while long wavelength geoid with low degree (n=6, 7) provides information about the mass anomalies in the crust mantle zone [Bowin, 1985]. Their isostatic model gravity/geoid admittance were used to determine the mass anomalies in the crust mantle zone. The gravity admittance map (Fig. 9, a) from the Pratt Hayford isostatic model shows a uniform mass density distribution that corresponds with the result obtained from the mass density structure of EGM 2008 gravity admittance (Fig. 9, b). The mass anomalies from EGM 2008 geoid admittance corroborate with the gravity admittance structures from the isostatic model. The mass anomalies vary between 60 and 82 m/km in the northeast while it varies between 38 m/km and 60 m/km in the southeast of the project area. Consequently, it could be inferred that the northeast zone has high density mass deposits while the southeast has low density mass deposits. The northeast is therefore a favourable location for solid minerals while the southeast is a favourable location for hydrocarbon deposits. The gravity/geoid admittance structures corrobo-rate the findings established in the geoid/residual geoid and seismic reflection studies.

Conclusion
Gongola basin geoid is primarily affected by the lateral density variation from the crust/mantle discontinuities. Based on this, the lateral density variation indirectly affected the accuracy of the geoid significantly by about 3 cm.
The geoid and residual geoid undulations have shown in this research to corroborate the findings of the reflection seismic data in the determination of the lateral extent of known geological structures. This has also shown that the geoid can be used as a complementary method in the definition of the location and radial distribution of the geological structures for mineral and hydrocarbon exploration. The application of Stokes' formula to create geoid undulation requires no masses outside the geoid. Usually, a constant density of 2.67 g/cm 3 is used in the determination of the geoid which introduces error in the reduced gravity anomalies (Helmert's condensation) and consequently the geoid. In this paper, isostatic models of Airy-Heiskanen and Pratt-Hayford were utilized in the determination of the geoid by considering the planar and spherical approximation models the indirect effect of the topographic lateral density variation on the geoid was computed as additive correction for the improvement of the accuracy of the computed geoid. Additional density information deduced from seismic and well log information were considered for the variable density computation. The geopotential geoid undulations were computed from the EGM 2008 model. The residual geoid was obtained by subtracting the local isostatic geoid from the geopential geoid. Geoid and gravity admittance studies were also carried out to complement the results from the residual geoid.
The planar and spherical approximation results showed similar characteristics; but a change in magnitude in both models. Our results suggest that the effects of topographic lateral density variations in geoid determination are significant and should be considered in rift basins. The geophysical analysis of the geoid results show that the north-east has positive residual geoid which indicates the presence of high density intrusive igneous rocks, while the south-east has negative residual geoid which indicates the dominant presence of low density sedimentary rocks. The results also show that the radial distribution of the anomalous mass obtained using the geoid/residual geoid anomaly uniquely matched that obtained using the seismic reflection data which inferred the presence of hydrocarbon accumulation in the southeast zone of the project area. The gravity and geoid admittance studies corroborated the residual geoid and seismic reflection results.