Lineament characteristics using gravity data in the Garoua Zone, North Cameroon: Natural risks implications

The Garoua Zone in North Cameroon, the subject of this study, is known to have undergone tectonic movements during the Cretaceous, but the zone's structural data remain poorly known. This study exploits the Bouguer anomaly to improve knowledge of Garoua tectonics structures. In order to characterize these structures, two methods are used: Euler's deconvolution method and the method of the horizontal gradient of the vertical derivative. Superposition of the Euler's solutions map for index N=1 with the map from the horizontal gradient of the vertical derivative method allows determination of gravimetric lineaments, interpreted as faults or as linear contacts, from which we deduce a structural map of the study area. Based on this map, we identify sixteen lineaments, of which we count eight as linear contacts and eight as faults. Among the faults, we denote one of depth between 4 and 8 km, five faults of depth ranging between 8 and 13 km, and two faults of depths between 13 and 36 km. Analysis of these faults shows that the seven deepest faults might present a natural risk in our study area. For purposes of civil protection, such deep faults should be monitored and taken into consideration in the implementation of large public works. The structural map, established herein from data on the in‐depth extension of each fault, thus increases scientific knowledge in the area that can be used to site public works in ways that reduce risk.


Introduction
The study area is in North-Cameroon. It lies between latitudes 9°N-10°N and longitudes 13°20'E-14°30'E. Numbers of studies have been carried out in this zone. The southern part of our study area was the subject of a structural study; three faults ( Figure 1) were detected. Two faults come from the study of Mouzong et al. (2014) while one results from the study of Déruelle et al. (2007). Stuart et al. (1985) suggest a thin crust beneath the Garoua rift by using eleven seismic stations. This last result has further been confirmed by Dorbath et al. (1986), Baudin (1991), Kamguia et al. (2005) and Tokam et al. (2010) who underlined an uplift of the upper mantle and associated the thinning of the crust to the opening of South Atlantic Ocean that would have been responsible for the formation of the Benue trough in the Cretaceous.
None of these studies gives complete information on the layout of the various faults of the study area, nor their extensions in depth. However, such structural details are important in that they will not only increase scientific knowledge of the area but also will permit better choices of sites to search for useful substances, such as minerals, or for implementation of large public works, which should not be placed above a deep fault sensitive to the stresses of friction between, or collisions of, tectonic plates, especially since the risk of a fault's destabilization increases with its depth.
The North Cameroon region is relatively aseismic by world standards. Even though it has not experienced large scale earthquakerelated disasters in the past, tectonics activities have still affected our study area, evidenced by the number of tectonic lines existing in this zone ( Figure 1)-thus the importance of locating the various faults in the study area and determining their depth extensions. To achieve this goal, we use the multi-scale analysis of the maxima of the horizontal derivative of the vertical gradient and the Euler deconvolution method on the gravity data.

Geological and Tectonic Setting
The North Cameroon region belongs to the mobile zone of Central Africa. This mobile zone is situated between the Western African craton and the Congo craton. In fact, the work of Wright et al. (1985) showed that the Benue trough is one of the most interesting sedimentary basins of West Africa because of the tectonic movements that account for the marine and continental sediments found there, and also because of the presence of volcanism and the intrusion of plutonic rocks. This trough is directed NE-SW and extends a distance of about 1000 km, with a width of 50-150 km. It contains a great part of the sedimentary basin of North Cameroon, including the Garoua basin. Its origin is related to the opening of the South Atlantic during the Cretaceous, which led to the separation of the African and South American continents (Guiraud and Maurin, 1992). It is made of the sedimentary marine, continental series and is divided into three parts: the low, middle and high Benue. The high Benue includes the Gongola rift and the Yola-Garoua rift. The Yola-Garoua rift extends into Cameroonian territory and is directed E-W while the Gongola rift goes to Niger and is directed N-S. The area of this study (Figure 1) is situated in the Cameroonian territory of the Benue rift (Yola-Garoua branch) and especially in the formation known as Garoua sandstone. Its base is made up of two great lithological sets which underwent a metamorphic and tectonic evolution: the polycyclic unit made of amphibolites, gneiss with pyroxenes, garnet gneiss and gneiss; and the monocyclic whole made of amphibolites, metagabbros, orthogneiss granodioritic, cordierite gneiss and spinel, granites with hypersthene, granites post-tectonics and dolerites (Nzenti, 1998). It comprises granites, gneisses, and micaschists intrusions (Hossere Badjouma, and the East area of Lagdo and Adoumri). Kamguia et al. (2005) have shown that sandstones are accumulated in the Garoua basin and the thickness of these sediments is greater than 4 km. They also suggest that the continental crust below the basin is thinner (about 24 km) than the normal crust, but may be a little thicker to the east. Tokam et al. (2010) used the joint inversion of Rayleigh wave group velocities and receiver functions to study the crustal thickness of the Garoua rift. These studies have found that the crustal thickness is variable across this part of Cameroon. The Moho has a depth of 26-33 km in the Garoua rift. Structurally, Mouzong et al. (2014) determined that the predominant trending of the fault is NW-SE. According to these authors, the potential hydrocarbon field area is structurally controlled and the depth of the basement in the area varies between 4.4 and 8.9 km.
On the topographic plane (Figure 2), the natural site is a plain of about 300 m height. The Garoua Zone is dominated by the hill situated at the NW of the study area whose height ranges from 400 to 597 m, and by the lesser hill (lower than 400 m) situated at the center and South of the basin.

Data
The gravity data used in this work were collected during gravity surveys of Central Africa by ORSTOM and referenced in Collignon (1968), supplemented by data from surveys of Princeton University in 1968, University of Leeds in 1982, and IRGM and University of Leeds between 1984and 1988, referenced in Poudjom-Djomani et al. (1996. The data acquisition campaigns were carried out using cars, along roads or tracks suitable for vehicles. The ORSTOM measurements were carried out at intervals of 3 km while the other measurements were carried out at intervals ranging from 4 to 10 km, and sometimes more. The coordinates were determined from topographic maps and compass routes while the altitudes were obtained by barometric leveling using the Wallace & Tiernan or Thommen (Type 3B4) altimeters and GPS (Global Positioning System). The gravimeters used were the Lacoste & Rombert (Model G, N° 471 and 823), Worden (N° 69, 135, 313, 1153), North American (N° 124 and 165), World Wide (N° 36) and the Canadian Scintrex (N° 305G) gravimeters. The measurements were connected to the gravitational bases of the ORSTOM base network in Africa, known as the Martin network (Duclaux et al., 1954). The values of gravity were measured with an average precision of 0.2 mGal and the error on the coordinates of a station varied from 0.1 to 1 minute, i.e., approximately 200 to 2000 m. This error was less than 200 m for coordinates measured by GPS. The reference system chosen was the IGSN71 (International Gravity Standardization Network, 1971). The data acquired by ORSTOM were attached to the reference system of Postdam 1930. In order to maintain homogeneity all the compiled data were converted into the IGSN71 system. The distribution of gravity measurement points is shown in Figure 3.  The Bouguer gravity anomaly was calculated at each point, assuming an average density of the earth's crust equal to 2.67 g/cm 3 in the terrain and plate correction. The Bouguer gravity anomaly is calculated by taking into account the drift, free-air, Bouguer plate corrections, and terrain correction. However, the data obtained from the IRD (Institut de Recherche pour le Développement) have already undergone drift, free-air, and Bouguer plate corrections. We thus had to determine the terrain correction values by using the topographic data presented in Figure 2. Two steps are necessary to calculate the terrain correction: the first step is to use the Digital Terrain Model SRTM.tem (Shuttle Radar Topography Mission) to create a regional terrain correction grid in Oasis montaj  software; the second step consists of using this grid to correct the terrain effects in the data. After executing these two steps, we obtain the terrain correction values. As these values were positive, they have been added to the gravity data in order to obtain the complete Bouguer gravity anomaly.

Processing of Gravity Data
The gravity data obtained are interpolated using the Kriging method; the spacing of the grid points is 3 km along each axis. The resulting grid is contoured at 5 mGal intervals using Arcgis software to produce the Bouguer gravity contour map of the study area ( Figure 3). Kamguia et al. (2007) has done a comparative study with different interpolation methods applied to the gravity data of Cameroon. In this work, they show that the Kriging method is one of the most flexible methods because it generates the best overall interpretation in most data sets. The filtering of gravity data allowed more information about the basement structure. Thus, in order to best exploit the gravity data, different complementary techniques were applied.

Determination of the Upward Continued Field at Different Heights
The upward continuation is an operator that acts like a filter by attenuating short wave lengths and revealing anomalies related to deeper structures with respect to the continuation height (when continuation height z increases, the depth of the top of the anomaly source increases). Let g(x, y, z) be a function defined in the three-dimensional spatial domain; its Fourier transformation in two dimensions, denoted as G(k x , k y , z), decomposes the function g(x, y, z) in terms of λ x and λ y wavelengths that it contains, where k x and k y are the wave numbers following the x and y directions, respectively. If we know the value of the spectrum for z=0, then we can determine its value for z=h by the relation: We can then perform the inverse transformation of G(k x , k y , h) to obtain g(x, y, h). We have used the Fourpot program to calculate the upward continuation (Pirttijärvi, 2009).

Determination of the Optimum Upward Continuation Height
The optimum upward continuation height h o of the gravimetric field is determined by the empirical method of Zeng HL et al. (2007). It consists of determining the Bouguer upward continuation height where the correlation between the upward continued fields at the successive heights presents a maximum deflection. The steps of data processing carried out for the determination of the optimum height in our study area are as follows.
(1) Upward continuation of Bouguer anomaly at heights from 5 to 120 km, by 5 km intervals.
(2) Calculation of correlation factors between the upward continued fields (g 1 and g 2 ) at two successive heights using the relation (2) proposed by Abdelrahman et al. (1989): where M and N are the number of sampling data along x and y directions, respectively. After applying the entire step shown previously for the determination of optimum upward continuation height, the correlation factor is plotted as a function of increasing continuation height by making each correlation factor correspond to the lower of the two successive heights (Figure 4). The deflection denoted C on Figure 4 at each height is given by the gap between the correlation factor curve and the line joining the two ends of the curve. Figure 5 presents the variation of the deflection with respect to the continuation height. This curve attains a maximum at the height h o =25 km, corresponding to the optimum upward continuation height in our region.

Calculation of the Vertical Derivative at Different Heights
The first-order vertical derivative of the gravity field at each height is calculated in the space domain using the method of finite differences proposed by Florio et al. (2006). This method is more stable than the continuation in the frequency domain (Gunn, 1975), which in some cases enhances data errors, depending on the signal/noise ratio. It also has the advantage of allowing the calculation of the vertical derivative at several heights, using a stable operator like upward continuation (Jacobsen, 1987). Given g DV to be the vertical derivative of potential field g at height h, it is calculated by formula (3).
is the field upward continued at the height h and is the field upward continued at a slightly higher level h+Δh. Δh is a small height difference lying between 1/10 and 1/100 of the data sampling interval (Noutchogwe, 2010). We have taken h = 0.5 km.

Determination of the Horizontal Gradient of the Vertical Derivative and Its Local Maxima
Let g DHDV be the horizontal gradient of the vertical derivative at a height h; its value is calculated in the spatial domain using the formula (4), where g DV is the vertical derivative of potential field g at height h as defined in equation (3). The positions of local maxima are determined by the Blakely and Simpson (1986) method. The horizontal gradients and their maxima were calculated using the BOUNDARY program of the Fortran 77 package of the United States Geological Survey (Phillips, 1997).
This program allows calculation, on a regular grid, of the maxima of the total horizontal gradient of the vertical derivative by comparing the value at a center point of a 3×3 grid window to the surrounding eight in the four directions: horizontal, vertical, and both diagonals. If the center point is found bigger than the other two surrounding points, its index value N increases and, for each direction, a second order polynomial is fitted and the peak position of the parabola is computed. The position of the largest peak value that is located inside the area of the central grid cell (square box in Figure 6) is then used as the maximum position of the horizontal derivative. The index N gives an indication of the linearity of a maximum: if N=1, the anomaly is linear along one of the four directions; if N= 4, the maximum is a local peak.
In fact, the works of Cordell (1979), Cordell and Grauch (1985) showed that the maxima of the horizontal gradient of gravity anomalies help to locate contacts associated with abrupt changes in density, interpreted either as lineaments or intrusive formations. Lineaments are expressed by a quasi-linear disposition of many maxima where the number of maxima showing one lineament must be greater than three and horizontal limits of intrusive bodies are shown by quasi-circular disposition of many maxima. The multi-scale analysis of the horizontal gradient involves an upward continuation of the gravity field to different heights with a view to characterize the importance of vertical extension of anomalous structures and to determine the sense of dip of identified contacts/faults. As a matter of fact, when the lineament presents a dip, the maxima of horizontal gradients are displaced progressively down-dip as the height of upward continuation is increased, but these maxima remain practically at the edge of the lineament in the case of a sub-vertical lineament.

Euler's Deconvolution Method
The Standard Euler 3D method is based on Euler's homogeneity equation. It is an equation that relates the potential field and its gradient components to the location of the source, with the degree of homogeneity N, which may be interpreted as a structural index (Thompson, 1982). This system uses a least squares method to solve Euler's equation (equation (7)) simultaneously for each grid position within a sub-grid (window). The method involves setting an appropriate structural index (SI) to solve the equation for an optimum (x 0 , y 0 , z 0 ) and B. As well, a square window size must be specified, consisting of the number of cells in the gridded dataset to use in the inversion at each selected solution location. The window is centered on each of the solution locations. All points in the window are used to solve Euler's equation for solution depth, inversely weighted by distance from the center of the window. The window should be large enough to include each solution anomaly of interest in the total field gravity grid, but ideally not large enough to include any adjacent anomalies.
In order to define the Euler's equation, Thompson (1982) gives a brief method to describe it. Let any real t, a three-dimensional function f is said to be homogeneous of degree n if it obeys to the expression: From this, he shows that the following (known as Euler's equation) is also satisfied: x ∂g ∂x + y ∂g ∂y By considering potential field data, the final Euler's homogeneous equation is then written as follows: 1 2 3 4 Figure 6. Determination of the field maxima from gridded values using 3×3 point windows (after Blakely and Simpson, 1986). Curves represent contours of the horizontal gradient of the gravity. A parabola is fitted through the four triplets along the four main directions. The number of parabolas that reach a maximum inside the area of center grid cell gives index N.
where (x 0 , y 0 , z 0 ) are the gravity source coordinates, g is the field intensity measured at position (x, y, z), N is the structural index and refers to the geometry of the source, and T, the total field determined at position (x, y, z) is the sum of regional field B and anomaly ΔT due to the disturbance source given by relation: Thompson (1982) and Reid et al. (1990) signaled that the choice of structural index is very important. They define the structural index as a measure of the rate of change of a field with distance. They also established a structural index N for a certain number of structures. It is proved that the structural index N, for a gravity field, may take a value between 0 and 2. It is estimated that a structural index varying from N=0 to N=1 is better adapted for faulted structures, and N=2 for a sphere.
The Euler deconvolution process is an appropriate and more interesting method because it allows determination of the tectonic structures existing in any zone. The success of this method in this application has been well demonstrated: Milligan et al. (2004)

Analysis of Bouguer Anomaly Map
Bouguer gravity map (Figure 3) consists of two major domains of gravity values: the domains of high anomalies sources (anomalies higher than the average, which is of -37 mGal), and the domain of low anomalies sources (anomalies lower than the average), sometimes separated from the first by more or less significant horizontal gradients characterized by narrowing of the isoanomaleis lines. The high anomalies sources can be due to intrusions of dense rocks in the basement. This supports the works of Mouzong et al. (2014) which associate these intrusive rocks to sedimentary rocks or the Precambrian rocks (magma, gneiss and old granites) in the basement. The zones of high anomalies sources are situated at Garoua, Gaschiga, Lago, and in the North-East of Douloumi according to a direction SW-NE to SSW-NNE; the zones of low anomalies sources are observed around Djidda, such as in the South and East of Guider and Bossoum. When we compare the zones showing low anomalies sources with the geology of the study area, the first thought that comes to mind is that of collapse of the basement. This idea is supported by the findings of Dorbath et al. (1986), Baudin (1991), Kamguia et al. (2005) and Tokam et al. (2010) who underlined an uplift of the upper mantle associated with thinning of the crust in the study area.

Analysis of the Horizontal Gradient of the Vertical Derivative at Selected Continued Height
The horizontal gradient of the vertical derivative maps with their corresponding maxima at selected continued heights of 5, 15 and 25 km are shown on Figures 7, 8 and 9. These maps present the zone of horizontal gradients with variable forms and amplitudes dominated by the gradients zones corresponding to the structures of type contact/fault and those corresponding to the presence of an intrusive body: (1) For example, in Figures 7, 8 and 9 respectively, at the East of the study area and its surrounding (Bibemi, Djida, and Guider), we observe that, the more we use the value of the upward continued height, the more circular are the contours lying in the SSE-NNW direction. This correlation suggests that there may exist, in this area, a dense geological structure that varies from shallow to great depth. This hypothesis is supported by Mouzong et al. (2014), who underline the presence of a positive anomaly, situated in the North of Bibemi, that can be interpreted as due to the shallowness of the basement rocks, which are composed mainly of granite, gneiss, schist and micaschist.
(2) In Figures 7 and 8  by Déruelle et al. (2007). Here we also observe that the North of the Gaschiga town is marked by an anomaly oriented in the E-W to SW-NE direction.
(3) In Figures 7 and 8, a fault in the SW-NE direction can be seen to exist in the Pitoa and Badjouma axis. These anomalies, which characterize the gradients zones, may be due to the presence of contact or faults in these zones. Figure 10 shows the superposition of the local maxima of horizontal gradient determined on the vertical derivative of the Bouguer anomaly upward continued at the heights 1, 3, 5, 10, 15, 20 and 25 km.

Multi-Scale Analysis of Maxima and Lineaments Cartography
According to Jacobsen (1987) and Koumetio et al. (2012), the upward continuation at the height h eliminates the effect of the sources situated above the depth h/2. Therefore, the depth of the top of the anomaly sources, obtained after an upward continuation at the height h, is superior or equal to h/2. In Figure 10, black maxima indicate that the effects of the sources situated above the depth 0.5 km are eliminated. This depth becomes 1.5 km for gray maxima, 2.5 km for blue maxima, 5 km for green maxima, 7.5 km for yellow maxima, 10 km for magenta maxima, and 12.5 km for red maxima.
Therefore, the presence of black, gray, blue, green, yellow, magenta, and red maxima on the same lineament indicates that the depth of its bottom is greater than 12.5 km. The presence of black, gray, blue, green, yellow and magenta maxima on the same lineament with the absence of red maxima indicates that the depth of its bottom is between 10 and 12.5 km.
In Figure 10, the presence of blue and green maxima on the lineament F2 with the absence of yellow, magenta and red maxima suggests that the depths of their bottoms are between 5 and 7.5 km. The presence of blue, green, and yellow maxima on the lineaments L4, F6, L8 and F8 with the absence of magenta and red maxima suggests that the depths of their bottoms are between 7.5 and 10 km. The absence of red maxima on the lineaments L1 and F4 suggests that their bottom depths are between 10 and 12.5 km. The presence of red maxima on the lineaments F1, L2, F3, F5, L5, L6, F7 and L7 indicates bottom depths greater than 12.5 km.
Figure 10 also shows that the depth of the bottom of the intrusive body (C) is not uniform: it is between 5 and 7.5 km for the presence of blue and green maxima; between 7.5 and 10 km for the presence of blue, green and yellow maxima; between 10 and 12.5 km for the presence of magenta maxima and the absence of red maxima; greater than 12.5 km for the presence of red maxima.

Analysis of the Euler's Solutions and Structural Interpretative Map
The Euler method has been applied to the Bouguer anomaly maps using a moving window of 3 km×3 km with structural index of 0, 0.5 and 1. The superposition of the maxima map with the Euler's solution allowed us to obtain Figure 11 for index N=0, Figure 12 for index N=0.5 and Figure 13 for index N=1.  Figure 10. Superposition of the maxima of the gradients computed at various upward continuation heights in the study area. ent faults. For example, the presence of a cross (see legend of Figure 13) on the fault F2 indicates that its depth is situated between 4 and 8 km. In the same manner, the presence of a square on F3 shows that the depth of this fault varies from 8 to 36 km. Table 1 shows directions and depths of different faults. In or-der to choose an index that better defines the structure of our study area, we compared the depths determined from the maxima of gradients field with depths from the Euler's solutions (Table 1). From this analysis we observe that the solutions are compatible only for index N=1, which is thus chosen to be the op-  timum structural index in our study area.
The analysis of faults ( Figure 13) allows us to note one fault (F2) which has depth situated between 4 and 8 km, five faults (F1, F4, F5, F6, and F8) with depths between 8 and 13 km, and two faults (F3 and F7) with depths situated between 13 and 36 km. From the map of Figure 13, we deduce the structural map of the study area presented in Figure 14. Here, we denote sixteen lineaments; the quasi-circular contact (denoted C) corresponds to the horizontal limit of the intrusive body.

Discussion
When we apply the multi-scale analysis of the horizontal gradient of the vertical derivative to the Bouguer anomaly, the maxima of the horizontal gradient are found above abrupt changes of dens-ity (Figures 7, 8 and 9). The calculation of the upward continuation at several heights (from 5 to 120 km, by 5 km intervals) leads us to the determination of optimum upward continuation height, which is 25 km. This altitude, which fixes the depth of investigation at h o /2=12.5 km in our study area, is similar to that determined by Noutchogwe (2010) in the Adamaoua zone. The results of the superposition of the maxima from the horizontal gradient of the vertical derivative have allowed us to detect two faults already obtained in previous studies, F1 by Déruelle et al. (2007) and F5 by Mouzong et al. (2014).
Other potentials faults (F2, F3, F4, F6, F7 and F8) are also highlighted and contribute to the structuring of our study area. Superposition of the geological map (Figure 1) with the structural map ( Figure 14) reveals that a majority of faults repose on gneiss and granite formations. This correlates with the work of Mouzong et al. (2014) who associate these faults to the presence in the basement of the Rrecambrian rocks resulting from local tectonic forces in our study area. The superposition of the Euler's maps for index N=1 with the maxima of gradients map ( Figure 13) shows that the depths of the faults F1, F2, F3, F4, F5, F6, F7, and F8 identified starting from the Euler's solutions are in agreement with those determined by the method of the horizontal gradient of the vertical derivative coupled with the upward continuation in the Garoua Zone.
In spite of the presence of all these faults, previous studies in this zone show that its seismicity is low. This implies that the study area is a non-active zone. However, the deeper a fault, the larger its destabilized area is expected to be; also the more easily it is activated by further tectonic events. To prevent serious consequences of such possible tectonic events, more focus is needed on deep-reaching faults. In particular, large construction projects such as major roads, railroads, and buildings are inadvisable in such areas, which can be identified by comparing information obtained from the horizontal gradient of the vertical derivative method and Euler's method. Our findings suggest that the faults F1, F3, F4, F5, F6, F7 and F8 (their depths are between 8-36 km) may require monitoring to insure adequate civil protection.
Although the Euler method is a good and appreciated method, the lineaments situated at the extreme zone of the study area need further investigation. The Euler's distribution solutions at this level are uncertain because of the effect of the survey edge on upward continuation, which may affect calculation of maxima. In order to improve understanding of the lineaments depths in these zones, the study area can be extended to verify that these lineaments at the edges of the studied area are or are not faults.
Consequently, the extreme lineaments of the study area, such as L2, L3, L4, and L7, should be better characterized by a regional study.

Conclusion
The analysis of the Bouguer anomaly map allows us to distinguish the anomalies from structures of different densities. The upward continuation at different gravity field heights leads us to determination of the optimum upward continuation height, which is 25 km. The Euler's deconvolution method and the horizontal gradient of the vertical derivative method have allowed is to highlight sixteen lineaments, among which we identify eight linear contacts and eight faults. The structural interpretation map resulting from the horizontal gradient of the vertical derivative method has been drawn. We suggest that this map presents important information critical to the work of civil engineers engaged in the construction of great works (public works) such as major buildings, roads, railroads, and so on.