# Review of Environmental and Geological Microgravity Applications and Feasibility of Its Employment at Archaeological Sites in Israel

Authors: Lev V. Eppelbaum
Journal: International Journal of Geophysics
 (2011)

---

## Abstract

Microgravity investigations are widely applied at present for solving various environmental and geological problems. Unfortunately, microgravity survey is comparatively rarely used for searching for hidden ancient targets. It is caused mainly by small geometric size of the desired archaeological objects and various types of noise complicating the observed useful signal. At the same time, development of modern generation of field gravimetric equipment allows to register promptly and digitally microGal (10−8 m/s2) anomalies that offer a new challenge in this direction.  An advanced methodology of gravity anomalies analysis and modern 3D modeling, intended for ancient targets delineation, is briefly presented. It is supposed to apply in archaeological microgravity the developed original methods for the surrounding terrain relief computing. Calculating second and third derivatives of gravity potential are useful for revealing some closed peculiarities of the different Physical-Archaeological Models (PAMs). It is underlined that physical measurement of vertical gravity derivatives in archaeological studying has a significant importance and cannot be replaced by any transformation methods. Archaeological targets in Israel have been ranged by their density/geometrical characteristics in several groups. The performed model computations indicate that microgravity investigations might be successfully applied at least in 20–25% of archaeological sites in Israel.


---

## Body

## 1. Introduction

Development of a new modern gravimetric and variometric (gradientometric) equipment (permitting registering earlier inaccessible small anomalies and improving the observation methodology) and creation of new methodologies for gravity data processing and interpretation have triggered an arising of microgravity methodology application in environmental and economic minerals geophysics.Microgravity is recognized now as a effective tool for analysis of various geological inhomogeneities in subsurface, monitoring of volcanic activity, and prospecting of useful minerals (e.g., [1–30]).  Butler [6, 7] carried out a careful investigation of the gravity and gravity-gradient determination concepts and corresponding interpretative procedures of microgravity. The types of noise (disturbances) arising in the microgravity investigations are studied in detail in Debeglia and Dupont [31]. Styles et al. [17] discussed a few actual problems suggested to removing of noise components arising in microgravity under complicated environments.


## 2. A Brief Review of Employment of Archaeological Sites: A World Experience

Obviously, the history of gravity (microgravity) application at archaeological sites is beginning from the work of Linnington [32]. The Linnington’s attempt to examine the Etruscan chambered tombs at Cerveteri (Italy) was not very successful, but next microgravity applications in archaeology have shown the prospectiveness of this method.Fajklewicz [3] investigation was suggested to examination of the gravity vertical gradient (Wzz) over underground galleries. Apparently, he was first indicating the significant difference between the physically measured Wzz and this value obtained by the way of any transformation. Interesting examples of microgravity anomalies examination from some archaeological targets are presented in Blížkovský [4].Kerisel [33] has been used microgravity method for estimation of the Cheops pyramid (Egypt) weight and peculiarities of it structure. Lakshmanan [34] microgravity investigations inside, above, and around the Cheops pyramid led to an evaluation of the structure’s overall density and of density changes in the structure.Slepak [35] obtained gravity anomalies of 30–80 microGals reflecting the remains of ancient buildings in the territory of Kazan Kremlin (Russia). A careful examination of microgravity anomaly at the Boden Vean (Cornwall, UK) allowed to suggesting existence of new archaeological remains [36].Pasteka and Zahorec [37] examined a complex microgravity anomaly in the area of a former church of St. Catherine (Dechtice, Slovakia) by the use of stripping technology. Microgravity anomalies of 15–20 microGals were registered over the remains of the Late Byzantine church walls in Umm er-Rasad (Jordan) [38].A small number of microgravity measurements were performed in the Bedem archaeological locality in Yugoslavia [39]. Abad et al. [23] have carried out an assessment of a buried rainwater cistern in the Carthusian Monastery (Valensia, Spain) on the basis of 2D microgravity modeling. Microgravity survey has been successfully employed at the Roman Amphitheatre of Durres, Albania [40]. Microgravity examination performed by Sarris et al. [41] allowed identifying negative gravity anomalies above the ancient tombs in the Roman Cemetery at Kenchreai (Korinthia), Greece.Castiello et al. [42] have been delineated underground ancient cavity in the complex urban environment of Naples. Detailed gravity examination has been effectively performed in the Red Square (Moscow) several years ago with the aim to delineate ancient underground galleries (personal communications).


## 3. Analysis of Microgravity Anomalies

It is known that the trivial formulas of quantitative analysis (based on simple relationships between the gravity field semiamplitude and center of the disturbing body) are widely presented in the geophysical literature (e.g., [43, 44]). However, absence of reliable information about the normal gravity field in the studied areas strongly limits practical application of these methods.


### 3.1. Some Common Aspects between the Magnetic and Gravity Fields

Gravity field intensityF is expressed as(1)F=-gradW,
where W is the gravity potential.For anomalous magnetic fieldUa we can write (when magnetic susceptibility ≤0.1 SI unit) [45](2)Ua=-gradV,
where V represents the magnetic potential.Let us consider analytical expressions of some typical models employed in magnetic and gravity fields (Table1).Table 1
Comparison of some analytical expressions for magnetic and gravity fields.
FieldAnalytical expressionMagneticThin bed (TB)Point source (rod)Zv=2I2bzx2+z2  (*)Zv=mz(x2+z2)3/2  (†)GravityHorizontal Circular Cylinder (HCC)SphereΔg=2Gσzx2+z2  (**)Δg=GMz(x2+z2)3/2  (††)HereZv is the vertical magnetic field component at vertical magnetization, I is the magnetization, b is the horizontal semithickness of thin bed (TB), m is the elementary magnetic mass, z is the depth to a center of body (for HCC and sphere) and depth to the upper edge of TB and rod (point source), and M is the mass of sphere.It is clear that expressions (1) and (2) are analogical ones and (*) and (**), (†), and (††), respectively, are proportional ones.Taking into account all above mentioned, we can apply for the gravity field analysis the advanced interpreting methodologies developed in magnetic prospecting for complicated environments [45]. For instance, we can interpret anomaly from the gravity HCC by the use of formulas applied in the magnetic TB (see Table 1).We can also calculate the “gravity moment,” which could be used for classification and ranging gravity anomalies from various types of targets. The “gravity moment” of HCC may be calculated by the use of corresponding formula for the magnetic TB [45]:(3)MΔg=12Δgah,
where MΔg is the gravity moment, Δga is the amplitude of gravity anomaly (in mGal), and h is depth of HCC occurrence (in meters).


### 3.2. Calculation of Inclined Terrain Relief Influence

A significant number of archaeological sites occur under conditions of rugged terrain relief. Uneven observation lines are responsible for variations in the distance from the point of measurements to the source that can strongly complicate quantitative analysis of gravity anomalies [46]. If anomalies are observed on an inclined profile, then the obtained parameters characterize a certain fictitious body. The transition from fictitious body parameters to those of the real body is performed using the following expressions (the subscript “r” stands for a parameter of the real body) [45]:               (4){hr=h+x0⋅tanω0xr=-h⋅tanω0+x0},
where h is the depth of the upper edge occurrence, x0 is the location of the source’s projection to plan relative to the extremum having the greatest magnitude, and ω0 is the angle of the terrain relief inclination (ω0>0 when the inclination is toward the positive direction of the x-axis).A simple example of interpreting gravity anomaly from a buried cavity on inclined profile is presented in Figure1. It should be noted that if the well will be drilled on the projection of gravity anomaly minimum to the earth’s surface could not met this target (due to disturbing effect of inclined relief). Application of the improved tangent’s method (developed in the magnetic prospecting for complex environments [45]) together with (4) permits to obtain the position of the cavity center (bold red point in Figure 1) with the sufficient accuracy. MΔg in this example consists of 1/2·10.5m·0.0175mGal≅0.092mGal·m.Figure 1
Quantitative analysis of gravity anomaly on inclined profile from a buried sphere.


### 3.3. Calculation of Second and Third Derivatives of Gravitational Potential

Second and third derivatives of gravity potential could be very useful for defining some important peculiarities of archaeological targets location in different physical-archaeological models (PAMs). It is well known that gravity field is a function of mass and derivatives of gravity field is a function of density. Taking into account small depth of archaeological targets and their not large geometrical size, observation of both vertical and horizontal derivatives of gravity field undoubtedly will permit to obtain new important information about the desired targets. Obviously, integrated analysis of gravity field and vertical and horizontal derivatives will significantly extend the possibilities of geophysical investigations at archaeological sites. It is necessary to underline that physical measurement of vertical gravity derivatives cannot be replaced by computing of this parameter obtained by any transformation procedures: the Wzz values computed from the field ΔgB, as rule, show decreasing values compared with the Wzz obtained from physical measurements (or computation at two levels).


## 3.1. Some Common Aspects between the Magnetic and Gravity Fields

Gravity field intensityF is expressed as(1)F=-gradW,
where W is the gravity potential.For anomalous magnetic fieldUa we can write (when magnetic susceptibility ≤0.1 SI unit) [45](2)Ua=-gradV,
where V represents the magnetic potential.Let us consider analytical expressions of some typical models employed in magnetic and gravity fields (Table1).Table 1
Comparison of some analytical expressions for magnetic and gravity fields.
FieldAnalytical expressionMagneticThin bed (TB)Point source (rod)Zv=2I2bzx2+z2  (*)Zv=mz(x2+z2)3/2  (†)GravityHorizontal Circular Cylinder (HCC)SphereΔg=2Gσzx2+z2  (**)Δg=GMz(x2+z2)3/2  (††)HereZv is the vertical magnetic field component at vertical magnetization, I is the magnetization, b is the horizontal semithickness of thin bed (TB), m is the elementary magnetic mass, z is the depth to a center of body (for HCC and sphere) and depth to the upper edge of TB and rod (point source), and M is the mass of sphere.It is clear that expressions (1) and (2) are analogical ones and (*) and (**), (†), and (††), respectively, are proportional ones.Taking into account all above mentioned, we can apply for the gravity field analysis the advanced interpreting methodologies developed in magnetic prospecting for complicated environments [45]. For instance, we can interpret anomaly from the gravity HCC by the use of formulas applied in the magnetic TB (see Table 1).We can also calculate the “gravity moment,” which could be used for classification and ranging gravity anomalies from various types of targets. The “gravity moment” of HCC may be calculated by the use of corresponding formula for the magnetic TB [45]:(3)MΔg=12Δgah,
where MΔg is the gravity moment, Δga is the amplitude of gravity anomaly (in mGal), and h is depth of HCC occurrence (in meters).


## 3.2. Calculation of Inclined Terrain Relief Influence

A significant number of archaeological sites occur under conditions of rugged terrain relief. Uneven observation lines are responsible for variations in the distance from the point of measurements to the source that can strongly complicate quantitative analysis of gravity anomalies [46]. If anomalies are observed on an inclined profile, then the obtained parameters characterize a certain fictitious body. The transition from fictitious body parameters to those of the real body is performed using the following expressions (the subscript “r” stands for a parameter of the real body) [45]:               (4){hr=h+x0⋅tanω0xr=-h⋅tanω0+x0},
where h is the depth of the upper edge occurrence, x0 is the location of the source’s projection to plan relative to the extremum having the greatest magnitude, and ω0 is the angle of the terrain relief inclination (ω0>0 when the inclination is toward the positive direction of the x-axis).A simple example of interpreting gravity anomaly from a buried cavity on inclined profile is presented in Figure1. It should be noted that if the well will be drilled on the projection of gravity anomaly minimum to the earth’s surface could not met this target (due to disturbing effect of inclined relief). Application of the improved tangent’s method (developed in the magnetic prospecting for complex environments [45]) together with (4) permits to obtain the position of the cavity center (bold red point in Figure 1) with the sufficient accuracy. MΔg in this example consists of 1/2·10.5m·0.0175mGal≅0.092mGal·m.Figure 1
Quantitative analysis of gravity anomaly on inclined profile from a buried sphere.


## 3.3. Calculation of Second and Third Derivatives of Gravitational Potential

Second and third derivatives of gravity potential could be very useful for defining some important peculiarities of archaeological targets location in different physical-archaeological models (PAMs). It is well known that gravity field is a function of mass and derivatives of gravity field is a function of density. Taking into account small depth of archaeological targets and their not large geometrical size, observation of both vertical and horizontal derivatives of gravity field undoubtedly will permit to obtain new important information about the desired targets. Obviously, integrated analysis of gravity field and vertical and horizontal derivatives will significantly extend the possibilities of geophysical investigations at archaeological sites. It is necessary to underline that physical measurement of vertical gravity derivatives cannot be replaced by computing of this parameter obtained by any transformation procedures: the Wzz values computed from the field ΔgB, as rule, show decreasing values compared with the Wzz obtained from physical measurements (or computation at two levels).


## 4. The Developed Algorithm for 3D Modeling of Gravity Field

TheGSFC (Geological Space Field Calculation) program was developed for solving a direct 3D gravity and magnetic prospecting problem under complex geological conditions [16, 45]. This program has been designed for computing the field of Δg (Bouguer, free-air, or observed value anomalies), ΔZ, ΔX, ΔY, ΔT, as well as second derivatives of the gravitational potential under conditions of rugged relief and inclined magnetization. The geological space can be approximated by (1) three-dimensional, (2) semi-infinite bodies and (3) those infinite along the strike closed, left-hand non-closed, right-hand non-closed and open (Figure 2). Geological bodies, are approximated by horizontal polygonal prisms (Figure 3).Figure 2
Types of geological bodies used in modeling.Figure 3
Computing derivatives of gravity potential for a horizontal polygonal prism.The program has the following main advantages (besides abovementioned ones): (1) simultaneous computing of gravity and magnetic fields; (2) description of the terrain relief by irregularly placed characteristic points;  (3) computation of the effect of the earth-air boundary by the method of selection directly in the process of interpretation; (4) modeling of the selected profiles flowing over rugged relief or at various arbitrary levels (using characteristic points); (5) simultaneous modeling of several profiles; (6) description of a large number of geological bodies and fragments. The basic algorithm realized in theGSFC program is the solution of the direct 3D problem of gravity and magnetic prospecting for horizontal polygonal prism limited in the strike direction (Figure 3). In the developed algorithm integration over a volume is realized on the surface limiting the anomalous body.Analytical expression for the first vertical derivative of gravity potential of (m-1) angle horizontal prism (Figure 3) has been obtained by integrating a common analytical expression:(5)Wz′=-∫sz(R+y)Rdxdzz(R+y)R|y1y2,
where R=x2+y2+z2, Sis the area of normal section of the prism by the plane of xOz(6)Wz′T={-y1ln⁡R11j+U1jR21j+U2j)-Gσ∑j=1m-1[Vjsinαj(ln⁡R12j+y2R22j+y2-ln⁡R11j+y1R21j+y1)+Vjcos⁡αj(Vj2R12jR22j+U1jU2jy22r1jr2j(y22+Vj2)sgn⁡(y2Vj)arccosVj2R12jR22j+U1jU2jy22r1jr2j(y22+Vj2))-(Vj2R11jR21j+U1jU2jy12r1jr2j(y22+Vj2)sgn⁡(y1Vj)arccosVj2R11jR21j+U1jU2jy12r1jr2j(y12+Vj2))+cos⁡αj(y2ln⁡R12j+U1jR22j+U2j-y1ln⁡R11j+U1jR21j+U2j)]∑j=1m-1[Vjsinαj(ln⁡R12j+y2R22j+y2-ln⁡R11j+y1R21j+y1)},
where G is the gravitational constant, σ is the density of the  body, and αj is the angle of the prism’s side inclination.Detailed description of analytical expressions for the first and second derivatives of gravity potential of the approximation model of the horizontal polygonal prism and their connection with magnetic field is presented in Khesin et al. [45].


## 5. Terrain Correction Applied for Nonconventional Approaches

It is well known that accuracy of microgravity investigations substantially depend on the accuracy of terrain correction (TC) computing. Two approaches presented below were applied for exact TC calculation for the detailed Bouguer gravity observed at ore deposits occurring in the Lesser and Greater Caucasus.


### 5.1. First Approach

A first method was applied in the Kyzylbulakh gold-pyrite deposit situated in the Mekhmana ore region of the Nagorny Karabakh (Lesser Caucasus) under condition of rugged relief. This deposit is well investigated by mining and drilling operations and therefore was used as a reference field polygon for testing this approach. A special scheme for obtaining the Bouguer anomalies has been employed to suppress the terrain relief effects dampening the anomaly effects from the objects of prospecting. The scheme is based on calculating the difference between the free-air anomaly and the gravity field determined from a 3D model of a uniform medium with a real topography. 3D terrain relief model with an interval of its description of 80 km (the investigated 6 profiles of 800 m length are in the center of this interval) was employed to compute (by the use ofGSFC software) the gravitational effect of the medium (σ = 2670 kg/m3). With applying such a scheme the Bouguer anomalies were obtained with accuracy in two times higher than that of TC received by the conventional methods. As a result, on the basis of the improved Bouguer gravity with the precise TC data, the geological structure of the deposit was defined [9].


### 5.2. Second Approach

Second approach was employed at the complex Katekh pyrite-polymetallic deposit, which is located at the southern slope of the Greater Caucasus (northern Azerbaijan). The main peculiarities of this area are very rugged topography of SW-NE trend, complex geology and severe tectonics. Despite the availability of conventionalΔgB (TC far zones were computed up to 200 km), for the enhanced calculation of surrounding terrain topography, a digital terrain relief model was created [16]. The SW-NE regional topography trend in the area of the Katekh deposit occurrence was computed as a rectangular digital terrain relief model (DTRM) of 20 km long and 600 m wide (our interpretation profile with a length of 800 m was located in the geometrical center of the DTRM). As a whole, about 1000 characteristic points were used to describe the DTRM (most frequently points were focused in the center of the DTRM and more rarely—on the margins). Thus, in the interactive 3D ΔgB modeling (by the use of GSFC software) was computed effect not only from geological bodies occurring in this area, but also from surrounding DTRM. In the issue of this scheme application, two new ore bodies were discovered [16].


## 5.1. First Approach

A first method was applied in the Kyzylbulakh gold-pyrite deposit situated in the Mekhmana ore region of the Nagorny Karabakh (Lesser Caucasus) under condition of rugged relief. This deposit is well investigated by mining and drilling operations and therefore was used as a reference field polygon for testing this approach. A special scheme for obtaining the Bouguer anomalies has been employed to suppress the terrain relief effects dampening the anomaly effects from the objects of prospecting. The scheme is based on calculating the difference between the free-air anomaly and the gravity field determined from a 3D model of a uniform medium with a real topography. 3D terrain relief model with an interval of its description of 80 km (the investigated 6 profiles of 800 m length are in the center of this interval) was employed to compute (by the use ofGSFC software) the gravitational effect of the medium (σ = 2670 kg/m3). With applying such a scheme the Bouguer anomalies were obtained with accuracy in two times higher than that of TC received by the conventional methods. As a result, on the basis of the improved Bouguer gravity with the precise TC data, the geological structure of the deposit was defined [9].


## 5.2. Second Approach

Second approach was employed at the complex Katekh pyrite-polymetallic deposit, which is located at the southern slope of the Greater Caucasus (northern Azerbaijan). The main peculiarities of this area are very rugged topography of SW-NE trend, complex geology and severe tectonics. Despite the availability of conventionalΔgB (TC far zones were computed up to 200 km), for the enhanced calculation of surrounding terrain topography, a digital terrain relief model was created [16]. The SW-NE regional topography trend in the area of the Katekh deposit occurrence was computed as a rectangular digital terrain relief model (DTRM) of 20 km long and 600 m wide (our interpretation profile with a length of 800 m was located in the geometrical center of the DTRM). As a whole, about 1000 characteristic points were used to describe the DTRM (most frequently points were focused in the center of the DTRM and more rarely—on the margins). Thus, in the interactive 3D ΔgB modeling (by the use of GSFC software) was computed effect not only from geological bodies occurring in this area, but also from surrounding DTRM. In the issue of this scheme application, two new ore bodies were discovered [16].


## 6. Feasibility of Microgravity Application on archaeological Sites in Israel

Analysis of the numerous archaeological and geological publications as well as the author’s investigation (e.g., [29, 47]) indicates that the ancient objects supposed for examination by the use of microgravity survey may be classified (in the order of decreasing) by the following way: (1)
underground ancient cavities and galleries,(2)
walls, remains of temples, churches, and various massive constructions,(3)
pavements and tombs,(4)
roman aqueducts (under favorable physical-geological environments),(5)
areas of ancient primitive metallurgical activity, including furnaces, (under favorable physical-geological environments).Examining the different archaeological targets in Israel, it was supposed that microgravity method might be effectively applied at least on 20–25% of ancient sites [29, 47].A simplified model example of buried pavement delineation is presented in Figure4. A buried pavement having the positive density contrast of 400 kg/m3and occurring at a depth of 1.8 m in uniform medium (PAM of one of the Megiddo sites was selected as a basis) could be easily recognized by a microgravity survey (Figure 4(a), anomalous effect from two bodies). Let us assume a low-density layer (2100 kg/m3) over the pavement. It makes the delineation of the pavement practically impossible in field conditions (registered anomaly is oscillating about 1 microGal) (Figure 4(a), anomalous effect from three bodies). At the same time values of the second vertical derivative of gravity potential Wzz (computed for the levels of 0.3 and 1.5 m) with a measurable accuracy testify to a presence of disturbing body (Figure 4(b)). The similar graphs were developed for the models of buried ancient caves of various radius, ancient iron furnaces, Roman aqueducts, and so forth.Comparison of Bouguer gravity and vertical gradient anomalies. (a) Bouguer gravity, (b) vertical gradientgz (Wzz) computed for the base of 1.2 m, and (c) archaeological sequence.
(a)(b)(c)Calculations ofWxxz and Wzzz might be also useful for delineation of anomalies from closely disposed objects and removing of regional background, respectively. Taking into account that many archaeological sites in Israel are located in the vicinity of natural (sea and lake basins, mountains, water channels, forests, etc.) and industrial (plants, buried pipelines, tunnels, railways, etc.) disturbing objects, this fact might have an essential significance. For instance, Figure 5 shows that computing of Wxxz (gxx) makes it possible to recognize reliably gravity effects from two closely located underground caves. Ancient Prehistoric caves corresponding to this PAM are investigated in several archaeological sites disposed closely to the Golan Heights.Computing of horizontal derivatives from models of two closely disposed caves. (a) computed gravity curve, (b) calculated first horizontal derivative of gravity fieldΔgx, (c) calculated second horizontal derivative of gravity field Δgxx, and (d) physical-geological model.
(a)(b)(c)(d)The planning microgravity investigations must be accompanied by development of preliminary physical archaeological models (PAM) [47, 48]. The PAM presented in Figure 6 reflects a real archaeological site located in the vicinity of Beit-Shemesh town (central Israel). The developed PAMwas used for the estimation of expected gravity anomalies amplitudes, calculation of the most optimal step of observations along profiles, and distances between profiles.Figure 6
Physical-geological model of buried prehistoric cave and computed 3D gravity anomalies. (a) location of projected profiles and disposition of buried cave (view over), (b) 3D computed gravity effects along profiles 1–5, (c) geological-archaeological sequence.


## 7. Conclusions

Archaeological remains in the Israeli territory are classified by the degree of microgravity method applicability. The described characteristics of the developed software for combined 3D gravity-magnetic modeling under complicated environments indicate that it is a powerful tool for microgravity examination at archaeological sites. Two earlier applied nonconventional schemes for computing of surrounding terrain relief in ore deposits may be successfully adopted for obtainingΔgB+TC in archaeological sites situated in the mountainous (submountainous) regions of Israel. A proposed scheme of gravity anomalies quantitative interpretation is briefly discussed and explained on a model example. It is shown the necessity of gravity field derivatives employment in some specific situations. The developed PAM of the typical prehistoric cave (central Israel) testifies to the usefulness of such models utilization.



---

*Source: 927080-2011-06-30.xml*