Depth Estimation of Geothermal Heat Structures by Euler Deconvolution of Gravity Data at Eburru Area, Kenya

Show more

1. Introduction

Gravity is an indirect geophysical method used in geothermal heat sources exploration to determine subsurface dense bodies (Georgsson, 2009). The basic principle in gravity surveying is changes in the Earth’s gravitational field arising from the differences in density of rocks within the subsurface. A region within the subsurface that has a different density from the surrounding masses results into a perturbed gravitational field referred to as a gravity anomaly (Kearey, Brooks & Hill, 2002). Volcanic centers, where geothermal activity is common, are signs of cooling magma or hot rock beneath these areas, as evidenced by the recent volcanic flows, volcanic domes and hydrothermal activities occurring in the form of fumaroles and hot springs (Mariita, 2011). These volcanic centers are usually mapped with gravity highs. Gravity studies in volcanic areas have effectively demonstrated that this method provides a good evidence of shallow subsurface density variations, associated with the structural and magmatic history of a volcano (Ndombi, 1981). Gravity highs correspond well to volcanic centers, faults and geothermal activities. For example, Olkaria Domes and Suswa geothermal centers are located on the crest of a high gravity (Mariita, 2011).

Euler deconvolution method provides location and depth estimates of anomaly sources within the Earth’s subsurface (Nyakundi, Githiri, & Ambusso, 2017). Therefore, Euler deconvolution locates the causative body and determines its depth from the observation level. It generates a map that exhibits the positions and depths of the observed gravity anomaly sources. Euler deconvolution does not adopt any geologic model (Pawan, Ramprasad, Ramana, Desa, & Shailaja, 2007), and thus, it can be used for gravity interpretation without considering the geology of the study area.

2. Methodology

2.1. Gravity Survey

Gravity data was collected from 195 observation points covering an area of approximately 144 km^{2} using a CG-5 gravimeter. The daily instrumental drift was addressed by repeating the measurements at the base station in the beginning and at the end of any day data acquisition. For every observation point, the time, northing, easting, altitude and gravity value in milligals were recorded. Drift, Latitude, Free air, Bouguer and Terrain corrections were applied to the collected gravity data to obtain complete bouguer anomaly (CBA). The CBA information was moved to Oasis Montaj Programme for Euler deconvolution analysis.

2.2. Euler Deconvolution

Euler deconvolution is a potential field data analysis technique for estimating the depth and position of a causative body (Nyakundi et al., 2017). The technique combines the potential field and its gradient components to locate the potential anomalous source, with the strength of homogeneity implied as a structural index and it is a suitable method for figuring out anomalies resulting from isolated and multiple sources (Dawi, Liu, Shi, & Luo, 2004). The technique relates the potential field, for example, gravity or magnetic field and its gradient components to the location of the source of an anomaly with a degree of homogeneity expressed as a structural index.

Euler deconvolution is based on the Euler equation of homogeneity (Reid & Thurston, 2014) expressed in Equation (1);

$\left(x-{x}_{0}\right){T}_{zx}+\left(y-{y}_{0}\right){T}_{zy}+\left(z-{z}_{0}\right){T}_{zz}=n\left({B}_{z}-{T}_{z}\right)$ (1)

where ${T}_{z}$ is a vertical component of gravity anomaly source with degree of homogeneity n, $\left({x}_{0},{y}_{0},{z}_{0}\right)$ is the coordinate of the gravity anomaly source in the Earth’s crust to be determined while $\left(x,y,z\right)$ is the measured co-ordinate. Parameters $\left({T}_{zx},{T}_{zy},{T}_{zz}\right)$ are the determined gradients in the x-, y- and z-directions. n is the structural index and ${B}_{z}$ is the regional gravity value to be approximated.

The new located Euler deconvolution method first determines the analytic signal, finds peaks in the analytic signal then uses these peak locations for Euler deconvolution. The analytic signal grid was calculated and displayed from derivative grids in this study. Analytic signal grid expression shown in Equation (2) is the square root of the sum of the squares of the derivatives in the x, y and z directions (Thompson, 1982).

$A=\sqrt{\left(\partial x\times \partial x\right)+\left(\partial y\times \partial y\right)+\left(\partial z\times \partial z\right)}$ (2)

where A is the analytic signal grid.

Located Euler deconvolution is favored as solutions are only determined over identified analytic signal peaks, the window size varies according to anomaly size and the final solution involves only a few more accurate depth estimates.

2.3. Structural Index

A structural index is applied when performing Euler deconvolution analysis. In regional interpretation of gravity data, structural indices of 0.5, 1.0 and 2.0 are common for fault, contact, sill and dyke location (Felipe & Valeria, 2017). A structural index is a measure of the rate of change of potential field with distance. The basic principle of the structural index is Euler’s equation of homogeneity (Reid, Allsop, Granser, Millett, & Somerton, 1990) expressed in Equation (3).

$x\frac{\partial f}{\partial x}+y\frac{\partial f}{\partial y}+z\frac{\partial f}{\partial z}=nf$ (3)

For gravity field data analysis, the expression in Equation (3) is rewritten as Equation (4).

$\left(x-{x}_{0}\right)\frac{\partial T}{\partial x}+\left(y-{y}_{0}\right)\frac{\partial T}{\partial y}+\left(z-{z}_{0}\right)\frac{\partial T}{\partial z}=N\left(B-T\right)$ (4)

where $\left({x}_{0},{y}_{0},{z}_{0}\right)$ is the position of anomaly source whose vertical component of gravity field T is measured at $\left(x,y,z\right)$. N is the structural index, and B is the regional value of gravity field (Reid & Thurston, 2014).

A structural index is calculated by determining how many infinite dimensions are present in a geologic representation (Thompson, 1982). The representation structural index is the infinite dimension number taken away from the maximum structural index for the field. Maximum structural index for gravity field is two because the gravity field from a point source dies off as $1/{r}^{2}$ (Kearey et al., 2002). A structural index of zero suggests that the potential field does not vary with distance from the anomaly source which is not true in real experience.

There are possible geologic interpretations for a particular structural index in gravity field data. The structural index of 0.5 is associated with the geologic model of a sill, dyke or step while a structural index of 1.0 is associated with a pipe-like model (Thompson, 1982). The structural index of 2.0 in a gravity field generates relatively deeper solutions compared to structural indices of 0.5 and 1.0 and is associated with a geological model of a sphere.

During this study, complete Bouguer anomaly (CBA) data was gridded, processed, starting with a run for a structural index of 0.5, results were then plotted, windowed, visualized on screen and then final map plotted whereby a 15% depth tolerance was used to determine the accepted solutions. The window size varied according to anomaly size. This procedure was repeated for 1.0 and 2.0 structural indices.

3. Elevation Map of Eburru Area

Figure 1 shows the terrain of the study area in heights above the sea level which ranges from the lowest at approximately 1810 m to the highest point at approximately 2519 m. It reveals the highest point to the south at approximately 2519 m and the lowest point to the North of the study area at approximately 1810 m. Generally, the elevation map exhibits a highland at the southern part of the Eburru crater and altitude decreases northwards towards Lake Elmentaita. The volcanic activity in the area probably resulted in increased elevation to the south due to mantle materials pushing within the Earth’s crust. High gravity values in this area could be due to hot denser materials from the mantle that penetrated towards the crust.

4. Complete Bouguer Anomaly (CBA)

Complete Bouguer anomaly (CBA) was obtained after conducting gravity corrections to the collected raw gravity data. The average density of rocks within the Kenya’s rift valley is 2.67 g/cm^{3} (Nyakundi et al., 2017) and was used during gravity reductions. If the density of rocks in the Earth’s crust was constant, the complete Bouguer anomaly would be zero hence would not reveal any anomaly.

Figure 1. Topography map for the Eburru area.

Figure 2 shows the complete Bouguer anomaly map of the Eburru study area. The CBA values range from a −272 mGal to −286 mGal. The variation of density in the underlying rocks within the Earth’s crust results in a very small change in the Earth’s gravitational field. In this study area, there is a 14 mGal variation in the Earth’s gravitational field between the highest density material and lowest density material.

Qualitative interpretation was done using the CBA map in Figure 2. These entails visual inspection of the map to identify areas of gravity highs and those with gravity lows. This reveals the variation of the density of underlying structures. Gravity high implies denser underlying bodies while gravity low implies the presence of less dense bodies than host materials in the Earth’s crust. Pink color on the map shows a region with a high gravity anomaly which decreases to blue color which shows region with a low gravity anomaly.

Northern parts of the study area towards Lake Elmentaita shows gravity high region with small sections of gravity lows. Towards the central part of the study area, there is a gravity low region trending from northeast towards the south and then to the west with an average CBA of −280 mGal. Eastern parts of the map show a negative anomaly and reveal a positive anomaly to the west. To the

Figure 2. Complete Bouguer anomaly map for the Eburru area.

southern part, there is generally gravity low region with some few spots of gravity highs. These gravity highs could be dykes in the form of high-density materials intruding from the mantle. To the Northeast, the map reveals gravity highs with few gravity lows. An intruding mass from the mantle has a density range of 2.7 g/cm^{3} - 3.2 g/cm^{3} which is a positive anomaly (Kearey et al., 2002). Heat source is linked with a positive anomaly as masses intruding are denser. Temperature of the Earth increases from the Earth’s surface towards its interior and therefore intruding materials from the mantle have elevated temperature. Due to temperature gradient, there is continuous heat flow towards the Earth’s crust through the intruding materials. If these materials are trapped by impermeable rock cap within the Earth’s subsurface, they form a geothermal reservoir. Therefore, regions with high gravity anomaly imply denser masses in the Earth’s crust, which could be heat sources.

5. Residual Bouguer Anomaly

Residual Bouguer anomaly is the outcome of the effect of density variation of near-surface bodies which are of interest. Geothermal wells are done to the depth of approximately 4000 m (World Energy Council, 2013) hence deep heat sources beyond 10,000 m may not be of interest. Shallow anomalous bodies of less than 500 m depth increase the amplitude of gravity anomaly due to their low wavelength thus generating noise to the result. Shallow structures are not of interest as geothermal heat cannot be trapped at these depths because elevated temperatures and pressure will erupt to the Earth’s surface. Therefore, in this study, the effect of shallow masses less than 500 m depth and deep masses beyond 10,000 m depth were removed to enhance features of interest between the bandpass.

Figure 3 shows a residual anomaly map where the effect of structures deeper than 10,000 m depth and the effect of shallow masses less than 500 m depth were removed from the complete bouguer anomaly. Figure 3 displays an enhanced effect of Bouguer anomaly distribution between 500 m and 10,000 m depth bandpass within the Earth’s subsurface. It displays a series of high Bouguer anomalies represented by pink color and low Bouguer anomalies represented by blue color. The residual gravity anomaly amplitude ranges between the highest at 3.4 mGal and lowest at −3.0 mGal. This shows that the study area has alternating underlying structures of high and low densities. This could be due to a

Figure 3. Residual gravity anomaly map for the Eburru area.

series of intruding bodies from the mantle into the subsurface. This also indicates a system of underground faults and fractures that probably controls the movement of hydrothermal fluids in this study area.

6. Euler Deconvolution Results and Discussion

The structural index of 0.5 generated five solutions at the depth range of 433 m - 2269 m occurring on gravity highs, as displayed in Figure 4. The figure exhibits the deepest Euler solution occurring to the north of the study area at a depth of 2269 m on a high gravity amplitude of 3.4 mGal. These could be hot dense materials from the mantle trapped at this depth in the Earth’s crust. The shallowest Euler solution occurs to the south of the study area at a depth of 433 m with a high gravity amplitude of 3.4 mGal. The depth is shallow because the solution occurs towards Eburru crater, where magmatic materials were extruded to the surface in the form of lava flow. The structural index of 0.5 in Figure 4 was observed to trend in NE-SW indicating ribbon-like structures in this direction.

The structural index of 1.0 generated five solutions at the depth range of 801 m - 1433 m occurring on gravity highs as displayed in Figure 5. The figure also exhibits the deepest Euler solution occurring to the north of the study area at a depth of 1433 m on a high gravity amplitude of 2.3 mGal. These could be hot dense materials in the geometrical shape of a pipe from the mantle trapped at

Figure 4. Located Euler deconvolution results for a structural index of 0.5.

Figure 5. Located Euler deconvolution results for a structural index of 1.0.

this location in the Earth’s crust. The shallowest Euler solution occurs to the south at a depth of 801 m on a high gravity amplitude of 3.4 mGal. It is observed that the trend related to structural index of 1.0 in Figure 5, is in the NE-SW direction which depicts a pipe-like structure.

The structural index of 2.0 generated five solutions at a depth range of 1170 m - 2246 m occurring on gravity highs as displayed in Figure 6. The deepest solution occurring to the north of the study area at a depth of 2246 m is shown in Figure 6 with a high gravity amplitude of 1.8 mGal. In this case, gravity anomaly amplitude for this solution is lower than that of structural index of 1.0 and 0.5 because it is deeper and has a geometrical shape of a sphere. This increases the wavelength and reduces the frequency and amplitude of the signal. Deeper anomaly sources have a large wavelength with small signal amplitude indicating reduced energies. These could be hot dense materials in the geometrical shape of a sphere intruding into the Earth’s crust from the mantle trapped at this depth in the Earth’s crust. The shallowest solution occurs to the south of the study area at a depth of 1170 m on a high gravity amplitude of 3.4 mGal. The structural index of 2.0 in Figure 6 trending in NE-SW indicates a spherical like structure.

Figure 6. Located Euler deconvolution results for a structural index of 2.0.

7. Conclusion

Located Euler deconvolution can be used to locate and estimate depth to gravity anomaly sources within the Earth’s subsurface. Located Euler deconvolution results for structural indices of 0.5, 1.0, and 2.0 shown in Figures 4-6 respectively produced five solutions on high gravity amplitudes indicating gravity anomaly sources at the same locations. There are deeper solutions towards the northern part of the study area which becomes shallower towards the south. This was interpreted to be due to the volcanic crater at the south of the study area where dense materials emanated from deeper depths to the surface in the form of lava flow. Euler solutions on the southern part of the study area are shallow because the dense bodies from the mantle pushed to the near-surface. These near-surface hot bodies could have contributed to the existence of hot springs and fumaroles in Eburru.

Gravity solutions interpreted to be intruding hot dense materials from the mantle that were trapped towards the north occur at depths range of 1433 m - 2269 m. These intruding hot dense materials could be possible sources of heat in the system. These gravity sources could be occurring beneath an impermeable sedimentary cap rock that prevents them from reaching shallow depths. Water from Lake Elementaita possibly seeps down through subsurface fractures and faults to the geothermal heat sources. This water gets heated and migrates through geological structures and forms geothermal reservoirs if trapped underground or manifests on the surface in the form of fumaroles and hot springs.

Acknowledgements

Special thanks to the Kenya Electricity Generating Company (KenGen) for gravity data acquisition.

References

[1] Dawi, E. L., Liu, T. Y., Shi, H., & Luo, D. P. (2004). Depth Estimation of 2-D Magnetic Anomalous Sources by Using Euler Deconvolution Method. American Journal of Applied Sciences, 1, 209-214.

[2] Felipe, F., & Valeria, C. (2017). What to Expect from Euler Deconvolution Estimates for Isolated Sources. In 15th International Congress of the Brazilian Geophysical Society. Academic Scientific Conference.

[3] Georgsson, L. S. (2009). Geophysical Methods Used in Geothermal Exploration. In Short Course IV on Exploration for Geothermal Resources. Scientific Workshop.

[4] Kearey, P., Brooks, M., & Hill, I. (2002). An Introduction to Geophysical Exploration. Ames, IA: Iowa State University Press.

[5] Mariita, O. N. (2011). Application of Geophysical Methods to Geothermal Energy Exploration in Kenya. In Short Course VI on Exploration for Geothermal Resources. Scientific Industrial Workshop.

[6] Ndombi, J. M. (1981). The Structure of the Shallow Crust beneath the Olkaria Geothermal Field, Kenya, Deduced from Gravity Studies. Journal of Volcanology and Geothermal Research, 9, 237-251.

https://doi.org/10.1016/0377-0273(81)90006-8

[7] Nyakundi, E. R., Githiri, J. G., & Ambusso, W. J. (2017). Geophysical Investigation of Geothermal Potential of the Gilgil Area, Nakuru County, Kenya Using Gravity. Journal of Geology and Geophysics, 6, 2.

[8] Pawan, D., Ramprasad, T., Ramana, M. V., Desa, M., & Shailaja, B. (2007). Automatic Interpretation of Magnetic Data Using Euler Deconvolution with Nonlinear Background. Pure and Applied Geophysics, 164, 2359-2372.

https://doi.org/10.1007/s00024-007-0264-x

[9] Reid, A. B., Allsop, J. M., Granser, H., Millett, A. J., & Somerton, I. W. (1990). Magnetic Interpretation in Three Dimensions Using Euler Deconvolution. Geophysics, 55, 80-91.

https://doi.org/10.1190/1.1442774

[10] Reid, A. B., & Thurston, J. B. (2014). The Structural Index in Gravity and Magnetic Interpretation: Errors, Uses, and Abuses. Geophysics, 79, J61-J66.

https://doi.org/10.1190/geo2013-0235.1

[11] Thompson, D. T. (1982). EULDPH: A New Technique for Making Computer-Assisted Depth Estimates from Magnetic Data. Geophysics, 47, 31-37.

https://doi.org/10.1190/1.1441278

[12] World Energy Council (2013). Geothermal, World Energy Resources.