Constrained 3D gravity interface inversion for layer structures: implications for assessment of hydrocarbon sources in the Ziway-Shala Lakes basin, Central Main Ethiopian rift.

Multi layer 3D gravity inversion for layered structures and density interfaces are performed in the Central Main Ethiopian rift bounded between 38000'-39030' E and 7000'-8030' N. The inversion is carried out in wave number domain using Parker-Oldenburg algorithm and is constrained with initial model information. The previous studies in the region focused on mapping crustal structures and Moho depths and least is known about the shallow earth. This study thus targets on mapping layers relief of shallow earth origin. Stacked horizons with depth to tops of density contrast are obtained from well log data and previous geophysical studies. These stacked grids represent major geological boundaries where density contrast exists. The model utilizes observed residual gravity anomaly and generates the structural relief maps of the respective layers with their corresponding gravity anomaly responses and the associated errors. Successive structural inversions are performed on three layers with their corresponding acceptable mean misfits' errors. The iteration process converges successively for each layer in each structural inversion and the result is validated against a priori information. In addition to the topography/thickness of each layers, this study for the first time identified a new Mesozoic horizon laying between a Tertiary ignimbrite layer and the crystalline basement at depths between -2499 m and -3060 m and having estimated maximum thickness of 561 m. The identified Mesozoic sediment formation underlies a thick volcanic cover of 2.5 km which might be a suitable geologic setting for the growth of hydrocarbon reserves in the area and could probably be the source of CO2 degassing.

Introduction

The geology of an area is characterized by varied rock types and geologic structures (Mammo, 2010) (Figure 2 and Figure 3). Geophysical methods are the main tools to investigate the hidden structures indirectly based on characteristics of physical properties of rocks (Williams, 2008) and geometries of interfaces (Feng et al., 2016). The studies of crustal structures and Moho depth have been made in the Main Ethiopian rift (Mahatsente et al., 1999; Tessema and Antoine, 2004; Tiberi et al., 2005; Maguire et al., 2006; Braitenberg et al., 2006; Mickus et al, 2007; Salimi and Motlagh, 2014; Pallero et al., 2015; Ekinci et al., 2021) The methods used to map these structures include controlled-source seismic survey and the 3D gravity forwarding modelling and inversion.

Figure 2

Geology map of the Ziway-Shala lakes Basin, Central Main Ethiopian rift modified from Tefera et al. (1996) and Molin and Corti (2015).

Figure 3

Mesozoic sediment surface outcrop, distribution and orientation of the outcrop surrounding the Ziway-Shala Lakes Basin: (a) Fault map of the Ziway-Shala Lakes basin modified from Agostini et al. (2011) and Molin and Corti (2015) (b).

The Parker-Oldenburg algorithm is a multi-layer 3D gravity inversion designed to estimate mainly the variations of interfaces of density contrast within the lith

osphere, such as Moho depths or crustal structures thicknesses (Tiberi et al., 2005; Zhao et al., 2020). This method is also used in different area to estimate interface topographic variations that are responsible for the observed gravity anomaly (Salimi and Motlagh, 2012; Zhao et al., 2020; Pallero et al., 2015; Feng et al., 2016; Sabah & Al-rahim, 2018; Altinoglu et al., 2018; Kebede et al., 2021a, Kebede et al., 2021b). This algorithm allows us to calculate the three-dimensional geometry of the interfaces (layers) iteratively using an apriori initial model input information that could reduce the non-uniqueness of inversion. This constraining information is obtained from previous geological and geophysical studies in the region under study. In addition to depths to the various gravity sources, calculated and misfit gravity anomalies are the results of the inversion. The convergence of iteration and the adjusted tolerable error limit should be guaranteed in each of the inversions run.

Seismic velocity model (Maguire et al., 2006) and 2.5 D gravity model (Mickus et al., 2007) were the two prominent surveys conducted under the Ethiopia Afar Geo-scientific Lithospheric Experiment (EAGLE) project. These studies were designed to map crustal structure and Moho depths along the rift axis and over the rift floor. Considering that the depth of the seismic and gravity model is 50 km and 120 km, respectively and limited vertical resolution, the fine-scaled structure of the shallow crust remains unclear.

In this study, ground based gravity data in the region are analyzed to understand layers structures/geometries (density contrast interfaces) of the shallow earth origin in the Ziway-Shala Lakes basin. To accomplish this task, 3D iterative structural inversion algorithm of Parker-Oldenburg (Oldenburg, 1974; Parker, 1973) is used. Furthermore, significance (if any) of the mapped structures for exploration of hydrocarbon resources in the region is considered.

Location of the study area

The Central Main Ethiopian rift (CMER) is part of Main Ethiopian rift bounded between 38000′E–39030′ E and 7000′N–8030′ N. The study area is characterized by the near NE–SW trending rifted basin, flanked by the uplifted blocks on both sides.

The Ziway-Shala Lakes Basin (CMER) encompasses four main residual rift floor lakes (Le Turdu et al., 1999) respectively from north to south listed as Ziway, Abiyata, Langano and Shala lakes. This basin is hydrologically closed (Chernet et al., 2001) with its elevation varying from shallowest 1342 m above sea level to highest 4168 m above sea level (Figure 1).

Figure 1

Location (UTM) and topographic map of the Central Main Ethiopian Rift and its surroundings; the 3D view shows the geomorphologic features of the region.

Geological and structural settings

The Main Ethiopian Rift (MER) where the study area lays was developed over a span of time from the Oligocene to the Quaternary (Le Turdu et al., 1999; Woldegabriel et al., 2000). The current geologic and geomorphic features observed in the MER are as a result of geologic processes which include tectonism, volcanism and sedimentation activities. These processes in the MER in general and CMER in particular resulted in observable geologic structures like faults (Woldegabriel et al., 2000), joints, contacts and fractures. These structures have surface expressions shown in the geologic map and structural map. These surface structures mainly include the Wonji Fault Belt (WFB) (Mohor, 1962), Silti Debrezit Fault zone (SDFZ) and their associated boundary faults (Boccaletti et al., 1998). They generally have N–S, NNE–SSW, NE–SW, E–W and NW–SE (Korme et al., 2004) orientations (Figure 3(b)). These structures are considered as young and active fault systems and are cross-cut by pre-existing NW–SE Mesozoic Ogaden rift fault (Korme et al., 2004; Kebede et al., 2021a, Kebede et al., 2021b) (Figure 3(a)).

Geologic sections within the CMER are poorly exposed (Woldegabriel et al., 2000). Woldegabriel et al. (2000), however, identified localized crystalline basement, Mesozoic sedimentary rocks, and Miocene and Pliocene mafic and silicic lavas and tephra are exposed on the western rift margin of the CMER. Woldegabriel et al. (2000) further identified crystalline basement, Mesozoic sandstone, limestone and mudstone on the north-eastern part of the CMER (along the southern Afar margin). The identified rocks have contributed much in the accumulation of sedimentary successions over the rift floor (Figure 2).

Mesozoic succession is present in the Blue Nile Basin which formed a NW-trending rift (Gani et al., 2008), over the eastern margin of the Ethiopian Rift adjacent to the Ogaden Basin, over the southern border of the Ethiopian Rift around Dire Dawa town (Bosellini et al., 2001) and at Kella along the western margin of the Central Main Ethiopian Rift (MER) (Abebe et al., 2010).

The Mesozoic outcrop is observed to surround the study region of the Ziway-Shala Lakes Basin (or CMER) (Figure 3(a)). Therefore, one could not deny the existence of Pre-Tertiary units (Mesozoic and Precambrian) at specified depth in the Central Main Ethiopian rift. Consequently, the conceptual model referring to these strata was constructed by Woldegabriel et al. (2000) and is schematically shown in Figure 4.

Figure 4

Schematic representation of stratified earth in the CMER adopted from (Woldegabriel et al., 2000).

Though volcanic, tectonic and sedimentation processes distort the strata and faulting characterizes the region under consideration, the strata corresponding to each fault block from top to bottom were approximated by Woldegabriel et al. (2000) as volcano-classics, volcanic, sedimentary and crystalline basement formations (Figure 4).

Similarly, a West-East oriented geological section was constructed over the Aluto-Langano geothermal field using Well-log data (core samples) information (Abera and Mizunaga, 2018; Wilks et al., 2017; Hutchison et al., 2016). The resulting layers’ stratigraphy (Figure 5) from top to bottom were identified as pyroclastic, Lacustrine sediments, Bofa basalt, tuff and breccias, and Tertiary Ignimbrite formation layers.

Figure 5

West-east geological cross-sections derived from core samples drilled in the Aluto-Langano geothermal field adopted from (Cherkose and Mizunaga, 2018). Depths of the eight drilled wells range from 1300 to 2500 m below ground surface.

The stratified grids defining the layers, the type of rocks in each layer for estimating layer density, depths to interfaces are extracted from the geologic section produced using the well-log data (Figure 5). These parameters are systematically put together to define an initial model information for gravity structural inversion given in section 2.2.

Methodology and data

Governing equation

The basis on which the gravity method depends (Telford et al., 1990) is summarized in two laws derived by Newton second law of motionand universal Law of gravitation

In Cartesian coordinates, , describes the mutual force between two point masses, m and at two distinct points.

Equating Eqs. (1) and (2), we have

Extending the point masses to interfaces of constant density, the gravitational potential resulting from these layer of constant density contrast, bounded below by the horizontal plane and above by the surface at a position, , is written as (Eq. 4) (Nagendra et al., 1996)Where,

z is positive upwards.

G is universal gravitational constant, and.

is density contrast.

The gravity anomaly caused by an interface with a constant density contrast is calculated by classic Parker’s formula, shown in equ. 5 (Parker, 1973).

Oldenburg (1974) rearranged this equation to compute the depth to the undulating interface from the gravity anomaly by means of one dimensional Fourier transform in an iterative equation (Eq. 6).Where:

is the Fourier transform of the gravity anomaly.

G is Newton’s gravitational constant, 6.67 × 10−11Nm2kg.−2

r is a horizontal plane.

is the density contrast across the interface (two media),

k is the wave number,

n is an integer

h(r) is the topography of the interface

z0 is the mean (reference) depth of the horizontal interface or observation plane.

is Fourier transform of the interface topography.

Parker–Oldenburg iterative expression (Eq. 6) allows us to calculate three-dimensional geometry of the interface (layer) iteratively using interface density, , and mean depth, z0 The first term of Eq. (6) is computed by assigning h(r) = 0. This value, h(r), is then used in the equation to evaluate a new estimate of h(r+1). This convergent process is stopped until a reasonable solution is achieved.

The inversion operation (Eq. 6) is unstable at high frequencies and thus it requires a high-cut filter (low-pass) in the frequency domain for the convergence of iteration (Eq. 6) to be guaranteed. Oldenburg (1974) proposed a low-pass filter of the form (Eq. 7)

The inversion is then governed by three parameters, density contrast, , datum level, z0, and the filter B(k). Existence of these free parameters results in non-uniqueness of the gravity inversion. With the adjustment of the filter parameters, WH and SH, however, we are guaranteed the convergence of the iterative process. According to Oldenburg (1974), the error norms and convergence criterion respectively given asWhere, E is some sufficiently small predefined value.

The convergence criterion (Eq. 8) set by Oldenburg (1974) for such fixed point iterative problems (Eq. 6) is different from the convergence criterion that exists on any Numerical methods text book.

Constraints for gravity inversion

A priori information on the subsurface structure is required when using Parker-Oldenburg algorithm. This extra information includes the layer numbers, their mean depths and density contrast. The identified parameters are used as a starting model for determining the calculated and misfit gravity anomalies through Forward modeling.

The average crystalline basement top estimated from gravity data analysis is 3.0 km deep (Kebede et al., 2020). Following the evidence documented in section 1.2, Mesozoic layer is identified as fifth layer below Tertiary ignimbrite and above crystalline basement layer. Here, however, there is a need to approximate the depth to the top of the Mesozoic layer.

Multi-Layer 3D Gravity forward calculation is performed by choosing different depths to top of the Mesozoic layer (Table 1, Column 6) and sensitivity of the changes to the top depth of this layer is documented. The best-fitting depth to Mesozoic layer is 2.9 km by the error values in Table 1.

Table 1

Sensitivity analysis conducted to determine the depth to top of Mesozoic layer through Multi-Layer 3D Gravity forward calculation. Start the invasion at a shallow depth of 2.6 km then follow 2.7 km, 2.8 km and 2.9 km. The least mean square error was found to be at a depth of 2.9 km.

SN	Minimum	Maximum	Mean	Std. dev	Depth to top of Mesozoic layer (km)	Remark
	in m
1	-16.4676	38.79607	6.00792	6.42845	2.6
2	-16.8031	38.46064	5.67248	6.42845	2.7
3	-17.1385	38.12521	5.33705	6.42845	2.8
4	-17.4739	37.78977	5.00162	6.42845	2.9	Chosen depth

Based on well-log data in the geological section (Cherkose and Mizunaga, 2018) and gravity forward modeling conducted in this research, six layers with depths to top and density are approximated. The initial model is defined by a number of stacked surface layers (Figure 6) with different density contrasts listed in Table 2.

Figure 6

Initial model created from the available geological and geophysical data with different colors representing stacked horizons where density contrasts occur. The top most layers represent the elevation (DEMs) of the study area (CMER).

Table 2

The depth and density from top to bottom estimated from Well-log data, seismic study, gravity interpretation result and previous geological studies.

Layers	Layer Geology	Sub Layers Geology	Depth to top of the layer (m)	Initial Density for the interfaces (g/cm3)	The source of the information obtained from
1	Pyroclastic and Lava flows (Silicic products)	Per alkaline rhyolite lava (47–201)_ 2.6 g/cm3	0	2.477	(Alemu, 1992)
		silicic breccias and weakly welded tuffs (201–604)_ 2.43 g/cm3
2	Lacustrine sediments		1000	2.34
3	Bofa basalt, tuff and breccias		1500	2.81	(Alemu, 1992)
4	Tertiary Ignimbrite		2500	2.58	(Alemu, 1992)
5	Mesozoic sediments		2900	2.5	(Mammo, 2010)
6	Crystalline basement		3000	2.74	(Kebede et al., 2020)

These stacked grids represent major geological boundaries where density contrasts occur.

The overall 3D forward modeling (Parker, 1973) and 3D interface inversion (Oldenburg, 1974) procedures are summarized in the following flow chart (Figure 7).

Figure 7

Summary of the 3D forward modeling and Inversion process.

The decision symbol of the flowchart (Figure 7) helps to check whether the initial model perfectly suits our predefined error tolerance (Eg. 0.01 mGal). If the initial model doesn’t fulfill the criterion, we work iteratively on setting up a new initial value of the parameters. Layer geometries, density and depth (topography) of the interfaces are some of the parameters that need to be improved through repeated refinement and forward modeling. The second criterion written at bottom in the decision symbol (Figure 7) helps to decide whether the 3D gravity interface structural inversion result is acceptable or not. This is based on a misfit (error) anomaly map computed by taking the difference of observed and calculated anomalies. Based on governing equation (Eq. 6), 3D interface inversion of residual gravity anomalies over the Ziway-Shala lakes basin are performed to estimate layers’ undulations of crystalline basement, Mesozoic sediment and Tertiary Ignimbrite formation layers.

In addition to structural inversion, single density interface inversions are performed to estimate the density of layers after structural inversion run. This inversion calculates the error grid cell by cell and does an estimate of the density for all cells, then averages them to get an overall density change.

Gravity data

About 3013 ground based secondary gravity data (Figure 8(a)) were reprocessed, homogenized to the International Gravity Standardization Network 1971 (IGSN71) and corrected (reduced) to produce the Bouguer anomaly (Figure 8(b)). The upward continuation height of 6.0 km (Kebede et al., 2020) is used to estimate the regional anomaly (Figure 8(c)) and this anomaly is subtracted from Bouguer anomaly to have the residual anomalies of the region under study (Figure 8(d)).

Figure 8

Gravity stations distribution map (a) Bouguer anomaly map (b) regional anomaly map (c) and residual anomaly map (d).

Jacobsen (1987) stated that when the potential field data is upward continued to a height 2''; it maps the sources found at and below the depth ''. The residual gravity anomalies depicted in Figure 8(d) has revealed gravity source bodies located to depths of 3 km. It is these anomalies that are subjected to structural inversion for layers geometries in the region under study.

Results

Once the initial model is compiled as documented in section 2.2, the multi layer, surface based frequency domain forward and inverse modeling are carried out using Parker-Oldenburg algorithm. A top most layer in Figure 9 is the residual gravity anomalies inverted for interface geometries.

Figure 9

Six geological horizons stacked and configured at depths to top of 0 m, 1000 m, 1500 m, 2500 m, 2900 m and 3000 m depth and corresponding layers densities 2.477 g/cm, 2.34 g/cm, 2.81 g/cm, 2.58 g/cm, 2.5 g/cmand 2.74 g/cm are identified from well-log data, previous geological and geophysical analysis in central Main Ethiopian rift.

Inversion on crystalline basement

We start the inversion with the deepest interface (crystalline basement) formation layer with the aim of fitting the calculated and the observed residual anomalies. Correspondingly, after an inversion the original layer geometry is altered. The gravity data analysis in the study area shows that the crystalline basement horizon is found at an estimated mean depth of 3.0 km (Kebede et al., 2020). The reference (initial) depth of the horizon is then set at –3 km and layer density set to 2.74 g/cm3. Top of this layer is Mesozoic sediment with a density of 2.5 g/cm3. The average density contrast is thus 0.24 g/cm3. The smallest and greatest cut-off frequency parameters are chosen as 0.5 and 0.7 km−1, respectively. The maximum inversion iteration number is set to 20 and convergence limit is set as 0.01mGal. Since the objective of the study focuses on shallow earth, the residual gravity anomaly is inverted for layers undulation and thickness. The result of which given in (Figure 10(a–d)).

Figure 10

Maps of the residual gravity anomaly (a) calculated gravity anomaly due to topography of crystalline basement (b) and the difference between observed and calculated gravity anomaly (c). The basement undulation map (crystalline basement-Mesozoic layer interface) obtained through application of structural inversion on crystalline basement layer and using residual gravity anomaly (d). The coloring system is Histogram equalization.

The structural inversion on this layer modified the initial constant layer and generated crystalline basement topography of the area (Figure 10(d)). The inversion converges after 12 maximum iterations run. In this modeling, the minimum error is –12.05387 m, maximum error is 37.12428 m, and means square error is 5.74523 m and standard deviation 5.46933 m.

The calculated crystalline basement topography (Figure 10(d)) shows that the layer’s shallowest depth starts at a depth of 2899 m and extends down from this depth. The basement topography gets deeper at some specific area near the northern water divide (Meki area and its surroundings), floor of Gademota caldera and N to NNW of Abiyata Lake. Here only the anomalies within the rift floor get interpreted to avoid edge effects caused by the calculation process.

Structural inversion on mesozoic surface

Initially before structural inversion carried on the Mesozoic surface, the depth of this horizon was set at –2.9 km and the layer density was set to 2.5 g/cm3. Top of this layer is Tertiary ignimbrite layer with a density of 2.58 g/cm3. The average density contrast is thus 0.08 g/cm3. The smallest and greatest cut-off frequency parameters are chosen as 0.5 and 0.7 km−1, respectively. These high-cut filter parameters should guarantee the convergence of the iterative scheme. The convergence limit is set at 0.01 mGal. The maximum inversion iteration number is set to 20. The iteration converges after an iteration run of 2 with root mean square error of 5.80728 m and standard deviation of 5.43950 m. The inversion result (Figure 11a, b) shows that the shallowest of Mesozoic sediments was found at around 2499m which is –401m from the reference depth of 2900m and extends down to –3060m. The maximum thickness of this layer is thus 561m. The sediment is found at the rift floor near the Meki area extending south-west to Gademota Caldera and to north of Abiyata lake.

Figure 11

Mesozoic layer topography in 2D view (a) and 3D view (b).

Structural inversion on tertiary ignimbrite layer

To run the structural inversion on the Tertiary ignimbrite layer, the depth of the horizon was set at –2.5 km and the layer density set to 2.58 g/cm3. Top of this layer is the Bofa Basalt layer with a density of 2.81 g/cm3. The average density difference put into calculations is considered to be 0.23 g/cm3. The cut-off frequency, the maximum iteration number and convergence limit is similar to the previous inversion setup. After structural inversion, the layer extends up –1608 m from the reference depth –2500 m and moves down –2900 m from this depth. The iteration converges at an iteration number 9. The RMS error is 3.82638 m and standard deviation is 4.65044 m. The maximum thickness of this layer is 1292 m. The shallowest depth is found at the rift floor, near Meki extending south-west to Gademota Caldera and to north of Abiyata lake (Figure 12). Thick Tertiary Ignimbrite is correspondingly found at the mentioned locations.

Figure 12

Relief map of Tertiary Ignimbrite Formation layer.

Discussions

Layers geometry (structures) through structural inversion

The constraining parameters used in gravity interface inversion are compiled from previous geological and geophysical studies carried out in and around Central main Ethiopian Rift. This information helps to generate a better inversion result within the set tolerable error limit. A priori information defined in section 2.2 helps us only to define constant density for each layer with corresponding single density values. Through the cases, the layers’ geometries are changed as an inversion iteration progresses. Inversion is carried out first on the crystalline basement layer, secondly on the Mesozoic layer and finally on the Tertiary ignimbrite layer. Summary of the inversion results are described in Table 3.

Table 3

Successive Structural Inversions statistics on stated layers.

Inv. N	Minimum	Maximum	RMS error	Std. dev	Inversion surface/layer	Total _ iteration for Convergence	Minimum value (m)	depth to top of the layer (m)	Shallowest depth (m)	Estimated max. layer thickness (m)
	in meter
1	-12.05	37.12	5.75	5.47	Crystalline basement	12	-	-3000	-2899
2	-11.87	37.14	5.81	5.44	Mesozoic layer	2	-3060	-2900	-2499	-561
3	-12.61	33.76	3.83	4.65	Tertiary Ignimbrite	9	-2900	-2500	-1608	-1292

In all inversion runs the iteration converges and mean square error ranges from 3.83 to 5.75 m. The estimated layers thickness (Table 3) for each of the three inversion layers and their undulations give a good first approximation to the topographic map of the layers.

The inversion on crystalline basement layer modifies the geometry of the layer (grid) and produces the undulation of the surface with estimated shallow depth of 2899 m. This depth is measured from the surface and extends down from this depth. The inversion on the Mesozoic layer produces the shallow topographic depth of 2499 m and the deepest depth of 3060 m. The estimated maximum layer thickness for the Mesozoic sedimentary formation is 561 m. The inversion carried on the Ignimbrite formation layer resulted in a shallow top of 1608 m and deepest depth of 2900 m. The thickness of the Ignimbrite formation layer is then 1292 m. Comparison of this layer thickness with Tertiary ignimbrite layer geologic section thickness (Figure 5) are made and compared. The comparison is found to be in a good agreement.

The inversion on the remaining top three layers such as Pyroclastic and Lava flows (Silicic products), Lacustrine sediments and Bofa basalt, tuff and breccias formation layers gives noised result and should be excluded from interpretation.

Implications for hydrocarbon potential

Characterizing the geology of an area is an important task to understand the role that geology plays in the growth of petroleum reserves. Gravity method is one of the geophysical methods used to delineate oil and gas-bearing sedimentary basins and locate main trap structures (Zhang et al., 2019). These generally include locating source rocks, reservoir rocks; trap (seal) and the wider basin structures that are responsible for accumulation of hydrocarbon resources.

The hydrocarbon source rocks potential investigations are the main parameter needed to be studied for oil and gas resources. These types of rocks are composed of very small mineral fragments in between them are remains of organic material. This study identified a Mesozoic layer which lay in between the crystalline basement layer and the Tertiary ignimbrite (sedimentary) layer. The identified Mesozoic interface was believed to be formed as a result of transgressive-regressive processes in the Mesozoic era.

The estimated thickness of volcanic overburden including Tertiary sedimentary is 2.5 km (Figure 5). Growth of petroleum reserves needs such a kind of geology where the sedimentary basin is overlain by volcanic rocks (Mammo, 2010). For example, in Yerer Tulu Welel Volcanic Lineament zone the presence of volcanic activity and low enthalpy geothermal systems lead to the heating of the underlying Mesozoic sediments which then liberate a considerable amount of carbon dioxide gas and other gases (nitrogen, sulphur and argon) (S. Kebede, 2013).

The identified Mesozoic and Tertiary sedimentary formation thickness is estimated to –1853 m. The sedimentary rocks in Mesozoic formation possibly help the growth of hydrocarbon resources in the region considered. It is observed that temperature in the region increases with depth (Cherkose and Mizunaga, 2018) (Figure 5), when this temperature exceeds 120 °C the organic remains within the rocks begin to be “cooked” and as a result of which oil and natural gas are formed. These suggest the presence of organically rich sedimentary rocks (source rocks).

In the study region, the volcanic overburden claimed to turn the sediments into hard rock strata with further burial guaranteeing temperature increase. Furthermore, the region is dominated by active volcanism and tectonics which is manifested by hydrothermal fluids and fumaroles. These fluids are important in mobilizing and moving hydrocarbons (Schutter, 2003) to reservoirs (Figure 13). The hydrocarbons matured in Mesozoic sedimentary rocks and can migrate vertically or laterally into structurally higher igneous rocks (tertiary ignimbrite). These rocks are considered to be porous and permeable and might act as good hydrocarbon reservoirs.

Figure 13

Theoretical hydrocarbon system model based on elements commonly found in rift basins taken from (Schutter, 2003).

A Lacustrine sediment layer exists above the Tertiary ignimbrite layer at top depth of 1000 m (Figure 5) which might act as a seal or trapping structures.

Hutchison et al. (2015) estimated total carbon dioxide (CO2) emitted from Aluto volcanic center using CO2 degassing to be in between 250–500 t d−1. Seismicity study in the same area was conducted to map magmatism, faulting and hydrothermal circulation (Wilks et al., 2017). Both studies have shown existence of magma storage at about 5 km depth. In contrast, Magneto-telluric (MT) study conducted at Aluto volcano (Samrock et al., 2015) showed that there are no magmatic sources directly under Aluto. The other study that confirmed Samrock et al. (2015) result was conducted by Hübert et al. (2018), which evidenced the non-existence of conductive regions at 5 km depth.

The previous two studies (Wilks et al., 2017) (Hutchison et al., 2015) are against the later two studies (Samrock et al., 2015) (Hübert et al., 2018) on the location of magma storage. Magnetic data analysis (Kebede and Alemu, 2020) conducted in the region found to support the results of the latter two findings. If Aluto volcano is deprived of heat (magma) source beneath, one can ask the source of CO2 degassing.

The possible source of this gas is claimed to be as a result of combustion of fossil fuels (coal, oil, and natural gas) or magma storage. In North and North-West of the Central Main Ethiopian Rift for example the CO2 gas emanating from the Mesozoic sediment de-carbonation (S. Kebede, 2013). This has led to the formation of numerous naturally sparkling springs, notably around Ambo, Woliso, Diddessa Valley, South of Lake Tana, and Filwuha thermal springs around Addis Ababa (S. Kebede, 2013).

This study for the first time identifies Mesozoic formation in the region considered which could possibly be the source of carbon dioxide gas. Furthermore, this study suggests that the plausible source is not related to magma chamber.

In addition, from gravity data analysis conducted in this study, one can observe nearly circular shaped geometries (lineaments) in each structural inversion run which might indicate remains of an old caldera structure nesting smaller ones. The presence of such circular gravity lineaments around positive anomalies indicate shape of concentric faults (Naouali et al., 2016) edging Tertiary ignimbrite and Mesozoic sediments in the region considered. The identified old and big caldera-like structural features hosts (nests) the smaller and younger calderas called Aluto and Gadamota calderas.

The Mesozoic layer structure/geometries mapped in this study are based only on the input constraining information (Table 2) and doesn’t reflect the Mesozoic source at any other layers and outcrops at the surface. The resource exploration however needs detailed investigations on how the geology of the area significantly controls the growth of hydrocarbon (petroleum) reserves with more refined data sources. The existence of a Mesozoic layer in the basin leads the researcher to suspect hydrocarbon resources though detail investigation is required.

Constant density inversion after geometric change of layers

Constant density inversion calculates the error grid cell by cell and does an estimate of the density for all cells, then averages them to get an overall density change. The procedure is based on Parker’s algorithm and with an initial layer density for the interfaces (Table 2) and layers configuration (Figure 9). To compute a single density for an interface we first estimate layers’ horizons relief (topography) as shown in Figure 10(d), Figure 11(a, b) and Figure 12 using structural inversion. The computed and optimized constant density of each selected model layer is given in Table 4.

Table 4

Interface geology, depth to top, Initial density and optimized density columns. Single density inversion is run after structural inversion geometric modifications.

Layers	Layer Geology	Sub Layers Geology	Depth to top of the layer (m)	Initial Density for the interfaces (g/cm3)	single density inversion after structural inversion (g/cm3)
1	Pyroclastic and Lava flows (Silicic products)	Per alkaline rhyolite lava (47–201)_2.6 g/cm3	0	2.477
		silicic breccias and weakly welded tuffs (201–604)_ 2.43 g/cm3			2.4983898748
2	Lacustrine sediments		1000	2.34	2.3861097941
3	Basalt, tuff and breccias		1500	2.81	2.8280985006
4	Tertiary Ignimbrite		2500	2.58	2.5816582477
5	Mesozoic sediments		2900	2.5	2.4892385754
6	Crystalline basement		3000	2.74	2.7538995576

The recalculated gravity anomalies and the misfit (error) gravity anomalies after constant density inversion are shown respectively in Figure 14(a) and (b).

Figure 14

Multi-Layer 3D Gravity Forward Calculation of calculated gravity anomaly (a) and the misfit gravity anomaly (b) using optimized constant density of each stacked layer (Table 4).

Before geometric and constant density inversion the ranges of observed residual gravity anomalies lies between –11.11mGal to 11.29 mGal. The density and geometric change optimizes the calculated gravity anomalies to lie in between –7.32 mGal and 8.61mGal and the misfit anomaly lies in between –5.04 mGal and 3.2 mGal.

Conclusion and recommendation

Shallow crust structures of the Ziway-Shala lakes basin are investigated using 3D structural inversion of Parker-Oldenburg iterative program. The algorithm calculates the subsurface topography resulting from density contrasts between interfaces. In each step the errors between the observed and calculated anomalies are computed at all points by iteratively improving the misfit. The study initially approximated the layers as stratified strata though tectonic, volcanic and sedimentation processes distort the strata and faulting are the main characteristics properties of the region under consideration. 3D structural inversion of residual gravity anomaly is performed with parameters of layers of information obtained from previous geological and geophysical studies. Six layers are identified, from top to bottom includes Pyroclastic and volcanic rocks formation, Lacustrine sediment formation; Bofa baslt, tuff and breccias formation; Tertiary ignimbrite formation; Mesozoic sediment formation and crystalline basement formation layers. The density contrast is defined by defining the density to each six horizons. The geometries of the last three layers (relief maps) are produced in each inversion run (Figure 9(d), Figure 10 and Figure 11a,b). The layers thicknesses are validated using well-log data drilled for the purpose of geothermal study. This study reveals the presence of a Mesozoic layer starting at a shallow depth of –2499 m extending to a deeper depth of –3060 m and having a maximum thickness of –561 m. Beneath the Mesozoic formations, a crystalline basement layer of shallowest depth starts at an approximate depth of –2899 m and extends downwards. The inversion results give fair first approximations to the topographic map (lithostratigraphy) of the layers in the region considered. The fact that the basin is a closed sedimentation process characterizes the region under consideration. The research has also identified the presence of Mesozoic sediments underlying a thick volcanic cover of thickness 2.5 km. This mapped stratified geology underlying the sedimentary basin and overlain by volcanic rocks might be favorable for the growth of hydrocarbon deposits in the area. The researcher thus claims the possible existence of hydrocarbon resources in the area. From previous observations made in area the source of CO2 gas could not be related to a magma chamber as there is no magma under the Aluto area. Therefore, the source of CO2 could probably be due to metamorphic de-carbonation of carbonate rocks of Mesozoic origin.

Parker–Oldenburg algorithms used in this study assumed either a constant or a laterally varying density for layers involved. In reality, densities vary both vertically and laterally for a given strata. The results obtained are approximation and crude generalization of interfaces/layers topography. Thus, it is much required to work on gridded potential field data either air borne or ground based for a better geologic and structural mapping. Furthermore, it is required to update the results obtained by incorporating more refined initial models whenever information about the subsurface is required.

Declarations

Author contribution statement

Hailemichael Kebede: Conceived and designed the experiments; Performed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Abera Alemu: Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Funding statement

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Data availability statement

Data will be made available on request.

Declaration of interest’s statement

The authors declare no conflict of interest.

Additional information

No additional information is available for this paper.