Tracing of paleo-shear zones by self-potential data inversion: case studies from the KTB, Rittsteig, and Grossensees graphite-bearing fault planes

Abstract

This paper describes a new method for tracing paleo-shear zones of the continental crust by self-potential (SP) data inversion. The method falls within the deterministic inversion framework, and it is exclusively applicable for the interpretation of the SP anomalies measured along a profile over sheet-type structures such as conductive thin films of interconnected graphite precipitations formed on shear planes. The inverse method fits a residual SP anomaly by a single thin sheet and recovers the characteristic parameters (depth to the top h, extension in depth a, amplitude coefficient k, and amount and direction of dip θ) of the sheet. This method minimizes an objective functional in the space of the logarithmed and non-logarithmed model parameters (log(h), log(a), log(k), and θ) successively by the steepest descent (SD) and Gauss-Newton (GN) techniques in order to essentially maintain the stability and convergence of this inverse method. Prior to applying the method to real data, its accuracy, convergence, and stability are successfully verified on numerical examples with and without noise. The method is then applied to SP profiles from the German Continental Deep Drilling Program (Kontinentales Tiefbohrprogramm der Bundesrepublik Deutschla - KTB), Rittsteig, and Grossensees sites in Germany for tracing paleo-shear planes coated with graphitic deposits. The comparisons of geologic sections constructed in this paper (based on the proposed deterministic approach) against the existing published interpretations (obtained based on trial-and-error modeling) for the SP data of the KTB and Rittsteig sites have revealed that the deterministic approach suggests some new details that are of some geological significance. The findings of the proposed inverse scheme are supported by available drilling and other geophysical data. Furthermore, the real SP data of the Grossensees site have been interpreted (apparently for the first time ever) by the deterministic inverse scheme from which interpretive geologic cross sections are suggested. The computational efficiency, analysis of the numerical examples investigated, and comparisons of the real data inverted here have demonstrated that the developed deterministic approach is advantageous to the existing interpretation methods, and it is suitable for meaningful interpretation of SP data acquired elsewhere over graphitic occurrences on fault planes.

Background

Graphite precipitation on the planes of shear zones can play an important role in understanding the tectonic evolution of an area (e.g., Glover and Vine 1995; Haak 1989; Oohashia et al. 2012). This precipitation produces an anomalous self-potential (SP) response that can be inverted to retrieve the characteristic parameters (e.g., depth, extension in depth, and amount and direction of dip) of the graphite sheet formed on these shear planes (Bigalke et al. 2004; Dmitriev 2009, 2012; ELEKTB Group 1997; Stoll et al. 1995).

Kontny et al. (1997) discussed mechanisms of the graphite occurrence on these planes and reported that graphite in shear zones is predominantly associated with chlorite, plagioclase, and pyrite. In order to explain the SP anomaly measured over a graphite precipitation, Stoll et al. (1995) proposed an electrochemical model that considers the redox conditions and the electrode kinetics at the graphite surface. This model therefore focuses on the redox potential that varies with depth rather than on the groundwater flow and gradients of the chemical potential (Bigalke and Grabner 1992).

The ELEKTB Group (1997) studied the electrical conductivity of the continental crust in the German Continental Deep Drilling Program (Kontinentales Tiefbohrprogramm der Bundesrepublik Deutschla - KTB). They concluded that the high electrical conductivities in the upper crust are primarily caused by graphite accumulations rather than by fluids and that these anomalous conductivities are related to shearing stress regimes. However, it is essential to note that for both graphite and saline fluids, it is their connectivity rather than their pure presence that is the most crucial factor controlling conductivity (ELEKTB Group 1997; Emmermann and Lauterjung 1997; Haak et al. 1991; Stoll et al. 1995). The ELEKTB Group (1997) also found and concluded that the original (primary) graphite existing in the gneisses cannot be interconnected. This primary graphite represents original organic material in the protoliths and is so finely dispersed that it does not contribute to the overall electrical conductivity (Emmermann and Lauterjung 1997). The most important aspect for electrical conductivity is the secondary graphite, which is ubiquitous in the cataclastic shear zones, where it is always associated with iron sulfides and chlorite. This secondary graphite often forms a quasi-continuous (connective) and abundant coating along shear planes, and it is locally concentrated in millimeter-thick layers that constitute good electrical conductors over hundreds or thousands of meters in depth along the fault planes (Emmermann and Lauterjung 1997; Zulauf 1990).

Emmermann and Lauterjung (1997) reported that all the evidences suggest that this graphite was precipitated from hydrocarbon-bearing fluids at about 400°C and 2 kbar. The shearing process is the reason beyond the interconnection of the graphite on the shear planes (ELEKTB Group 1997). Nover et al. (2005) carried out a laboratory experiment to understand better the mechanisms that lead to the abundant formation of graphite in overthrust structures. They concluded that graphite is found on shear planes and attributed this graphite to shear movement, shear strain, and strain energy (Bustin et al. 1995; Mathez et al. 1995).

In light of these comprehensive findings in the literature, it is evident that SP methods can be used to map and trace paleo-shear zones associated with graphite sheets (films) precipitated on the shear planes of these zones. The measured SP signatures generated by these graphite sheets (Bigalke et al. 2004; Lehmann et al. 1998; Srivastava and Agarwal 2009; Stoll et al. 1995) can be inverted to extract information of some geologic significance such as depths, extensions in depth, and inclinations of the sheets. Here, this paper is exclusively devoted to the purpose of interpreting the SP data measured (observed) over the shear zones whose fault planes are coated by graphite films. For the diverse applications of the SP methods, the reader is referred to, for example, Essa et al. (2008), Srivastava and Agarwal (2009), Soueid Ahmeda et al. (2013), and Mehanee (2014a) and the references therein.

A few methods have been developed for the particular purpose mentioned above. Ram Babu and Atchuta Rao (1988) developed an iterative inversion approach, based on the Marquardt method, to interpret a SP profile by a two-dimensional (2D) sheet, a sphere, or a 2D horizontal cylinder model. The initial values of the model parameters (initial guess) are calculated automatically, prior to performing the SP iterative inversion, by first solving a linear inverse problem. They then applied the method to noise-free numerical examples in the restricted class of a sphere and a horizontal cylinder model. Ram Babu and Atchuta Rao (1988) also applied the method to a field data example with a sheet-like structure. They then tabulated and compared their results against those obtained by Murthy and Haricharan (1985) and noted some variations in the inverse parameters of the amplitude coefficient, extension in depth, and depth to the center.

It appears that the influence of the initial guesses on the convergence, stability, and misfit of the Marquardt method was not addressed or discussed in Ram Babu and Atchuta Rao (1988). Ram Babu and Atchuta Rao (1988) carried out the Marquardt inversion in the space of the model parameters themselves, which does not guarantee the positivity of these parameters. As an example, a negative value for the depth to the top of the sheet means that the sheet virtually extends into the air, which is geologically erroneous and unacceptable in this context (Sharma and Biswas 2013). It appears that the occasional negativity of the model parameters evolved from the Marquardt scheme was not addressed or treated in Ram Babu and Atchuta Rao (1988). Note that both the occasional negativity of the inverse model parameters and the effect of the initial guess on the convergence are satisfactorily addressed and treated here as will be seen in the “Formulation of the inverse problem and the proposed steps of its solving” and “Model 1” sections.

Jagannadha Rao et al. (1993) developed an iterative inversion technique that is also based on the Marquardt method for the interpretation of a SP profile by a 2D inclined thin-sheet body. The main objective of the interpretation is to determine the characteristic parameters of the sheet, which are as follows: h (depth to the top), a (extension in depth), θ (amount of dip), and k (amplitude coefficient or scale factor) (see the “Formulation of the forward modeling (direct) solution” section and Figure 1 of the present paper). Jagannadha Rao et al. (1993) performed the Marquardt inversion in the space of the model parameters themselves (h, a, θ, and k), which does not force the positivity of these parameters as will be seen and discussed in the “Formulation of the inverse problem and the proposed steps of its solving” and “Model 1” sections below. Again, the occasional negativity of the model parameters evolved from the Marquardt scheme was not addressed or treated in Jagannadha Rao et al. (1993).

Figure 1

A sketch showing a cross-sectional view, geometry, and the parameters of the sheet structure.

According to Jagannadha Rao et al. (1993), their approach differs from conventional approaches in that the initial values (initial guess) of the model parameters are calculated automatically within the computation program itself prior to carrying out the Marquardt iterative inversion. Precisely, Jagannadha Rao et al. (1993) determined a tentative initial guess for the three inverse parameters h, a, and θ of the sheet by measuring the maximum and minimum anomalies of the SP profile subject to the interpretation and by scaling a few distances (the so-called characteristic distances) out of this profile and using some characteristic ratios that are defined by these characteristic distances. As for the determination of the initial guess of the amplitude coefficient (k) inverse parameter, Jagannadha Rao et al. (1993) approximately estimated the initial value of k with knowledge of the predetermined initial values of the parameters h, a, and θ and by solving a separate least squares minimization problem. Before employing the aforementioned approximate initial guess of h, a, and θ in the Marquardt inversion approach, Jagannadha Rao et al. (1993) applied some additional computational procedures (also based on some characteristic distances) to further improve this initial guess.

Note that Jagannadha Rao et al. (1993) used a very good initial guess (nearly the true solution) in the Marquardt method when inverting the numerical example presented in their paper. It is also relevant to note that the same numerical example will be inverted and analyzed further in the present paper (see the “Model 1” section) in terms of the positivity of the model parameters and convergence using two different initial guesses. Jagannadha Rao et al. (1993) also used very good initial guesses (nearly the actual model parameters) when inverting and interpreting the two field data examples shown in their paper.

It appears that the use of a very good initial guess was essential for the convergence of the Marquardt method formulated in Jagannadha Rao et al. (1993) and hence for solving this particular inverse problem. It appears also that the influence of a relatively far initial guess on the convergence and stability of the Marquardt method was not addressed or discussed in Jagannadha Rao et al. (1993).

Stoll et al. (1995) introduced a new approach based on a trial-and-error modeling method for the interpretation of a SP profile measured over 2D sheet structures. This approach has been applied to a number of SP profiles measured over drilling-confirmed graphite precipitations on shear planes from the KTB and Rittsteig regions, Germany (Bigalke et al. 2004; Stoll et al. 1995). Note that these aforementioned SP profiles as well as the SP data set of Grossensees, Germany, are also inverted and analyzed here (“Field data inversion” section).

Cooper (1997) developed a computer code “SPINV” based on an iterative least square method for the interpretation of a SP profile by some geometrically simple model such as a sphere, a cylinder, a rod, and a 2D thin-sheet body. The SPINV code solves for the model parameters in the linear space, that is, solving for the model parameters themselves and not their logarithms. However, Cooper (1997) did not show numerical nor field data inversion results. It appears that the influence of the initial guess on the convergence and stability of the solver as well as the imposition of the positivity of the model parameters evolved from the inversion were not discussed either.

Radhakrishna Murthy et al. (2005) developed a technique that falls into the category of the Werner deconvolution, in which all the individual anomaly values and their positions along the profile are not considered during the calculations. Instead, the positions of the pairs of points (whose anomaly values differ from each other by a constant value and whose distances are measured from a reference point), which are selected from the profile, are used to construct linear equations whose coefficients are related to the body parameters. The body parameters are then obtained from these coefficients. Radhakrishna Murthy et al. (2005) applied the developed method to sheet-like mineralization data from the Surda area of the Rakha mines in India and the results obtained were found to be in good agreement with those reported in the literature.

Srivastava and Agarwal (2009) introduced a method based on the concept of enhanced local wave (ELW) number for the interpretation of SP anomalies over 2D and three-dimensional (3D) sources including sheet structures. It appears, however, that the method was applied to recover only one parameter (depth of the structure). It is relevant to note that mapping a paleo-shear zone requires knowledge of the depth, extension in depth, and amount and direction of inclination (Bigalke et al. 2004; Stoll et al. 1995) of the structure. Monteiro Santos (2010) described a method based on particle swarm optimization (PSO) to interpret SP data over 2D sheet-type structures. However, it appears that the developed method was applied to a real data set that is not sheet-type. Mehanee et al. (2011) presented a graphical method for determining only the width and depth of a thick sheet structure using the formula of the geometrical shape of a dike. Though the method was successfully applied to a number of case studies from mineral exploration, it is applicable to the total SP data (that is, the sum of the residual and regional SP signatures).

Dmitriev (2012) developed a powerful computer code “SPI-SV” (self-potential inversion, Siberian version) that uses a high-performance random search algorithm (with combined linear and nonlinear tactics) for the interpretation of a SP profile by up to 29 polarized bodies at a time. Dmitriev (2012) reported that this number of bodies is sufficient to image the details of a geologic section and that the code can simultaneously carry out automated fitting for a given number of bodies with respect to their characteristic parameters (e.g., depths, extensions in depth), which vary within user-specified intervals. The code may be used for forward modeling or inverse calculations or for both applications involving various geometrical bodies such as thick prisms, trihedral prisms, dipping thin sheets, and dipping layers with variable dip.

According to Dmitriev (2012), SPI-SV runs in two steps: in step 1, the positions and geometries of the bodies are determined in an automated mode, and in step 2 (the interactive mode), the faulty inverse parameters are corrected. Dmitriev (2012) successfully applied SPI-SV to three real data examples from mineral exploration, as well as the KTB SP profile in which 12 conductors were used to perfectly fit its data.

Sharma and Biswas (2013) introduced a very fast simulated-annealing (VFSA) global optimization technique for the interpretation of a SP anomaly measured over a 2D inclined sheet-type structure. They simultaneously inverted for the scale factor (amplitude coefficient, k), depth to the center of the body (z), depth extent (a), and amount of dip (θ) of the structure. The actual computation (not CPU) time of the developed technique is about 52 s on a personal computer (PC). They insightfully analyzed the method and reported that the presented approach can occasionally produce geologically erroneous inversion results. In such occasions, Sharma and Biswas (2013) precisely found that the depth to the top (h) of the sheet (computed from the inverse parameters z, a, and θ evolved from the VFSA inversion technique) is negative, which subsequently means that the sheet virtually extends into the air, which is geologically unacceptable and inappropriate in this context. As indicated earlier, the matter of the aforementioned unacceptable negative depth has been addressed in detail and satisfactorily treated here, as will be seen in the “Formulation of the inverse problem and the proposed steps of its solving” section.

Recent advances to solve the inverse problem of SP data are based on the use of deterministic approaches that utilize the regularized least squares techniques and the differentiability of the objective function subject to minimization (e.g., Mehanee 2014a; Soueid Ahmeda et al. 2013). This paper builds on the existing literature and describes a new and very fast method capable of interpreting the observed SP anomalies (caused primarily by graphite precipitation formed on a shear plane(s) within the continental crust) by a single 2D conductive thin sheet to trace paleo-shear zones. The proposed iterative approach solves the SP inverse problem in the framework of the regularized least squares techniques (Tarantola 2005; Tikhonov and Arsenin 1977; Wannamaker and Doerner 2002) using a mixed solver of the steepest descent (SD) and Gauss-Newton (GN) methods.

This newly proposed algorithm simultaneously determines the depth to the top (h), depth extent (a), amount of dip (θ), and amplitude coefficient (k) of the sheet from the SP data measured along a profile. It solves for h, a, and k in the logarithmic space (log(h), log(a), and log(k)) and for the dip angle (θ) in the linear space to impose the positivity of the inverse parameters. The aforementioned mixed solver of the SD and GN methods is essential in order to maintain the stability and convergence of this iterative approach in which arbitrary but reliable initial guesses can be used. Furthermore, this scheme alternatively formulates the forward problem (which is defined by a closed-form solution) in terms of the depth to the top of the sheet rather than in terms of the depth to its center. This assures that the inversion scheme always produces geologically acceptable and meaningful inverse results (that is, the body never virtually extends in the air), as will be discussed in the “Formulation of the inverse problem and the proposed steps of its solving” section. In addition, the presented inversion approach uses the entire SP observed data of a profile rather than just a few characteristic data points and distances out of the profile to simultaneously retrieve all inverse parameters of the causative conductive sheet.

To the best knowledge of the writer, deterministic inversion of SP data by sheet-type structures was not developed before. Furthermore, the real SP data sets of the KTB, Rittsteig, and Grossensees sites, Germany (which have been carefully analyzed here (see the “Field data inversion” section) and generated primarily by graphite occurrences on fault planes of shear zones) have not been interpreted before by deterministic inversion.

This paper is structured as follows. First, the forward modeling problem is described. Second, the formulation of this particular inverse problem in the logarithmed and non-logarithmed model parameters and its solving by the SD and GN methods is discussed. Third, the developed method is verified and analyzed on numerical models without noise and is then applied to noisy data. Finally, the method is carefully analyzed on seven field data examples from Germany.

Methods

Formulation of the forward modeling (direct) solution

The SP anomaly (mV) due to an inclined 2D sheet-like structure (Figure 1) at a point P(x) along a profile normal to the strike of the structure is given by the following formula (Jagannadha Rao et al. 1993; Rao and Babu 1983; Sharma and Biswas 2013):

$$ V\left(r_{1}, r_{2}, k\right) = k \; \log \frac{{r_{1}^{2}}}{{r_{2}^{2}}}, $$

((1))

where r ₁ and r ₂ are the distances from the top and bottom edges of the sheet to the observation point (m), respectively, and k (k=I ρ/2π, where I is the current density and ρ is the resistivity of the surrounding medium) is the amplitude coefficient (mV) (Jagannadha Rao et al. 1993; Rao and Babu 1983).

Sharma and Biswas (2013) formulated the forward solution (1) in terms of the depth to the center of the sheet (z) as

$$ V\left(x, z, a, k, \theta\right) = k \; \log \frac{\left[x^{2} + \left(z - \frac{a}{2}\sin\theta\right)^{2} \right]}{\left[ (x-a\cos\theta)^{2} + \left(z + \frac{a}{2}\sin\theta\right)^{2} \right]}, $$

((2))

where x is the coordinate of the measurement station (m), a is the extension (m) of the sheet in the z-direction (which is chosen as positive downward), and θ is its dip angle measured clockwise from the horizontal axis (Figure 1).

Formula (2) has been employed in the inversion scheme of Sharma and Biswas (2013) to obtain z, k, a, and θ subject to the condition $\left (z - \frac {a}{2}\sin \theta \right) \ge 0$. It is essential that this condition be satisfied as $\left (z - \frac {a}{2}\sin \theta \right)$ amounts to the depth to the top of the sheet (Figure 1). Consequently, the violation of this condition means that the sheet is virtually extending into the air, which is geologically unacceptable and invalid. Sharma and Biswas (2013) insightfully discussed this matter and discarded the inverse models that virtually extend into the air when the simulated-annealing global optimization inversion method is used (see the “Background” section).

In order to completely eliminate the aforementioned problem (that is, the virtual presence of the conductive sheet and its associated fault plane in the air), an alternative formulation for expression (2) is described here. Formula (2) can be re-written equivalently in terms of the depth to the top (h) of the sheet rather than the depth to its center (z) as (Figure 1)

$$ V\left(x, h, a, k, \theta\right) = k \; \log \frac{\left[x^{2} + h^{2} \right]}{\left[ (x-a\cos\theta)^{2} + \left(h + a\sin\theta\right)^{2} \right]}. $$

((3))

Although the SP modeling results obtained from formulas (2) and (3) are numerically identical, it was found (as will be seen in the next section) that inversion based on the usage of the forward modeling formula (3) is far better than that based on the usage of formula (2). It can be readily seen that formula (3) never produces inverse models extending into the air; the rationale behind this argument will be explained and clarified further in the next section. Thus, the inverse scheme that is based on (3) can be used for meaningful interpretation of the SP data measured over thin sheet-type structures (such as graphite films coating shear planes) buried in the continental crust.

Formulation of the inverse problem and the proposed steps of its solving

The conventional method of solving ill-posed inverse problems is based on the minimization of the Tikhonov parametric functional (e.g., Mehanee 2014a; Menke 1990, 2012; Tarantola 2005; van den Doel et al. 2013; Wannamaker and Doerner 2002; Zhdanov 2002):

$$ \Phi^{\alpha}\left(m,d_{\circ}\right) = \left\Vert A(m) - d_{\circ} \right\Vert^{2} + \alpha \; \left\Vert m \right\Vert^{2} = \text{min}, $$

((4))

where the first term is the data misfit functional determined as a square norm of the difference between the observed (measured) and predicted (computed) data, the second term is a stabilizing functional (the stabilizer), α is the regularization parameter, m is a column vector of the model parameters; m=[z,k,a,θ]^T, A(m) is the forward modeling operator that acts on m to produce the predicted data, d _∘ is a column vector of the observed SP data; $d_{\circ }=\left [d_{{\circ }_{1}}, d_{{\circ }_{2}}, d_{{\circ }_{3}},.....d_{{\circ }_{N}}\right ]^{T}$, T is the transposition operator, and N is the number of data points of the SP profile subject to inversion.

In order to guarantee and impose the positivity of the model parameters that we seek to invert for, we solve (4) in the mixed logarithmed and non-logarithmed space of the model parameters (that is, log(z), log(k), log(a),θ) rather than in the non-logarithmed space (that is, the space of the parameters themselves; z,k,a,θ) (Mehanee 2014a). The new logarithmed-space objective functional takes the form:

$$ \Psi^{\alpha}(\widetilde m,d_{\circ}) = \left\Vert \widetilde A(\widetilde m) - d_{\circ} \right\Vert^{2} + \alpha \; \left\Vert \widetilde m \right\Vert^{2} = \text{min}, $$

((5))

where $\widetilde m = \left [\log (z), \log (k), \log (a), \theta \right ]^{T}$.

The nonlinear minimization problem (5) is solved iteratively using a serial hybrid technique that automatically and successively combines both the SD and GN methods (e.g., Tarantola 1987; Zhdanov 2002). First, the SD method is used until a normalized misfit (defined as $\frac {\left \Vert \widetilde A(\widetilde m)-d_{\circ }\right \Vert }{\left \Vert d_{\circ }\right \Vert } \times 100\%$) below about 20% is reached. After that, the GN method, whose initial guess is the final inverse solution produced from the SD method, is used and applied to the same observed data. The iterative process of the GN method stops when the misfit becomes quasi-steady (that is, the update of the model parameters becomes negligible and insignificant).

It is found, as will be seen in the “Numerical models” section, that the inverse scheme can occasionally converge slowly or stagnate when the SD method is solely used. Here, the GN method is essential and advantageous as it is capable of overcoming the aforementioned stagnation. Furthermore, the GN method converges faster, but it requires a good initial guess (Zhdanov 2002) in order to converge. That is why the SD method, which does not demand a good initial guess to converge, is employed first in the developed scheme in order to generate and prepare an adequate initial guess for the GN method. Overall, solving this particular inverse problem by the abovementioned hybrid approach is very fast; it takes on average about 20 s of computation time on a 3.3 GHz PC.

The regularization parameter (α) is chosen such that it makes some balance between the misfit functional term and stabilizer (e.g., Mehanee 2014b; van den Doel et al. 2013). The conventional approach in solving an inverse problem is to run inversion several times using different starting values for the regularization parameter (Ramlau 2005; Reginska 1996), which is quite acceptable for the inverse problem under discussion as the approach developed here for solving this particular inverse problem is very fast; thus, the computation time is not really a concern. Ramlau (2005) examined and assessed the influence of the largeness and smallness of the regularization parameter on convergence. The selection of a range of regularization parameter values for solving a specific inverse problem is specific to that inverse problem.

Since the inverse scheme described above solves for the depth to the center (z) and extension in depth (a) (Figure 1) in the space of their logarithms [log(z), log(a)] when the minimization problem (5) is supplemented with the forward modeling formula (2) (which is essentially dependent upon the depth to the center z, among other parameters), then the model parameters z and a recovered from the inversion are always positive real numbers (e.g., Mehanee 2014a). However, there is a possibility that the condition $\left (z - \frac {a}{2}\sin \theta \right) \ge 0$ is violated in some inversion cases (Sharma and Biswas 2013), which consequently means that the sheet structure extends into the air - this is unacceptable and geologically erroneous as discussed earlier.

If we were to use (2) in (5), then it is essential we incorporate in (5) an additional term with a particular mathematical operator to impose the condition $\left (z - \frac {a}{2}\sin \theta \right) \ge 0$ in order to guarantee that the structure is always buried inside the subsurface. This term is considered as an additional regularizing term whose addition to the minimization scheme (5) can be computationally expensive and complicated.

Rather, we can alternatively and readily supplement the minimization problem (5) with the forward modeling formula (3) (which is a function of the depth to the top h, among other parameters), in which case the inverse algorithm solves for h and the extension in depth (a) (Figure 1) in the space of their logarithms [ log(h), log(a)]. This newly introduced inverse scheme (which combines (5) and (3)) is by far advantageous as it guarantees that the anomalous sheet structure is buried inside the earth. This is because the depth of expression (3) is taken to the top (h) of the body (in Figure 1, note that the z-axis is chosen as positive downward) and that this depth (h) and the extent (a) of the body (retrieved from inversion) are always positive real numbers being computed in the space of their logarithms.

Therefore, the inverse scheme formulated based on both the objective functional (5) and the forward modeling formula (3) is capable of extracting meaningful inversion results when applied to SP data sets. It is noted that all the inversion results shown in the “Numerical models” and “Field data inversion” sections are solely obtained from the minimization problem (5) and the modeling expression (3).

This inverse problem is nonlinear, and therefore, as indicated earlier, it must be solved iteratively. Thus, it requires the calculation of the Jacobian (Fr$\acute {\mathrm {e}}$chet) matrix, which is evaluated analytically by differentiating the forward modeling formula (3) with respect to the model parameters we invert for.

Results and discussion

Numerical models

In order to illustrate the accuracy and effectiveness of the proposed SP inversion algorithm, it is first applied to and analyzed on a published numerical example (Model 1). Following this, noise-free and noisy data sets of realistic models that can resemble actual geologic settings in the continental crust are inverted to carefully assess the developed algorithm. The synthetic (numerical) data sets were produced using formula (2) or (3). As indicated in the “Formulation of the forward modeling (direct) solution” section, these two SP formulas generate identical forward modeling results.

Model 1

Jagannadha Rao et al. (1993) illustrated the Marquardt inversion approach on a SP profile generated by a sheet anomaly that was 1 m deep (h=1 m), 3 m long (a=3 m), and had a 90° dip (θ=90°). It appears that the magnitude of the scale factor (amplitude coefficient) parameter (k) of the sheet was not provided in their paper. However, based upon simulation experiments, it can be readily inferred that this parameter has a magnitude of about 100 mV (k=100 mV). Figure 2 shows the SP response of this model, which will be subject to inversion here.

Figure 2

Model 1: noise-free data (after Jagannadha Rao et al. [ 1993 ]). Inversion results obtained by the Marquardt method (in which the inverse parameters are formulated in the non-logarithmed space; h,a,k, and θ) from “Initial guess 1” that is shown in the figure. See the text for details. The observed and predicted SP data are coincident with each other.

Jagannadha Rao et al. (1993) used an initial guess of h=0.95 m, a=2.48 m, k=100 mV, and θ=90° when inverting this SP profile. It can be readily seen that this initial guess is very close to the actual model (h=1 m, a=3 m, k=100 mV, and θ=90°). According to Jagannadha Rao et al. (1993), the iterative process of the Marquardt inversion terminated when the objective function subject to minimization reached a value of 0.02 mV ².

In order to examine the influence of the initial guess on the convergence and stability of both of the method presented here and the Marquardt method, two different initial guesses (namely, “Initial guess 1” and “Initial guess 2”) are used in both methods to invert the SP anomaly shown in Figure 2. These two initial guesses have been chosen such that “Initial guess 1” is very close to the actual model, while “Initial guess 2” is slightly far from the actual model.

Figure 2 shows the inversion results obtained using “Initial guess 1” by the Marquardt method (e.g., Jagannadha Rao et al. 1993; Marquardt 1963; Ram Babu and Atchuta Rao 1988) in which the inverse parameters are formulated in the non-logarithmed space (that is, h, k, a, and θ). As is seen, the Marquardt method has successfully recovered the true parameters from this initial guess in just eight iterations with a final normalized misfit of about 10⁻⁸. Behavior of the normalized misfit and objective functional of the Marquardt method is presented at the middle and bottom panels of Figure 2. Figure 3 illustrates the corresponding inverse solution recovered from “Initial guess 1” by the Marquardt method as well in which the inverse parameters this time are expressed in the logarithmed and non-logarithmed space (log(h), log(k), log(a), and θ) in order to maintain their positivity. The panels of this figure indicate that the Marquardt method has also successfully converged and retrieved the actual parameters when applied in the aforementioned space.

Figure 3

Model 1: noise-free data (after Jagannadha Rao et al. [ 1993 ]). Inversion results obtained by the Marquardt method (in which the inverse parameters are formulated in the mixed logarithmed and non-logarithmed space; log(h), log(k), log(a),and θ) from “Initial guess 1” that is shown in the figure. See the text for details. The observed and predicted SP data are coincident with each other.

The inversion results obtained from “Initial guess 1” by the method proposed here (in which the inverse parameters are formulated in the logarithmed and non-logarithmed space; log(h), log(k), log(a), and θ) are presented in Figure 4. The figure shows that the true parameters are recovered in 10 iterations with a normalized misfit of about 10⁻⁸. Evolution of the normalized misfit of the SD and GN methods is presented at the middle panel of Figure 4. Evolution of the objective functional is rendered at the bottom panel of this figure. Note that, as has been discussed in the “Formulation of the inverse problem and the proposed steps of its solving” section, the SD method was applied to the abovementioned SP profile. Following that, the GN method was applied to the same data to avoid possible occasional stagnation with the SD method. It is worthy to note that the preliminary inverse solution evolved from the SD method is automatically used as the initial guess for the GN method. Thus, both the Marquardt method and the method developed here have successfully converged and retrieved the actual model parameters when separately applied to this SP profile using “Initial guess 1.”

Figure 4

Model 1: noise-free data (after Jagannadha Rao et al. [ 1993 ]). Inversion results obtained by the method developed here (in which the inverse parameters are formulated in the mixed logarithmed and non-logarithmed space; log(h), log(k), log(a), and θ) from “Initial guess 1” that is shown in the figure. “SD” and “GN” stand for the steepest descent and Gauss-Newton methods. See the text for details. The observed and predicted SP data are coincident with each other.

The model parameters of “Initial guess 2” are: h=10 m, a=20 m, k=1,000 mV, and θ=10°. The Marquardt inversion method (formulated in the space of the model parameters themselves (h, k, a, and θ)) produced negative values to the inverse model parameters k and h (which are positive real numbers in theory) when applied to the same SP data using “Initial guess 2”. Therefore, the code terminated after the first iteration.

Figure 5 depicts the inverse solution retrieved by the Marquardt method in which the inverse parameters are formulated in the logarithmed and non-logarithmed space (log(h), log(k), log(a), and θ), from “Initial guess 2” using a varying (that is a non-fixed) α during the iterative process. It can be seen from this figure that the Marquardt method suffered some divergence and could not adequately fit the observed data; it therefore could not recover the actual model parameters. The inversion results obtained by the method developed here from “Initial guess 2” are illustrated in Figure 6. The figure shows that the developed approach is stable, convergent, and it successfully recovered the true model parameters.

Figure 5

Model 1: noise-free data (after Jagannadha Rao et al. [ 1993 ]). Inversion results obtained by the Marquardt method (in which the inverse parameters are formulated in the mixed logarithmed and non-logarithmed space; log(h), log(k), log(a), and θ) using “Initial guess 2” that is shown in the figure. The observed and predicted SP data are marked by open and solid circles. The values of the regularization parameter (α) are also shown; see the text for details.

Figure 6

Model 1: noise-free data (after Jagannadha Rao et al. [ 1993 ]). Inversion results obtained by the method developed here (in which the inverse parameters are formulated in the mixed logarithmed and non-logarithmed space; log(h), log(k), log(a), and θ) from “Initial guess 2” that is shown in the figure. See the text for details. The observed and predicted SP data are coincident with each other.

On the basis of the aforementioned analysis, it can be concluded that the inverse scheme developed here is more advantageous than the Marquardt method in terms of convergence, stability, and flexibility in the choice and usage of initial guesses. Note that all inversion results shown hereinafter are based on the method developed in this paper.

Model 2

The SP numerical response of Model 2 (h=77 unit, a = 110 unit, k = 125 mV, θ = 75°, N = 121) was generated (Figure 7, top panel). This response was used as the observed data in inversion. It is common practice in this class of research to use the dimension “unit” for h, a, and x (e.g., Mehanee 2014a and the references therein). In numerical examples, one can realistically choose the aforementioned unit to equal, for example, 10, 20, or 25 m. In real data inversions, this unit dimension is replaced with the actual dimension (e.g., m) of the SP profile, as will be seen in the “Field data inversion” section.

Figure 7

Model 2: noise-free data. Obtained inversion results (top panel). The observed data and predicted response are coincident with each other. Evolution of the normalized misfit of the SD and GN methods (middle panel). Evolution of the objective functional of the SD and GN methods (bottom panel).

The initial guess used, the true model parameters, and the model parameters retrieved from inversion are shown in the top panel of Figure 7. Evolution of the normalized misfit and objective functional of the SD and GN methods is presented at the middle and bottom panels of Figure 7. It can be seen that this algorithm has converged and successfully recovered the true value of all four inverse parameters of the sheet body.

Model 3

The numerical data set of this model (h=37.5 unit, a=50 unit, k=50 mV, θ=30°, N=226) has been contaminated with about 7% noise. This set was inverted in an analogous manner to Model 2. The top panel of Figure 8 shows the approximative solution recovered from inversion. It can be seen that the obtained inverse parameters are acceptably accurate. Behavior of the normalized misfit and objective functional of the SD and GN methods is, respectively, rendered at the middle and bottom panels of this figure.

Figure 8

Model 3: data contaminated with 7% noise. Inversion results (top panel). Behavior of the normalized misfit and objective functional of the SD method (middle panel) and the GN method (bottom panel). The observed and predicted data are shown in open and solid circles, respectively.

Model 4

This model represents a large scale shear zone, and it is essentially different from the three aforementioned models in the sense that it is defined by a much greater extension in depth (a) and a much larger dip angle (h=100 unit, a=2,000 unit, k=30 mV, θ=130°, N=226).

Analysis of the noise effect

The data were corrupted by about 3% noise and then inverted as described earlier. Figure 9 shows that the approximative solution evolved from inversion is stable and reasonably accurate. The middle and bottom panels of this figure, respectively, present the behavior of the normalized misfit and objective functional of the SD and GN methods.

Figure 9

Model 4: data contaminated with 3% noise. Inversion results (top panel). Behavior of the normalized misfit and objective functional of the SD method (middle panel) and the GN method (bottom panel). The observed and predicted data are shown in open and solid circles, respectively.

In order to examine the effect of the noise further, this profile was contaminated by 15% noise. Figure 10 illustrates the corresponding inversion results. This figure suggests that the embedded noise outliers do not have a significant effect on the accuracy of the yielded results.

Figure 10

Model 4: data contaminated with 15% noise. As in Figure 9.

On the basis of the analysis presented in the “Model 3” and “Analysis of the noise effect” sections, it can be concluded that the developed technique can extract reasonably accurate inversion results even when the observed data are distorted by noise artifacts similar to those reported above.

Analysis of SP profile truncation (shortening)

In order to assess and understand better the influence of truncation of the SP profiles on the accuracy and stability of the inversion results, various truncation levels (from the left, right, or from both flanks of the profile) are applied to the noise free data of this model (Model 4). The short profiles evolved are separately inverted and discussed.

Figure 11 shows the inversion results obtained from the SP profile that has been truncated from the left side. The corresponding inversion results retrieved from the SP profile that has been overtruncated from the left side as well are presented in Figure 12. Both figures show that the developed method is capable of recovering the actual model parameters of the model.

Figure 11

Model 4 truncated from the left : noise free. Inversion results (top panel). Behavior of the normalized misfit and objective functional of the SD method (middle panel) and the GN method (bottom panel). The observed and predicted data shown in open and solid circles, respectively, are coincident with each other.

Figure 12

Model 4 over truncated from the left : noise free. Inversion results (top panel). Behavior of the normalized misfit and objective functional of the SD method (middle panel) and the GN method (bottom panel). The observed and predicted data shown in open and solid circles, respectively, are coincident with each other.

The results of the SP profiles that have been truncated and over-truncated from the right side are depicted in Figures 13 and 14, respectively. Figure 15 illustrates the inversion results of the SP profile over truncated from both the left and right sides. It can be seen from these three figures that the developed inverse scheme has successfully rendered the true model parameters.

Figure 13

Model 4 truncated from the right : noise free. Inversion results (top panel). Behavior of the normalized misfit and objective functional of the SD method (middle panel) and the GN method (bottom panel). The observed and predicted data shown in open and solid circles, respectively, are coincident with each other.

Figure 14

Model 4 over truncated from the right : noise free. Inversion results (top panel). Behavior of the normalized misfit and objective functional of the SD method (middle panel) and the GN method (bottom panel). The observed and predicted data shown in open and solid circles, respectively, are coincident with each other.

Figure 15

Model 4 over truncated from both the left and right : noise free. Inversion results (top panel). Behavior of the normalized misfit and objective functional of the SD method (middle panel) and the GN method (bottom panel). The observed and predicted data shown in open and solid circles, respectively, are coincident with each other.

In light of the analysis illustrated above, it can be concluded that the developed approach can recover accurate inverse results from truncated SP profiles. It is relevant to note that this analysis was essential prior to inverting the short SP profiles of the field examples investigated in the “Field data inversion” section.

Field data inversion

In order to examine and assess the applicability of the presented technique, seven SP profiles from the KTB, Rittsteig, and Grossensees sites, Germany, were inverted and interpreted. These profiles were acquired over graphite sheets coating shear zones within the continental crust, and they show the corresponding residual SP data sets of these zones. These data sets do not require residual-regional separation as they are the residual anomalies. That is why they have been inverted and interpreted without a regional correction in numerous published papers (e.g., Bigalke et al. 2004; Dmitriev 2012; Mehanee 2014a; Revil et al. 2001; Srivastava and Agarwal 2009; Stoll et al. 1995). Regional anomalies are due to what are generally called telluric currents, and these may amount to 100 mV per kilometer or more (Meiser 1962; Yüngül 1950). Usually, the regional component is corrected for in the stage pertinent to processing the raw SP data by removing some constant part (Cooper 1997; Jagannadha Rao et al. 1993; Mehanee 2014a; Ram Babu and Atchuta Rao 1988; Yüngül 1950).

KTB-borehole anomaly, Germany

Two research boreholes were drilled into the Variscan basement at Oberpfalz, Bavaria, in the Zone of Erbendorf-VohenstrauB (ZEV) (Franke 1989; Kontny et al. 1997). A pilot hole, the KTB-VB, went to about 4 km deep, and a subsequent superdeep main hole, the KTB-HB, went to about 9.1 km deep; the two boreholes were about 200 m apart (Emmermann and Lauterjung 1997).

On the basis of these two boreholes, it is evident that the graphite is associated with a number of steeply inclined shear planes (ELEKTB Group 1997; Nover et al. 2005; Stoll et al. 1995). As indicated in the “Background” section, this graphite forms an interconnected coating along shear planes over a distance of hundreds or thousands of meters, and it is locally concentrated in millimeter-thick layers that constitute good electrical conductors (Emmermann and Lauterjung 1997).

Kontny et al. (1997) reported that the ZEV zone is characterized by a medium-pressure, high-temperature metamorphism. The rocks of the ZEV are separated from the Permo-Mesozoic sedimentary basin by the Franconian fault system. This fault system has a NW-SE strike, dips steeply towards the NE direction, and crosscuts the KTB borehole at a depth of about 7 km.

A SP profile (redrawn from Figure two(a) of Stoll et al. 1995) was recorded in the vicinity of the KTB boreholes (Siwczyk 1990; Stoll et al. 1995). This profile (Figure 16) was digitized at 5-m intervals, and it is characterized by two distinct anomalies (Anomaly I and Anomaly II) with magnitudes of approximately −500 mV and −600 mV.

Figure 16

The KTB SP profile (redrawn from Figure two(a) of Stoll et al. 1995 ).

Various initial guesses have been used in the developed scheme when inverting Anomaly I, all of which converged to a normalized misfit of about 9.26% and yielded quasi-similar inverse solutions. The corresponding results are shown in the panels of Figure 17, which suggest that the inverse parameters of the causative conductive sheet range from (h=27 m, a=2,168 m, k=59 mV, θ=120°) to (h=27 m, a=2,213 m, k=59 mV, θ=120°). The suggested interpretation conforms well to, and is closely supported by, the redox potential measurements carried out in one of the KTB boreholes to a total depth of 1,500 m (see Figure three(b) in Stoll et al.1995). These redox potential measurements exhibited continued variations (that is, a non-zero gradient) along the entire 1,500 m investigated. In an analogous manner, Anomaly II has been inverted. Figure 18 shows the results obtained for which the misfit reached about 6.82%, which suggest that the inverse parameters of this sheet are (h=26 m, a=745 m, k=84 mV, θ=65°).

Figure 17

The KTB Anomaly I: inversion results retrieved from various initial guesses. The measured data (redrawn from Figure two(a) of Stoll et al.1995) and the predicted response are shown in open and solid circles, respectively.

Figure 18

The KTB Anomaly II: inversion results retrieved from various initial guesses. The measured data (redrawn from Figure two(a) of Stoll et al.1995) and the predicted response are shown in open and solid circles, respectively.

Both the amount and direction of dip of the two interpretive conductive sheets recovered by the developed algorithm (top panel of Figure 19) from the inversion of the SP data of Anomaly I and Anomaly II are in very good agreement with those of the shear planes confirmed from drilling and other geophysical research on which graphitization has occurred (Figure 20).

Figure 19

The KTB site. Sketch showing the approximative results (see Figures 17 and 18 for the precise numerical values obtained from inversion) recovered by the developed inversion method (top panel) and the results yielded by the trial-and-error modeling method (Stoll et al.1995) (bottom panel) for the KTB SP profile. The bottom panel is redrawn from Figure five of Stoll et al.1995.

Figure 20

The KTB site. Sketch showing the geologic setting in the vicinity of the KTB-HB borehole (redrawn From Figure six of Stoll et al.1995; for the legend of the various rock types shown, see that figure). Note that the fault system F2 extends further in depth to about 4 km (see Figure four in Emmermann and Lauterjung1997).

Stoll et al. (1995) interpreted and fitted these two SP anomalies by 2D trial-and-error modeling with two sheet-like electric conductors (bottom panel of Figure 19): (h=50 m, a=490 m, θ=60°) and (h=30 m, a=490 m, θ=60°). Stoll et al. (1995) and Revil et al. (2001) pointed out that the electric field is established by the gradient of the redox potential (Figure three(b) in Stoll et al.1995), a local quantity dependent upon the local conditions and on the local pressure and temperature, which is interconnected by an electronically conducting material (graphite). This gradient is an important factor as it regulates the amplitude of the measured SP anomaly (ELEKTB Group1997).

Revil et al. (2001) applied a 2D tomography technique to this profile and interpreted the aforementioned SP anomalies as dipolar sources extending from 100 to 200 m deep. Srivastava and Agarwal (2009) inverted each of these two anomalies using the ELW number technique and obtained a depth of 135.5 and 80 m for Anomaly I and Anomaly II, respectively. On the basis of the resulting inverse parameters, they interpreted these anomalies as quasi-horizontal cylinders or quasi-sheet-like structures, which suggests that the abovementioned depths are in fact measured to the center of the anomalous bodies (Abdelrahman et al.1998). These analyses and comparisons show that the method developed here can provide realistic and reliable inversion results when applied to observed SP data generated by a sheet-type structure such as those inverted here.

Rittsteig anomaly, Bavaria, Germany

Surface SP and surface-downhole induced polarization (IP) measurements have been acquired in a 583-m deep borehole drilled in Rittsteig to help trace paleo-shear zones in this area (Bigalke and Junge2004; Bigalke et al.1999). Figure 21 (redrawn from Figure four of Bigalke et al.2004) illustrates a SSE-NNW SP profile that shows two distinct SP anomalies (Anomaly I and Anomaly II). This profile has been sampled into 5-m intervals.

Figure 21

The Rittsteig SP anomaly (redrawn from Figure four of Bigalke et al. 2004 ).

Bigalke et al. (2004) reported that the aforementioned borehole intersected several graphite shear zones (see Figure three in Bigalke et al.2004) and that out of all of these zones, the two shear zones that are of major importance are the ones located at depths of 320 and 460 m. The first shear zone cuts through Moldanubian biotite-muscovite schists at a depth of 320 m. The second shear zone separates Tepla-Barrandian amphibolites from Moldanubian biotite-muscovite schists at 460 m depth. These two main shear zones will be discussed further below.

Both Anomaly I and Anomaly II were inverted in a manner similar to that described earlier using three different initial guesses. The inversion results of Anomaly I, all of which converged to a normalized misfit of about 11.77%, are presented in Figure 22. Anomaly II has two negative peaks at x = 0 and 20 m, for which the position of the top of the interpretive sheet was fixed at x = 0 m in the inversion. Figure 23 shows the approximative inverse solutions of Anomaly II, all of which reached a normalized misfit of about 22.38%. The relatively large misfit of Anomaly II could be attributed to the presence of the two negative peaks located at 0 and 20 m. These two peaks were probably due to two or more separate graphitic layers formed on shear planes. The fitting of this profile could be improved further by employing simultaneous inversion of multi-inclined sheets (Dmitriev2012), which will be the subject of future research as indicated in the “Conclusions” section.

Figure 22

The Rittsteig Anomaly I: inversion results retrieved from various initial guesses. The measured data and the predicted response are shown in open and solid circles, respectively.

Figure 23

The Rittsteig Anomaly II: as in Figure 22 .

The top panel of Figure 24 illustrates a sketch for the recovered inverse solutions; F1 and F2 are the interpretive graphitic sheet models for Anomaly I and Anomaly II, respectively. As can be seen, this interpretation suggests two paleo-shear zones dipping to the SSE (F1) and NNW (F2). It is worthy to note that the left side of the top, middle, and bottom panels also shows the abovementioned two main shear zones intersected by drilling; both of which dip in the SSE direction as confirmed from drilling (Bigalke et al.2004).

Figure 24

The Rittsteig site. Sketch showing the interpretive sections obtained from the SP data analysis (top and middle panels) and from the induced polarization (IP) data analysis (bottom panel). Top panel corresponds to the approximative average results (for the corresponding precise numerical values obtained from inversion, see Figures 22 and 23) recovered by the developed inversion method. Middle panel corresponds to the results obtained by the trial-and-error modeling method (Bigalke et al.2004). The middle and bottom panels are redrawn, respectively, from Figure five of Bigalke et al. (2004) and Figure four of Bigalke and Junge (1999). A borehole and the faults confirmed from drilling are shown at the left side of each panel.

The conductive sheet F1 of Figure 24 can probably be considered as a collective interpretive model for both of the two shear zones (shown at the left sides of the panels of the figure) that were encountered in drilling.

The middle panel of Figure 24 shows the interpretive model for this same SP profile proposed by Bigalke et al. (2004) based on the 2D trial-and-error SP modeling method. The bottom panel of Figure 24 shows the interpretive model suggested by Bigalke and Junge (1999) based on the surface-downhole IP measurements carried out in the aforementioned borehole. It appears that the IP research mainly focused on the conductive sheet F1.

Grossensees anomaly, Germany

Figure 25 presents the SP contour map that was taken over a graphitic shear zone (number 4 in Figure one of Stoll et al.1995) in Grossensees. The northern part of this map shows a prominent anomaly with NW-SE strike that appears to be associated with three closures. Having said that, it would probably be wiser to take a profile across each of these three closures rather than just one profile across the center of the middle closure. The three profiles, $A\bar {A}$, $B\bar {B}$, and $C\bar {C}$, which were taken normal to the abovementioned strike, are illustrated in Figure 25. Each profile has been discretized into 20-m intervals. In fact, to the best knowledge of the writer, the literature does not appear to contain any available geological or geophysical information about this particular site.

Figure 25

The Grossensees self-potential anomaly map (provided by Stoll, 2013 via personal communication).

Figure 26 renders the inversion results of profile $A\bar {A}$ obtained from three various initial guesses, all of which converged to a normalized misfit of about 6.86% and yielded nearly the same inverse solution. Profiles $B\bar {B}$ and $C\bar {C}$ were inverted in an analogous manner. The normalized misfits reached for these two profiles were about 12.02% and 17.54%, respectively, and the corresponding results are presented in Figures 27 and 28. Figure 29 depicts the suggested geologic cross sections for the three profiles that were inverted.

Figure 26

The Grossensees anomaly (profile $\boldsymbol {A \; \widetilde A}$ shown in Figure 25 ): inversion results retrieved from various initial guesses. The measured data and the predicted response are shown in open and solid circles, respectively.

Figure 27

The Grossensees anomaly (profile $\boldsymbol {B \; \widetilde B}$ shown in Figure 25 ): inversion results retrieved from various initial guesses. The measured data and the predicted response are shown in open and solid circles, respectively.

Figure 28

The Grossensees anomaly (profile $\boldsymbol {C \; \widetilde C}$ shown in Figure 25 ): inversion results retrieved from various initial guesses. The measured data and the predicted response are shown in open and solid circles, respectively.

Figure 29

The Grossensees site: sketch showing the approximative results shown in Figures 26, 27, and 28 . See these figures for the precise numerical values recovered from inversion.

Based on the inverse results of profiles $A\bar {A}$, $B\bar {B}$, and $C\bar {C}$ presented in Figure 29, a final interpretive geologic cross section of a conductive sheet of (h=50 m, a=7,500 m, θ=142°) is suggested for the three anomalous SP closures located at the northern side of the Grossensees site. The parameters of this sheet are the average values of the corresponding inverse parameters of profiles $A\bar {A}$, $B\bar {B}$, and $C\bar {C}$ that are shown in Figure 29. As stated earlier, it appears that the presented Grossensees data have been interpreted here for the first time ever, and unfortunately, it appears that the literature does not have any available information (e.g., drilling and/or geologic) about this site that can be used to assess the accuracy of the interpretive geologic cross section suggested here.

It is very important and relevant to note for the reader that the SP profile shown in Figure two(b) of Stoll et al. (1995) was taken on the strike (Stoll, 2013 via personal communication) of an SP anomaly at the Grossensees site. In other words, that particular profile was not actually normal to the strike of the SP anomaly. Since the SP profile subject to inversion should be normal to the strike of the anomaly (e.g., Mehanee 2014a; Sharma and Biswas 2013), then that profile is not suitable for inversion, and accordingly, it has not been considered here (nor in Stoll et al. (1995)) for analysis. Thus, it has been completely excluded from the interpretation herein. Consequently, alternative profiles ($A\bar {A}$, $B\bar {B}$, and $C\bar {C}$) normal to the strike of the anomaly were obtained from the Grossensees SP map shown in Figure 25, and these were inverted and interpreted here as discussed above.

Extended discussion

A newly developed rigorous deterministic inversion approach was employed to trace paleo-shear planes coated by graphite precipitations from the SP data taken along a profile. The approach fits each anomaly of the observed (measured) SP profile by a single conducting sheet. Such a sheet represents and resembles the thin graphite film(s) present on the buried shear plane(s); see, for example, the interpretive conducting sheets suggested by the presented approach (F1 and F2 in the top panel of Figure 19) and their corresponding actual shear fault systems, which were confirmed from drilling and other geophysical data at the KTB site (F1 and F2 in Figure 20).

Comparisons of the interpretive section suggested here based on the deterministic approach (top panel of Figure 19) for the KTB site against the published interpretive section of Stoll et al. (1995) that is based on the trial-and-error modeling method (bottom panel of Figure 19) exhibit some reasonable agreement in the right side of the profile. However, the left side of the profile exhibited some variations in terms of the depth extension and direction of dip, as can be seen in the comparisons. It appears that the results of both the right and left sides obtained here by the developed scheme are supported by actual shear fault systems shown in Figure 20.

As for the Rittsteig SP profile, Figure 24 shows the interpretive sections obtained from the SP data analysis (top and middle panels) and from the induced polarization (IP) data analysis (bottom panel). One can see some significant dissimilarities between the interpretive model suggested by the proposed scheme (top panel) and the published model of Bigalke et al. (2004) (middle panel). The interpretation yielded from the IP data (Bigalke and Junge 1999) (F1 in the bottom panel) and drilling information (Figure 24 left side of each panel) appear to support the conductive sheet recovered by the present method, which dips SSE (F1 in the top panel). As for the right SP anomaly (Anomaly II of Figure 21), the algorithm proposed here suggests a paleo-shear zone with an amount of dip of 40° NNW (F2 in the top panel).

It has been found, based on the field data inversions, that the conjugate gradient (CG) method is occasionally unstable and divergent (corresponding plots not shown here). Based on the same field data experiments, it has been found that the SD method is stable and convergent for this particular SP inverse problem as has been seen. That is why the SD method was used in the inversion approach developed here.

Conclusions

A new method for the interpretation of the SP data measured over sheet-type structures has been developed. The algorithm of this method lies in the framework of the deterministic (gradient type) inversion, in which a parametric (objective) functional is minimized using both the steepest descent and Gauss-Newton methods in the mixed space of the logarithmed and non-logarithmed inverse model parameters. It was found essential that the algorithm inverts for the depth to the top (h) of the structure rather than the depth to its center in order to assure that the algorithm produces meaningful and realistic geologic results. Furthermore, the method is capable of retrieving the extension in depth (a), amount of dip (θ), and the amplitude coefficient (k) of the sheet structure.

The presented approach is very fast - on average it takes about 20 s of computation time. It is suitable for inverting real SP data generated by graphitic layers with a few millimeters of thickness occurring on shear faults. The evaluation of such data offers the chance to construct a reliable interpretive cross section(s) in a minimal amount of time.

The accuracy and robustness of the proposed method were successfully verified on synthetic examples with and without noise. After that, the method was applied to seven SP profiles from the KTB, Rittsteig, and Grossensees sites in Germany for tracing graphite-bearing shear planes. Comparisons of the interpretive geologic cross sections (obtained here based on the developed deterministic approach) against the existing published interpretations of the SP data of the KTB and Rittsteig sites have shown that the deterministic approach is capable of suggesting and extracting some new information that can more accurately describe the underlying shear zones. Moreover, the real SP data of the Grossensees site have been interpreted and geologic cross sections illustrating the shear faults proposed for this site are suggested.

The computational efficiency of the developed algorithm and analyses and comparisons of the numerical and real data inverted here have demonstrated that the newly developed deterministic SP method is advantageous to the existing SP interpretation methods that are based on the trial-and-error technique and on gradient-type approaches. Therefore, the newly described method is suitable for meaningful interpretation of SP data acquired elsewhere over interconnected graphite occurrences on fault planes within the continental crust.

The gradient-type inversion scheme described here does not account for the coupling effects that involve multiple sources (that is, the superposition response of several adjacent thin sheets), which is a limitation. Nevertheless, this approach produced some encouraging results as has been seen in the KTB and Rittsteig profiles that are characterized by multi-shear planes. These effects can be accounted for by employing simultaneous inversion of multi-inclined sheets, which will be the subject of future research.