Modern Perspectives of the Microgravity Application in Archaeology: A Review

1. Introduction

Development of a new modern gravimetric and variometric (gradientometric) equipment (permitting registering earlier inaccessible minor anomalies and improving the observation methodology) and creation of new methods for gravity data processing and interpretation have triggered the application of microgravity methodology in environmental and economic minerals geophysics.

Microgravity is recognized now as an effective tool for analysis of various geological inhomogeneties in subsurface, including archaeological purposes (Linnington, 1966; Arzi, 1975; Fajklewicz, 1976; Blížkovský, 1979; Lakshmanan and Montlucon, 1987; Slepak, 1999; Di Filippo et al., 2005; Yuan et al., 2006; Abad et al., 2007; Eppelbaum, 2009; Eppelbaum et al., 2010; Panisova et al., 2013; Fais et al., 2015; Gołebiowski et al., 2018; Tianyao et al., 2000; Pašteka et al., 2020, 2022; Eppelbaum, 2022), karst localization (e.g., Colley, 1963; Beres et al., 2001; Rybakov et al., 2001; Debeglia et al., 2006; Eppelbaum, 2007; Eppelbaum et al., 2008; Whitelaw et al., 2008; Al-Zoubi et al., 2013; Ezersky et al., 2023), engineering problems (e.g., Yule et al., 1998; Butler, 2001; Branston and Styles, 2006; Bradley et al., 2007; Khosravi et al., 2019; Porzucek and Loj, 2021), mining and hydrocarbon geophysics (e.g., Sjostrom and Butler, 1996; Karshenbaum et al. 1997; Eppelbaum and Khesin, 2004; Styles et al., 2006; Nind and Seigel, 2007; Slonka, 2011; Gadirov and Eppelbaum, 2012), and volcanology (e.g., Watanabe et al., 1998; Deroussi et al., 2009; Yagi et al., 2024).

Bichara et al. (1981) first identified the specific problems encountered in microgravity investigations. Butler (1984a, 1984b) conducted a thorough investigation of the concepts and corresponding interpretative procedures for gravity and gravity-gradient determination in microgravity. The types of noise (disturbances) arising in the microgravity investigations are studied in detail in Debeglia and Dupont (2002). Styles et al. (2006) discussed several actual problems suggested for removing noise components arising in microgravity under complex environments. The advanced procedures developed, for example, in mining or engineering in microgravity, can be effectively transferred to archaeological investigations and vice versa.

The planning of microgravity investigations must be accompanied by the development of preliminary physical-archaeological models (PAM) (Eppelbaum, 2010). Such PAMs enable the estimation of the expected gravity effects from archaeological targets, allowing for the determination of the area under study, the observation step, GPS observation accuracy requirements, and the necessary accuracy of observations. Then, at subsequent stages, more complex PAMs can be developed that reflect the generalized properties of the archaeological objects being studied.

2. Different Kinds of Noise Appearing in Microgravity Investigations at Archaeological Sites

2.1. Typical kinds of Noise in Microgravity

Microgravity surveys in archaeological sites are complicated by several kinds of noise (Eppelbaum et al. 2010, 2011) (Figure 1). These disturbances are briefly considered below.

Artificial (Man-Made) Disturbances

The industrial component of noise mainly arises from surrounding buildings, engineering constructions, and various underground and transportation systems. The instrumental component is associated with the ‘shift zero’ and other technical properties of gravimeters (variometers) and their spatial locations. Human errors, notably, can accompany geophysical observations at any time. Finally, poorly documented (or undocumented) results from previous surveys can distort the preliminary development of the PAM.

Natural Disturbances

The Natural noise (disturbances) can be subdivided into Nonstationary and Stationary. Nonstationary noise includes the tidal variations, mainly from the Moon and the Sun’s attraction, and some variations associated with the geodynamic events at depth. Meteorological and atmospheric factors reflect changing pressure, temperature, humidity, and electrical fields. Such meteorological conditions as rain, lightning, snow, and hurricanes hamper the microgravity observations, too. Soil-vegetation factors associated with specific soil types (e.g., swampy soil or loose ground in deserts) and dense vegetation, which can sometimes hinder movement along the profile, also need to be considered. Finally, the movement of water in basins and rivers can also affect the readings.

Stationary disturbances consist of two main components. First are Geological-geophysical and archaeological factors. The application of any geophysical method depends primarily upon the existence of physical property differences between the objects under study and the surrounding medium. The physical limitations of the methods reflect measurable density contrast properties between the archaeological targets and the geo-environmental sequence. The thermal-density dependence of sedimentary rocks, as discovered by Gadirov and Eppelbaum (2015), may necessitate reinterpreting the results of some detailed gravity surveys in archaeological sites (where sedimentary thickness exceeds 4-5 meters). The spatial coordinate determination plays a significant role in microgravity, as even a slight inaccuracy in the coordinates yields a value compatible with the expected anomalies from the archaeological targets.

The variety of anomalous sources is composed of two factors: the variable surrounding medium and the variety of archaeological targets; both aspects are crucial and greatly complicate the interpretation of microgravity data. Additional disturbances can introduce the variable subsurface and various gravity trends.

The separate group of disturbances creates obvious Soil-vegetation factors, as land microgravity measurements cannot be carried out in unfavorable environmental conditions.

Uneven terrain relief has three branches: (1) Rugged surrounding terrain relief (which can be removed or eliminated by Terrain Correction), (2) Hardship of observations (the author of this paper has performed the gravity observations at the Greater Caucasus on the inclined relief with an angle of almost 40^o), and (3) Topographic mass effect which is generally twofold and involves the form and physical properties of the topographic features of the terrain relief as well as the effect of variations in the distance from the point of measurement to the hidden target (Eppelbaum, 2022).

The possibility of reducing topographic effects by grouping points of additional gravimetric observations around the central point located on the survey network was demonstrated in Khesin et al. (1996). A group of 4 to 8 extra points is situated above and below along the relief, approximately symmetrically and equidistant from the central point. The topographic effect is reduced to the obtained difference between the gravity field in the center of the group and its mean value for the whole group.

Let us consider the topographic effect, which is very important in microgravity studies, in detail. The correct removal (elimination) of regional gravity trends is not a trivial task (e.g., Telford et al., 1993; Jacoby and Peter, 2009). Figure 2 shows two cases of non-horizontal gravity observations with the presence of an anomalous body occurring close to the Earth’s surface.

The distorting effect of a non-horizontal observation line occurs when the anomalous target differs from the surrounding medium by density, which produces an anomalous vertical gradient. Analysis of the Δg_B anomalies observed on the horizontal and inclined rugged relief indicates that the gravity effects in these situations are cardinal different (Figure 2).

Although the Bouguer reductions were applied to the observations on the inclined relief, the computed Bouguer anomaly is characterized by the small negative values in the downward direction of the relief, whereas the anomaly on the horizontal profile has no negative values. Theoretically, it is not possible since the gravity effect from a body with a positive contrast density (750 kg/m³) occurring in the homogenous medium must be only positive and tend to zero at far distances from the target.

Thus, applying the conventional corrections does not eliminate this trend because the observation points for the anomalous object are significantly different, and the large anomalous body occurs close to the observation points (Eppelbaum, 2011). Hence, a special methodology is required for the gravimetric quantitative interpretation of anomalies observed in the rugged relief conditions (see next section).

2.2. Terrain Correction Applied for High-Accuracy Gravity Investigations: Two Nonconventional Approaches

How can we enhance the effectiveness and reliability of interpretation? Undoubtedly, it must be a multi-stage process. We must begin with nonconventional methodologies for reducing the surrounding topographic effect and terrain correction computation.

It is well known that the accuracy of microgravity investigations substantially depends on the accuracy of terrain correction (TC) computing. The two approaches presented below were applied to calculate the exact TC for the detailed Bouguer gravity observations at ore deposits in the Lesser and Greater Caucasus.

2.2.1. First Approach

The first method was applied in the Gyzylbulakh gold-pyrite deposit, situated in the Mekhmana ore region of Azerbaijan (Garabagh, Lesser Caucasus), under conditions of rugged relief. This deposit has been thoroughly investigated through mining and drilling operations; therefore, it 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, thereby dampening the anomaly effects from the objects being prospected. 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. A 3-D terrain relief model with an interval of 80 km (six investigated profiles of 800 m length are in the center of this interval) was employed to compute (using GSFC software) the gravitational effect of the medium (σ = 2670 kg/m³). By applying such a scheme, the Bouguer anomalies were obtained with an accuracy two times higher than that of TC received by the conventional methods. As a result, based on the improved Bouguer gravity with the precise TC data, the geological structure of the deposit was defined (Eppelbaum, 2019).

2.2.2. Second Approach

The 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 the very rugged topography of a SW-NE trend, complex geology, and severe tectonics. Despite the availability of conventional Δg_B (TC of far zones were computed up to 200 km), a digital terrain relief model was created to enhance the calculation of surrounding terrain topography (Eppelbaum and Khesin, 2004). The SW-NE regional topography trend in 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 in the geometrical center of the DTRM). About 1,000 characteristic points were used to describe the DTRM (most frequently, these points were focused on the center of the DTRM, and more rarely, on the margins). Thus, in the interactive 3D Δg_B modeling (using GSFC software), the effect was computed not only from geological bodies occurring in this area but also from the surrounding DTRM. In the issue of this scheme application, two new ore bodies were discovered (Eppelbaum and Khesin, 2004).

3. Feasibility of Microgravity Application at Different Types of Archaeological Sites in the World

Let us consider the application of the microgravity technique to various types of archaeological targets worldwide.

3.1. Egyptian Pyramids

Egyptian pyramids are significant (and not fully understood to this day) archaeological objects where microgravity studies can introduce crucial information. Kerisel (1988) employed the microgravity method to estimate the weight and peculiarities of the Cheops pyramid (Egypt). Lakshmanan and Montlucon (1987), Lakshmanan (1991), and Issawy and Radwan (2012) conducted microgravity investigations inside, above, and around the Cheops pyramid, leading to an evaluation of the structure’s overall density and density changes within the structure. Saleh et al. (2022) have detected two buried chambers and an unrecognized negative gravity trench at the Pyramid of Senusret II, Lahun, Fayoum, Egypt.

The potential possibilities of microgravity investigations have not been fully explored. An integrated analysis of Δg, Δg_x, and Δg_z carried out at different levels will enable a more precise interpretation of the complex gravity patterns from these enigmatic objects. A comprehensive estimation of the expected gravity effects was demonstrated in Pašteka et al. (2022).

Ivashov et al. (2023) proposed generalizations of previous investigations and a perspective on using a UAV to deliver a gravimeter to the desired point on the step surface of the Khufu Pyramid.

3.2. Ancient Underground Caves

A careful examination of microgravity anomaly at the Boden Vean (Cornwall, UK) allowed for the suggestion of new archaeological remains occurring (Linford, 1998). Castiello et al. (2010) have delineated an ancient underground cavity in the complex urban environment of Naples (Italy). The residual gravity map, with the vertical gravity map highlighted, showed the presence of empty cavities or cavities filled with low-density sedimentation in an archaeological site in Sardinia (Italy) (Fais et al., 2015).Chromčák et al. (2018) demonstrated the detection of ancient cavities in the Saint Mary Church (Kláštor, Slovakia), with the obtained results subsequently confirmed by GPR data analysis. Nicolas et al. (2024) study enabled the reliable detection of the ancient pits and voids at the Ouels Mine (Castel-Minier, France).

3.3. Ancient Pavements and Roads

Styles et al. (2006) studied the major collapse that occurred within the ancient Roman road connecting Dover to London. Here, the corrections for the effects of the main open-cut workings and for the influence of surrounding buildings and topography were fruitfully introduced.

A simplified model example of buried pavement delineation is presented in Figure 3. A buried pavement having the positive density contrast of 400 kg/m³ and occurring at a depth of 1.8 m in a uniform medium (PAM of one of the Megiddo sites (Finkelstein and Martin, 2022) was selected as a basis) could be easily recognized by a microgravity survey (Figure 3A, anomalous effect from two bodies). Let us assume a low-density layer (2100 kg/m³) over the pavement. It makes the delineation of the pavement practically impossible in the field conditions (registered anomaly is oscillating about one microGal) (Figure 3A, anomalous effect from three bodies). At the same time, the values of the second vertical derivative of gravity potential W_zz (computed for the levels of 0.3 and 1.5 m) with a measurable accuracy testify to the presence of a disturbing body by anomalies appearing over the pavement’s edges (Figure 3B).

3.4. Ancient Water Reservoirs

Abad et al. (2007) carried out an assessment of a buried rainwater cistern in the Carthusian Monastery (Valencia, Spain) based on 2D microgravity modeling. Ebrahimi et al. (2019) conducted a successful gravity investigation in conjunction with an electric resistivity analysis over the ancient aqueduct (qanat) in northeastern Iran; the qanat was identified by microgravity anomalies of –(10-20) µGal. Bárta et al. (2020) successfully applied detailed gravity surveys to detect old water reservoirs in the historical center of Prague.

A PAM of the buried Roman aqueduct in Israel has been developed on the example of the excavated Caesarea aqueduct (Porath and ‘Ad, 2015). The following parameters were assumed: depth of the HCC center – 2.5-2.8m (σ = 0), radius – 1.2 m, thickness of the aqueduct’s walls – 0.08 m (σ = 2600 kg/m³), the density of the surrounding medium σ = 2100 kg/m³, and length – 200 m. Despite the disturbing effects of the aqueduct walls, the calculated gravity effects were estimated as 10-15 µGal, i.e., quite detectable by the modern gravimeters.

3.5. Ancient Underground Galleries and Tunnels

Fajklewicz’s (1976) investigation was suggested for the examination of the vertical gravity gradient (W_zz) over the underground galleries. He was first indicating the significant difference between the physically measured W_zz and this value obtained via any transformation. Microgravity surveys have been successfully employed at the Roman Amphitheatre of Durres, Albania, both over galleries and walls (Di Filippo et al., 2005). A detailed gravity examination was effectively performed in the Red Square (Moscow) several years ago to delineate the ancient underground galleries (personal communications).

The microgravity survey successfully produced a map of the 150-year Williamson tunnels beneath inner-city Liverpool (England), nicely demonstrating the validity of applying the method to industrial archaeology even in ‘brownfield’ sites (Cuss and Styles, 1999). The precise gravity observations at the National Library court in Moscow enabled the revelation of an earlier unknown tunnel at a depth of 1.5 m (Chernov, 2009). Orfanos and Apostolopoulos (2011) concluded, based on their detection of an ancient tunnel in the Lavrion area (Greece), that a combination of microgravity and resistivity is a powerful interpretive tool. Madej et al. (2018) applied microgravimetry for the delineation of the buried tunnels in Walbrzych, Poland, for searching the ‘Gold Train’ of the Second World War.

3.6. Ancient Building Architecture

Slepak (1999) obtained gravity anomalies of 30-80 microGals reflecting the remains of ancient buildings in the territory of Kazan Kremlin (Russia). Pašteka and Zahorec (2000) examined a complex microgravity anomaly in a former church of St. Catherine (Dechtice, Slovakia) using stripping technology. Yuan et al. (2006) have conducted a comprehensive gravity field modeling over the Emperor Qin Shi Huang Mausoleum (China), where data from other geophysical fields were used to agree on the geometric and density properties. Microgravity anomalies of 15-20 microGals were registered over the remains of the Late Byzantine church walls in Umm er-Rasad (Jordan) (Bataynen et al., 2007). Rabbel et al. (2018) applied a set of geophysical methods to study the Roman church in the city of Iznik (southern Turkey); microgravity studies enabled the localization of a low-density zone in the near-surface, which was interpreted as a grave filled with loose sediment.

Palgiara and di Filippo’s (2009) microgravity survey has led to the detection of remains of walls, foundations, fillings, crypts, and passages in the Basilica of San Paolo fuori le Mura (Rome). The microgravity measurements in the Church of St. George, Svätý Jur (Slovakia) significantly supplemented GPR constructions (Panisova et al., 2016). Detailed gravity-magnetic surveys at Tepe Hissar-Damghan archaeological site, northeastern Iran (Sarlak and Aghajani, 2017) enabled the recognition of earlier unknown chambers and walls. Gołebiowski et al. (2018) have effectively utilized the microgravity technique in conjunction with other geophysical methods to delineate the medieval underground salt chambers in the village of Wiślica, Poland.

Blížkovský (1979) has identified, using the stripping technology, the crypt under the floor of the St. Venceslas Church in the town of Továčov (Czech Republic). Tianyao et al. (2000) outlined the underground palace in the Mao Ling mausoleum (Ming Dynasty) in Shaanxi (China) using a combined analysis of microgravity and vertical gravity gradient data. The microgravimetric prospection inside the Don Church (Valencia, Spain) enabled the localization of the crypt, which was found in the entrance of the church (Padín et al., 2012).

Detailed investigations in the historic centre of Kazan enabled the recognition of remains of ancient walls (Slepak and Nugmanova, 2004). Combined geophysical investigations with application of microgravity in the Late Byzantine ‘Lion Church’ in Umm er-Rasas (Jordan) enabled the identification of floors of a building with the remains of walls, rooms, paths, and foundations(Bataynen et al., 2007). Non-conventional integration of gravity and its horizontal gradient with the radon gas anomalies was performed in the Alkhantepe archaeological monument (Southern Azerbaijan) (Safarov et al., 2019). Calculation of the second derivative of the gravity potential enabled the contour of ancient walls in the Cot Sidi Abdallah site (Morocco) (Sari et al., 2020).

It is well known that in urban or industrial development conditions, methods for calculating gravitational effects from nearby buildings are of great importance. New techniques for solving this problem are described in Pánisová et al. (2012), Yu (2014), Zahorec et al. (2014), and Loj and Porzucek (2019).

3.7. Ancient Cemeteries and Tombs

The Linnington’s (1966) attempt to examine the Etruscan chambered tombs at Cerveteri (Italy) was not very successful but following microgravity applications in archaeology have shown the potential of this method. Detailed gravity survey in the Kenchreai cemetery (southern Greece) enabled the reliable localization of two ancient tombs (Sarris et al., 2006). Pašteka et al. (2020) suggested methodical recommendations for the detection of ancient crypts and tombs and presented practical examples of crypt localizations in various churches in Slovakia and the Czech Republic. Microgravity studies by Cesnek et al. (2019) enabled the recognition of a cavity placed 0.60m under the floor in the chapel of St. Mary’s Assumption, city of Púchov (Slovakia). Zvara et al. (2021) demonstrated microgravity employment with the application of the minimum support focusing stabilizers at St. Nicholas church in Trnava (Slovakia); using this approach demonstrated the robustness of the results obtained. AbouAly et al. (2023) successfully applied the microgravity technique in the animal cemetery at Saqqara, Egypt, to detect the ancient crypts.

4. Quantitative Analysis of Gravity Anomalies in Complex Physical-Environmental Conditions

4.1. Simple Conventional Formulas for Quantitative Analysis

It is known that the trivial formulas of quantitative analysis (based on the simple relationships between the gravity field semi-amplitude and the center of the disturbing body) are widely presented in the geophysical literature (e.g., Telford et al., 1993; Parasnis, 1997). However, the absence of reliable information about the normal gravity field (superposition of gravity anomalies of different orders) in the studied areas strongly limits the practical application of these methods.

4.2. Some Common Aspects Between the Magnetic and Gravity Fields

Gravity field intensity F is expressed as

$$F=-gradW,$$

(1)

where W is the gravitational potential.

For anomalous magnetic field U_a, we can write (when magnetic susceptibility ≤ 0.1 SI unit) (Khesin et al., 1996):

$$U_a=-gradV,$$

(2)

where V represents the magnetic potential.

Let us consider analytical expressions of some typical models employed in magnetic and gravity fields (Table 1).

Here Z_v is the vertical magnetic field component at vertical magnetization, I is the magnetization, b is the horizontal semi-thickness of 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.

Expressions (1) and (2) are analogical ones, and equations (3) and (5), (4) and (6), respectively, are proportional ones. Considering all the above-mentioned factors, we can apply the gravity field analysis and advanced interpreting methodologies developed in magnetic prospecting for complex environments (rugged terrain relief and superposition of different rank anomalies) (Khesin et al., 1996; Eppelbaum et al., 2001; Eppelbaum, 2019). For instance, we can interpret an anomaly from the gravity HCC using 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 using the corresponding formula for the magnetic TB (Eppelbaum, 2019):

$$M_{\Delta g}={\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\Delta g_{a}h,$$

(7)

where $M_{\Delta g}$is the gravity moment, $\Delta g_a$is the amplitude of the gravity anomaly (in mGal), and h is the depth of the HCC occurrence (in meters). Calculating the “gravity moment” will allow for the classification and ranging of the microgravity anomalies.

4.3. Calculation of Inclined Terrain Relief Influence

A significant number of archaeological sites are in rugged terrains. 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 (Alexeyev et al., 1996). If anomalies are observed on an inclined profile, then the obtained parameters characterize a particular 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) (Eppelbaum, 2019):

$$\left\{\begin{array}{l}{h}_r=h+x\cdot \tan {\omega }_0\\ x_r=-h\cdot \tan {\omega }_0+x\end{array}\right\},$$

(8)

where h is the depth of the upper edge occurrence, x_o is the location of the source’s projection to plan relative to the extremum having the most significant magnitude, and ω₀ is the angle of the terrain relief inclination (ω₀ > 0 when the inclination is toward the positive direction of the x-axis).

A simple example of interpreting a gravity anomaly from a buried cavity on an inclined profile is presented in Figure 4. It should be noted that if the borehole is drilled on the projection of the gravity anomaly minimum to the earth’s surface, this target could not be met (due to the disturbing effect of inclined relief). Application of the improved tangent method (developed in the magnetic prospecting for complex environments (Khesin et al., 1996; Eppelbaum et al., 2001)) together with Eq. (8) permits the position of the cavity center (bold red point in Figure 4) to be obtained with sufficient accuracy. In this example, $M_{\Delta g}$ consists of ${\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\cdot 10.5\mathrm{m}\cdot 0.0175\mathrm{mGal}\cong 0.092\mathrm{mGal}\cdot \mathrm{m}.$

4.4. Calculation of Second and Third Derivatives of Gravitational Potential

The second and third derivatives of the gravity potential could be beneficial for defining certain crucial peculiarities of archaeological target locations in various physical-archaeological models (PAMs). It is well known that the gravity field is a function of mass, and derivatives of the gravity field are a function of density. Considering the small depth of archaeological targets and their relatively large geometrical size, the observation of both vertical and horizontal derivatives of the gravity field will undoubtedly permit the obtaining of new and essential information about the desired targets. An integrated analysis of the gravity field and its vertical and horizontal derivatives will significantly expand 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 this parameter obtained by any transformation procedures. The W_zz values calculated from the field Δg_B, as a rule, show decreasing values compared with the W_zz obtained from physical measurements (or computation at two levels).

Calculations of W_xxz and W_zzz may also help delineate anomalies from closely disposed objects and remove the regional background, respectively. Considering that many archaeological sites in the world are situated near natural features (such as sea and lake basins, mountains, water channels, and forests) and industrial sites (including plants, buried pipelines, tunnels, and railways), this fact may have significant implications. For instance, Figure 5 shows that computing W_xxz (g_xx) makes it possible to reliably recognize gravity effects from the two closely located underground caves. Ancient Prehistoric caves corresponding to this PAM have been investigated in several archaeological sites located near the Golan Heights (northern Israel).

5. Examples of Quantitative Analysis

5.1. Model Examples

A model of archaeological cavity presented in Figure 6 occurs at the boundary of two layers with different densities. A quantitative examination of this gravity anomaly using conventional methods (e.g., Telford et al., 1993; Parasnis, 1997) yields significant errors. At the same time, employment of quantitative interpretation methods developed in magnetic prospecting (improved modifications of the tangents, characteristic points, and areal method) (e.g., Eppelbaum et al., 2001; Eppelbaum, 2019) enabled the determination of the center of this anomalous body with a high accuracy (Figure 6). Calculated parameter$M_{\Delta g}$ in this example consists of ${\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\cdot 7\mathrm{m}\cdot 0.23\mathrm{mGal}\cong 0.805\mathrm{mGal}\cdot \mathrm{m}.$

The following example (Figure 7) illustrates the potential of quantitative gravity analysis for identifying anomalies in the buried remains of the ancient Roman road. It is known that the ancient Roman roads in Israel were developed from the dense limestone, which is much denser than the surrounding deposits. The typical density of limestone is 2600 kg/m³, and the surrounding density was assumed as 2200 kg/m³ (according to the combined archaeological-geophysical data from the Beit Gouvrin site, central Israel (e.g., Magness and Avni, 1998; Eppelbaum, 2010)). The position of the ancient road along the strike was chosen as +1.5 m and -1.5 m. It can be seen from this figure that over the ancient Roman remains, a gravity anomaly of about 14 µGal is computed. Quantitative analysis of this gravity anomaly from a model of a horizontal thin bed (horizontal plate) is complicated. At the same time, after the calculation of W_zx (Δg_x), these two anomalies can be interpreted by the methodologies developed in magnetic prospecting (Eppelbaum, 2019). As a result, we obtain for these two anomalies that indicate the positions of the edges of the ancient road remain. Quantitative analysis of these anomalies reveals the sources located in the middle of the vertical thickness, which aligns nicely with the initial PAM.

5.2. Field Examples

Padín et al. (2012) discovered a buried ancient crypt in the Don Church located in the town of Alfafar (Valencia, Spain). For the quantitative analysis, an interpretative model of the horizontal circular cylinder (HCC) was applied, and the improved tangent method was used as an interpretation technique (Figure 8). As can be seen from the figure, the obtained center of the HCC is slightly shifted down the dip. This can be explained by some discrepancy between the real tunnel model and the interpretation model of the HCC. Gravity moment here is 0.0364 mGal·m.

The following example (Figure 9) demonstrates the employment of the quantitative analysis of microgravity anomaly using improved tangent and characteristic point methods for microgravity anomaly examination in the Roman Amphitheater (initial data after Di Filippo et al. (2005)). The HCC center obtained for the tunnel is slightly shifted upwards, which can be explained by the influence of adjacent walls with a relatively high density. The calculated gravity moment consists of 0.0624 mGal·m.

The last example of quantitative analysis shows the interpretation of a negative gravity anomaly from the well-studied buried chamber in the Boden Vean, St. Anthony Meneage (Cornwall, UK). Here, the improved tangent and characteristic point methods were applied (the HCC interpreting model was used). The interpretation results nicely agree with the available archaeological data (Figure 10). The gravity moment value is 0.0215 mGal·m.

6. The Developed Algorithm and Some Examples of Gravity Field 3D Modeling

6.1. The Developed Algorithm

The GSFC (Geological Space Field Calculation) program was developed to solve direct 3-D gravity and magnetic prospecting problems under complex geological conditions (Eppelbaum and Khesin, 2004). This program has been designed to compute the field of g (Bouguer, free-air, or observed value anomalies), ΔZ, ΔX, ΔY, ΔT, and the 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, L.H. non-closed, R.H. on-closed, and open. Geological bodies are approximated by horizontal polygonal prisms (Figure 11).

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 the rugged relief or at various arbitrary levels, using characteristic points (it is essential by the low-flying remotely operated vehicle application); (5) Simultaneous modeling of several profiles; (6) Description of a large number of geological bodies and fragments. The basic algorithm realized in the GSFC program is the solution of the direct 3-D problem of gravity and magnetic prospecting for a horizontal polygonal prism limited in the strike direction (Figure 11). In the developed algorithm, integration over a volume is realized on the surface, limiting the anomalous body.

An analytical expression for the first vertical derivative of the gravity potential of (m-1) angle horizontal prism (Figure 11) has been obtained by integrating a common analytical expression:

$$W_{z\text{'}}=-{\displaystyle \underset{s}{\int }\frac{z}{\left(R+y\right)R}}dxdz{|}_{y_1}^{y_2},$$

(9)

where $R=\sqrt{x^2+y^2+z^2},$and S is the area of the normal section of the prism by the plane of xOz.

$$\begin{array}{l}{W}_{z\text{'}T}=\left\{\begin{array}{c}-G\sigma {\displaystyle \sum _{j=1}^{m-1}\left[\right.V_j\sin {\alpha }_j\left(\ln \frac{R_{12j}+y_2}{R_{22j}+y_2}-\ln \frac{R_{11j}+y_1}{R_{21j}+y_1}\right)}+\\ V_j\cos {\alpha }_j\left(\mathrm{sgn}\left(y_2{V}_j\right)\arccos {\displaystyle \frac{V_j^2{R}_{12j}{R}_{22j}+U_{1j}{U}_{2j}{y}_2^2}{r_{1j}{r}_{2j}(y_2^2+V_j^2)}}\right)-\\ \left(\mathrm{sgn}(y_1{V}_j)\arccos {\displaystyle \frac{V_j^2{R}_{11j}{R}_{21j}+U_{1j}{U}_{2j}{y}_1^2}{r_{1j}{r}_{2j}(y_2^2+V_j^2)}}\right)+\\ \cos {\alpha }_j\left(y_2\ln {\displaystyle \frac{R_{12j}+U_{1j}}{R_{22j}+U_{2j}}}-y_1\ln {\displaystyle \frac{R_{11j}+U_{1j}}{R_{21j}+U_{2j}}}\right)\left.\right]\end{array}\right\},\\ \end{array}$$

(10)

where G is the gravitational constant, σ is the density of the body, and α_j is the angle of the prism’s side inclination (see Figure 11).

A detailed description of analytical expressions for the first and second derivatives of the gravity potential of the approximation model of the horizontal polygonal prism and their connection with the magnetic field is presented in Eppelbaum (2019).

6.2. Examples of 3D Modeling

The PAM presented in Figure 12 reflects a real archaeological site located in the vicinity of Beit-Shemesh town (central Israel). The developed PAM was used to estimate expected gravity anomaly amplitudes, calculate the optimal step size of observations along the gravity profiles, and determine the distances between profiles.

The following figure demonstrates 3D gravity effects computed from a model of the buried furnace (Figure 13). The parameters of computation are given in the figure, for the initial model was assumed to be an ancient furnace discovered in the Tel Hara Hadid ancient metallurgical site (southern Israel) (e.g., Gilat et al., 1993). The Bouguer gravity anomaly at the level of 0.3 m above the earth is about 65 µGal, and at the level of 1.5 m, it is less than 40 µGal. The calculation of the second vertical derivative of the gravity potential $\partial g/\partial z$(W_zz) indicates that the graph is more localized, which makes this parameter more promising for the ancient furnace localization.

7. Discussion

7.1. The Ways of Microgravity Investigation Improvement

How can we improve the microgravity investigations in archaeological sites? Several perspective approaches can be suggested (see below).

7.1.1. Improved Procedures for the Surrounding Terrain Relief Correction

Considering that microgravity studies in archaeological sites typically yield very small values (ranging from several microGal to 200 microGal in some favorable cases), the precise calculation of terrain correction becomes critical. Some possible ways to improve the surrounding terrain relief corrections are considered in subsections 3.2.1 and 3.2.2. A comprehensive review of the available terrain correction methodologies is presented by Nowell (1999). An ultrahigh resolution global model of the gravimetric terrain corrections presented by Hirt et al. (2019).

7.1.2. Effective Gravity Data Transformation

Various transformations of gravity fields were tested: informational approach and third horizontal derivatives of gravity potential (Eppelbaum, 2009, 2011; Sari et al., 2020), entropy, wavelet and Morlet computations, and self-adjusting, information indicator, and adaptive filtering (Al-Zoubi et al., 2013; Eppelbaum, 2019). Hajian et al. (2012) proposed a local linear model tree for estimating cavity depths from the microgravity data.

The selection of the concrete transformation method strongly depends on the PAM of the artifacts under study, measurement accuracy, and the level of geological and environmental noise. Klokočník et al. (2014, 2017) successfully applied such non-conventional procedures as “strike angle”, “gravity invariants”, “virtual deformations (compression and dilatation)”, and “radial derivatives” in the regional gravimetry. There is no reason not to use these methods in the precise gravity investigations in archaeological sites.

7.1.3. New Procedures for Download Continuation

Undoubtedly, the procedure of the gravity download continuation has critical importance in archaeological geophysics. It is necessary to note several essential works in this field: Pašteka et al. (2012, 2020), Zhang et al. (2018), Cheng and Yang (2022), Li et al. (2023). Wang et al. (2023) carried out a helpful comparison of continuation method for the potential field data.

In the case of a borehole presence on the archaeological site (or in the vicinity of it), the borehole gravity gradiometry can be applied (e.g., Rim and Li, 2012). In turn, Zhdanov et al. (2020) have demonstrated the usefulness of the joint iterative migration of surface and borehole gravity gradiometric data.

7.1.4. Application of the Wavelet Approach and Diffusion Maps

The application of the wavelet approach and diffusion maps for advanced integrated microgravity anomalies, together with magnetic data, began in 2011 (Eppelbaum et al., 2011; Eppelbaum, 2015). An archaeological cave was chosen as an example, of which there are many in various regions of the world. The complex host geological media and an unfavorable signal-to-noise (S/N) ratio often prevent us from revealing the geophysical anomalies from the desired targets. Our approach to recognizing distinct geophysical field characteristics and integrating geophysical data involves applying modern developments in wavelet theory and data mining. The wavelet approach is applied to derive enhanced (e.g., coherence portraits) and combined images of geophysical fields simulated for typical PAM of archaeological cavities. The methodology is based on a matching pursuit with wavelet packet dictionaries (e.g., Averbuch et al., 2014, 2016, 2019), enabling us to extract desired signals even from strongly noisy data. The recently developed technique of diffusion clustering, combined with the wavelet methods, is utilized to integrate geophysical data and detect existing irregularities in the subsurface structure. The most crucial factor is that these obtained results can be applied to reveal the archaeological cave (based on a “machine learning” approach) from field geophysical observations. The combination of the above approaches enables the creation of a new methodology that enhances the reliability and confidence in the application of any individual geophysical method, as well as the integration of geophysical methods.

Wavelet methodology has a significant impact on cardinal problems of desired signal (in our case – gravity) processing, such as denoising, signal enhancement, and distinguishing signals with closely related characteristics. Based on the abovementioned, a three-phase approach was developed for the integrated gravity–magnetic (GM) localization of archaeological cavities. The first phase involves 3D computer modeling of geophysical effects from a few initial typical PAMs of caves. The second phase involves the development of signal processing approaches for analyzing the large number of various PAMs (a total of ninety types were utilized). In the third phase, a new methodology for combining different geophysical data interpretations (based on modern developments in the wavelet technique of signal and image processing) is presented.

Figure 14 illustrates an initial model of the archaeological cave occurring in a trivial two-layer medium. In this case, both gravitational and magnetic anomalies are clearly visible.

The following calculations (Figure 15) illustrate the geological section without an archaeological cavity, incorporating the effects of five bodies. Figure 16 presents PAM with the archaeological cavity occurring in a complex geological medium. It is the 32nd PAM from the 90s developed ones.

The developed wavelet approach to recognizing archaeological caves through the examination of geophysical method integration involves advanced processing of each geophysical method and the non-conventional integration of different geophysical methods with one another. Modern developments in wavelet theory and data mining are utilized to analyze integrated data. The methodology, based on matching pursuit with wavelet packet dictionaries, enables the extraction of desired signals even from strongly noisy data.

For analyzing the GM data, we employed a technique based on an algorithm that characterizes a GM image using a limited number of parameters (Eppelbaum, 2015). This set of parameters serves as a signature of the image and is used for discriminating between images containing a cavity and those that do not. The constructed algorithm consists of the following main phases: (a) collection of the database, (b) characterization of GM images, and (c) dimensionality reduction. Then, each image is characterized by the histogram of the coherent directions (Averbuch et al., 2019). As a result of the previous steps, we obtain two sets of A and N values for the signature vectors of images from sections containing archaeological cavity and its absence, respectively. The obtained 3D set of the data representatives can be used as a reference set to classify newly arriving data.

The interim goal is to extract the characteristic features of these two sets in a way that the maps of these sets in the feature space are well separated. Once this goal is achieved, these maps can be used as reference sets for the construction of classifiers, which automatically determine whether newly arrived data indicates the presence of a cavity or not. The goal is achieved through several steps. Either of the GM data vectors comprises 40÷60 computation points (step of computation was 0.5 m). Then, we couple these two vectors into one combined data vector (Averbuch et al., 2016). For further processing, we need the data vectors, whose lengths are a power of 2. Thus, we apply to each pair of vectors the following operations. Each vector is symmetrically extended via its left and right ends to a length of 64. The window is applied to the extended vector to avoid discontinuities when merging (Eppelbaum, 2015).

And, at last, the extended and mirrored gravity and magnetic data vectors are merged into one combined vector of length 128. The window is applied to the extended vector to avoid jumps when merging. Thus, it was observed that the curves differ by subtle geometrical details. To capture these differences and differentiate between N and A sets, wavelet packet transforms and dimensionality reduction using Diffusion Maps (Averbuch et al., 2014; Eppelbaum, 2015) were employed. At the end, the scattered projections of the data curves onto the diffusion eigenvectors. As a result of the above operations, the original data was embedded into a 3-dimensional space where data related to the A (archaeological cave) subsurface are well separated from the N (non-cave) data. This 3D set of data representatives can be used as a reference set for classifying newly arriving data. Geophysically, it refers to a reliable division of the studied areas into those containing archaeological caves and those devoid of caves.

7.1.5. Development of High-Precision Gravity Measurements on the Remotely Operated Vehicles

It is worth noting that an autonomous robot for detecting subsurface voids and tunnels using microgravity (a beneficial development) was developed approximately twenty years ago (Wilson et al., 2006). This research was continued in Wilson et al. (2009), but unfortunately, it has not been further developed or implemented. Stray et al. (2022) have reported on developments in the field of quantum sensing for gravity cartography, enabling the obtaining of high-resolution gravity signals with a low noise level in moving vehicles.

A significant number of archaeological sites are situated in coastal areas, where high-precision measurements of the gravitational field and gravity gradient tensor are feasible using autonomous underwater vehicles (e.g., Wu and Tian, 2010; Vu et al., 2024; Qiao et al., 2025).

The gravity observations on the unmanned aerial vehicles (UAV) (e.g., Kaub et al., 2018; Tang et al., 2019; Passey et al., 2020; Luo et al., 2022) have not yet achieved the accuracy necessary for localizing archaeological targets.

The recently published research by Peng et al. (2025) provides a comprehensive comparison of various instrumental possibilities for precise gravity field measurements.

7.2. Feasibility of Microgravity Application in the Archaeological Sites Worldwide

Analysis of the numerous archaeological-geophysical sources (e.g., Lakshmanan, 1991; Cuss & Styles, 1999; Tianyao et al., 2000; Di Filippo et al., 2005; Yuan et al., 2006; Abad et al., 2007; Batayneh et al., 2007; Chernov, 2009; Eppelbaum, 2009; Fais et al., 2015; Rabbel et al., 2018; Cesnek et al., 2019; Pašteka et al., 2020; Eppelbaum, 2022; AbouAly et al., 2023; Nicolas et al., 2024), as well as the author’s investigation (e.g., Eppelbaum, 2009, 2022), indicates that the ancient objects supposed for examination using microgravity survey may be classified (in order of decreasing) in the following way: (1) underground ancient cavities and galleries, (2) walls, remains of temples, churches, and various massive constructions, (3) crypts and tombs, (4) ancient pavements and aqueducts, (5) areas of ancient primitive metallurgical activity (including furnaces) [under favorable physical-geological environments]. Examining the presence of archaeological targets worldwide, it was assumed that the microgravity method could be effectively applied to at least 20-25% of ancient sites. Unfortunately, the review of literature indicates that the number of microgravity investigations is limited and does not reach the 20% mentioned. This can be attributed to an underestimation of the increased possibilities of microgravity, both instrumentally and interpretatively.

Considering the small depth of archaeological targets and their relatively minor geometrical size, observing both vertical and horizontal derivatives of the gravity field will undoubtedly provide new and essential information about the desired targets. An integrated analysis of the gravity field and its vertical and horizontal derivatives will significantly expand 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 this parameter obtained by any transformation procedures: the Δg_z values calculated from the field Δg, as a rule, show decreasing values compared with the Δg_z obtained from physical measurements.

7.3. Integrating Microgravity with Other Geophysical Methods

It was proven that, in some cases, a simple model of geophysical method integration can be used. Assume that q is the risk of the erroneous solution. Then the anomaly detection reliability γ is (Khesin et al., 1996)

$$\gamma =1-q.$$

(11)

If a set of methods is focused on investigating some independent indicators of equal value, the anomaly detection reliability γ can be described by an error function (probability integral) as:

$$\gamma =erf\left({\displaystyle \frac{{\displaystyle \sum _i\sqrt{{\nu }_i}}}{2}}\right),$$

(12)

where ν is the ratio of the anomaly square to the noise dispersion for each i-th geophysical field.

Now, let us assume that three points indicate the anomaly and that the mean square of the anomaly for each field is equal to the noise dispersion. For a single method, the reliability of the detection of an anomaly of a known form and intensity by Kotelnikov’s criterion is expressed by $erf\left(\sqrt{{\upsilon }_i}/2\right)$. Then, the value of reliability for individual methods is 0.61, and those of a set of two and three methods are 0.77 and 0.87, respectively (according to equation (12)). This means that the q value (risk of an erroneous solution) at i-th integration of two or three methods (see equation (11)) decreases by the factors of 1.7 and 3, respectively. A comparison of the risk with the expenditure enables one to find an optimum set of methods.

Microgravity is often combined with other geophysical methods, including GPR, resistivity, and magnetics (see Section 3 and subsection 7.1.4). Safarov et al. (2019) have integrated precise gravity studies in Azerbaijan with the examination of the radon gas anomalies. Pánisová et al. (2012) employed microgravity, combined with close-range photogrammetry. Finally, it seems promising to combine microgravity with the multi-remote sensing analysis and elements of information theory (e.g., Eppelbaum et al., 2024).

8. Conclusions

The typical kinds of noise appearing in the microgravity investigations are shown. The archaeological remains in the world are classified by the degree to which the microgravity method is applicable. A proposed scheme for the quantitative interpretation of gravity anomalies (based on the advanced developments in the magnetic field analysis) is briefly discussed, along with the various models and field examples. The new suggested parameter – ‘gravity moment’ – will help classify the different microgravity anomalies. The characteristics developed by the software for combined 3D gravity-magnetic modeling in complex environments indicate that it is a powerful tool for the microgravity examination at archaeological sites. Two earlier non-conventional schemes for computing the surrounding terrain relief in precise gravity studies of ore deposits can be successfully applied to obtain the exact values Δg_B+TC in archaeological sites. The possibility of employing gravity field derivatives in some specific geological-archaeological situations is underlined. The developed physical-archaeological models of some typical ancient remains demonstrate the effectiveness of utilizing such models. The definite perspectives in the automatic recognition and arrangement of microgravity anomalies are associated with the perspective development of modern machine learning procedures, such as “diffusion maps”. The integration of microgravity with other methods is briefly discussed.