A new approach automatic separation of the Bouguer gravity anomaly, using a new concept for 2D-semi-inversion of the sphere-shaped model

ABSTRACT The present research represents a new, fast and easy-to-apply semi-inversion technique for the direct separation of the Bouguer gravity anomaly to its corresponding depths or thicknesses of the prior known rock formations from a borehole (as a control point). As well as the possibility of tracing these formations on the profile points to the Bouguer gravity map according to their density contrasts with basement rock. The proposed method is based on the fact that the sum of the gravitational effect at a point on the Earth’s surface is equal to the sum of the subsurface gravity effects arising from the points of masses of the causal bodies along the vertical line between the Earth’s surface and the basement rocks. In this method, two sedimentary basin models were built based on the distribution of densities from a prior known borehole. Assuming that the depths of the pointed tops of these rock formations can be replaced by the centres of sphere bodies (or point masses) as causative sources of gravitational effect and by the use of a simple algorithm, it was possible to calculate the depths of the rock formations and trace them on the profile, using a new concept called the zero-offset gravity measurement. In this concept, the measuring gravity is vertically above the centre of the sphere-shaped body as it represents the causative source, and its depths are equal to the radius of the sphere body. The present method was assessed on hypothetical models and synthetic data and applied to two real data in field cases which vary in geological and lithological aspects. The first assessed location is Abu Roash Dome Area, southwest Cairo, Egypt, and the second location is the Mors Salt Dome, North Jutland, Denmark.


Introduction
The gravity method depends mainly on the measuring of variations in the Earth's gravitational field resulting from the variance in density of the subsurface rock formations and/or their structure deformations, where all gravity anomalies come from horizontal variations in density (Nettleton 1971). Therefore, there would be no gravity anomalies if the Earth's materials were layered or deposited horizontally in uniform density, whatever may be the vertical variation in density. Any geologic condition that results in a horizontal variation in density will cause a horizontal variation in gravity or gravity 'anomaly', as shown in Figure 1, where the flat layers are disturbed by a structural uplift (Sharma 1976). Therefore, the density contrasts between subsurface rock formations should be estimated before postulating their structures and/or lithologic distributions.
The main aim of interpretation of gravity data is essentially the estimation of the location, shape and depth of the causative sources (Hajian et al., 2012). Therefore, the present research is based on the idea that the observed value of the Bouguer gravity anomalies at each isolated point of the profile's line can be replaced with their equivalent sum of the calculated gravity effects of the point masses (as causative point sources), representing top points of rock formations that prior know their depths, thicknesses and densities from the borehole (as a control point). But, the density contrasts used in this calculation are the vertical density contrasts between each rock formation and the basement rocks, not with the surrounding rocks as in the classical calculation of the gravity effect at the point mass (sphere centre) model. The new technique proved its high sensitivity to the relative vertical changes in the densities or the vertical density contrasts of subsurface rocks' formations and the basement rocks or (of the Earth's deepest layer of crust).
Inversion of gravity data is non-unique in the sense that the observed gravity anomalies in the plane of observation can be explained by a variety of density distributions (Essa 2012). This arise the most challenging problem of ambiguity, for interpreting the gravitational field data, which is still facing the researchers, where the modelling of potential field data is considered to be a non-linear problem. However, a unique solution may be found, when assigning a simple geometrical shape to the causative body (Salem et al., 2002). Fortunately, almost most of the geological structures can be approximated, by one or more of the available simple geometrical shape models, to represent the causative sources for gravity anomalies.
The present research has arrived at a unique solution for Bouguer gravity anomalies and contributed to reducing as possible the ambiguity problem in interpreting the gravity data. This was carried out by building two deposition basin models using prior known rock formations' depths and densities in a borehole (as a controlling point) penetrating the being investigated area. Along the profile of Bouguer's gravity points, data was inverted to their equivalent rock formations' depths, as well as tracing them from point to point, using a simple algorithm of a sphere-shaped model.

Materials and methodology
The end product of a gravity field survey (after applying all the necessary corrections) is usually a contoured anomaly map called Bouguer map (Sharma 1976). 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 higher-density medium should result in a negative Bouguer anomaly (Reynolds 1997).
In the present research, Microsoft Excel software, Surfer-15 Golden software, and a simple algorithm script written in MATLAB environment were used to implement the technique of automatic separation of the Bouguer gravity anomaly map. The gravity data of the map was in the form of a profile. The point values of this profile then were separated into the equivalent of the point masses of the gravity effects, where those point masses were assumed to represent causative sources of corresponding top points of the surfaces of rock formations.

Theory equations of sphere-shaped body model
The sphere-shaped body (point mass) model was used by several researchers and authors as a conventional method for estimation inversion parameters such as the depth of some buried geological structural features, where the individual gravity anomaly can be inverted as a single isolated body (e.g. Siegel et al., 1957;Nettleton,1962Nettleton, & 1976Gupta 1983;Abdelrahman et al.,1993;Reynolds 1997;Salem et al., 2004;Asfahani et al., 2008;Essa 2012;Mehanee 2014).
The gravity effects resulting from the sphere-shaped body (point masses) model more-or-less fit the gravity effects resulting from the causative point sources, which represent the rock formations' top surface points. In this case, the inversion process is possible since the rock formations' top surfaces are of prior known depths as also their densities from the borehole (as control point). Thus, in the present semi-inversion technique, it is alsopossible to determine and trace the depths of the point mass sources representing the rock formations' tops according to their physical properties (densities).

The gravity effect of sphere-shaped body (or point mass)
The gravitational effect of buried sediment of mass (M) in a sphere shape (Figure 2), with density distribution (ρ), is a function of radius (R) alone, is identical with that due to a point mass (M) situated at the centre of the sphere (O). The mass of the sphere (its volume times its density) is given as follows: Thus, the force at the point (C) on the Earth's surface, according to Newton's Law of attraction, is given by equation as follows: where (r) is the distance between the centre of the sphere and the observed point on the Earth's surface (datum 1) and given by the following equation: and (G) is the universal gravitational constant. Since most commonly used gravimeters are restricted to detecting only the vertical component of acceleration (gz), which is given as follows: where (θ) is the angle of gravity effect with vertical axis (depth) and given as: by substituting equations (1)-(3) in equation (4), the gravity effect of the sphere at point (C) due to vertical component (where θ ¼ 90 0 ) is given as follows: where (Δρ) in this research represents the density contrast of the isolated sphere (or point mass) and the basement rocks (∆ρ = ρ sphere À ρ basement ), without respecting density contrast with the surrounding rocks, since the measurements are dealing with vertical components.

The concept of zero-offset gravity measurement (ZOGM)
Assuming (z) as the depth from the observed point to the centre of the sphere or the point mass equals to the radius of the sphere (R = z, as first constraint), and the gravity effect is measured directly at a point vertically above the centre of the sphere (x = 0, as second constraint). So, if the previous two constraints are satisfied at any observed measured point, then it is called ZOGM, at which it gives the maximum gravity effect value. Therefore, assuming that each of the measured, observed gravity effects is considered to be a single isolated maximum value, hence equation (6) can be modified and written as follows: Equation (7) can be also inverted to obtain the vertical depth from the surface of the measured datum to the centre of the sphere, by prior known density contrast and the measured gravity anomaly at that point as follows: The absolute value is taken for equation (8), to avoid the negative values of either density contrasts (Δρ) or maximum gravity value (g max ).
Equation (8) as seen has three variables which are the measured depth to the centre of the sphere (z); the observed gravity, which is equivalent to the maximum gravity effect (g_max) and the density contrast with basement rock (∆ρ). Thus, the depth from the Earth's surface to the centre of the sphere, where the latter datum representing a point on the assumed top surface of any rock formation layer, can be calculated by prior known gravity effect (from gravity anomaly's profile) and its density contrast (from a borehole, etc.).

The semi-inversion technique
The mathematical treatments of the inversion process that were previously carried out using the concept of ZOGM are called in the current research 'semiinversion'. Because of neglecting the entire shape of the sources causing the gravitational effects. In addition, the application of both forward modelling and inverse modelling concepts of the present technique is used for calculating and tracking the depth of rock formations for the profile points, according to the density contrasts between rock formations and the basement rock, available from a prior known borehole. As a result of using the density distribution model for rock formations, it is obvious that the constraints in equation (8) are satisfied for implementing the semi-inversion process within the range of density contrasts of the formations' rocks with the base rocks.

Deposition basin models of rock density distributions
To implement the new semi-inversion method, two hypothetical deposition basin models are being built up and used. Whereas the rock formation depths and density distributions inside the deposition basin models are constrained by the borehole data of the area being considered to be as follows:

Heterogeneous densities distribution model (model 1)
According to borehole data, the rock formations in the deposition basin model are distributed heterogeneously, based on the Walther's Law (Figure 3).

Homogeneous densities distribution model (model 2)
According to the same previous borehole data, the rock formations in the deposition basin model are distributed homogeneously, based on the Steno's Superposition Law or Deposition History Principle ( Figure 4). This model also agrees with well-known Bouguer Slab Model, of correcting gravitational data.

Testing of ZOGM Concept of hypothetical sphere-shaped body model
The reason for creating two deposition basin models is to prove that the calculated gravitational effects of the sphere-shaped body model using the ZOGM concept at any point on the Earth's surface give the same value when using either of the two models for the hypothetical depositional distribution, as long as using the same restrictions (x = 0 and R = z). Since the calculated gravitational effects depend only on the vertical density contrast of each isolated rock formation and the basement rock, it is possible to perform the proposed semi-inversion process for each of the profile points (x, gBt) to separate the Bouguer gravity anomaly into rock formation depths. Using an algorithm of the previous solution of the sphere-shaped body model equation, the calculated values of the model are used as initial values for the iteration process in the semiinversion technique that applied for each of the profile points.

Testing of the forward synthetic model of three-layer case
A three-layer synthetic structure model was created, as seen in Figure 5. Assuming the number of deposited rock formations is N = 3 for the three-layer case. These  rock formations are depositing as stacked layered horizontally, at depths from top to bottom, as z 1 = 8 km, z 2 = 10 km and z 3 = 12 km (decreasing by 2 km in high areas and increasing by 3 km in low areas) and their densities, ρ 1 = 3.000 g/cm 3 , ρ 2 = 3.100 g/cm 3 and ρ 3 = 3.200 g/cm 3 , respectively, as they were deposited and overlying a basement rock of density ρ b = 2.670 g/cm 3 , and it tops hypothetically at depth z b = 20 km ( Figure 5). The vertical density contrasts calculated to be as follows: which are ∆ρ 1 ¼ 0.33 g/cm 3 , ∆ρ 2 ¼ 0.43 g/cm 3 and ∆ρ 3 ¼ 0.53 g/cm 3 , respectively, and the average density contrast is also calculated to be: which is Δρ = (0.33 +0.43 +0.53)/3 = 0.43 g/cm 3 . In the synthesis model, it was intended to increase and decrease some depths of the rock formations to create irregular surfaces in different places, also simulating the grabens and horst shapes and testing the stability of the present semiinversion method.

Creation of synthetic profile' line.
To simplify the concept of ZOGM, for the previously threelayer case ( Figure 5), three points are assumed representing three-layer tops (O 1 , O 2 and O 3 ) on a vertically coaxial (AC). Such that (C) is a point on the top of the basement rock ( Figure 6a). And, the three points are at the same time representing the centres of three sphereshaped bodies. Thus, their radii (R 1 , R 2 and R 3 ) correspond and are equal to that of rock formations' depths (Z 1 , Z 2 and Z 3 ). According to the ZOGM concept, the calculated gravity effect value of the observed point (A) at the Earth's surface is obtained for each isolated sphere by using equation (7) as follows: = ((4/3) * π * (0.43) * 0.00667 * (10)) = 0.14808 m.Gal.

Forward model calculation.
The forward model of the previous three-layer case was calculated at hypothetical 51 points on the Earth's surface (x(i), g Bt i ð Þ), and the results are shown in Table 1.

Bouguer gravity anomaly calculation, for two models (1) & (2)
In model (1) of the heterogeneity density distributions, where the rock formation densities are distributed, according to Walther's Law of deposition, and by applying equations (15) and (16), for calculating the gravity effect and their corresponding depths, the obtained results are shown in Table 2. Then, jointly inverting these data of model by incorporating them through the solution algorithm for each point of the total gravity anomalies of synthetic profile line, which was created from 51 ZOGM points, and the results are summarised in Table 3.
Whereas in model (2) of the homogeneity density distributions, where the rock formation densities  (15) and (16), for calculating the gravity effect and their corresponding depths, the obtained results are shown in Table 4. Then, jointly inverting these data of model by incorporating them through the solution algorithm for each point of the total gravity anomalies of synthetic profile line, which was created from 51 ZOGM points, and the results are summarised in Table 5.

Inversion for the forward synthetic model of three-layers case
To invert the gravity anomaly of the synthetic model of the three-layer case (Figure 6b), equation (15) is rearranged to calculate depths as follows: where (g Bt i ð ÞÞ is total gravity anomalies along the vertical co-axial of centres of sphere-shaped bodies at observed points of indices (i), and (Δρ i ð Þ) represents the density contrasts at (i = 1, 2, and 3). Hence, the (g Bt i ð Þ) can inverse into corresponding rock formation depths (z (i)), by solving equation (16) iteratively, as well as tracing the rock formation depths of profile line points ( Figure 6c). Also, the depth of the basement is obtaining as follows: This means the basement's rock depth at any point is twice the average vertical depths of overlaying (N) layers at the point. Therefore, the inversion of the previous synthetic model was performed at 51 points on the Earth's surface (j = 51 points), where the depths by using the two deposition models and basement depths are given, respectively, as follows: where g Bt j ð Þ is the gravity effects calculated for a number of points (j), g Bt i ð Þ is the model (1 or 2) calculated gravity effect, the average vertical density (AVD) and the effective density (D effective ).
The average vertical density is given as follows: where i = 1, 2, and N = 3 for three-layer case. While the effective density is obtained by summation of equation (9) and multiplied with thicknesses, the result is divided by the summation of thicknesses as follows: The simultaneous tracking of rock formation depths (z j ð Þ) profile points (51 points) is done by using equation (18), which depends on the vertical density contrasts of formations with basement rocks and the correlated profile's gravity (g Bt j ð Þ), with the modelling gravity effects calculated for hypothetical deposition basin (g Bt i ð Þ). The results of calculating gravity effects and depths for the hypothetical three-layer cases are summed in Table 6.

Testing noisy data of the synthetic three-layers model
For testing the influence of the noise on the present method, 10% random noise has been added to the synthetic gravity anomaly data (Figure 7a), using the following equation: where g Bt err i ð Þ is the output contaminated gravity anomaly of the input gravity anomaly data g Bt i ð Þ, corresponding to each observed i ð Þ, and then added 10% of Random Noises RAND i ð Þ � 0:1 ð Þ, that are used in equation (22).
The corresponding g Bt err i ð Þ values are gotten from equation (22) and used in depth inversion by applying equations (18) and (19), where the results of rock formation depths and also basement depths are obtained in Table 7, in comparison with both freenoisy data and noisy data, respectively.
The inversion and tracking depth of original synthetic gravity anomaly data (Figure 7b) are comparable with the inversion and tracking depth for data being contaminated with 10% random error or noise (Figure 7c), which means the new method is stable because it depends only on the vertical density contrasts of formation between rock formations and basement rock.
In equation (18), the calculation of the depth is based on automatic recovering of density contrasts of rock formation top boundaries with basement rock from used hypothetical deposition basin model (1 or 2).
It is noted from Tables 2-5 that the resulted depths are the same for the two models (1 & 2), which means the depths calculated by the present method do not depend on the density distributions and even the difference in gravity effects that was calculated for each of two models, but means that the calculated depths only depends on the vertical density-density calculation (or average density contrasts, effective densities and average vertical density contrasts).

Real-field application cases
The 2D semi-inversion method has been assessed for two field cases of different lithological and geological structural aspects, which are the Abu Roash Dome Area, southwest Cairo Egypt, and the Mors Salt Dome, North Jutland, Denmark. The semi-inversion procedures for the gravity data investigating real areas are, generally, summarised in the following steps (1) digitised and re-contouring the available Bouguer gravity anomaly map with a proper equal contour interval (m. Gal).
(2) building two hypothetical deposition models (heterogeneous and homogeneous) of rock formation depths and their densities from available borehole data (controlling point).
(3) digitising a gravity anomaly profile from the location of the borehole data to any of being investigated geological structures and/or lithologic features in the area. (4) interpretation by inverting the profile' gravity anomaly values at each point automatically into their equivalent rocks' formation depths values of at the same points. (5) by implementation of the proposed MATLAB algorithm that was previously prepared, and (6) comparing the results with other researchers' and authors' methods.

Abu Roash Dome Area, southwest Cairo, Egypt
The famous Abu Roash Area is located 10 km to the southwest of Cairo (Abdel Khalek et al., 1989). The Abu Roash Dome Area constitutes a complex Cretaceous sedimentary succession with outstanding tectonic features, as shown in Figure 8.  The depth of the basement rocks of the Abu Roash Dome Area has been estimated by gravitational methods of many researchers such as Abdelrahman et al. (1985), Abdelrahman et al. (1993), Essa (2007), and other authors.

Semi-inversion of the Abu Roash Bouguer gravity anomalies
The steps of semi-inversion approach as previously described were carried out upon two digitised profiles AA' and BB' for the available Bouguer gravity anomaly  map that covers the Abu Roash Dome Area (Figure 9). The original Bouguer gravity anomaly map was compiled by the Egyptian General Petroleum Corp. (G.P.C., 1984), using a Bouguer density of 2.35 g/cm 3 . In addition, the Abu Roash-1, borehole data (El-Malky 1985) was used for choosing the interesting rock formations' depths and their corresponding densities (Table 8). The borehole data was used to build the two hypothetical models (1 & 2), as shown in Figure 10.
The gravity effect values were theoretically calculated for each point of the selected five formations' model (1) of the Abu Roash Dome Area (Table 9). Hence, their corresponding depths were calculated by semi-inversion of the gravity effects using the ZOGM concept (Table 10). It is obvious that the inverted depths are the same as were estimated before in the drilled borehole (Abu Roash-1). By carrying out the same calculations  for model (2), the results are obtained (Tables 11, and 12). It is obvious that the inverted depths for both of the two models are the same as were estimated before in the drilled borehole (Abu Roash-1). And, Figure 11 and Figure 12 represent the graphics of the inverted depth results of the profiles AA' and BB' each with the same interval of 2.09 km.

Interpretation of the Abu Roash' inversion results and comparing results
As seen from Table 10, the basement rock depth values are 1.94 km, at point (xc = 37.49482 km, gB = −6.08805 m. Gal), and 1.92 km, at point (xc = 241.6333 km, gB = −6.0403 m. Gal), at the profile AA'. As seen also from Table 12, the basement rocks depth values are 2.09 km, at point (xc = 175.6856 km, Figure 11. Bouguer gravity anomaly of profile AA' and its semi-inversion to depths of Abu Roash Dome Area. gB = −6.55445 m. Gal), and 1.95 km, at point (xc = 227.1058, gB = −6.12654 m. Gal), at the profile BB'. Moreover, the average deepest basement depth is 5.07 km, which was obtained by the present method.
The inverted basement rock depth values of the two profiles AA' and BB' of the new method (1.94, 1.92, 2.09 and 1.95) are comparable (Table 13), with depths estimated in the Abu Roash-1 drilled borehole (1.902 km) or calculated by Abdelrahman et al. (1995a) (1.620 km) and Essa (2007) (1.91 km). The small differences in the depth estimated values of basement rock formations in the present method from the drilled value are due to the direction of the profile (dip or strike), and this from of the geological point of view is expected.

Mors Salt Dome, North Jutland, Denmark
The Mors Salt Dome area is confined to the island of Mors, which is situated in the northwestern part of Jutland, Denmark, as shown circled with red ( Figure 13a). The Mors island covers an area of about 360 km 2 . It is 10-15 km wide and about 35 km long with the SSW-NNE trend (Jorgensen et al.,2005). The dome-like structure on central Mors is made of chalk which covers the top of the Erslev salt diaper. In the bedrock map of Mors, the structural contour lines at 25 m intervals show the elevation of the pre-Quaternary surface (Figure 13b).

Semi-inversion of the Mors Bouguer gravity anomalies
The steps of semi-inversion approach as previously described were carried out upon two digitised profiles AA' and BB', respectively, for the available Bouguer gravity anomaly map and the Mors Salt Dome (Figure 14). The Bouguer's gravity anomaly map (Saxov 1956;Madirazza 1980;Sharma 1997), covers the Mors Salt Dome Area. In addition, the Erslv-2 and the Mors-1, available boreholes data (Gomm 1982), were used for choosing the interesting formations' depths and  their corresponding densities (Tables 13 and 14), respectively. The borehole data (Table 15) were used to build two models (1 & 2) to distribute the densities of the rock formations, one heterogeneous and the other homogeneous, respectively ( Figure 15).
The gravity effect values were theoretically calculated for each point of the selected six formations' model (1) of the Mors Salt Dome Area (Table 16). Hence, their corresponding depths were calculated by semi-inversion of the gravity effects using the ZOGM concept (Table 17). It is obvious that the inverted depths are the same as were estimated before in the drilled borehole (Erslv-2 and/or the Mors-1). By carrying out the same calculations for model (2), the results are obtained (Tables 18,and 19). It is obvious that the inverted depths for both of the two models are the same as were estimated before in the drilled borehole (Erslv-2 and/or the Mors-1). And, Figure 16 and Figure 17 represent the graphical of the inverted depths result of the profiles AA' and BB' each with intervals of 2.01 km and 2.7 km, respectively.

Interpretation of the Mors inversion and comparing results
The interpretation of profile AA' (model 1) restricts inside the closed contour value of 17 m, Gal or less where the highest area of gravity anomaly, nearest to the top of salt dome and borehole Erslv-2. From Table 16 basement is 4.86-5.11 km. Whereas, the interpretation of profile BB' (model 2) also restricts inside the closed contour value of 17 m. Gal or less, as seen from Table 18, the range of inverted depth of the Mors Salt Dome is 4.03-4.24 km, and the range of inverted depth of the basement is 4.84-5.10 km. Table 20 illustrates the range of inverted depth value of the salt dome in the present method. To some extent it is comparable with the value of 4.8 km calculated by the sphere model (Sharma 1997), the value of 4.70 km calculated by Normal Full Gradient (Aghajani et al., 2009) and the value of 4.82 km calculated by General regression neural network (Hajian et al. 2012).

Discussion
The new method is being considered to be the first real stage of an attempt to overcome the problem of ambiguity in gravity interpretation, and it was based on the accumulated work of the previous researchers, using the sphere-shaped body (the point masses) models to represent rock formation tops. The other methods such as Abdelrahman et al. (1995a) were used to the least-squares method based on the analytical expression of simple numerical horizontal gravity gradient anomalies to estimate depth and shape of a buried body. Essa (2007) utilised the numerical fourth horizontal derivatives computed from the observed gravity anomaly, using filters of successive window lengths to estimate the depth and shape of a buried structure, and Sharma (1997) used the half-width method for depth estimation of buried spherical body. Aghajani et al. (2009) used the normalised full gradient of gravity anomaly, in estimating the depth of the anomalous spherical bodies, and Hajian et al. (2012) used the feed-forward back-propagation neural networks to simultaneously estimate shape factor and depth of gravity anomalies. It is worth mentioning that the depths of basement rock were estimated by the other researchers' methods, using the sphere-shaped body model. But, they had to separate the Bouguer's gravity anomalies into residual and regional anomaly components which have to be carried out before only the basement rock depth estimation techniques were used. As well as, all techniques were concentrated only to estimate the depth of boundary      between the basement and overlying sediments. Whereas, in the present method, the sphere-shaped body (or the point masses) model was used to separate Bouguer gravity anomalies directly into the depths of the rock formations and the underlying basement rocks depths, where the models were constrained through the introduced the new idea for the concept of ZOGM.

Conclusions
The present research embraces a new way of thinking in treating the ambiguity in interpreting the Bouguer gravity anomaly data. The method applications on real areas showed comparable values with drilled boreholes, besides it is to some extent the stability of the method, which will encourage optimising it in the future. As well, the present technique is similar to some extent tracing formations from borehole data to the seismic cross section in the seismic interpretation process. Therefore, the method can be recommended as a reconnaissance tool in a preliminary petroleum bid-round evaluation procedure and before the seismic surveys, where it will reduce the cost of price for geophysical exploration.  and valuable constructive criticism discussion. The author would also like to thank Prof. Dr. Ahmed Said Salim for his always encouraging and also his family, especially his wife, for providing a good environment for carrying out this research.

Disclosure statement
No potential conflict of interest was reported by the author(s).