The Self-Potential Anomaly Produced by a Subsurface Flow at the Contact of Two Horizontal Layers and Its Quantitative Interpretation

In the present paper the problem of a polarized cylinder with a small cross-section, which is located at the contact of two horizontal layers with different resistivities, is studied. Such a polarization geometry simulates the self-potential (SP) field produced by a horizontal flow at the contact between the two layers. First, the expression of the self potential at the space domain is derived, applying the image technique. Then, the expression for the Fourier transform of the SP anomaly is found and the behavior of the amplitude spectrum is studied. Based on this study, a direct interpretation method at the spatial frequency domain is proposed, in order to calculate the depth of the flow and the reflection coefficient of the stratified medium. Experimentation with a synthetic model shows that the method works well (small deviations between true and calculated values). When the SP curve contains noise, deviations between calculated and true depths are smaller than those between calculated and true reflection coefficients. The proposed method, which is also applied on SP data from a geothermal system (Mauri et al., 2010), may be useful in detecting underground water or heat flows.


Introduction
The self-potential (SP) method has been extensively used to detect groundwater and hydrothermal flows.During the last decade, a variety of applications of the SP method have been developed, such as detecting ground water and hydrothermal flows at volcanic areas [1][2][3][4], mapping infiltration in a karstic environment [5], monitoring evolving plumes at subsoil [6] and studying landslides [7].Recent experimental work has also been done to study the electrokinetic and hydraulic parameters, which influence the SP anomaly [8][9][10].A thorough review on the applications of the self-potential method in the detection and monitoring of subsurface flows is made by Jouniaux et al. [11].
According to Schiavone and Quarto [12], a horizontal groundwater flow at the contact of two layers is expected to produce a self-potential anomaly of electrokinetic origin, at ground surface.The geometry of the subsurface flow in the stratified medium is presented in Figure 1.The two layers have different electrokinetic coupling coefficients C 1 (upper layer) and C 2 (lower layer) and electrical resistivities ρ 1 and ρ 2 , respectively.In Figure 2 a profile of the SP anomaly at ground surface is presented.The direction of the profile is perpendicular to that of the subsurface flow.The SP anomaly is expected to be positive for (C 1 − C 2 ) P > 0 and negative for (C 1 − C 2 ) P < 0. P is the pressure difference between the contact of the layers, where the flow takes place, and the ambient.Electric potential drops are expected to be produced along the direction of the water flow, but they do not affect the polarity of the SP anomaly across the water flow.The polarity could be affected by the presence of a vertical component of the flow, which does not exist in the model under consideration.
Considerable research has been made to model the electrokinetic self-potential anomaly in terms of its production mechanism (water flow in a porous medium) and its International Journal of Geophysics Representation of a self-potential anomaly that is produced by a horizontal subsurface flow at the contact of two layers (Schiavone and Quarto [12]).dependence on the resistivity distribution of the ground [13][14][15].On the other hand, in order to invert the measured SP anomaly, which means to calculate the depth of the horizontal flow, it is important to have in mind that the physical system of the subsurface flow at the contact of the two layers is equivalent to a system of dipoles oriented perpendicularly to the contact, as Ogilvy et al. [16] have pointed out.Therefore, the inverse problem consists in determining the location of the underground dipole current sources.

SP (mV)
If the width of the underground current flow is small relatively to its depth from ground surface, the source of the SP anomaly may be represented by a system of two elongated horizontal lines with an opposite charge.This elongated dipole is usually called polarized cylinder of a small radius R and depth b (see Figure 3).A first approach to the inverse problem is to assume an electrically homogeneous ground and calculate the depth of the cylinder by one of the several direct quantitative interpretation methods which have been developed, such as those of Bhattacharya and Roy [17] or Murty and Haricharan [18].A serious drawback of these approaches is that if the electrical inhomogeneity of the ground is not taken into account, considerable errors in the depth calculation of the polarized body may be made [19,20].
An alternative approach on the inversion problem is to fit the SP curve measured at ground surface by an automated trial and error process, assuming different distributions of electric current sources at a nonhomogeneous ground for each iteration, till a satisfactory fit is achieved [21,22].Such a process may work well if the resistivity distribution of the ground is known a priori.
The self-potential tomography, which has been first developed by Patella [23,24], solves the inverse problem without any a priori knowledge of the resistivity distribution, using correlation integrals and defining probability of tomography.However, Gibert and Sailhac [25] demonstrated that the so-called probability of tomography defines images of SP data in the wavelet domain that must not be interpreted as underground images of SP sources.They pointed out that an appropriate inversion is necessary to achieve underground image.An alternative methodology of self-potential tomography, based on the wavelet transform, has been developed by Gibert and Pessel [26].
In the present paper, a direct interpretation method of SP anomalies is proposed in order to calculate the depth of the horizontal subsurface flow at the contact between two horizontal layers, as well as the resistivity contrast of the ground, which is expressed as reflection coefficient.The flow is simulated by a horizontal polarized cylinder.The quantitative interpretation is carried out at the spatial frequency domain of the self-potential field.The potential advantage of the proposed method, which may be useful in detecting horizontal groundwater or geothermal flows at a stratified ground, is that although it takes into account the ground electrical inhomogeneity, it works with no a priori information about the resistivity distribution.

The SP Anomaly at Space Domain
In Figure 3 the geometry of the problem of the horizontal polarized cylinder at the contact of two horizontal layers is presented.The polarized cylinder is the equivalent of a horizontal subsurface flow at the same location, which produces a self-potential anomaly V(x) at ground surface.The intersection of the cylinder (see solid circle in Figure 3) has a small radius R and its center is at depth b below ground surface.b is also the thickness of the upper layer.
For a vertically polarized horizontal cylinder at depth b in an electrically homogeneous medium of resistivity ρ, the SP anomaly V(x) is given by [17,18] M is the polarization of the body and it is given by [17,18] International Journal of Geophysics 3 I is the current intensity per unit length, at the poles of the cylinder.It can be positive or negative, according to the polarity of the SP anomaly (positive or negative, resp.).
If the resistivity ρ 1 of the upper layer is different from the resistivity ρ 2 of the bedrock and the cylinder is located at depth h greater than b (see the lower dashed circle in Figure 3), then the expression for the self-potential anomaly V 2 (x) can be found by the image technique [27], which gives Polarization M 2 is given by k is the reflection coefficient, which is defined by If the cylinder is located at depth h smaller than b (see the upper dashed circle in Figure 3), then the expression for the self-potential anomaly V 1 (x) is given by Polarization M 1 is given by The SP anomaly V(x), which is produced by a polarized cylinder at the contact between the two layers, may be seen as the limit of the potentials V 1 and V 2 for h → b − and h → b + , respectively.Taking into account the relations (3), ( 4), ( 6) and ( 7), it can be easily proved that the limits of V 1 and V 2 are equal and V(x) is The self-potential field produced by a horizontal flow at the contact of the two layers is described by relation (8).Taking the Fourier transform of this expression, it is possible to study the SP anomaly at spatial frequency domain.

Derivation of the Fourier Transform of the SP Anomaly
The Fourier transform U(u) of the SP anomaly V(x) is defined by u is the spatial frequency.Every u is related to a certain x value by u = 2π/x.The following tabulated integral [28] may be used in order to find the Fourier transform of Based on relation (10), it can be easily shown that (11) Combining ( 8), (9), and (11) gives Summing the terms of the series of relation (12), the following expression for U(u) is obtained: K is given by The amplitude spectrum A(u) of U is defined as the absolute value of U, so Using tabulated integrals, it can be also proved that the function A(u) is continuous at u = 0.In Figure 4 the amplitude spectrum of the SP anomaly is presented.Studying the behavior of A(u), certain simple mathematical expressions may be derived, in order to develop a method to calculate the quantities b and k.

Quantitative Interpretation of the SP Anomaly
In Figure 5, the variation of the natural logarithm of amplitude A against u is presented.It can be seen that as long as u increases, the quantity ln A tends to present a linear variation with u.
The mathematical expression for this linear variation can be derived by relation (15) if u is big enough so that k• exp(−2bu) may be neglected, meanwhile exp(−bu) remains more than zero.In such a case, the natural logarithm of A(u) becomes ln A(u) ≈ ln|K| − bu.(16) This means that, for relatively high u values, for which k• exp(−2bu) ≈ 0, the absolute value of the slope of the straight line ln A(u) is equal to the depth b of the cylinder.The minimum value u for which relation ( 16) is valid depends on b and k.In practice, the interpreter has to define the linear part of ln A against u, by which b may be calculated.Now let u t be the u value for which A(u t ) is equal to one tenth of A at null spatial frequency.Such a condition is expressed by According to relation (15), A(0) is given by Combining relations (15), (17), and ( 18) and solving for k give The quantity exp(−2bu t ) is much smaller than 10•exp(−bu t ) (ten times smaller for bu t = 0 and even more as long as bu t increases).Neglecting exp(−2bu t ), relation (19) takes the following simplified form: Therefore, determining the u t value for which the amplitude A is equal to A(0)/10, it is possible to calculate the reflection coefficient k, by relation (20).
According to the mathematical study of the behavior of the amplitude spectrum, a direct quantitative interpretation method is proposed, which may be executed by the following steps.
(1) Calculate the amplitude spectrum A(u) of the SP anomaly measured in the field.
(2) Take the natural logarithm ln A of A and calculate the slope of the linear part of ln A against u.The absolute value of the slope is equal to the depth b.
(3) Find the u t value for which the amplitude A(u t ) is equal to A(0)/10.The reflection coefficient k between the two layers can be calculated by relation (19) or (20).
With this interpretation method it is possible to calculate the depth of the horizontal subsurface flow and the reflection coefficient k between the two layers.

Experimentation with a Synthetic Model
In order to test the efficiency of the proposed quantitative interpretation method, a synthetic model was used with 4Iρ 1 ρ 2 R/[π(ρ 1 + ρ 2 )] = 5000 mV.m, b = 50 m, and k = 0.6.The SP anomaly V(x) for 256 values of x and a 10 m interval was calculated.The Fourier analysis of V(x) was done, applying the well-known FFT algorithm.
In Figure 6, the amplitude spectrum A(u) of the selfpotential anomaly is presented, for u between 0 and 0.07 rad/m.The value u t for which A becomes equal to A(0)/10 may also be noted.
In Figure 7 one can see how ln A varies with u.The range of u is from zero to the Nyquist frequency.For very low and for very high values of u the variation is not linear.The nonlinear behavior for low spatial frequencies is in agreement with the theoretical predictions, as it can be seen in Figure 5.The scattering of ln A values for very high spatial frequencies is due to numerical errors, which may seriously affect the calculations when A(u) is very small.In a broad part of the spectrum, ln A varies linearly with u.
In Figure 8, the linear part of the spectrum is presented.The least squares straight line was calculated, and the slope of the line was found to be equal to −50.02.Thus, the calculated b is equal to 50.02 m, which is very close to the true value of 50.00 m (the relative error is 0.04%).
As it can be seen in Figure 6, u t = 0.030 rad/m.Putting b = 50.02and u t = 0.030 at the right part of relation (20), a value of k equal to 0.55 is calculated.Thus, there is an 8% relative error between the calculated and the true value of k, which is equal to 0.600.If k is calculated by (19), which is an accurate expression, it takes the value k = 0.56 (relative error equal to 7%).Therefore, the approximate formula (20) gives a k value which is very close to that of the exact expression (19).
The experimentation with the synthetic model produced encouraging results as long as the reliability of the proposed quantitative interpretation method is concerned, since the calculated b and k values are very close to the true ones.

Application on Real-Field Data
The proposed method was tested on an SP profile along Waita volcano (Japan), where hydrothermal activity is developed (Mauri et al. [4]).The self-potential anomaly is presented in Figure 9. From the mathematical point of view the SP anomaly may be simulated by a polarized cylinder with an upper negative pole and a lower positive pole.The two negative peaks at the centre of the anomaly indicate the possible existence of at least two polarized bodies, but according to the proposed method only one body may be detected.The horizontal projection of the cylinder is estimated to be located at the centre of the anomaly, which means at distance x between −500 m and −300 m (mean value of x = −400 m).
The SP anomaly was manually digitized and 128 values were stored (sampling interval equal to 31.25 m).These values were analyzed by an FFT.In Figure 10 the amplitude spectrum A(u) of the anomaly is presented.
It can be observed that because of the "wavy" behaviour of the curve A(u), the value u t cannot be exactly defined and it is estimated to be between a little more than 0.0045 rad/m and a little less than 0.0090 rad/m.u t takes the value of the midpoint of this interval, which is u t = 0.0068 rad/m.The International Journal of Geophysics error σ(u t ) can be taken equal to half of the uncertainty range between 0.005 and 0.009, which means σ(u t ) = 0.002 rad/m.
In Figure 11, the linear part of ln A against u is presented.The inclination of the least squares straight line is equal to −194; therefore the depth of the SP source, which is also the depth of the horizontal contact between the two layers, is b = 194 m.The standard error in the calculation of the inclination of the least squares line was found to be σ(b) = 27.
The reflection coefficient k may be calculated by relation (20), putting b = 194 and u t = 0.0068, which gives k = 0.63.If k is calculated according to relation (19), then k = 0.64.The error in k may be calculated by the square root of the sum of the squares of (∂k/∂u t )•σ(u t ) and (∂k/∂b)•σ(b), for b = 194, u t = 0.0068, σ(u t ) = 0.002, and σ(b) = 27.This error was found to be σ(k) = 0.16.
The results of the quantitative interpretation may be summarized as follows According to relation (5), a reflectance coefficient equal to 0.6 gives a resistivity ratio equal to 4/1 (lower to upper layer).Mauri et al. [4] interpreted the SP anomaly by a tomography methodology using the wavelet transform.They detected five SP sources all along the self-potential profile, at depths varying between 30 m and 400 m.More precisely, they found sources at x = −1850 m, −1230 m, −400 m, +365 m and +1220 m, with respective depths b = 295 m, 150 m, 400 m, 100 m, and 30 m.The methodology which is proposed in the present paper cannot detect more than one source and the calculated depth has an intermediate value between the five source depths.
On the other hand, Mauri et al. [4] do not report anything about the resistivity distribution of the ground.It would be interesting to compare our prediction about a horizontal layer boundary at 190 m depth and a 4/1 resistivity ratio with data obtained by vertical electric sounding or electrical tomography, if there is any.

A Study of the Sensitivity of the Method
The efficiency of the proposed interpretation method is controlled by the ground resistivity distribution and the noise.
The role of ground resistivity may be understood by taking into account relation (8).High values of ρ 1 and ρ 2 produce strong SP anomalies, which help to make a reliable calculation of depth and reflection coefficient.On the other hand, a very low value of ρ 1 or ρ 2 may produce a very weak SP anomaly, which cannot be reliably interpreted.This may happen, for example, if the upper layer is a clay formation.Since clay has an electronic conductivity, its resistivity is very low and this does not favor the development of a strong SP anomaly.
In order to study the sensitivity of the quantitative interpretation method in the presence of noise, the synthetic model of paragraph 5 was used (model parameters: 4Iρ 1 ρ 2 R/[π(ρ 1 + ρ 2 )] = 5000 mV•m, b = 50 m, and k = 0.6).Random noise was introduced, with a zero mean value and amplitudes between 0% and 40% of the amplitude of the SP anomaly of the synthetic model.Each SP curve, which contained the signal of the polarized body and random noise, was interpreted according to the proposed method.The absolute relative error values |Δb/b| and |Δk/k| were calculated according to the relations: b c and k c are the calculated values of depth and reflection coefficient, respectively.Figure 12 presents the relative absolute error |Δb/b| against the noise percentage (ratio of noise amplitude to signal amplitude).It can be observed that up to a 5% noise, the error is null and keeps being small (less than 20%), when the noise percentage is not more than 10%.For higher noise percentages, |Δb/b| increases and exceeds 50% for a noise percentage more than 30%.
Figure 13 presents the relative absolute error |Δk/k| against the noise percentage.Compared to |Δb/b|, the error in the calculation of the reflection coefficient is considerably International Journal of Geophysics bigger.For a 5% noise content, it is 15% and for a 30% noise percentage or more the deviations between real and calculated k are so high that it is meaningless to calculate k.
In Figure 12, it can be observed that there is a decrease in the error |Δb/b| between 20 and 30% noise.Actually, |Δb/b| depends on the exact values of the random noise and, more important, on the exact points of the amplitude spectrum which are considered to follow a straight line.These points were chosen by a visual criterion, according to which the beginning and the end of the line was defined.Such a criterion is inevitably subjective and may produce deviations in the ascending tendency of the error with the noise content.The decrease of |Δk/k| in between 15 and 20% noise (see Figure 13) may be partly attributed to the deviation of b and partly to the "wavy" behavior of the curve A(u), because of which there is an uncertainty in the value of u t .

Conclusions
Self-potential measurements at ground surface are often considered as a preliminary exploration method, which has to be followed by other techniques, such as geoelectrical sounding or electrical tomography, in order to form a more precise picture of the subsurface water flow and the geological formations in which this flow occurs.
On the other hand, it is well known that if ground electrical inhomogeneities are not taken into account, considerable errors in the calculation of the parameters of the SP source may be made.
Applying the quantitative interpretation method which is proposed in the present paper, the depth of the SP source (which is a horizontal subsurface flow) may be reliably calculated since the electrical inhomogeneity is not ignored, although it is not known a priori.Furthermore, it is possible to estimate the ground resistivity contrast by calculating the reflection coefficient, according to relation (19) or to relation (20).The latter is an approximate formula to calculate k; however, it gives a value which is very close to that of the exact expression (19).Therefore, the SP method gains reliability and provides information about the geological stratification of the ground, which may be taken into account during the application of other exploration techniques.
Experimentation with a synthetic model provided encouraging results long as the reliability of the proposed quantitative interpretation method is concerned.Good estimates of depth b may be done, for an up to 10% noise percentage.The calculation of the reflection coefficient k is considerably more sensitive to noise than the one of b.If the noise content is not more than 5%, the relative error in the calculation of k is less than 20%.
The ground resistivity distribution controls the amplitude of the self-potential anomaly.Weak SP anomalies are expected to be produced over a layer with a very low resistivity, such as a clay formation.
The proposed method was applied to field data obtained at a hydrothermal field with more than one SP sources.The method could calculate only one depth value, which was found to be between the depth values of the sources.Further application on other self-potential anomalies produced by only one source is necessary in order to assess the efficiency of the method.This could be the subject of a future research.

C 1 C 2 Figure 2 :
Figure2: A profile of the self-potential anomaly produced by a horizontal flow perpendicular to the paper.The location of the flow is indicated by the symbols XXX (Schiavone and Quarto[12]).

Figure 3 :
Figure 3: Representation of the SP anomaly produced by a polarized cylinder at the contact of two horizontal layers with different resistivities.

Figure 6 :Figure 7 :
Figure 6: The amplitude spectrum of the SP curve of the synthetic model.

Figure 8 :Figure 9 :
Figure 8: The linear part of the variation of ln A against u, fitted by the least squares line.

Figure 10 :Figure 11 :
Figure 10: The amplitude spectrum of the SP anomaly.

Figure 12 :Figure 13 :
Figure 12: Absolute relative error in the calculation of depth against noise percentage.