Gilgil area is located between Naivasha and Nakuru in Nakuru County, Kenya. It lies 121 km north of Nairobi. Gilgil area is in the Kenyan rift where a number of geothermal fields lie. Preliminary surface investigations have been carried out in Suswa, Longonot, Olkaria, Eburru, Menengai, Bogoria, Baringo, Korosi, Silali and Emurangogolak geothermal fields [1]. Drilling has been done in Eburru and Olkaria. The present power station is in Olkaria. Thus this study was carried out to establish the potential of Gilgil area as a geothermal reservoir.

Gravity surveying has been done to gain information on geothermal potential areas. Gravity technique in geophysical exploration deals with measurements of changes in the Earth’s gravitational field strength [2]. Gravity measurements and observations are done on the earth’s surface. The gravimeter is an instrument used to measure changes in the Earth’s gravitational field on the Earth’s surface and records its values in milligals. It helps to find bodies within the subsurface of the earth which have greater or lesser density than the surrounding host rocks. Gravity can also constrain data during interpretation of other geophysical techniques such as seismic and magnetic.

Gravitational field is natural on the earth’s surface similar to magnetic and radioactivity. It is a natural field technique that uses gravitational field of the earth. There is no energy required to be put into the subsurface to gain information [3]. It reveals change in these natural gravitational field that is attributed to economic feature of concern within the subsurface. This feature portrays a subsurface area of anomalous mass and causes localized change in gravity referred to as gravity anomaly.

The geology of Gilgil area is as a result of volcanism and tectonic activities of the rift valley. The volcanism of the rift preceded and accompanied the rift tectonism. Gilgil area is dominated by quaternary volcanic ash and diatomaceous silts in the plain areas and some volcanic tuff, lava flow and diatomite deposits in the higher escarpments [4]. Alkaline volcanism composed of pumiceous pyroclastics, ashes, trachytes, ignimbrites, phonolites and phonolitic trachytes, tuffs, agglomerates and acid lava dominates Gilgil area. Also volcanic soil and diatomite deposits dominate the area with trona impregnated silts bordering Lake Elmenteita [4]. The area is also characterized by repeated volcanicity followed by movement. The eruptives in each episode start with basalt [5]. The southern part of Gilgil is within the Olkaria volcanic complex. Craters, fumaroles, hot springs and steam vents are found in several places within the Olkaria and Eburru area [5]. The earlier tectonic geology is reflected in the step-faults of Satima and Kinangop generating Kinangop plateau. Grid faulting generated Gilgil plateau while the Mau escarpment is as a result of fault flexures. The major fault escarpments influence topography of the rift floor that influences the drainage flow pattern [6].

In gravity technique, the geology is examined on the foundations of changes in the Earth’s gravitational field emerging from deviations of mass within the underlying rocks. The fundamental concept is the idea of a causative body, which is a rock of unusual density from the host masses (Figure 1). This causative body portrays a subsurface region of abnormal density and results in change in the Earth’s gravitational field called gravity anomaly [7]. An area of approximately 68 km² was covered during this study. Data was gathered from 147 measurement points using CG-5 gravimeter. The base stations were formed for the purpose of drift corrections. Stations were spaced at 500 m apart. At each station the time, northing, easting, altitude and gravity value in milligals was recorded.

Figure 1: Map showing the geology of the study area.

Regional density

The average density of rocks in the study area ρ_awas taken as 2.67g/cm³ [4]. Density of an intruding body ρ_b ranges from 2.70 g/cm³-3.20 g/cm³ [7]. Density contrast is given by

ρ =ρ_b - ρ_a

Density contrast range was found to be 0.03 g/cm³-0.53 g/cm³ and was employed for modelling. Body of density contrast 0.03 g/cm³- 0.53 g/cm³ is associated with heat source at its basin because it best forms at plumes and hotspots below the continent. Mostly forms as an extrusive rock such as lava flow but can also form as intrusive bodies like dike or sill [4]. During gravity forward modelling of this study, a density contrast of 0.25 cm³ and 0.28 cm³ produced the best fit between observed gravity anomaly and computed gravity anomaly.

The Bouguer anomaly map

Gravity anomalies are obtained after reductions have been done to the observed gravity data. If there were no mass distribution within the Earth’s subsurface, the gravity anomaly would be zero.

The contour map in Figure 2 was generated from processed gravity data. This map shows contour intervals of 1 mgal with the highest value at −181 mgal and the least value at −211 mgal. To the Northeast, the map reveals gravity highs with few gravity lows. To the southeast, the map reveals gravity lows with a few gravity highs. To the northwest, the map reveals gravity lows with a small part of gravity high. An intruding rock has density ranging from 2.70 g/cm³-3.20 g/cm³ which is a gravity high [7]. Geothermal reservoir is associated with a gravity high because materials coming from the Earth’s mantle are of higher density than materials found in the Earth’s crust.

Figure 2: Contour map profiles for the Gilgil area.

Euler deconvolution

Euler deconvolution technique provided automatic approximations of a causative body location and its depth within the Earth’s subsurface. Therefore, Euler deconvolution located the boundary of the said resource and its depth from the surface. The most important outcome of Euler deconvolution is the description of trends and depths [8]. In this study, a structural index of 1.0 was used as it best delineates fractures and intruding dykes in the subsurface which are associated with heat sources.

Euler solutions along profile PP’ as shown in Figure 3 suggest a causative body which occurs at maximum depth of 2053.74 m. It reveals a fault at 1000 m and 2000 m along the profile [9]. It also shows a causative body at 1000 m and 3000 m along the profile which has a material of higher density than the host rock.

Figure 3: Euler solutions along profile PP’.

Euler solutions along profile QQ’ as shown in Figure 4 reveals a causative body which occurs at a maximum depth of 792.74 m. It has imaged a body of higher density than the surrounding rock at 1000 m along the profile [10]. At 200 m along the profile, there is a shallow causative body.

Figure 4: Euler solutions obtained along profile QQ’.

Euler solutions along profile RR’ as shown in Figure 5 shows a causative body which occurs at a maximum depth of 1194.21 m. The body occurs between 1200 m and 2000 m along the profile and has a higher density than the host rock. It also reveals a fault between 1200 m and 2000 m along the profile.

Figure 5: Euler solutions obtained along profile RR’.

Euler solutions along profile SS’ as shown in Figure 6 shows an intrusive body which occurs at maximum depth of 2837.15 m. It has also imaged faults at about 2000 m and 4000 m along the profile which is filled by a material of higher density than the surrounding host rock.

Figure 6: Euler solutions obtained along profile SS’.

Euler solutions along profile TT’ as shown in Figure 7 reveals a causative body which occurs at a maximum depth of 4331.38 m. It has also imaged a fault at about 2000 m to 3000 m along the profile which has a higher density than the surrounding host rock.

Figure 7: Euler solutions obtained along profile TT’.

Forward modeling

This was done using GRAV2DC software in surfer 8 computer programme. Modelling entailed construction of an appropriate model based on geological information of the study area [11]. The cross section data was transferred to GRAV2DC software for forward modelling. The parameters determined by Euler deconvolution acted as start-up parameters for the model bodies. The model’s gravity anomaly was computed and compared to the observed anomaly. Features of the model were altered to increase the correspondence of observed anomaly and computed anomaly. In this interpretation, the depth and density contrast of a causative body was determined [12]. Models constructed are as shown in Figures 8-12.

Figure 8: Model fit on residual bouguer anomaly profile PP’.

Figure 9: Model fit on residual bouguer anomaly profile QQ’.

Figure 10: Model fit on residual bouguer anomaly profile RR’.

Figure 11: Model fit on residual bouguer anomaly profile SS’.

Figure 12: Model fit on residual bouguer anomaly profile TT’.

Profile PP’ is on the northern part of the study area as shown in Figure 2 and it cuts across a gravity high anomaly region trending in a NW-SE direction. Models on profile PP’ as shown in Figure 8 reveals two subsurface intrusive bodies [13]. The first body has a density of 2.92 g/cm³ and imaged at a depth of 169.48 m while the second body has the same density of 2.92 g/cm³ and imaged at a depth of 159.51 m. This gravity high could be due to hot intrusive bodies of high density from the mantle which are probably feeding the hot spring in the area.

Profile QQ’ is on the southern part of the study area and it cuts across a gravity high anomaly region trending in a NW-SE direction. Models on profile QQ’ as shown in Figure 9 reveals an intrusive body of density 2.95 g/cm³ and imaged at a depth of 50.33 m. Presence of recent volcanic soil shows there was volcanic activity which deposited high density materials close to the surface hence the imaged body could be a cooling dyke injection [14].

Profile RR’ is on the north western part of the study area and it cuts across a gravity high anomaly region trending in a NW-SE direction. Models on profile RR’ as shown in Figure 10 shows an intrusive body of density 2.92 g/cm³ and imaged at a depth of 370.22 m. This was presumed to be a dense body imaged under a volcanoe which is probably a hot intruding dyke hence a heat source at the basin [15].

This model shown in Figure 11 was oriented in a SW-NE direction to constrain the density contrast and depth of profile RR’ and PP’. It generated the same values as obtained in Figure 8 and Figure 10. The density contrast of the intruding bodies was found to be 0.25 g/cm³. The depth for body 1 was 365.99 m and body 2 was 169.48 m. It has imaged faults responsible for underground thermal movement and massive intrusions which could be heat sources [16].

The model fit shown in Figure 12 was oriented in a nearly N-S direction to constrain the density contrast and depth of profile PP’ and QQ’. It generated the same results as in Figures 8 and 9. The density contrast for body 1 was 0.25 g/cm³and body 2 was 0.28 g/cm³. The depth was 129.16 m for body 1 and 50.61 m for body 2. These bodies were interpreted to be dense intruding dykes into the subsurface which could be heat sources.

Profile PP’ is on the northern part of the study area as shown in Figure 2. It cuts across a gravity high anomaly region trending in a NW-SE direction. Models on profile PP’ as shown in Figure 8 reveals two subsurface intrusive bodies. The first body has a density of 2.92 g/ cm³ and imaged at a depth of 169.48 m while the second body has the same density of 2.92 g/cm³ and imaged at a depth of 159.51 m. This positive gravity anomaly could be a result of hot intrusive bodies of high density from the mantle under the volcanic complexes which are probably feeding the hot spring in the area hence there could be a heat source at the basin.

Profile QQ’ is on the southern part of the study area as shown in Figure 2. It cuts across a gravity high anomaly region trending in a NWSE direction. Models on profile QQ’ as shown in Figure 9 reveals an intrusive body of density 2.95 g/cm³ and imaged at a depth of 50.33 m. This was presumed to be due to phonolitic trachytes during the lower Pleistocene period in the study area. Presence of recent volcanic soil shows there was volcanic activity which deposited high density materials close to the surface. Probablythere was a volcanic activity in the area which stopped hence the imaged body could be a cooling dyke injection.

Profile RR’ is on the north western part of the study area as shown in Figure 2. It cuts across a gravity high anomaly region trending in a NW-SE direction. Models on profile RR’ as shown in Figure 10 shows an intrusive body of density 2.92 g/cm³ and imaged at a depth of 370.22 m. This was presumed to be a dense body imaged under a volcanoe which is probably a hot intruding dyke hence a heat source at the basin.

Profile SS’ shown in Figure 2 was drawn to cut across profile RR’ and profile PP’ for the purpose of constraining the density contrast and depth of imaged body. Also profile TT’ in Figure 2 was drawn to cut across profile PP’ and QQ’ for the purpose of constraining the density contrast and depth of imaged body. Profile SS’ and profile TT’ gave the same density contrast and depth as profile PP’, profile QQ’ and profile RR’ as shown in Figure 11 and Figure 12. The two profiles have imaged faults and massive intrusions which could be heat sources. Due to hydrothermal activity and imaged fractures in the area, probably these hot intrusive bodies for profile PP’ and RR’ are responsible for the hot spring in the area.

The Gilgil prospect area is located in the Kenyan rift where a number of geothermal fields lie. It is characterized by major fracture lines, quiet volcanic craters, fumaroles and hot springs. Fractures resulting from extensional tectonics of continental rifting probably provides a good structural set up that allows water from the rift scarps to penetrate deep into the crust, towards the hot magmatic bodies as modelled under the volcanoes and normal faults conducting hot fluids from deep into possible geothermal reservoirs at shallower depth. Major fracture lines were imaged in this prospect area by Euler deconvolution. These faults are at different depths from the surface as shown by Figure 3 through Figure 7. There are those at the deep basement while others at the shallow subsurface. The deep faults transport thermal fluids from deep parts of the Earth to the subsurface. Also the shallow faults in the Earth’s subsurface direct the flow of thermal fluids on the upper part of the basement. The top faults direct the flow of water from the rift scarps to the hot masses underground. This probably led to the trapping of a heat source in the area as evidenced by hot springs. The deep water circulation would therefore collect heat from the bodies and discharge it through hot springs along faults and fractures as observed in the study area.

Geothermal occurs along major fracture lines, inactive volcanic craters and where there are hot springs. Therefore the modelled bodies across the selected gravity profiles lie relatively at shallower depths as shown from the models in Figure 8 through Figure 12. Thus, the high heat flow observed in the area as evidenced by hot springs could be due to shallow dense intruding bodies within the rift floor faults. Gravity technique was able to locate these dense bodies within the Earth’s subsurface as positive gravity anomalies. It is postulated that intrusives, in the form of dykes would be tapping heat from large magma bodies at few kilometres from the surface. Therefore this study was able to detect gravity highs that show evidence of a buried dense body compared to the surrounding rocks. The buried dense bodies were interpreted as intruding dyke injections within the subsurface which could be heat sources.

This gravity study was done to gather information on the possibility of geothermal occurrence in Gilgil area. It has provided information which will be used as a start point for future detailed geophysical study. Gravity technique is ambiguous and this implies that any anomaly could be a result of many possible sources. To reduce this ambiguity during interpretation, this study recommends the application of other geophysical techniques like seismic, magnetotelluric (MT) and magnetic for outcome comparisons. This ensures confirmation of this gravity results before drilling is done which is expensive.

Special thanks go to research supervisors, Dr. Willis Ambusso of Kenyatta University and Dr. John Githiri of Jomo Kenyatta University of Agriculture and Technology for their technical guidance, suggestions and encouragement during this research work. Thanks to the staff physics department of Kenyatta University for their valuable suggestions.