Structural Mapping Inferred from Gravity Data to Image the Upper Lithospheric Structures and its Hydrocarbon Implication from Gwandu Formation NW, Nigeria

: : Presently a digitized composite satellite gravity data covering the Gwandu formation in Sokoto Basin were acquired and processed with a view to interpret the Bouguer anomalies as well to equally image the upper lithospheric structures beneath the Study area and its environs. The research work was aim to study the structural settings of crustal movement in the Gwandu formation. A least-square fitting polynomial surface of a third-degree order was applied in separating regional and residual gravity components from the Bouguer anomaly. The attributed low gravity sedimentary infill from the residual anomalies were tectonically trends NE -to- SW about the vicinities of Tambuwal, Goronyo, Gada and Argungu, Kolmalo and about Yauri, Koko and Jega, Kamba as well as Bagudo. Data enhancement techniques such as first vertical derivative, total horizontal derivative (THDR), analytic signal, spectral depth analysis, and the standard Euler deconvolution (SED) were applied to enhance deep-seated structures. Results from the Spectral Analysis revealed that the average thickness of the sediments varies from 1.679 km to 4.181 km, outsized enough for hydrocarbon prospect. The derivative maps revealed parallel to sub-parallel trending NW – to- SE, E -to- W fracture zones within the sedimentary infill underlying the study area, coinciding with the cretaceous zones. Hence, the identified lineaments (faults or lithologic contacts) and structures in the area can be attributed to the tectonic setting of the area and probable migratory routes for hydrocarbon migration. More detailed ground gravity and seismic studies may lead to discoveries of structural or stratigraphic traps.


INTRODUCTION
Commercial deposits of oil and gas have discovered and are being produced from contiguous structurally and stratigraphically rifted basins of Niger Republic, Chad Republic, and Sudan. The Gwandu formation is part of the formation in Sokoto Basin in Kebbi States of Nigeria has drawn the attention of geoscientists and stake holders in the petroleum industry largely because of intensified petroleum exploration efforts in the inland basins and the fact that the Nigeria's Department of Petroleum Resources has identified the Sokoto Basin as one the potential hydrocarbon-bearing inland basins. There is equally another compelling reason for the exploration efforts in the Sokoto Basin (Adamu & Likkason, 2022;Adamu et al., 2021;Kogbe, 1972Kogbe, , 1979Kogbe, , 1981Obaje et al., 2009Obaje et al., , 2020. Potential field gravity method has proved very effective for providing useful information known to guide various exploration campaigns, be it regional studies (Olawale et al., 2020;Moumoni et al., 2016;Zahra & Oweis, 2016), economic mineral or oil and gas exploration (Ugbor et al., 2020). In defining basin's tectonic framework from Gwandu formation in the southeastern part of Sokoto basin for delineation of outlining favourable region for hydrocarbon prospect, an important reconnaissance mapping of structural features has been generated by the analyses of satellite gravity data. Thus, gravity survey is the primary method in geophysical exploration as a regional structural mapping tool (Reid, 2003;Reid et al., 1990;Reid et al., 2003). The success of gravity survey be determined by the existence of a significant density contrast between altered rocks or structures and their host rocks (Reid, 2003). Therefore, this present study set out to take advantage from the acquired satellite gravity data to image the upper lithospheric structures based on structural mapping over the Gwandu formation and its adjoining environs in Sokoto Basin, Nigeria. Consequently, the goal and objectives of this research was to explore the use of satellite gravity data in prospects mapping towards the identification of prolific zones, depth determination and linear structures within the area of study that may be potential targets for hydrocarbon trap or related mineralization.

STUDY AREA: LOCATION AND GEOLOGY
The study area is fall within Gwandu formation in the Sokoto basin and it spans about 26,000 km 2 , situated between latitudes 04° 00ʹ N and 05° 30ʹ N and longitudes 12° 00ʹ E and 13°00ʹ E with an average elevation of 218 m above mean sea level (Fig. 1). There are well linked through foot path and fairly motorable road. These roads are mainly used by rural dwellers to access their Houses and settlement; it also provides linkage to some of the survey traverse.
The area is generally underlain by sedimentary rocks of the Gwandu formation. It contains a number of prominent ridges and groups of flat-topped, steep-side hills capped by ironstone. Other hills covered with ironstone debris occur in all stages of disintegration, rising out of the sandy plain over which the products of erosion have been distributed. Rock exposures are rare on the plain, buy numerous on the hillsides, however, they are small and obscured by rain-wash and ironstone scree (Fig. 2). The non-marine origin of Gwandu formation is certain, and the sediments can be correctly attributed to a continental environment, or more precisely to a lacustrine environment (Kogbe, 1976). The best outcrops of the Gwandu formation occur around Birnin Kebbi and Argungu. The sediments consist of massive while clays interbedded with coarse and medium-grained red sandstones and mudstones with occasional peat bands. The type section proposed for the formation by Kogbe (1979) shows the typical lithologic characteristic of the formation. Beneath the lateritic capping is a hard-ferruginous sandstone layer which is easily eroded into network of gullies. These are underlain by red sandy clays and white massive mudstones which are invariably, stained pale brown or pink. The mudstone with sandstone intercalations extends monotonously throughout the sections. Similar sections of the Gwandu formation occur on the slopes of the Gwandu outliers within the Kalambaina Formation on the outskirts of Sokoto Township near the cement factory. The sands at the surface are quite red in color, often showing color handing and poor stratification. The mudstones often show a nodular structure with nodules suggestive of local turbulence in the depositional environment. By correlation with palynomorphs from tropical Tertiary deposits earlier mentioned, the age of the Gwandu Formation was tentatively put as Eocene-Miocene (Kogbe, 1979).

METHOD
In this present study, satellite gravity data were acquired, analysed, and processed for the purposes of image the lithospheric structures beneath the Sokoto basin and its surroundings. The gravity method depends on the measurements of variations in the gravity field caused by horizontal variations of density in the subsurface (Ismail et al., 2001;Askari, 2014;Mammo, 2004). It is the essential method in a number of specific geological studies, as in mapping near surface voids, quantitative studies of metallic ore bodies, characterizing salt structures, and monitoring changes of fluid/gas content in volcanoes. The gravity method has also been used in regional characterization of the earth to determine the structures of the crust, identifying potentially favorable regions for resource exploration, and developing conceptual exploration models (Hinze et al., 2013;Tiberi et al., 2005).  Kogbe, 1976) Satellite Gravity data The Bouguer gravity anomaly data ( Fig. 6) were obtained from the International Gravity Bureau, which maintains the WGM 2012 (Balmino et al., 2011;Zelalem et al., 2018). The WGM 2012 Bouguer gravity anomalies were calculated using the Earth Gravity model, EGM 2008, satellite gravity and the Topography 1 arc-minute (ETOP1) database for topography. This is to provide harmonic coefficients for a spherical harmonic expansion up to an order of 10, 800 to produce the Earth's Topography derived gravity model at 1' x 1' resolution (ETOPG1) model (Balmino et al., 2011;Zelalem et al., 2018).

Edge detection techniques First vertical derivative (FVD)
This enhancement technique applied on the Bouguer anomaly map is to sharpen up anomalies and allow clearer imaging of the causative structures vertically. The first vertical derivative filtering action is applied to enhance the high frequency components at the expense of low frequency components which are deeply seated. The first vertical derivative can also be applied to the gravity data by using the amplitude spectra of the gravity field (Gunn et al., 1997;Reid, et al., 1990). This transformation magnifies the short wavelength features (Mammo, 2010), which could reflect the residual anomalies. The gridded gravity anomaly data in (Fig. 3) can be expressed as a function in Cartesian co-ordinate system, i.e. F = f(x, y, z). The vertical derivative which shows the change of field with depth (z) is expressed as: The residual anomalies obtained through subtraction of upward continued 5 km and regional of order one from Bouguer anomaly are compared against vertical derivative anomaly map to see whether the shallow earth anomalies are adequately picked or not.

Total Horizontal derivative (THDR)
The horizontal gradient which can be applied on the Bouguer anomaly map and this turns out to be a simplest approach to estimate contact locations of the bodies at depths. The biggest advantage of the horizontal gradient method is its low sensitivity to the noise in the data because it only requires calculations of the two first order horizontal derivatives of the field (Phillips et al., 2007).

Total gradient (analytic signal) method
The analytic signal method (Nabighian, 1972;Roest et al., 1992) assumes that the sources are isolated dipping contacts separating thick geologic units. Peaks in the analytic signal amplitude, which is derived from the first horizontal and vertical derivatives of the observed gravity field, are used to locate the contacts and estimate their strike directions. If the geologic units are not thick, the depth estimates from the analytic signal method will be too shallow. The analytic signal method is moderately sensitive to noise in the data and interference effects between nearby sources. This transformation is often useful at low latitudes (e.g. this study area), particularly for the magnetic data because of the inherent problems and also with the Bouguer anomaly. The analytic signal (AS) is the square root of the sum of the squares of the derivatives in the x, y and z directions.
Where ( ) 2 , ( ) 2 ( ) 2 are the squares of the first derivatives of the total gravity field in the x, y and z-directions respectively, it is very useful in locating the edges of gravity source bodies (Geosoft Oasis-montaj, 2007). The advantage of using AS technique to determine gravity dense parameters from Bouguer anomalies is the invariance of direction (inclination). The success of the AS method results comes from the fact that the location and depth of gravity sources are found with only a few assumptions about the nature of the source body, which is usually assumed as a source (for example, step, contact, horizontal cylinder or dike). For these geological models, the shape of the amplitude of the AS is a bell-shaped symmetric function located directly above the source body. We compute ASI as: ASI is a measure of the depth, D to the source and can be estimated as in (Nabighian, 1972;Roest et al., 1952); Where fv is the first vertical derivative of the total magnetic anomaly or gravity anomaly field, and D is the depth to the source body, N is known as a structural index and is related to the geometry of the gravity or magnetic source. For example, N = 4 for sphere, N = 3 for pipe, N = 2 for thin dike and N = 1 for magnetic contact (Reid et al., 1990;Reid, 2003).

Depth to basement estimation (Sedimentary thickness) Upward continuation technique
It is a mathematical technique used to separate the anomaly of the deeper geology from shallower geology (Jacobsen, 1987). It is a transformation of gravity anomaly computed at a point, Q (x0; y0; z0 = 0) on the mean sea level to a point P (x0; y0; z0 = -h) on some higher flat surface upward continued to, z = h < 0 (Fig. 3). The gravitational attraction per unit mass of an anomalous source Body of mass, dm; at mean see level, O (x0; y0; z0 = 0) at distance R from the source location point Q (x, y, and z) is computed by: The vertical component, dgz, gives the gravity anomaly, Δg0, of the source at the mean sea level surface point as: Thus, the gravity anomaly, Δgp, of the anomalous source at the upward continued surface point P is: Eq. (6) represents the gravitational attraction of dm at a height, 'h', above mean sea level surface. This filter attenuates the short wavelength anomaly components while enhancing the long wavelength ones. This equation therefore represents an algorithm for developing an upward Continuation filter.

The Standard Euler deconvolution
Euler deconvolution technique is used to estimate the source depth locations of the gravity or magnetic signature in a region; and applied to profile or gridded data to solve the Euler's homogeneity equation (Thompson, 1982;Fitz-Gerald et al., 2004) of the form: Where f is the observed field at location (x, y, z) and B is the base level of the field [regional value at the point (X, X0, Y, Y0, Y, Y0)] and N is the structural index or degree of homogeneity (Reid et al., 1990). Source depth has been estimated based on the Euler deconvolution technique applied to the gridded data set by solving the Euler's homogeneity equation. This depth estimation depends on the type of structural index N which varies from -1 to 2. In the case of gravity data set, N = -1 denotes the contact type bodies, N = 0 for thin sheet edge and thin sill dyke, N = 1 for cylinder, thin bed and fault, and N = 2 for point or spherical bodies. In this work, we have used N = 0, which is theoretically not accurate but provides approximation of contact type bodies for depth calculation using the first vertical derivative of gravity data. The Bouguer gravity anomaly across the south-western part to northwestern part shows large deviation from local stability. The mass deficiency characterizes the sedimentary infill of the Sokoto basin. It is understood that the Sokoto basin is supposed to be overcompensated by as much as 68 mGal (Fig. 6). This combination of mass deficiency over the Sokoto basin and the excess mass over some escarpment features is due to isostatic occurrences within the strong Pan-African lithosphere. The regional correction is applied to Bouguer gravity anomaly (Reid et al., 1990).

Analysis of Spectrum technique
Spectral analysis estimates the depths to sources by making use of a statistical model which assumes the ground as an ensemble of series of independent blocks (Spector & Grant, 1970). Various sources at different depths (deep, shallow and noise components) contribute to the cumulative response of an ensemble of sources across a given spatial frequency spectrum. The Bouguer anomaly grid was subdivided into nine (9) overlapping blocks. Each block has a dimension of about 30 by 30 km (Fig. 4). The radially averaged power spectrum of the subdivided blocks was calculated and plotted. The linear segment on the spectrum is directly proportional to the depth (Naidu, 1968;Spector & Grant, 1970). The depth, h of an ensemble of sources is given by: where s is the slope of the power spectrum and h is the depth. The higher wavenumber portion corresponds to the shallow sources while the low wavenumber end of the spectrum corresponds to the deeper sources. The high wavenumber end of the spectrum is a consequence of noise present in the data and terrain clearance effect (Tadjou et al., 2009). The anomalous value of the gravity field at a point is the sum of the gravity effects of widespread and deep-seated mass distributions and smaller, localized mass distributions near the observation point. The interpretation of Bouguer gravity anomalies often involves isolating anomalies of interest (residual gravity anomalies) (Mickus et al., 1991;Linsser, 1967). The observed Bouguer gravity anomaly field consists of two components: a regional and residual gravity anomaly field. Regional/residual separation process was applied to gravity data-set in order to estimate the amplitude of the regional background. Upward continuation can be used to separate a regional gravity anomaly resulting from deep sources from the observed gravity (Kebede et al., 2020;Mammo, 2004;Linsser, 1967). This is an operation that shifts the data by a constant height level above the surface of the earth (or the plane of measurements). It is used to estimate the large scale or regional (low frequency or long wave length) trends of the data. Since the target depth is the sedimentary infill which is approximately undulating 1 km -3 km, the data is upward continued at 5 km to remove the short wavelength anomalies. Jacobsen (1987) demonstrated that if a potential field is upward continued to a certain height z, then it is possible to focus on sources situated at a depth greater than z = 2 (Mammo, 2004). The residual gravity anomaly (Fig. 5) is computed by removing the regional (Fig. 6) effect from the complete Bouguer gravity anomaly (Fig. 7).  The Bouguer anomaly map (Fig. 7) indicates lateral changes in the earth's gravity field and has a maximum anomaly value of about -13.9 mGal at the northern and northwestern parts and minimum anomaly value of about -66.3 mGal at the southern, southeastern and at the southwestern parts of the study area. High gravity fewer negative anomalies are concentrated in the southern, southwestern and nearly in the northwestern parts of the area of study. In general, these gravity anomaly zones may be due to the presence of deep or thicker sediments rock. Low gravity anomalies are concentrated in the northeastern, and in the northwestern parts. Generally, the negative gravity anomalies are probably due to presence of thick sedimentary sections in these parts of the study area (Kogbe, 1972;1979;Obaje et al., 2009;.

Correlation between Free-air versus Bouguer anomaly of the Study Area
The Linear Regression Analysis of Free-air Gravity against Bouguer anomaly were carried out from the acquired gravity data. The best density fit for Bouguer gravity reduction was obtained by dividing the slope factor of the equation Y= 0.08537618 * X-27.836420 with the Bouguer slab constant of 2πG (approximately x = 0.042, where G is the gravitational constant known as 6.67 x 10 -11 Kgm 2 /s 2 ). This approximately results are 2.06 gr/cc for the value of the best Bouguer reduction density, and was subsequently applied for recalculating the Bouguer slab corrections and the corrections for the regional terrain effects to produce the Extended Bouguer Gravity map shown in (Fig 7). By using the elastic thickness of the crust of 20 km and the assumed densities of the topography 2.23 gr/cc (Fig 8), the crust 2.67 gr/cc and the mantle 3.07 gr/cc, suggest that the image of upper lithospheric fall within the Northeastern Region is, at present, compensated only by about 20 percent, implying the depth of the crustal root of less than 31 km (Kearey & Vine, 1990). A total compensation may be achieved when the mean value of the free-air gravity is very close to 0 mGal at which, the depth of the compensating crustal root reaches 4 km below sea level (Keary & Vine, 1990). Free-air gravity maps are sensitive to elevation and tend to masked the deep geological bodies. The analysis of deeply seated geological features was carried out using the Bouguer gravity in which, effects of surface geology which masked the target of interest at depth were removed. First Vertical derivative (FVD) The first vertical derivative (FVD) map ( Fig. 9) revealed sharpened edges of anomalies and shallow features, trending NW-SE. Sedimentary infill having high or low gravity anomalies are evenly distributed within the map. Similar to the residual map, the FVD map reveals shallow features with less noise and which are prominent at the northern and southern part of the study area. Structural variations are similar to the structures observed on the residual map and dissimilar to the smoothened and broadened structures discernible on the regional map (Nabighian, 1984). Horizontal derivative transformation was utilized to decipher the contact positions of anomalous bodies at various depths. It employs the use of the first-order horizontal derivatives of the regional field. With this, a contour map displaying the depths to causative sources which was generated, analyzed and used to infer the hydrocarbon potential of the study terrain. The Horizontal derivative (Fig.10) allows estimation of depth, in kilometers, of a group of geologic dense body structures that vary in width, depth, and thickness, it reveals the maximum depth as occurring at the eastern end of the study area. The HDR map of the study area shows values ranging 0.145 mGal/km to 0.219 mGal/km (Fig. 10). The maxima of HDR values are mainly located in the northeastern part of the study area where prominent ridges of Gwandu, Kalambaina, Gamba, and Gundumi formation are located and in southern part (Kogbe, 1972;1979;1981). The maxima in HDR also represent the geological contacts, the interpreted contacts from HDR match known geological contacts. Figure 10. Horizontal derivative of gravity data of (total gradient) (mGal) Total gradient (analytic signal) The analytic signal (AS) map is independent of the strike and dip of perturbing gravity anomalies as well as the direction of interest and represents the envelope of both the horizontal and vertical derivatives over all possible directions of the earth's gravitational field. The AS transformation ensures that anomalies are positioned directly above their respective causative bodies. It does that by differentiating the total regional field gradient at each measurement point in three perpendicular directions. The mathematical basis of this transformation technique is detailed in the works of Nabighian (1972Nabighian ( , 1984 and Roest et al., (1992). High analytical amplitudes were observed at the southern and eastern sections of the map (Fig.11). Locations with such large analytical amplitudes often coincide with areas of substantial limestone deposits. Presently, red molted massive clay deposits with sandstone intercalation and ironstone are being extracted from sedimentary geological formations in Gwandu, Dange, and Dukamaje (Kogbe, 1972;1979;1981;Obaje et al., 2020). On the other hand, low analytical signals are seen mainly at the western portion of the map, as well as at the south-western section (Fig.11). These low values are probably anomaly signatures from sedimentary in-fills. Figure 11. Analytic signal (AS) map of the study area.

Depth estimation The Standard Euler deconvolution (SED)
In this investigation, Euler deconvolution technique has been carried out on the gravity data. However, the available Bouguer gravity anomaly map was subjected to this method using the structural indices (0, -1, 0.5 and 1.0) and window size equal to 10, as well as tolerance value of 7. An Euler map (Fig. 12) was derived which shows clustering of circles in linear shape indicating the nature of probable contacts between the rock units. The linear clustering circles are suggestive of faults and or contacts with depth values ranging between 500 m and 5000 m for contact, and 500 m to 8000 m for dyke for all the lineaments. This gives an insight of approximate depth range of all the lineaments/ fractures. These solutions are trending in NW-SE, ENE-WSW, NE-SW, E-W and NNW-SSE directions (Fig.12). The use of Euler deconvolution has emerged as a powerful tool for the direct determination of depth and probable source geometry in gravity data interpretations. The method can locate or outline the confined sources, dykes and contacts with remarkable accuracy. Euler deconvolution has been widely used in the automatic interpretation, because it requires no prior knowledge of the source direction and assumes no particular interpretation model (Thompson, 1982;Reid et al., 1990Reid et al., , 2003Stavrev, 1997). In this work, we propose to estimate both the source location and N using the Euler deconvolution, assuming nonlinear background.
We approximate the regional field using a rational function, in which both the numerator and denominator are linear.
Step-fault locations and their depths deduced from Bouguer gravity anomaly using Euler deconvolution for (a) Contact_0-7-10 and (b) Dyke_0.5-7-10 Spectral depth analysis Figures 13 (a-i) show the radial power spectrum from nine (9) of the subdivided blocks (Block 1 -Block 9). The high and low wavenumber anomalies reflect sources that are close to the surface and deep seated, respectively. The low wavenumber portion (layer1) which corresponds to the deep-seated sources comprises the layer of interest. The graph shows two distinct cut-off wave number (cyc/km) of the power spectrum trends, along with the indicated depth to the corresponding infill boundaries. The first layer consists of deep spectrum values with the best line fit as di = 3.908 km, 5.352 km, 4.123 km, 3.07 km, 2.894 km, 5.016 km, 5.366 km, 3.686 km and 4.247 km corresponds to the depth of top of intrusive structures. The second layer indicates the shallow depth values of the spectrum with the best line fit di = 1.540 km, 2.564 km, 2.413 km, 0.770 km, 0.769 km, 0.1696 km, 1.986 km, 1.629 km and 1.742 km with correspond to the depth of metalliferous ore bodies ( Fig. 13 (a-i)). These depths estimate values were subsequently used for constraint to the correlation of interpretation of Euler deconvolution. The average depths estimated from the subdivided blocks were compiled to generate 2D and 3D depth maps (Fig. 14  (a-b)) of the study area. The 2D and 3D depth maps both show a relatively deep and shallow region (curved area) with sedimentary thicknesses/depth to the basement ranging between 2.0 and 7.5 km in the south-eastern part of the study area. Away from this zone, lower average sedimentary thickness range (< 3.25 km) was observed in the study area. The 3D surface map of the basement comprises a series of depressions and uplifts (Fig.14b), a typical representation of fold features. These folds are thought to have been formed during the tectonic event that affected the study area in Cenomanian and Santonian times, which led to the (slight brittle and) intense ductile deformation of the basement surface (Kogbe, 1972(Kogbe, , 1979(Kogbe, , 1981Petters & Ekweozor, 1982). The computed depth values can have an implication on the hydrocarbon potential of the study area. An average sedimentary depth of about 4.18 km compares satisfactorily well with previous literature and implies that the area is favorable for the formation and accumulation of hydrocarbon given other factory for a petroleum system. According to Gluyas & Swarbrick (2005), an average sediment thickness of about 3 -4 km is required for the accumulation and entrapment of hydrocarbon. Only the southwestern end and northeastern part of the study area revealed sediment thickness up to 4 km. In addition, the presence of numerous intrusive structures at varying depths signifies a high-temperature environment, which reduces the possibility of hydrocarbon presence. The potential source rocks within the study area would have been broiled beyond the hydrocarbon-generating window, although, the study area possesses immense potential for metalliferous ore deposits.

CONCLUSION
In this study, the acquired satellite gravity data have been analysed and interpreted using combination of enhanced mathematical techniques so that updated information about structural features for concealment of Hydrocarbon accumulation or related mineralizations, edges, lineaments, trends and depth of the gravity source responsible for potential lithological structures in some parts of the Sokoto Basin of northwestern Nigeria were revealed. The separation of regional residual anomaly map reveals variations and distributions of anomalies at various sizes concealed as shear and weak zones characterized by positive and negative (high and low) gravity anomalies values-trending majorly in NW-SE and E-W direction. The interpretation of the results in the study area also reveals the varying amplitudes of the anomaly signature which implies that the gravity source body is not evenly distributed across the study area. The region of low gravity anomaly values from Bouguer anomaly map suggest possible discontinuity or faults or fractures, whereas the region with high gravity anomaly, less negative indicate that the area underlain with iron deposited mineral or escarpment features/ridges. The estimated depths vary between 2.0 km and 7.0 km which suggest that the gravity source body suspected to be 3.25 km is a near surface feature. The local geological mapping of the study area revealed that the area is underlain mainly by red molted iron, sandstone intercalation, ironstone phosphatic nodules, limestone and carbonaceous mudstone. Results showed that the application of potential fields in particular satellite gravity data were successful in imaging the lithospheric structures and estimating the depth to the causative bodies of the sedimentary infill in some parts of the Sokoto Basin.