Review of Processing and Interpretation of Self-Potential Anomalies: Transfer of Methodologies Developed in Magnetic Prospecting

The self-potential (SP) method is one of the most inexpensive and unsophisticated geophysical methods. However, its application is limited due to the absence of a reliable interpreting methodology for the complex geological-environmental conditions. To exclude disturbances appearing in the SP method, a few ways for their removal (elimination) before quantitative analysis are presented. A brief review of the available interpretation methods is included. For the magnetic method of geophysical prospecting, special quantitative procedures applicable under complex physical-geological environments (oblique polarization, uneven terrain relief and unknown level of the normal field), have been developed. The detected common peculiarities between the magnetic and SP fields make it possible to apply the advanced procedures developed in magnetic prospecting to the SP method. Besides the reliable determination of the depth of anomalous targets, these methodologies enable the calculation of the corrections for non-horizontal SP observations and to determine the orientation of the polarization vector. For the classification of SP anomalies, is proposed to use a new parameter: the ‘self-potential moment’. The quantitative procedures (improved modifications of characteristic point, tangent techniques and the areal method) including the determination of the SP vector and SP moment, have been successfully tested on models and employed in real situations in mining, archaeological, environmental and technogenic geophysics. The obtained results indicate the effectiveness of the presented methodologies.


Introduction
The Self-Potential (SP) method is based on the study of natural electric fields (in some cases, this method is named as 'spontaneous polarization'). The term "natural" here means that the field is not created by an external artificial source. Permanent electric fields arise in the course of redox, filtration, and diffusion-adsorption processes in the upper part of the geological section. The registration of these fields is the goal of the SP method, and the geophysical interpretation of the parameters generating this field is the main purpose of SP data examination. An oxidizing object (e.g., ore body, archaeological or environmental target) is a galvanic cell, the existence of which requires: (1) the contact of electric conductors with different types of conductivity (electronic and ionic), and (2) the difference in the redox conditions at different contact points of these conductors. An appearance of these conditions is usually impossible without underground water contact [1].
In the geological section, the conditions for the formation of a galvanic cell arise on the targets with electronic conductivity, if these bodies are in the water-saturated rocks with ionic conductivity. The change in the redox conditions at the contact point of the electronic conductor (anomalous target) and the surrounding medium is associated with a decrease in the oxygen content with a depth.
Fox's [2] SP observations at copper vein deposits in Cornwall (England) laid the foundation of the application of all electric methods in geophysics as a whole. SP is an effective, prompt and comparatively simple geophysical method. Equipment for SP method is one of the most inexpensive in the field of geophysics (Table 1). Table 1. Average prices for geophysical potential field equipment (on the basis of author's personal contacts with the various manufacturers of geophysical equipment (January 2021)).

Method
Gravity Magnetic Resistivity Self-Potential Price of equipment US $ 60,000-110,000 20,000-25,000 35,000-55,000  Conventional equipment employed in the SP method consists of a microvoltmeter, a pair of non-polarizable electrodes, cable and CuSO 4 solution (the latter is necessary for the better contact of employed electrodes with the environment).
Usually (not always, depending on geological and environmental peculiarities of studied deposit (cite)), ground penetration radar (GPR) and electric resistivity tomography (ERT) produce more detailed geophysical-geological information. However, they are much more expensive and, most importantly, water content in subsurface strongly complicates application of these methods. At the same time, the presence of water for the SP method is only positive factor, since it enables to increase SP anomaly intensity, e.g., [3,4].
The values of SP anomalies vary in wide range: from fractions of a millivolt in detailed studies to several Volts in regional surveys (the origin of the latter has not been determined until recent time) [3,8,12,13,15,23,26,30,32,34].

Self-Potential Observations: Analysis of Disturbances
The main kinds of noise appearing in the SP method are compiled in a block-scheme ( Figure 1). Some of these noise effects are considered below in detail.

Electrode Noise in SP Method
Although the fact that SP electrode is called as "non-polarizable", after some time it accomplishes some polarization effects from the surrounding media. However, taking into account that we measure the value ∆U (U 1 − U 2 ), it is important to keep not an absolute non-polarizability, but an equivalent polarization on both of the applied electrodes. For checking this equivalent, the following procedure can be employed in field conditions (of course, measurements in a physical laboratory are more precise). We can write a trivial equation for the first electrode: U 1 + e 1 (U 1 is the first "medium" signal, and e 1 is the noise accumulated in the first electrode). For the second electrode, correspondingly, we have U 2 + e 2 (U 2 is the second "medium" signal, and e 2 is the noise of accumulated in the second electrode). According to [3], we measure: where ∆U 1 is the self-potential field generated between the first and second SP electrodes. Let us change electrodes by their places. In this case, we will obtain: where ∆U 2 is the self-potential field generated between the second and first SP electrodes. After this, calculating a difference between ∆U 1 and ∆U 2 , we receive: If the value (e 1 − e 2 ) is significant (e.g., ≥3 mV), the noised electrodes must be replaced by new ones. A similar methodology for the electrode noise detection was suggested in [36].

Temporal Variations in SP Method
Parasnis [4] carried out SP measurements of the same profiles in the Akulla region (Sweden) seven times during the period of 1960-1967. These measurements indicate a good repeatability despite the fact they were conducted under different climatic conditions. Similar investigations performed by other researchers, e.g., [3,37], basically confirm the good repeatability of the different SP time observations. For the estimation of accuracy ε of SP field measurements, the following trivial formula is often employed in various geophysical methods: where N is the total number of SP observations, and δ = (∆U init − ∆U repeat ), where ∆U init means the ordinary (initial) measurements, and ∆U repeat -the repeated measurements.
The number of repeated measurements should be at least 8-10% of the total number of observations. If the value of ε exceeds some a priori assumed value (this value usually depend on the concrete spread of the SP amplitudes), results of SP survey can be rejected as non-reliable.

Terrain Relief Correction
In the SP method, terrain relief influence is two-fold. On the one hand, negative SP anomalies caused by electromotive force can be created over the positive topographic forms (this phenomena strongly depends on the peculiarities of underground water circulations). Comparison of the SP graphs with topographic data usually allows for the identification of anomalies of this type by the character mirror image of the terrain relief in SP anomalies [38].
From other side, as follows from the very detailed SP measurements of Ernstson and Schrerer [39], at the inclined topographic surface, the SP field increases directly with increases in the relief highs ( Figure 2). In the last case, for the elimination of the terrain relief influence, a correlation method developed in magnetic prospecting [38], and VLF studies [40] can be applied. The essence of this method is as follows. The method employs, for the removal of the terrain relief effect from the observed field, ∆U obser , a linear least-squares relation ∆U appr (application of more complex approximations (e.g., parabolic, etc.) is also possible): where h is the height of relief, b is the angle coefficient, and c is the free member.
Value ∆U appr approximates the observed field as a function of elevation h (anomalous zones usually do not include the correlation field), and then we obtain the corrected (residual) field ∆U corr , where the relief influence is essentially eliminated: It should be noted that this correction only eliminates the effect of the inclined topographic masses with certain electric properties. Special methods allow for the calculation of the difference in altitudes of SP observations to the anomalous target on an inclined profile, which are considered later in this paper.
Wang and Geng [41] studied the problem of terrain correction in the SP method in detail in several field examples. They noted that the mechanism of SP anomalies formed by terrains is rather complex, and therefore it is difficult to obtain the corresponding analytical formulas. Three types of relief fitting were applied: (1) linear, (2) quadratic and (3) exponential. After careful analysis, it was concluded that the linear fitting is more optimal since it does not create factious anomalies [41]. Thus, this investigation confirms the application of the aforementioned linear least-squares correlation method.

Calculation of SP Anomaly Distortion Due to Observations on Uneven Surface
SP anomalies (and anomalies of other potential fields) distort due to observations on uneven surfaces (and correspondingly, from different distances to anomalous objects). This disturbing effect usually is calculated at the end of the interpretation process (see Section 4.1).

Net Justification in Areal Observations
Net justification of SP data is conventionally performed by the use of procedure identical to justify the observations in gravity and magnetic prospecting, e.g., [42]. Other strategies for removing this effect are presented in [43].

Influence of Meteorological Factors
Many scientists note that rain increases the intensity of SP anomalies, e.g., [3,4,43]. Therefore, occasionally, an artificial irrigation of site intended for SP research is recommended.

Presence of Magmatic Associations
Obviously, the development of magmatic associations (or other kinds of hard geological rocks) in a site destined for field investigations does not allow for the grounding of SP electrodes. In such cases, it is impossible to carry out the SP measurements. The same reason may limit the water circulation at the subsurface, which can weaken, or even completely cancel a generation of the SP anomalies.

Some Environmental Factors
The SP anomaly level may affect some environmental factors. One of these factors is the shadowing of a part of the investigated area. For instance, Revil and Jardani [43] documented the fact that the difference between the SP electrodes placed in cold and warm media may exceed 10 mV. Another factor is a presence of some hygrophilous plants (e.g., hazel and almonds) whose roots can pick over a lot of moisture from the upper part of the geological section (thereby hindering the SP anomaly generation).
Ernstson and Schrerer [39] monitored SP and temperature anomalies during 15 months in 1980-1981 ( Figure 3). The correlation between SP and soil temperature is interpreted as result of the influence of the thermal diffusivity and convection processes in the subsurface [39]. Between these parameters, a correlation of r = 0.64 was established. It is a low relationship, but in any case, should be taken into account. Perrier and Morat [44] also detected a correlation between the amplitudes of the SP variations and soil temperatures. According to these authors, the state of the soil in the first 30 cm seems to play an important role. The authors concluded that the joint monitoring of electric potential and temperature appears to be a powerful tool to monitor the underlying soil processes [44].
An essential relationship between the SP data and temperature values recorded in boreholes at intermediate and large depths was established for the Cerro Prieto geothermal field (Baja California, NM, USA) [45]. An existence of this phenomenon under normal thermophysical conditions is not studied yet in detail.

Review of Available Quantitative Interpretation Methods
The calculation of theoretical anomalies due to SP has long been based primarily on Petrovsky's [46] solution, which was derived for a vertically polarized sphere [37]. Later on, some substantial solutions for sheet-like bodies and inclined plates were obtained [3]. The electric polarization vector was generally considered to be directed along the sheet-like body dipping (along the longer axis of the conductive body).
Other methods of SP anomaly quantitative interpretation includes anomalous body with a simple geometrical shape which approximates the anomaly source. Its parameters (i.e., the depth, the angle between the horizon and the direction of the polarization vector) are usually determined by: (1) graphically using characteristic points of the anomaly plot, or (2) by trial-and-error method consisting of visual comparison of the observed anomaly with a set of master curves [3].
Zaborovsky [37], Semenov [3], and Murty and Haricharan [47] applied the SP anomaly generated by plates to the following formula: where j is the current per unit length, ρ is the host medium resistivity, r 1 and r 2 are the distances from the plate left and right ends to the observation point.
Fitterman [48] gave the method of SP anomaly calculation for field sources of an arbitrary shape based on numerical integration using Green's function. This potentially promising approach is highly computer intensive and does not provide sufficient accuracy.
A few groups of authors performed spatial-frequency analysis of SP anomalies produced by polarized bodies of various geometrical shapes [49][50][51][52][53]. For instance, Rao et al. [49] analyzed the Fourier amplitude and phase spectra of SP anomalies to evaluate the parameters of the sheet; the application of this method on two anomalies (synthetic and field data) has given satisfactory results. 1D and 2D spatial frequency analysis of SP anomalies from the polarized sphere initiated in [52] was successfully evaluated in [53].
A series of publications [54][55][56][57][58] provided a large number of methodological approaches. However, these approaches (based mainly on application of the gradient analysis and calculation of derivatives) by their common usefulness, have not caused a quantitative jump in this field.
Gibert and Pessel [59] applied the continuous wavelet transform for the localization of SP anomalies. The wavelet analysis provided both an estimate of the location and of the nature of the target responsible for a given self-potential signal. The wavelet-based techniques as well as analytic signals were employed to interpret SP anomalies caused by subsurface fluid flow in Mt. Etna [60].
Patella [61] suggested an application a tomographic presentation of SP images. It consisted of scanning the geological-geophysical section through SP profiles, by the unit strength elementary charge, which is given a regular grid of coordinates within the section. At the each point of the section the charge occurrence probability function is calculated. The complete set of calculated grid values is employed to draw colored sections. This method was evaluated in [62].
Mendonça [63] employed the Green functions to simplify the evaluation of SP anomalies from buried conductors. This approach was used by this author to simulate geoelectric targets in mineral exploration and to obtain current source terms by inverting the SP data set.
Srivastava and Agarwal [64] developed the Enhanced Local Wave number technique wherein the nature of the causative source has been determined by computing structural indices based on its horizontal location and depth. This approach was tested on several mineral deposits with complex ore body distribution (among others, in the Rajasthan Copper belt (India) and sulphide deposit in Quebec (Canada)).
There are a number of recent interpretation techniques based on minimizing the difference between the observed and theoretical anomalies. The minimization is achieved by sequential optimization of the interpretation parameters through computer-aided iterations e.g., [65][66][67][68]. These techniques are usually effective under simplified physicalgeological conditions. An attempt to apply the neural network approach to compute the shape factor and depth of the causative target from SP anomaly was undertaken in [65]. The validity of the neural network inversion is demonstrated through the SP data taken from a graphite body in the southern Bavarian Woods (Germany). A comparable and acceptable agreement is shown between the results of inversion derived by the neural network and those from the real field data. Without denying the general promise of this approach, it should be noted that some inconclusiveness of examples given in this work.
Hristenko and Stepanov [66] demonstrated a system of 2D modeling of an SP field from several anomalous bodies with introduced Gaussian noise under conditions of uneven relief. Some applied transformations allowed for the exclusion of the effect of near-surface geological inhomogeneities from the total SP field and to enhance anomalies from the buried targets.
Gobashy et al. [67] proposed a method based on utilizing the optimization algorithm, which is an effective heuristic solution to the inverse problem of SP field due to a 2D inclined bed. The realization of this perspective algorithm under complex physicalgeological conditions should be continued. Rao et al. [68] proposed a global optimization methodology expressed in the development of inversion algorithm for 2D inclined plates. Absence of the geological sections in the mentioned work complicates the examination of this methodology.
The procedure based on the interpretation of potential anomalies due to simple geometrical structures using Fair function minimization [69] is of certain interest. Giannakis et al. [70] suggested a hybrid optimization scheme for SP measurements due to multiple sheet-like bodies. This procedure demands a wide verification on concrete field examples. Of special interest is the recent research of [31], who developed the least-square subspace preconditioned method to compute the known Tikhonov solution to reliably detect the depth and the shape of shallow electrical current density sources.
Fedi and Abbas [11] proposed the 'depth from extreme points method' and employed it on several models and field examples. The method yields estimates of the source horizontal location, depth (top or center), and geometry.
Kilty [71] published a paper that acknowledged the analogy between the current density of SP and magnetic induction. This author suggested interpreting SP anomalies based on the conventional methods developed for magnetic prospecting. However, the proposed methodologies are not effective for complex physical-geological conditions. A similar approach, but with the improved interpretation methodology (specially developed for conditions of the rugged terrain relief, oblique polarization and unknown level of the normal field) was proposed in [38]. Eppelbaum and Khesin [72] proposed a new elaboration of the interpretation process. The present paper shows a final generalization of this approach.

Some Common Aspects of Magnetic and SP Fields
The magnetic field is a potential one (when value of target's magnetization is not very high) and satisfies Poisson's equation: where U a is the anomalous magnetic field and V represents the magnetic potential. SP polarization is generated by the spontaneous manifestation of electric double layers on contacts of various geological (or environmental and artificial) objects. The electric fields E of the electric double layer l caused by natural electric polarization are defined as the gradient of a scalar potential Π i : The potential Π i satisfies Laplace's equation everywhere outside of the layer l [73]. Jardani et al. [74], in addition to Equation (10), proposed the following expression for explanation of SP pattern generating by geothermal flow: where σ is the conductivity (in S/m), J s is the current density vector (in A/m 2 ), and ζ is the volumetric current density (in A/m 3 ). Equation (11) also indicates a potential character of the SP field. Analytical expressions for some of the interpreting models for magnetic and SP fields are presented in Table 2 (after [38]). Table 2. Magnetic and SP fields: Comparison of analytical expressions for some interpreting models.

Field Analytical Expression
Magnetic

Self-Potential
Horizontal circular cylinder (HCC) (14) Sphere Here, Z v is the vertical magnetic field component at vertical magnetization; J is the magnetization; b is the horizontal semi-thickness of TB; m is the magnetic mass (point pole magnetic charge); ρ 1 is the host medium resistivity; ρ 2 is the anomalous object (HCC or sphere) resistivity; U 0 is the potential jump at the source body/host medium interface; r 0 is the polarized cylinder radius; R is the sphere radius; x is the current coordinate; z is the depth of the upper edge of TB (center of HCC or sphere) occurrence.
Formulas describing the potential character of magnetic (Equation (9)) and SP (Equation (10)) fields are identical. The proportionality of analytical expressions (12)- (15) for magnetic and SP fields is obvious. It allows us to employ advanced interpretation methods developed in magnetic prospecting for conditions of rugged terrain relief in SP data analysis, the superposition of anomalies of different orders and oblique polarization (we assume that the SP polarization vector is analogue of the magnetization vector). It is supposed that the majority of interpretation methodologies developed for the magnetic and gravity fields are applicable for the SP method.
In several works, the application of some advanced procedures from the most developed potential fields (gravity and magnetics) was demonstrated using self-potential data. Akgün [75] applied the Hilbert transform (usually employed in gravity and magnetic prospecting) for analysis of SP data. Sindirgi et al. [76] successfully tested the SP anomalies method of the total normalized gradient developed in gravity prospecting. This investigation was continued in [77] with the application of the Euler deconvolution. Extended Euler deconvolution techniques were successfully tested in some mineral deposits and observations in deep boreholes [78]. Biswas [79] suggested the study of SP anomalies using the 2D analytic signal developed in magnetic prospecting.
Of special interest is the recently published research [80] where the author proposed the employment of the posterior distribution model of the SP anomaly inversion. The advantages of this approach are that the SP data may contain single and multiples of SP sources and this method does not require prior assumptions over the shape of the anomaly source.

Quantitative Analysis of SP Anomalies by the Use of Advanced Methodologies Developed in Magnetic Prospecting
The improved methods for SP anomaly analysis include characteristic points, tangents and areal methods (these methods are described in detail in the publications suggested to magnetic anomaly interpretation: e.g., [38,40,72,81,82]. Formulas for interpreting SP anomalies by use of the characteristic point method are presented in Table 3. Several figures below display some peculiarities of the characteristic point and tangent methods application. Table 3. Formulas for quantitative interpretation of magnetic anomalies from bodies approximated by thin bed (TB) and horizontal circular cylinder (HCC) using the improved characteristic point method (after [40], with modifications).

Formulas for Calculation of Parameters Necessary for Quantitative Analysis
Thin Bed HCC Thin Bed HCC Indices "0" and "c" designate the thin bed (TB) and horizontal circular cylinder (HCC) models, respectively. Values h 0 and h c are the depths to upper edge of TB and center of the HCC, respectively. Parameter ∆h designates measurements of self-potential field at different depths of the electrodes' grounding.
Preliminary results for SP anomalies to estimating HCC radius and length of horizontal upper edge may be obtained from 3D magnetic field modeling (taking into account a common similarity of these fields) ( Table 4). Table 4. Nomenclature of variables applied for quantitative analysis of SP anomalies due to of thin bed and horizontal circular cylinder models (see Table 3).

Variable Description
θ Generalized angle reflecting the degree of SP anomaly asymmetry as a function relation of an anomalous body depth of occurrence, geometric form, value of polarization  The improved versions of tangent and areal methods are presented in detail in [40,81,82].

SP Observations on an Inclined Profile
When potential geophysical anomalies are observed on an inclined profile, the obtained parameters characterize some fictitious body [13]. The transition from the parameters of fictitious target to those of real target is realized using the following expressions: where h is the depth of the body upper edge occurrence (or HCC (sphere) center), x 0 is the shifting of the anomaly maximum from the projection of the center of the anomalous body to the earth's surface (produced by an oblique polarization), and ω 0 is the angle of the terrain relief inclination (ω 0 > 0 when the inclination is toward the positive direction of the x-axis), the subscripts "r" and "f " stand for parameters of real and fictitious bodies, respectively. The direction of the electric self-polarization vector ϕ p is calculated from the expression: and on an inclined relief: The performed calculation of the vector ϕ p direction on concrete field examples indicates that for the interpreting models closed to the model of a thin bed, the direction of this vector approximately coincides with the anomalous body-dipping. It enables us to obtain a supplementary interpretation parameter.
Besides the geometric parameters of the anomalous target, the self-potential moment can also be determined (see Table 3). For the models of thin bed and HCC, the self-potential moment can be calculated by the use of Equations (19a) and (19b), respectively (see also  Tables 3 and 4).
where ∆U a is the amplitude of SP anomaly (in mV), h 0 is the depth of the upper edge of thin bed (in meters), h 2 c is the squared depth to the center of the HCC (in m 2 ), and θ is the some generalized angle (see Tables 3 and 4). The self-potential moment, by analogy with the magnetic field analysis, can be used for classification of various SP anomalies (and, correspondingly, hidden targets). Comparison of these values between themselves (separately for interpreting models of thin bed and HCC) indicates size of anomalous bodies, their depth and space position as well as ratio of electrochemical properties of the 'anomalous target/surrounding medium'.
Initial methodologies for quantitative analysis of magnetic anomalies under complex physical-geological conditions for the model of the thick bed were presented in [38] and its significant evaluation (including the intermediate models between the thick and thin beds) was presented in [82].
For observation on an inclined profile, the real self-potential moment can be calculated as follows: M ∆U,r = M ∆U,f cosω 0 (20) Here, the subscripts "r" and "f " stand for a parameter of real and fictitious selfpotential moments, respectively.
Undoubtedly, the calculation of all aforementioned parameters from SP data should be connected to a unified computerized system with the slight participation of an interpreter.
For the testing of some SP anomalies, the developed software for the 3D computation of the magnetic field (e.g., GSFC [13]) was applied. In this case, magnetic vector orientation can be utilized as analogue of the self-potential vector.

Quantitative Analysis of SP Anomalies
The developed interpretation system in the SP method is applicable for complex physical-geological conditions (oblique polarization, inclined relief and superposition of SP anomalies of various orders).

Testing on Theoretical Models
First of all, the aforementioned interpretation methods were successfully tested on the theoretical SP anomalies for models presented in [3,68,83].

Chyragdere Sulfur Deposit (Central Azerbaijan)
It is interesting to compare SP studies carried out over the Ghyragdere sulfur deposit  This figure shows that the mining works in the underground shaft (1930, 1937 and 1938) strongly distort the SP field observed at the earth's surface (distance from the observation points to ore deposit consisted several tens of meters). This testifies to the tight correlation between the mining processes and SP anomalies. It would be fascinating to separately compare the volumes and contours of the mined ore with the SP isolines for the abovementioned years, but these documents have been lost over the past years.

Sariyer Sulphide-Pyrite Deposit (Near Istanbul, Turkey)
Yüngül [85] documented the results of the survey in the Sariyer area (Istanbul). The performed interpretation indicates that the obtained position of HCC center is in the line with geometrical and physical parameters of the sulphide-pyrite ore body ( Figure 5). Here and in some other figures, displayed parameters d 3 and d 4 relate to the improved tangent method (this method is described in detail, for instance, in [81]). Calculating the selfpotential moment by the use of Equations (19b) and (20), we obtain M ∆U = 31, 800 mV · m 2 . This value testifies an intermediate size of the ore body. The calculated direction of selfpotential vector by use of Equation (18) is a vertical one.  Figure 6 displays results of SP anomaly quantitative interpretation using characteristic points and tangent methods (areal method based on the calculation of the area occupied by the SP anomaly has also been applied). The interpretation results, as can easily see from Figure 6, have a good agreement with the location of the ore body. The self-potential moment (here, a model of thin bed was selected and Equations (19a) and (20) were applied) M ∆U = 1 2 60 mV · 6.5 m · 0.93 = 181 mV · m. The obtained value indicates the small size of the hidden ore body. The calculated direction of self-potential vector (Equation (18)) roughly corresponds to the dipping of the ore body ( Figure 6).

Katsdag Polymetallic Deposit (Azerbaijan)
Three SP anomalies were successfully interpreted in the Katsdag copper-polymetallic deposit (southern slope of the Greater Caucasus, Azerbaijan) under conditions of rugged terrain relief (Figure 7). Anomalies 1 and 2 are intensive ones, but anomaly 3 is comparatively small. It is important to underline here an essential difference between the quantitative results of SP anomalies analysis calculated without and with estimation of the rugged relief influence. The interpretation results are shown for anomalies 1-3 (Figure 7), but the SP moment is calculated only for anomaly 1 (after applying Equations (19a) and (20)) is M ∆U = 1 2 180 mV · 20 m · 0.984 = 3450 mV · m. This value indicates an average size of the buried ore body.

Filizchai Polymetallic Deposit (Azerbaijan)
A very intensive SP anomaly (almost 500 mV) was observed in the portion of Filizchai copper-polymetallic field (southern slope of the Greater Caucasus, Azerbaijan) under the conditions of highly complex terrain relief (Figure 8). The results of the interpretation (improved methods of characteristic points and tangents to the thin bed interpretation model were applied) also indicate significant difference of position of the anomalous body calculated without rugged terrain relief influence (blue circle) and after the calculation of this influence (red circle). The calculated SP moment (after Equation (19a)) is M ∆U = 1 2 440 mV · 90 m = 19, 800 mV · m. The employment of Equation (20) gives us some decreased value: 14,600 mV·m. It is a sufficiently high value of the SP moment, showing the presence of a large ore body. However, such large SP anomalies of ore origin are rarely observed. The direction of the self-potential vector was calculated by the use of Equation (18). The position of this vector agrees well with the dipping of this pyrite-polymetallic body (Figure 8).  Here, three different interpreting models were utilized (thin bed, HCC and thick bed) ( Figure 10). All three applied models are suitable ones. The calculated SP moment (for HCC model) is (after employing Equations (19b) and (20)) M ∆U = 21, 890 mV · m 2 . This value corresponds to the intermediate ore bodies. The calculated position of the polarization vector (HCC model) coincides with the dipping of the polymetallic body ( Figure 10).

Canyon Makhtesh Ramon (Negev Desert, Southern Israel)
The Makhtesh Ramon erosional-tectonic depression (canyon), 40 km long and approximately 8-km wide, is situated in the Negev Desert (southern Israel), 65 km southwest of the Dead Sea. On the basis of integrated geological-geophysical investigations in this area several tens of microdiamonds (largest sample is 1.35 mm) and a large amount of the mineral-satellites of diamond were detected [96]. Many geological-geophysical indicators showed that at least a part of the indigenous sources (kimberlites or lamproites) of the aforementioned minerals can occur in subsurface in this area. The compiled SP map (Figure 11a) displayed the presence of some anomalous zones. Quantitative analysis of the SP anomaly was carried out along profile A-B (Figure 11b) crossing one of the mentioned zones. Results of the performed interpretation indicated that the anomalous inclined body (having a geometrical form close to the thin bed model) occurred at the depth of 40 m, which agrees with the preliminary available geophysical data.
The territory of Israel contains more than 35,000 discovered archaeological sites of different age and origin. For SP observations, several typical archaeological sites located in different areas of the country were selected [18,99]. All SP measurements were performed using microvoltmeter with high input impedance and distinctive non-polarizable electrodes (Cu in CuSO 4 solution). The interpretation results obtained earlier at these sites were revised and generalized (the unified methodologies were employed) [100].

Roman Site of Banias (Northern Continuation) (Northern Israel)
The remains of the city of Banias are located in northern Israel at the foot of Mt. Hermon. Banias was the principal city of the Golan and Batanaea regions in the Roman period and occupied an area of more than 250 acres [101]. Here, different ancient remains of Roman and other historical periods were found. The area of the present study is located several km north of the well-investigated Banias site. In the nearest vicinity of the area of geophysical examination (SP and magnetic surveys), the remains of an ancient Roman cemetery and aqueduct were discovered [102]. Mineralogical and geochemical analyses of the excavated Roman chambers indicated that these objects were composed from the special type of hot-worked limestone.
SP observations were carried out by the grid of 1 × 1 m. The compiled SP map ( Figure 12) nicely indicates two anomalies. Interpretation profiles I-I and II-II cross the centers of these anomalies ( Figure 13). The Upper edges of the recognized anomalous targets occur at the depth of 1.1-1.3 m. Presence of these anomalous sources was confirmed by a comprehensive magnetic data analysis. Angle ϕ p for anomaly I consists of 75 • , but it is calculated from the opposite side (due to inversion of parameters d 3 and d 4 ). SP moment for the anomaly I (thin bed) is M ∆U = 1 2 16.5 mV · 1.2 m = 9.9 mV · m. Interestingly, this SP moment is almost 1500 times smaller than the same parameter calculated for the giant SP anomaly in the Filizchai deposit, Azerbaijan (see Figure 8). Such large differences in the aforementioned SP moments is caused by significant differences in the size of these anomalous targets and their electrochemical parameters. For anomaly II, taking into account that parameters d 3

Nabatean Site of Halutza (Southern Israel)
The Halutza site is located 20 km southwest of the city of Be'er-Sheva (southern Israel). It was the central city of southern Palestine in the Roman and Byzantine periods and was founded as a way-station for Nabatean (7th-2nd centuries BC) traders traveling between Petra (Jordan) and Gaza. This site was occupied mainly throughout the Byzantine period (4th-7th centuries AD) [103,104].
At this site, self-potential and magnetic measurements were carried out in a 20 × 10 m area with a 1 × 1 m grid [18]. The buried targets (ancient Roman limestone constructions) have produced negative anomalies in both geophysical potential fields. The results of the quantitative examination (here interpretation models of thin bed were utilized) are practically identical ( Figure 14A,B). Amplitude of SP anomaly reached 40 mV and is the largest from the anomalies considered in this site. The depth of both of these anomalous targets is about 0.85 m. The calculated moment for the self-potential anomaly is M ∆U = 20 mV · m. The ϕ p angle is calculated from the opposite side and consists of 70 • ( Figure 14B). It may be concluded that the recognized anomalous target approximated by thin bed in the SP method has no vertical dipping, but coinciding with the ϕ p angle. The obtained quantitative parameters of ancient constructions agree with the results of archaeological excavations performed in the vicinity of this site.

Christian Site of Emmaus-Nikopolis (Central Israel)
The Christian archaeological site, Emmaus-Nikopolis, is well known in ancient and Biblical history. The site is situated roughly halfway between Jerusalem and Tel Aviv (central Israel). The Crusaders rebuilt it on a smaller scale in the 12th century [101]. Nikopolis is displayed in almost all Christian Pilgrim texts from the fourth century onward; a majority of archaeological sources of this site are named as Emmaus-Nikopolis. Many scientists note that this site is characterized by a multilayer sequence, e.g., [104].
SP measurements in this site were performed by the grid of 1 × 1.5 m. In the compiled SP map (Figure 15), one local anomaly was selected for quantitative analysis (Profile A-B in Figure 15). The determined depth of the target upper edge is about 1.5 m (depth of the HCC center is about 2.1 m) ( Figure 16). The ϕ p angle consists here of 85 • and is calculated from the opposite side.    Figure 15). The black arrow shows the direction of selfpolarization vector.
The self-potential moment of this anomaly (HCC model is employed) consists of 4.5 mV·m 2 . This low value is mainly caused by the small depth and size of the cave. Fragments of some of the glass vessels discovered in this burial cave allowed for it to be attributed to the Byzantine period. Interestingly, to note that magnetic field examination allowed for the recognition in the same cave by the integrated effect from several small magnetic anomalies produced a glass vessel (made from fired clay) [13].

Buried Cavities in Dolomitic Limestone (Southern Italy)
Several impressive examples of SP application for the detection of underground cavities in southern Italy were displayed in [25]. Let us consider one of these field cases, where buried karst cavity exists in dolomitic limestones ( Figure 17). The cavity is horizontally extended and a significant SP anomaly (up to 100 mV) over it was observed. Quantitative examination along profile A-B crossing the center of this anomaly has been performed ( Figure 18). For interpretation models of thin bed (upper edge) and center of HCC were obtained depths of 6.0 and 9.5 m, respectively. The self-potential vector is oriented nearvertically. Self-potential moment of anomaly from this cave is about 3300 mV · m 2 .   Figure 17). The red cross designates the position of the middle of the thin bed upper edge, the red circle shows the position of the center of the horizontal circular cylinder (HCC). The black arrow shows the orientation of the polarization vector.

Cavities in the Djuanda Forest Park (Bandung, Indonesia)
The next example displays the results of SP and electric resistivity observations carried out above cavities (built during the WW II in early 1940s) at the Djuanda Forest Park, Bandung, Indonesia. Fascinatingly, the resistivity section ( Figure 19B) nicely shows two bright anomalies whereas the SP graph ( Figure 19A) indicates a significant anomaly (amplitude is more than 20 mV) over cave II, whereas the anomaly over cave I is only emerging (its amplitude is about 2 mV). Obviously, this fact is associated with the hydrogeological peculiarities of the geological section of the subsurface. Quantitative analysis of SP over cave II (here, an interpretation model of thin bed was applied) gave satisfactory results generally coinciding with the results of resistivity section. Calculated self-potential moment here is M ∆U = 52 mV · m.  Figure 20 displays the SP and resistivity graphs over the subvertical-fissured zone. SP quantitative examination showed significant disagreement between the results of interpretation and the available geological section. However, it is interesting to note that the performed quantitative analysis of the resistivity curve ρa (presented in the upper part of Figure 20) by the use of the same methodology gave similar results. Theoretical possibilities of such an analysis were reported in Eppelbaum [27] and evaluated in Eppelbaum [13]. Shevnin [115] also indicated a good correlation between the SP and resistivity methods.
Obviously, the disagreement can be explained by some erosion of the upper part of anomalous body (fissured zone) and corresponding changes of its physical properties (possibly appearing to be close to the physical properties of the host media). Nevertheless, the orientation of the self-potential vector coincides with the fissured zone dipping ( Figure 20). The value of the self-potential moment (model of thin bed was used) is estimated as M ∆U = 23 mV · m. The low value of the SP moment apparently reflects a small electrochemical ratio of the 'anomalous target/surrounding medium'.

Technogenic Geophysics
Let us designate technogenic geophysics as geophysical studies applied to the detection or determination of certain parameters of hidden modern industrial objects. The SP method fits well in such studies.
Onojasun and Takum [34] successfully applied SP investigations for the localization of an underground concrete water pipeline at the Kwinana industrial area of Southern Perth (Western Australia). A comprehensive examination of SP imaging over a metallic contamination plume was performed in [35]. The authors concluded that the SP method can be successfully used to monitor the underground metallic contaminants.
It was established that the examination of SP anomalies is significant for the localization of corrosion in buried oil, gas and water pipes, e.g., [31,32,106,116,117].
Underground Metallic Water-Pipe (Southern Russia) Fomenko [33] presented a case of a typical SP field distribution over the buried metallic water-pipe ( Figure 21). This anomaly has been interpreted by the use of tangent and characteristic point methods. The obtained position of the HCC center agrees (with some assumption) with a center of the hidden water-pipe. The calculated self-potential moment M ∆U = 79.8 mV · m 2 . Figure 21. Quantitative examination of the SP anomaly from the buried metallic pipe (southern Russia). Observed SP graph and near-surface section are taken from [33]. The red circle indicates the determined position of the center of HCC, and the red arrow shows the orientation of the self-potential vector.

Generalization of the Calculated Self-Potential Moments
The calculated self-potential moments for the variety of studied targets are compiled in Table 5. The values of self-potential moments presented in Table 5 demonstrate a wide range of the parameters calculated for the models of thin bed and horizontal circular cylinder. M ∆U values can be divided in three groups: (1) comparatively large values corresponding to comparatively big ore bodies, (2) middle values relating to targets studying in environmental and technogenic geophysics, and (3) relatively small values reflecting archaeological targets. However, even the smallest M ∆U value calculated, for instance, for the ancient cave in the site of Emmaus-Nikopolis, has independent importance.

Conclusions
The self-potential method is one of the oldest and simultaneously non-expensive geophysical methods. One of its main preferences is that the presence of water in the subsurface does not limit its method capabilities. The various disturbances complicated the SP observations under different physical-geological environments are analyzed. The available interpretation methodologies are briefly discussed. The proved common aspects between the magnetic and self-potential fields enable for the interpretation of SP anomalies in the modern interpretive procedures developed for complicated environments in magnetic prospecting (oblique magnetization (polarization), rugged topography and an unknown level of the normal field). An explanation of these procedures was presented in detail.
These interpretation procedures applied for SP anomalies enable reliably geometric parameters to be obtained of buried, anomalous targets occurring in complex physicalgeological environments. The suggested calculation of direction of the self-polarization vector allows to estimate, in many cases, the dipping of anomalous objects. The proposed calculation of the self-potential moment can be used for the classification of the observed SP anomalies. Testing these procedures in numerous cases of mining, environmental, archaeological and technogenic geophysics in various regions of the world indicates the high practical effectiveness of these developed methodologies.

Data Availability Statement:
The field examples presented in the paper are the author's data (Caucasus, Israel) or were reproduced from the published sources (with the corresponding citations) for the following analysis.