A fast imaging method for the interpretation of self-potential data with application to geothermal systems and mineral investigation

We describe a rapid imaging approach for the interpretation of self-potential data collected along profile by some geometrically simple model of cylinders and spheres. The approach calculates the correlation coefficient between the analytic signal (AS) of the observed self-potential measurements and the AS of the self-potential signature of the idealized model. The depth, electric dipole moment, polarization angle, and center are the inverse parameters we aim to extract from the imaging approach for the interpretative model, and they pertain to the highest value of the correlation coefficient. The approach is demonstrated on noise-free numerical experiments, and reproduced the true model parameters. The accuracy and stability of the proposed approach are examined on numerical experiments contaminated with realistic noise levels and regional fields prior to the interpretation of real data. Following that, five real field examples from geothermal systems and mineral exploration have been successfully analyzed. The results agree well with the published research.

Figure 1.Geometry and parameters of the assumed source models.Top, middle and bottom panels present the sphere, semi-infinite vertical cylinder, and infinitely long horizontal cylinder models.models was not developed before.It is relevant to note that rigorous 3D inversion of self-potential data is computationally expensive, and requires a priori information for the model parameters (3D electric conductivity distribution) we invert for 17 .
The paper is structurally described as follows.The "Self-potential direct solution" section presents the direct problem (forward modeling solution).In "The method" section, the formulation of the proposed imaging scheme is explained.The "Numerical examples" section validates the method using synthetic models contaminated with a wide range of noise, regional self-potential signatures, and interference anomalies.Real data examples are carefully analyzed and discussed in the "Field examples" section, and finally some findings are reported.

Self-potential direct solution
The self-potential signature (V) of some simple geometrical sources at an observation point ( x j , z) along profile (Fig. 1) is given by, e.g., Yungul 53 and Mehanee 66 : where x j (m) is the coordinate of the measurement station (Fig. 1), j is the index of the measurement station, x • (m) is the origin point of the self-potential profile, z • and z (m) are the coordinates of the buried body and the observation station, θ (degrees) is the polarization angle, q (dimensionless) is the shape factor (q = 1.5, 1 and 0.5 for sphere, horizontal cylinder and semi-infinite vertical cylinder), n is the number of data points and K is the electric dipole moment.It is noted that the unit of K (mV m 2q−1 ) is function of the shape factor (q) 62 , and that θ is measured clockwise (Fig. 1) and ranges from 0 to − 180 • in the above formula 8 .

The method
The AS expression 70 is: where ∂V ∂z and ∂V ∂x j are the vertical and horizontal derivatives of the self-potential anomaly.The analytic signal's amplitude ( |A s (x j , z)| ) of the self-potential anomaly is given by Nabighian 70 : Taking the vertical and horizontal derivatives of formula (1), and by substituting the results in expression (3), we get: The R-parameter (correlation coefficient) is dependent upon both of the amplitude analytic signal of the actual self-potential data (Aso) observed along profile and that of the analytic signal amplitude of the calculated (theoretical) self-potential data (Ast) produced by an assumed source (for example a sphere): The analytic signal [Aso] can be assessed numerically using expression (3), whereas the analytic signal [Ast] of an assumed source is calculated analytically from expression (4).
To calculate the parameters of a presumed source, the imaging parameter (R-parameter) is first plotted on a 2D map from which the depth z o can be readily read.It is worthy noting that when the parameters of the SP profile (Fig. 1) coincide with the buried anomalous source, the imaging parameter achieves its maximum value (so-called here R-max).Further pertinent detail is provided in the "Numerical examples" section.As is seen from expression (5), the computation of R(x o , z o , θ ) does not need knowledge of the electric dipole moment K, which is calculated from the maximum self-potential response.Using expression (1), the predicted self-potential response A s (x j ,z) =|K| Table 1.Model 1: horizontal cylinder.Shape index (q) and its corresponding maximum value of the R-parameter.

Numerical examples
The approach proposed here has been examined on synthetic self-potential data generated by various source models (for example, horizontal cylinder, sphere and vertical cylinder).The suggested scheme is first verified on numerical experiments without noise.After that, the data have been contaminated with realistic noise, and interpreted in order to assess the stability of the scheme e.g. 8,84.Second, in order to assess the stability of the scheme further, the effect of the regional background (that is embedded into the measured self-potential data) on the results is carefully investigated.
Model 1: horizontal cylinder.The self-potential response (Fig. 3a) of an idealized body of a horizontal www.nature.com/scientificreports/derivatives (horizontal and vertical) of the self-potential anomaly (Fig. 3a).The corresponding analytic signal amplitude (Fig. 3c) is then calculated from the spatial derivatives (Fig. 3b) using expression (3).The mosaic surface S (which was gridded into 1-m spaces in the x-and z-directions) extended, respectively, to 120 × 12 m in these directions (that is ( x o , z o ) ∈ S = (0, 120) × (1, 12)), and was used to compute and map the correlation coefficient (R-parameter).Expression ( 5) is employed to calculate the R-parameter for each possible sources (q = 0.5-1.5),where the largest value (R-max) of the R-parameter is attained at the true assumed source (that is q = 1 and R-max = 1.0) (Table 1).Figure 3d presents the image of the R-parameter composed using expression (5) assuming source of a horizontal cylinder model.The R-parameter's maximum value is marked by the black dot, which denotes the true model parameters of the buried structure (Fig. 3d).
We reiterate that the R map (Fig. 3d) shows the 2D distribution of the obtained R-parameter values.The R-parameter measures the goodness of fit between the observed and predicted self-potential data, and is not representative of geologic structures.An R value of 1 means that the observed and predicted self-potential data are in perfect fit.
To further assess the developed imaging scheme, a number of shape values have been investigated.The scheme is found stable and can retrieve the true values of the model parameters as can be seen from the results presented in Fig. 4 and Table 1.
We added 20% random noise into the self-potential data (Fig. 3a) of the above-mentioned synthetic model (Fig. 5a). Figure 5b,c show the spatial derivatives of the noisy self-potential anomalous signature (Fig. 5a) and their corresponding AS amplitude.Figure 5d shows the maximum R-parameter value (black dot) with a magnitude of 0.98.The estimated model parameters (K= 2906.30mV m, z o = 7.4 m, θ = −42.60• , and x o = 60 m for an assumed shape factor q of 1.0) (Fig. 5d) are in good agreement with the actual ones.
It can be concluded from the above analysis that the R-parameter imaging method can produce accurate model parameters when the self-potential data are contaminated with noise.
Model 2: sphere model.A synthetic self-potential anomaly for a sphere-shaped body (K = 1000 mV m 2 , z o = 4 m, θ = −25 • , x o = 60 m and profile length = 120 m) combined with a first-order regional anomaly generates the composite self-potential anomaly that is shown in Fig. 6a.The simulation formula of the composite anomaly has the form:  To assess the accuracy of the developed imaging scheme, the composite SP data profile (Fig. 6a) has been contaminated with 20% random noise (Fig. 7a). Figure 7b,c

Field examples
The scheme is analyzed on five published real self-potential data from geothermal systems and mineral exploration in the following sections.The Hi'iaka self-potential anomaly, the Kilauea volcano, Hawaii, the United States of America.
A number of self-potential geophysical surveys were carried out in 1973, 1995, 1997 and 2012 over a basaltic dike intrusion (referred to as the Hi'iaka dike, Hawaii) (Fig. 8).The Hi'iaka dike intruded into the upper part of the Kilauea volcano, which is associated with the Hi'iaka and Pauahi craters eruption alongside the Kilauea rift zone 87,88 .It made a 100-m long surface fracture that erupted magma southwest of the Hi'iaka crater.www.nature.com/scientificreports/Measurements of geophysical surveying recommended that the fracture is continued about 1.5 km under the subsurface in the southwest direction (Fig. 8).
The measurements of the self-potential data over the Hi'iaka dike intrusion commenced in 1973 by Zablocki 89 and Zablocki 90 , and continued in 1995, 1997 and 2012 83 .Localized fluid disruption is blamed for the SP anomaly 83,89,90 .Davis 83 stated that "geothermal reservoirs are found above intrusions of magma such as dikes or dike swarms, which set up hydrothermal circulation generating hot water and steam from which energy can be tapped".The self-potential anomaly profiles of 1973, 1995, 1997, and 2012 are digitized into 10-m intervals (Fig. 9(1)-( 4)).The Hi'iaka SP anomaly profiles have been interpreted by Davis 83 using the self-potential inversion approach of Sill 38 .Davis 83 interpreted the profiles by a trapezoidal source (approximated by a dike shaped model) located at different depths ranging from 50 to 190 m, and attributed the increase in depth to the magma cooling and heat loss at the top of the dike.
The aforementioned Hi'iaka self-potential anomaly profiles of 1973, 1995, 1997 and 2012 are interpreted using the scheme developed here (Fig. 9(1)-( 4)).The derivatives of the self-potential response of each profile (Fig. 9(1)a,(2)a,(3)a,(4)a), and the corresponding analytic signal amplitude are rendered in Fig. 9(1)b,(2)b,(3) b,(4)b, and in Fig. 9(1)c,(2)c,(3)c,(4)c.The R-parameter values are reported in Fig. 9 2 tabulates the recovered model parameters for each profile, and shows that the observed self-potential anomaly is fit by a dike-like model with a shape factor of 0.6-0.8(that is q = 0.6-0.8).Analysis shows that there is good match between the depths of the interpretive source (trapezoidal model) stated in the published literature and the depths obtained here (Table 3).The match between the observed and calculated self-potential data for each profile is depicted in Fig. 9(1)a-(4)a, which is quite good.
It is re-noted that Davis 83 reported that the SP anomaly remained strong throughout the measurement duration.However, the SP anomaly of 2012 is about 60% of that of 1973.Therefore, the variation in the recovered depths (54-177 m, Table 3) of the interpretated self-potential profiles (measured in 1973-2012) is not unexpected, and is attribuated to the magma cooling and heat loss at the top of the dike 83 .
The Osnabrück self-potential anomaly, Germany.A self-potential anomaly near the Osnabrü ck area (Fig. 10), Northwest Germany 6 has been carried out to trace a graphite anomaly that has a quasi-vertical form in the Lias-epsilon shales.Gurk et al. 6 found a significant single self-potential anomaly of roughly − 600 mV (Fig. 11a), which supports conductive graphite minerals.The 500-m long self-potential anomaly profile is meshed into 5-m intervals (Fig. 11a).4).According to the recovered R-parameters, the subsurface anomalous body is approximated by a horizontal cylinder like-structure with a horizontal location of 250 m and an estimated depth to the center of 23 m, which is in good agreement with the interpreted results of Gurk et al. 6 and Mehanee 8 (Table 5).The variation in the magnitude of the model parameter K is due to the inconsistent use of unit (Table 5) as the interpretive models are not quite identical; they range from thin sheet to quasi-horizontal cylinder.
In order to map the 2D electric conductivity (inverse of resistivity) distribution in the underground, Gurk et al. 6 measured a radio magnetotelluric data (apparent resistivities and phase) profile on the initial 400 m of the self-potential profile described earlier.The corresponding 2D inverse results of Gurk et al. 6 revealed a prominent conductive anomalous body (Fig. 12), the location and depth of which correlate well with the results inferred from the approach developed here (Fig. 11).
The Suleymankoy self-potential anomaly, Turkey.The SP anomaly of Suleymankoy 53 was carried out for copper deposits.The mine is characterized by alpine ophiolite containing several copper deposits.The anomaly is gridded at intervals of 2 m long (Fig. 13a).The self-potential anomaly of Suleymankoy is interpreted using the presented R-parameter imaging technique.Figure 13b-d depict the corresponding derivatives, the AS amplitude, and the imaging, which reveals an R-max value of 0.9985.The parameters revealed from interpretation are K = 1898 mV m 2q−1 , z o = 27 m, θ = −130 • , x o = 76 m and q = 0.8 (Fig. 13a-d and Table 6).The observed and calculated self-potential data are rendered in Fig. 13a.Table 7 presents a comparison between the obtained results and those mentioned in the published literature.
As can be seen from Table 7, the reported depths are in reasonable agreement, whereas the electric dipole moments (K) encountered some variations, which could be attributed to the various approximations employed in the interpretation schemes used in this table and how the parameter K is calculated from these schemes.
The Malachite mine self-potential anomaly, The United States of America.The Malachite Mine is an amphibolite belt bounded by schists and gneisses.The self-potential profile over the Malachite Mine is digitized at intervals equal to 2 m (Fig. 14a).
The Malachite SP anomaly is interpreted using the R-parameter imaging technique.Figure 14b,c render the gradients, and the AS amplitude.Upon applying the imaging technique, an R-max of 0.9841 (Fig. 14d) was obtained with the following interpretive parameters: K = 515 mV, z o = 15 m, θ = −112 • , x o = 88 m, and q = 0.67 (Fig. 14a,d and Table 8).The subsurface structure was approximated by a semi-infinite vertical cylindrical structure with a horizontal spatial location of 88 m.The estimated depth (15 m to the top of the structure) is in good match with the drilling information and previous interpreted works (Table 9).Table 9 shows that the parameter K encountered some variation; this is attributed to two main reasons.First, the inconsistent use of unit as the interpretive models are not quite identical; they range from vertical cylinder to quasi-vertical cylinder.Second, the nature of the approximations employed in interpretation schemes reported in this table, and how the parameter K is calculated from these schemes.
The Bavarian woods self-potential anomaly, Germany.Figure 15a depicts the self-potential anomaly collected over a graphite ore body at the southern Bavarian woods, Germany 95 .The self-potential anomaly profile is digitized using a 1-m sampling interval.Several authors have interpreted this anomaly profile.Al-Garani 96 interpreted the anomaly as a quasi-vertical cylinder using the neural network inversion with z o = 33 m (depth to the top).Mehanee 66 analyzed this SP profile with a horizontal cylinder of a depth to the center z o of 46 m using a regularized inversion.Gokturkler and Balkaya 5 described the anomaly by a horizontal cylinder model using a genetic algorithm ( z o = 45.03 m; to the center), simulated annealing ( z o = 47.59 m; to the center) and particle swarm optimization algorithm ( z o = 47.59 m).Di Maio et al. 97 fit the profile by a horizontal cylinder by applying a spectral analysis and tomographic approach with z o = 44.9m (to the center).
We interpreted this self-potential anomaly profile using the R-parameter imaging technique.The analytic signal amplitude anomaly calculated from the horizontal and vertical derivatives (Fig. 15b) is presented in Fig- ure 15c.The R-parameter values are rendered in Fig. 15d.An R-max of 0.9983 pertains to an interpretive model of K = 13821 mV m 2q−1 , z o = 41 m, θ = −62 • , x o = 263 m, and q = 0.94 (Fig. 15a,d and Table 10).The presented analysis shows that the depth and shape (which resembles a quasi-horizontal cylinder model) are generally in good agreement with the aforementioned results but the results of Al-Garani 96 , who interpreted the data by a quasi-vertical cylinder (q = 0.7) (Table 11).We do not expect all interpretation methods to yield the same results as each method has its own assumptions and limitation.

Discussion
As mentioned above, the R map shows the 2D distribution of the obtained R-parameter values.The R-parameter measures the goodness of fit between the observed self-potential data, and the theoretical self-potential data generated from the model parameters (z, K, θ ) of the interpretive idealized model.The R parameter value does not provide uncertainty estimate for the evolved model parameters.
The non-uniqueness is one of the most challenging issues in geophysical data interpretation e.g. 99, where multiple approximative solutions can equally fit the observed data.Joint inversion could help minimize this issue and provide better understanding e.g. 84.It is worthy noting that it is very rare to solely use/measure one data kind when it comes to a detailed geophysical prospecting program.In industry, multiple data sets are essential for comprehensive understanding and for maximizing the potential of the underlying exploration program.Usually multiple geophysical data along with geological information are used, inverted and interpreted in an integrated manner (the so-called joint interpretation) to hopefully recover and select a unique inverse model the data of which match the measured data sets, and that fits into the underlying geologic setting of the area under study.May be this is the best we can do in exploration geophysics in order to resolve the non-uniqueness issue of an inverse problem solution.

Conclusions
A rapid imaging scheme has been developed for the interpretation of self-potential data.In about 2 s on a simple PC, the scheme can estimate the parameters of the interpretive model (which is in the context of sphere, horizontal cylinder or vertical cylinder) that resembles the buried structure.The developed scheme uses the amplitude of the AS of the self-potential data undergoing interpretation and the amplitude of the AS of the self-potential data calculated by the assumed interpretive model to construct the corresponding 2D image of the so-called R-parameter.The scheme attains the largest value of the R-parameter when the recovered parameters coincide with the actual ones.It is noted that the R-parameter is independent of the electric dipole moment (K).The analyzed numerical examples demonstrated the stability of the developed scheme, and that its accuracy can be affected by the nearby geological structures.The five field data examples (from geothermal systems and mineral prospecting) analyzed here show that the scheme is capable of producing good results that agree well with those reported in other published research.The developed imaging scheme can have some potential in geothermal investigation and reconnaissance studies.Table 7.The Suleymankoy self-potential anomaly data, Ergani, Turkey.Estimated parameters.

Figure 2 .
Figure 2. Flowchart showing the workflow of the developed scheme.

Figure 3 .
Figure 3. Model 1: SP data (noise free).(a) SP anomaly produced from a horizontal cylinder model.(b) Horizontal (HD) and vertical (VD) derivatives of the SP data rendered in (a).(c) Analytic signal (AS) amplitude.(d) Image of the R-parameter (R).

Figure 4 .
Figure 4. Model 1: noise-free data.Relationship of the R-parameter, shape factor, and depth.

Figure 8 .
Figure 8.The Hi'iaka self-potential anomaly, the Kilauea volcano, Hawaii, USA.Profile of the self-potential measurements (dashed line).Inferred location of the Hi'iaka dike (solid line) (taken from Davis 83 with permission from Elsevier).

Figure 11b and c
Figure 11b and c present the derivatives, and the AS amplitude of the SP anomaly, respectively.The R-parameter maximum value (R-max = 0.9976, black dot, Fig. 11d) was determined with the corresponding best interpretive parameters: K = 46527 mV m 2q−1 , z o = 23 m, θ = −97 • , x o = 250 m, and q = 1.2 (Fig. 11a,d and Table4).According to the recovered R-parameters, the subsurface anomalous body is approximated by a horizontal cylinder like-structure with a horizontal location of 250 m and an estimated depth to the center of 23 m, which is in good agreement with the interpreted results of Gurk et al.6 and Mehanee 8 (Table5).The variation in the magnitude of the model parameter K is due to the inconsistent use of unit (Table5) as the interpretive models are not quite identical; they range from thin sheet to quasi-horizontal cylinder.In order to map the 2D electric conductivity (inverse of resistivity) distribution in the underground, Gurk et al.6 measured a radio magnetotelluric data (apparent resistivities and phase) profile on the initial 400 m of the self-potential profile described earlier.The corresponding 2D inverse results of Gurk et al.6 revealed a prominent conductive anomalous body (Fig.12), the location and depth of which correlate well with the results inferred from the approach developed here (Fig.11).

Figure 10 .Figure 11 .
Figure 10.The Osnabrü ck self-potential anomaly, Germany.Location of the survey area (star) north of Osnbrück, and isolines of Vitrinite Reflectance of the maturity map (taken from Gurk et al. 6 with permission from Elsevier).NL: Netherlands, B: Belgium.
is then calculated.Fig.2presents a flowchart showing the workflow of the developed scheme, which takes about 2 s on a simple PC to estimate the parameters of the interpretive model that resembles the buried anomaly.

Table 2 .
The Hi'iaka self-potential anomaly, the Kilauea volcano, Hawaii, USA.R-max versus q.See the text for details.

Table 8 .
The Malachite Mine self-potential anomaly, Colorado, USA.R-max versus q.

Table 9 .
The Malachite mine self-potential anomaly, Colorado, USA.Estimated parameters.

Table 10 .
The Bavarian woods self-potential anomaly, Germany.R-max versus q.

Table 11 .
The Bavarian woods self-potential anomaly, Germany.Estimated parameters.