Applications of Gravity Method to Map Groundwater Flow in the Ziway Langano Corridor, Central Main Ethiopian Rift

1. Introduction

The East African rift system consists of the Ethiopian, Western, and Eastern rifts, these are considered to be part of a large rift system because they share broad similarities in timing and are related in a geographic sense. For example, volcanism appears to have generally preceded or in some cases have been contemporaneous with major rift faulting [1]. The main Ethiopian rift (MER) is part of the East African rift system which is considered to be the classic example of a continental Rift. The MER is the northernmost of the above rift systems and is the southwestern arm of the Ethiopian–Red Sea–Gulf of Aden triple junction [2]. Gravity survey sets out to investigate the subsurface geology on the basis of anomalies causing gravity field to result from gravity properties of the underlying rocks. It is also used in mapping geological boundaries between contrasting lithologies including faults [3]. This paper deals with the application of gravity techniques to assess and map groundwater flow and their spatial distribution. The use of gravity survey method is to delineate the groundwater flow of Ziway Langano corridor. In the present study; gravity data have been compiled, processed and interpreted to examine the subsurface structures and groundwater potential in the Ziway-Langano corridor. Further, based on the results of gravity anomaly and profiles, this paper attempts to describe the tectonic and geologic implications of the study area.

1.1. Description of the Study Area

1.1.1. Location and Accessibility

The study area is located about 190 km South West of Addis Ababa, In the Central MER. The study area is situated within the Lakes district region between Lake Ziway to the north and Lake Langano to the south and forming some sort of corridor between the two lakes. Meki, Kator and Bulbla rivers are the major rivers which flow to Lake Ziway and Lake Langano. The study area covers an area of approximately 1,800 km². Geographically the study area is bounded between 7°45′N-8°0′N longitude and 38°30′E-39°00′E latitude. As shown in (Figure 1) the study area can be accessed through two roads passing from Adamitulu to the Aluto geothermal to Lake Langano, passing through the Aluto geothermal base camp and it joins again the asphalted road at Bulbula.

1.2. Climate and Vegetation Cover

The Ziway-Langano corridor appears to be a prominent topographic feature relative to the surrounding flat lands and the lake levels of Ziway (1685 m.a.s.l) and Langano (1535 m.a.s.l). The study area is characterized by warm to cool, semi-humid zone, where the annual rainfall is more than 1000 mm and it is a semi-humid to arid zone. Vegetation in the floor of the rift is dominated by the spiny acacia trees, eucalyptus and some short bushes.

1.3. Physiography and Geologic Setting of the Study Area

1.3.1. Physiography

The geomorphology of the area is a reflection of the volcanic and intensive tectonic activities in the geological time and modification through erosion processes. The MER which estimated to be about 80 km wide separates the uplifted Somalia and Ethiopian plateau and these plateaus rises more than 2000 m above mean sea level. The rift floor increases in elevation from the Turkana depression up to the main watershed between the Meki and Awash rivers north of the Lake Ziway to the north direction, generally the elevation of rift floor decreases toward the Afar depression [4].

1.3.2. Regional Geological Setting

The Ethiopian rift system being active rift type, up-doming was followed by volcanism and rifting. The earliest episodes took place in the Southern Ethiopia Rift (SER) and the Afar Rifts, before 25 Ma. Rift propagation proceeded both from the south and north directions towards the MER [5]. Rifting in the MER starts about 8 Ma [6], possibly due to the presence of transversal structures (GBL and YTVL to the south and to the north of MER, respectively), that stopped propagation of rifting for a long geological period. Several central volcanoes that form very thick pyroclastic deposits were erupted during the early stages of the main Ethiopian Rift. These volcanic products have about 5 Ma and are known as the Nazret group [7,8,9,10]. The formation of the rift margins of the Central main Ethiopian Rift, the Munesa and Guraghe escarpments were formed about 3.5 Ma [8]. These escarpments that are formed about 3.5 Ma are curved towards the rift and they are considered as the sources of the acidic pyroclastic deposits and large caldera collapses [8]. The important clue for such conclusion is that the curved escarpments of Munesa and Guraghe are related to the more than 700 m thick pyroclastic deposits that are found in both escarpments. Lithologies that are younger than 1.6 Ma are related to the NE-SW trending fracture systems together with their transversal structures that are NW-SE and ENE-WSW oriented. After 1.6 Ma the direction of the Somalian plate drifting changed from SE towards nearly E direction [5,9,11] causing the development of N-S/NNE-SSW trending fracture systems known as the Wonji Fault Belt (WFB).

2. Methodology

Relevant geophysical data were collected and analyzed for the study of ground water flow in the Ziway Langano corridor. From the results obtained appropriated maps, graphs and 2D models are constructed. This paper has utilized gravity survey methods, all the gravity data utilized in this paper were collected in a random distribution fashion. The selection of the survey patterns in effect has taken the accessibility of study area, scale of the study, time and financial constraints and others in to account. The field measurements have been recorded using a field recording format in at least three attributes (geographic coordinates in terms of Easting and Northing, elevation or height above sea mean level) and the anomaly reading associated with each geophysical method employed. After the necessary reductions and allowances for the effects of such as latitude, elevation, slab and terrain have been made to the raw gravity data. The complete Bouguer anomaly and its derivative anomalies have been computed using the Geosoft Oasis Montaj software (V 6.4.2). The complete Bouguer anomaly map represents the sum of the regional (deep seated) and local (residual) anomalies which arise from responses of the large, deep seated and local, shallow depth geologic materials and structures, respectively. Hence, there is a need to separate the responses of both regional and local features from the complete Bouguer anomaly data for effective and better interpretations of the intended objectives. Accordingly, following the gridding and complete Bouguer anomaly map generation step, the residual and regional anomaly separation process of the gravity data was accomplished using the low and high pass filtering techniques using geosoft oasis montaj Software (V 6.4.2).

3. Data Acquisition and Reduction

3.1. Gravity Data Acquisition

All the gravity data used in this paper amounting to a total of 333 (as depicted in the gravity data distribution plot map, Figure 2); were acquired from secondary sources, namely; the work of Alemu (1992) and from Geological Survey of Ethiopia (GSE).

3.2. Gravity Data Reduction

For most geophysical methods, it is necessary to apply a variety of corrections to the raw data obtained in the field in order to reduce or prepare the data for further enhancement and interpretation [12]. In gravity work, more than in any other branch of geophysics, large and in principle calculable effects are produced by sources that are not of direct geological interest [13]. The gravity data reduction processes was aimed at transforming the raw data sequentially in to a data table of final complete Bouguer anomaly values by correcting the variations in the earth’s gravitational field which don’t result from the differences of density in the underlying rocks. Accordingly, all the stations occupied in this paper have been referred to the International Gravity Standardization Net71 (IGSN 71) gravity datum. The theoretical gravity has been computed using the 1967 gravity formula (Geodetic reference system 1967). All the necessary standard correction types have been made, based on the standard gravity data correction procedures.

Figure 2. Distribution map of gravity data of the study area. Where (LZ) Lake Ziway, (AT)Adami Tulu (AL), AlutoLangano, and (LL) LakeLangano.

Figure 2. Distribution map of gravity data of the study area. Where (LZ) Lake Ziway, (AT)Adami Tulu (AL), AlutoLangano, and (LL) LakeLangano.

3.3. Gravity Anomaly

Gravity observations can be used to interpret changes in mass below different regions of the Earth. To see the mass differences, the broad changes in gravity from equator to pole must be subtracted from station observations by predicting the gravity value for a station’s latitude (theoretical gravity), then subtracting that value from the actual value at the station (Observed gravity), yielding a gravity anomaly [14]. A discrepancy between the corrected, measured gravity and the theoretical gravity is called a gravity anomaly [15]. It arises because the density of the Earth’s interior is not homogeneous. The most common types of gravity anomalies are the Bouguer anomaly and the free-air anomaly.

3.4. Data Processing and Presentations

3.4.1. Gravity Data Processing and Presentations

In this paper, the computed point gravity anomaly values were gridded using the minimum curvature gridding method so as to interpolate the field data on to a regular, square grid using an average grid cell size of 2000 m. A complete Bouguer anomaly map (Figure 3) was finally generated by using the complete Bouguer anomaly value column as gridding column, which was used as the final gravity anomaly map base for further filtering processes, enhancements and interpretations regarding the lateral subsurface density variations. The complete Bouguer anomaly map generated in (Figure 4) represents the sum of the regional (deep seated) and local (residual) anomalies which arise from responses of the large, deep seated and local, shallow depth geologic materials and structures, respectively. Hence, there is a need to separate the responses of both regional and local features from the complete Bouguer anomaly data for effective and better interpretations of the intended research objectives. Accordingly, following the gridding and complete Bouguer anomaly map generation step, the residual and regional anomaly separation process of the gravity data was accomplished using the low and high pass filtering methods using geosoft oasis montaj Software (V6.4.2).

Accordingly, the regional gravity anomaly map (Figure 4) was generated using the geosoft oasis montaj software by subtracting the residual gravity anomaly data from the complete Bouguer anomaly data. Following the separation of the residual and regional anomalies, different data enhancements and filtering techniques have been applied including respective maps such as horizontal gradient map (Figure 6) and Euler depth solutions for SI = 0.5 (Figure 7) are generated, discussed and interpreted.

4. Result and Interpretations

4.1. Complete Bouguer Anomaly Map

The Bouguer gravity anomaly reflects changes in mass distribution below the surface, except the Bouguer anomaly has had an additional correction, removing the effect of mass above sea level datum on land [14]. The variation of the Bouguer anomaly should reflect the lateral variation in density, such that a high density feature in a lower density medium should give rise to a positive Bouguer anomaly. Conversely, a low density feature in a high density medium should result in a negative Bouguer anomaly [16]. A first analysis of the complete Bouguer anomaly map (Figure 3) indicates that the gravity responses reflect the heterogeneity in the subsurface density distribution beneath the survey area. The observed gravity anomaly values range from a minimum of -198 mgal to a maximum of -177 mgal in the study area. The anomaly map further reveals three distinct anomaly zones, which include the medium to very high central anomaly zone which is oriented in a NNW to SSE direction and a medium to very low negative anomaly zones located in the Northeastern, eastern and south eastern part of the study area.

The positive peak (-177 mgal) which corresponds to the location of the Aluto volcanic center may thought to be caused by high density volcanic products. The positive peak (-177 mgal) which corresponds to the location of the Adami Tulu locality is thought to be caused by high density rhyolite intrusions. The relatively low anomaly (-198 mgal) observed over the study area is caused by low density sedimentary depositions occurring at the southern shore of Lake Ziway and northern shore of Lake Langano. The medium negative gravity anomaly (≈ -185 mgal) observed over the summit of the Aluto volcanic complex being superposed on the positive one is thought to be caused by the sedimentary deposit of the caldera and the volcanic ash and pyroclastic materials on the top of the mountain. Since the complete Bouguer anomaly includes effects from the sources of a local and shallow nature as well as regional and deep seated structures, these effects are described in the following section.

Figure 3. Complete Bouguer anomaly map. Where (LZ) Lake Ziway, (AT) Adami Tulu, (AL) Aluto Langano, and (LL) Lake Langano.

Figure 3. Complete Bouguer anomaly map. Where (LZ) Lake Ziway, (AT) Adami Tulu, (AL) Aluto Langano, and (LL) Lake Langano.

4.2. Regional Gravity Anomaly Map

The complete Bouguer anomaly map consists the sum of the effects of features of long wave length, deep seated (regional anomaly) and the near surface, short wave length (residual anomaly) components. The separation of these components is important to determine the effects of the deep seated and the shallow depth sources contributing for the complete Bouguer anomaly. Despite the availability of various methods of anomaly separation, a low pass filtering with a cutoff wavelength 10 km was applied to the complete Bouguer anomaly for generating the regional Bouguer anomaly map using geosoft oasis montaj (V 6.4.2) based on the following series of steps. During filtering, increasing the separation between the source and detector cause a decrease in amplitude and increase in wave length of the response, where the rate of attenuation of amplitude is wave-length dependent (i.e. Shorter wavelengths, high frequencies associated with near-surface sources attenuate more rapidly with height than the longer wavelengths and lower frequencies sources. Hence, shorter wave length, shallow-sourced responses plus any short wave length noises are suppressed and longer wave length responses dominate the filtered data [12]. The Regional gravity anomaly map (Figure 4) reveals three prominent anomaly zones which are similar to those of the complete Bouguer anomaly map (Figure 3). The only difference between them is that the anomaly features observed in regional anomaly map are smooth and elongated this indicates the effect of deep seated geologic structures.

Figure 4. Regional gravity anomaly map. Where (LZ) Lake ziway, (AT) Adami Tulu, (AL) AlutoLangano, and (LL) Lake Langano.

Figure 4. Regional gravity anomaly map. Where (LZ) Lake ziway, (AT) Adami Tulu, (AL) AlutoLangano, and (LL) Lake Langano.

4.3. Residual Gravity Anomaly Map

The residual gravity anomaly map was generated after the separation of the regional Bouguer gravity anomaly value from the complete Bouguer anomaly value following the use of high pass filtering methods with a cutoff wavelength 10 Km using geosoft Oasis montaj software (V 6.4.2). This map is expected to reveal the gravity anomalies resulting from shallow depth anomalous geologic structures in the subsurface. As shown in (Figure 5), the residual Bouguer anomaly map is dominated by a number of positive and negative short wave length anomalies occurring as patches of circular and elliptical features aligned almost in NE-SW and NNW-SSE directions. The positive highs in the study area coincide with the locations of high density volcanic products of the Aluto volcanic center and the Adami Tulu locality in the area are attributed to the presence of dense intrusions associated with these young silicic volcanoes. The circular lows which appear to be surrounded by the highs could be the expression of thick accumulation of the rift floor sediments. The central low which is surrounded by the highs of the Aluto high density volcanic products is interpreted as the effect of the sedimentary deposit of the Aluto caldera and the volcanic ash and the pyroclastic materials on the top of the mountain.

4.4. Horizontal Gradient Gravity Map

The horizontal gradient is a mathematical anomaly enhancement method that can be computed by:

$$HG=HG(x,y)=\sqrt{{({\displaystyle \frac{\partial \Delta g}{\partial x}})}^2+{({\displaystyle \frac{\partial \Delta g}{\partial y}})}^2}$$

Where $\Delta g$ is the residual anomaly,

The horizontal gradient method is used to locate the boundaries of regions of contrasting density from gravity data. The method is based on the premise that the maximum horizontal gradient of the observed gravity anomaly, which is caused by a planar anomalous body, tends to overlie the edges of the body (Figure 6). The edges mark locations of the contacts indicating that the HG method is designed to inspect fault and contact features. Horizontal gradient map can be produced from gridded residual gravity anomaly map data using the geosoft oasis montaj software. The compiled HG map is judged to be simple, interesting and a reputable end result which reveals the anomaly texture and highlight anomaly-pattern discontinuities. The HG maxima identified on such maps can be inspected to assist in the location of linear structures such as faults and lineaments. The horizontal gradient gravity map (Figure 6) of the survey is compiled from the gridded residual gravity anomaly data. The map reveals lineated maximum horizontal gradient anomalies which are dominantly oriented in the N-S direction. Locations of the maximum horizontal gradient magnitude are shown with black lines superposed along their orientation (Figure 6). The high horizontal magnitude plots outline lateral mass heterogeneities which are thought to be the expressions of lithological and structural boundaries. The high gradient magnitudes plots trace linear features trending N-S. The linear features are attributed to faults and fractures, and their continuities are observed in the study area. Moreover the N-S trending geological structures (faults and fractures) revealed by the horizontal gradient map can be considered to be the major geologic structures that favor the flow of ground water from Lake Ziway toward Lake Langano.

Figure 6. Horizontal gradientgravity map of the study area. Black lines represent faults traced in the study area. LZ (Lake Ziway), LL (Lake Langano), LA (Lake Shala), Law (Lake Awasa), MAB (Munesa – Asela border fault).

Figure 6. Horizontal gradientgravity map of the study area. Black lines represent faults traced in the study area. LZ (Lake Ziway), LL (Lake Langano), LA (Lake Shala), Law (Lake Awasa), MAB (Munesa – Asela border fault).

4.5. Euler Deconvolution Gravity Map

The Euler deconvolution gravity data enhancement method is used to estimate the depth and location of geologic sources that generate gravity anomalies. The Euler deconvolution equation which is thought to operate for gravity data in 3D was developed by [17] given by:

$$(x-x_o){\displaystyle \frac{\partial \Delta g}{\partial x}}+(y-y_0){\displaystyle \frac{\partial \Delta g}{\partial y}}+(z-z_0){\displaystyle \frac{\partial \Delta g}{\partial z}}=n(\beta -g)$$

Where $(x_o,y_{0,},z_0)$ is the position of a source (source location) whose complete Bouguer gravity anomaly is computed at (x, y, z), β is the regionl gravity anomaly $\Delta g_R$ value, and n is the structural index (SI), which can be defined as the rate of attenuation of the anomaly with distance. The SI must be chosen using prior knowledge of the source geometry. For example, SI = 2 for a sphere, SI = 1 for a horizontal cylinder, SI = 0 for a fault, and SI = –1 for a contact [18].

The Euler deconvolution method is applied to the complete Bouguer gravity data in order to estimate the depth and location of the gravity anomaly sources. Accordingly, the Euler deconvolution gravity map (Figure 7) of the study area is produced for a structural Index (SI = 0.5). The structural index (SI = 0.5) is chosen based on the objective of the survey, which is mapping the lateral and depth wise variation of linear structures (faults and fractures) found in the study area. Euler solutions obtained from the complete Bouguer gravity data show well clustering over the area of Ziway Langano Corridor. These clustering confirms the minor faults in the geologic map and infer new faults in different directions at different depths. Combining these results in one tectonic map could be giving a clear image of the subsurface structure of the area of Ziway-Langano corridor. The Euler gravity map (Figure 7) show that high clustering of the depth solutions for SI = 0.5 are values associated with regions of intense faulting, fracturing and extensional activities.

Figure 7. Standard Euler deconvoltiongravity depthmap for SI = 0.5. Where (LZ) Lake Ziway, (AT) Adami Tulu, (AL) Aluto Langano, (LL) Lake Langano.

Figure 7. Standard Euler deconvoltiongravity depthmap for SI = 0.5. Where (LZ) Lake Ziway, (AT) Adami Tulu, (AL) Aluto Langano, (LL) Lake Langano.

5. Conclussions

From the gravity survey, areas of low gravity observed on the Bouguer anomaly map is interpreted as highly water saturated region of the area because of the close proximity of the profile to the Langano and Ziway lakes and the absence of structures that act as a barriers for ground water movement. Based on the joint analysis gravity methods, the area has good ground water potential for further development from the central to the north and south part of the study area. It has also been possible to map the presence of weak zones/ fractures in the subsurface which could act as conduits for the circulation of fluids in the Ziway-Langano corridor. The lowering of gravity anomaly could suggest the circulation of ground water through the fractures/faults to result in ground water flow. The geologic and tectonic activities that occurred in Quaternary and even to recent years have directly affects the ground water flow of the study area.