Abstract
The geomagnetic deep sounding (GDS) method is one of electromagnetic (EM) methods in geophysics that allows the estimation of the subsurface electrical conductivity distribution. This paper presents the inversion modeling of GDS data employing Markov Chain Monte Carlo (MCMC) algorithm to evaluate the marginal posterior probability of the model parameters. We used thin-sheet model to represent quasi-3D conductivity variations in the heterogeneous subsurface. The algorithm was applied to invert field GDS data from the zone covering an area that spans from eastern margin of the Bohemian Massif to the West Carpathians in Europe. Conductivity anomalies obtained from this study confirm the well-known large-scale tectonic setting of the area.
1. Introduction
In geomagnetic deep sounding (GDS), we measure the natural Earth’s magnetic transient variations to infer large-scale subsurface conductivity distribution. Recent advances in magnetotelluric (MT) technique tend to put the attention to move towards local scale investigation of conductivity anomalies with more economic interests as in exploration for mineral, geothermal, or hydrocarbon. However, GDS is still considered as the most appropriate natural source electromagnetic (EM) method capable of imaging the Earth’s interior especially for tectonic study at the regional and continental scales (e.g., [1, 2]).
This paper describes the inversion modeling technique for GDS data in terms of conductivity distribution by using the Markov Chain Monte Carlo (MCMC) algorithm. The MCMC inversion algorithm has been applied for 1D modeling of MT [3, 4], vertical electrical sounding (VES) [5] and also controlled-source audio-frequency MT (CSAMT) [6] with satisfactory results. In the so-called thin-sheet modeling, the conductivity variations are assumed to be confined in a thin layer such that the model parameters are integrated conductivity over the thickness of the thin layer. This simplification suits well for modeling GDS data that are mostly sensitive to lateral variation of conductivity [7, 8].
We will first briefly review the concept and the formulation of the thin-sheet modeling and then describe the MCMC inversion algorithm. The result of inversion of real GDS data from Bohemian Massif-West Carpathians and also its interpretation will be discussed.
2. Thin-Sheet Electromagnetic Modeling
In electromagnetic (EM) geophysics, the thin-sheet modeling refers to an approximation of 3D conductivity variation by a thin layer with variable conductance, that is, integrated conductivity over the thickness of the thin layer. Such approximation is generally valid in the large-scale studies where the heterogeneities are confined in a layer with thickness much smaller than the penetration depth of EM fields. The thin-sheet modeling significantly simplifies the resolution of the Maxwell’s equations describing the EM fields in quasi-3D media. We used an algorithm employing integral equation method for thin-sheet EM modeling proposed by Vasseur and Weidelt [1]. By using the integral equation method, the computation domain covers only the anomalous zone where the conductance differs from that of normal (homogeneous or stratified) host medium (Figure 1).
The thin-sheet containing heterogeneities is discretized into rectangular uniform blocks in which the conductance is assumed constant. The total electric field in the th cell located at the heterogeneous layer is expressed by where is the normal electric field associated with the normal conductivity distribution, is the Green kernel representing the electric field at due to a unitary dipole at , and denotes the anomalous conductance of the th block, , . The () factors are imaginary part of complex numbers, angular frequency, and magnetic permeability of the free space, respectively. Equation (1) is in fact a system of linear equations of complex values of the orthogonal electric fields at each block associated with a certain conductance distribution, with being the total number of blocks [1].
For a large number of blocks in the anomalous domain, a direct matrix inversion to resolve (1) may be unstable and time consuming. Considering its form, (1) is more suitable to be resolved by the Gauss-Seidel method, where an initial value is iteratively updated to obtain the solution. At th iteration, the estimated electric field is calculated by
The normal electric field is commonly used as the initial value for , that is, . A relatively fast convergence can be reached since at the subsequent iterations all updated elements of the solution ( for ) are used to estimate the remaining elements. For low to moderate conductivity contrast, convergent solution is generally achieved after 10 to 20 iterations [9].
Resolving (1) results in electric fields that can then be used to calculate the corresponding orthogonal magnetic fields (). In the MT method the transfer function relating electric and magnetic fields is the complex impedance tensor, while in the GDS technique we are interested in the relationship between vertical and horizontal magnetic fields. The complex magnetic transfer function is usually presented as magnetic (real and imaginary) induction vectors using the following equations:
where and are unitary vector in the - and -axes, respectively. Note that in (4a) the direction of the real induction vector is inverted such that it conforms to the Parkinson’s induction arrows. With this convention, the real induction arrows point towards the conductive medium [10].
3. Markov Chain Monte Carlo (MCMC) Inversion
Assuming that the parameters for 1D host medium are known, the inverse problem consists in determining the conductance of the thin-sheet discretized in blocks (see Figure 1). Furthermore, possible conductances for each thin-sheet block are also available as prior information. These conductances are discrete values in a sufficiently large interval covering low to high conductance to minimize bias in the inverse problem resolution. Starting from a homogeneous thin-sheet, the conductance of the th block is updated by choosing from “a priori” values , , weighted by the probability where is the misfit related to a model in which while conductances of blocks other than th block are fixed at their current values. Assuming gaussian errors for observed data, more explicit expression for (5) is as follows: where and are the observed data and their standard deviations, while is calculated or theoretical response of the model and are data indices with being the number of data. In the case of GDS, data are the real and imaginary parts of the magnetic vector at each station for frequencies or periods used. is a normalizing constant for probabilities that involves a sum over possible conductance values.
Equation (6) represents the relative probability for , , as possible values for . We need to resolve the forward modeling times before a decision is made to select the best conductance value for each block, which is amenable for a reasonable number of , that is, 10 to 20. It is obvious that corresponding to smaller misfit will have greater probability to be selected. However, conductances with greater misfit will still have a chance to be selected; they are only less probable. This characteristic gives the MCMC inversion algorithm an ability to escape local minima frequently encountered in resolving a highly nonlinear inverse problem using linearized approach.
Updating model parameters sequentially using (6) generates a series of models (or states) obeying the fundamental Markov chain rule; that is, probability of a state depends on the past states only through the previous (or immediate past) state. Therefore, the probability in (6) serves as the transition probability of the Markov chain that governs the evolution in “time” or in iterations of the chain in a set of finite possible models. We assume without further details that key characteristics of the chain hold, for example, homogeneity or invariance and ergodicity that assure its convergence to meaningful results. After sufficiently long iterations, the Markov chain is independent of its initial state and exhibits a stationary state described by its invariant probability. It can be shown that the transition probability will be asymptotic to the invariant probability of the constructed Markov chain [11].
The iterative process described earlier is commonly called Gibbs sampler that may be employed to approximate a probability density function (PDF). It leads to the approximation of the marginal posterior probability for model parameter given the observed data, which is the solution of the inverse problems recast in the Bayesian inference context [11]. In general, empirical averages tend towards statistical averages such that statistical measures for model parameters (mean and variance) can be estimated from the evaluation of the estimated posterior quantities [12, 13].
The MCMC inversion algorithm was tested to invert synthetic data associated with thin-sheet models with satisfactory results [9]. In that test both resistive and conductive anomalies were relatively equally well resolved due to the addition of the Gauss-Seidel iterations in the misfit evaluation. The paper also underlined the importance of properly discretize the “a priori” conductance interval.
4. Inversion of GDS Field Data
The geology of the study area covering Czech, Slovakia, and parts of Poland is relatively complex since it is composed of two different geological units, that is, the Bohemian Massif and the western part of the East Carpathians. The Bohemian Massif and the Carpathians are part of the Hercynian and the Alpine chains, respectively, that are fundamental elements of the geology of Europe. Figure 2 shows a simplified tectonic map of the study area, extended to the East to cover the whole Carpathians and to the South to cover the Pannonian Basin. Much of the details about Bohemian Massif in the western part of the study area is not shown on such regional scale and simplified map. The most remarkable geological feature often related to the presence of the West Carpathian conductivity anomaly is the Pieniny Klippen belt. It is a narrow banded zone of extreme shortening and subvertical strike-slip fault zone of hundreds of kilometers long resulting from the convergence of inner Carpathians to the north area [14].
The study of the Carpathian conductivity anomaly has a very long history since its discovery deduced from geomagnetic soundings in the late 70s. Newly acquired data in recent studies still exhibit roughly similar regularities, in particular those related to the West Carpathian anomaly (WCA). Kováčiková et al. [15] employed thin-sheet inversion modeling based on the conjugate gradient minimization of misfit and a stabilizing functional. In order to describe better the sharp tectonic boundaries in this region, a minimum gradient functional was applied. This additional minimization focused the resulting image from their inversion. This approach produced blocky conductivity distribution in relatively narrow anomalous zones. Those conductive zones coincide with well-known conductivity anomalies present in the study area, mainly the Bohemian Massif anomaly (BMA) and the West Carpathian anomaly (WCA).
The GDS data set used in our study is fundamentally the same as previous studies (e.g., [15]), provided by the Institute of Geophysics, Academy of Sciences of the Czech Republic in cooperation with the Université Paris Sud, France (Menvielle, pers. comm.). The observation stations are distributed along several north to south and also north-west to south-east profiles crossing tectonic features of the area characterized by conductivity anomalies, in particular the WCA. Figure 3 shows the outline of the study area and the stations distribution along with country border and WCA outlines. The WGS84 geographical coordinate system used in Figure 3 and in all subsequent maps is for reference position only. For the thin-sheet modeling the area was then discretized into blocks with 35 by 35 km in size (see Figure 3).
The observed data are magnetic transfer function at 143 single stations for the periods of 20, 32, 64, and 98 minutes. The data processing is described in [16, 17]. Figure 4 shows the GDS data presented as real and imaginary part of the magnetic induction vectors at the period of 20 minutes. As shown in (4a) the real part vectors were inverted in accordance to Parkinson’s convention such that the vectors point towards more conductive zones. From Figure 4, qualitatively we can correlate the crossings of the real magnetic induction vectors with the prominent conductive anomalies present in the study area, especially the West Carpathian anomaly. The GDS data for other periods used in this study, that is, 32, 64, and 96 minutes (not shown), exhibit similar characteristics to those presented in Figure 4.
(a)
(b)
For the inversion, instead of interpolating the GDS data into more regular grid, for each block, we selected the datum located closest to the center of the block. In addition, we chose the most representative and coherent datum by evaluating the difference in a least-squared sense of the datum with its neighbors at the same block (if they exist). In this way we avoid introducing spurious effects that may be caused by interpolation process. However, there were blocks with no data at all that will be less constrained in the inversion process (see Figure 5).
From previous studies [11] and also from the results of Kováčiková et al. [15] we established a stratified (1D) medium for the host of a heterogeneous layer representing conductivity variation of the study area. The thin-sheet is 10 km thick at 20 km depth with a normal resistivity of 500 Ohm·m or equivalent to a conductance of 20 Siemens/m. The thin-sheet was discretized into 9 by 18 blocks of 35 by 35 km width (see Figure 3). The choice was dictated by the validity of the thin-sheet approximation; that is, the thickness of blocks should be less than their horizontal extension and also by the current available computation resources. Based on our experience with synthetic data [9], we used “a priori” conductance values from 10 up to 2000 Siemens/m discretized into quasi-logarithmic intervals to represent resistive to conductive anomalies in the region.
5. Results and Discussion
Starting from normal conductance (i.e., 20 Siemens/m) the inversion was performed up to 20 iterations. The posterior model was obtained by averaging the conductances of each block from the last 15 iterations. We present in Figure 5 the comparison between the calculated and observed data only for the real part of the induction arrows, since their physical correlation to the conductive anomaly is more obvious. In general, the fit is relatively good both for arrows’ direction and magnitude. In Figure 5, when both observed and calculated induction arrows coincide, although their magnitude might be different (in this case smaller), then it appears that only data are plotted.
The inverse model is presented in Figure 6. Despite the relatively rapid convergence of the inverse models, significant spatial variation of conductance is observed on the average model (see Figure 6(a)). The spatial discretization of the thin-sheet did not allow for a more regular conductance variation from block to block. In this case, additional constraint equivalent to the smoothness constraint as in 1D inversion [3–6] would be difficult to apply. It requires to discretize the thin-sheet into a large number of blocks resulting in blocks mostly with no data. On the other hand, the thin-sheet approximation imposes minimum block size compared to the thickness of the thin-sheet and the depth of penetration of the EM waves, such that the approximation is still valid.
(a)
(b)
In order to propose more meaningful interpretation, we rediscretized the thin-sheet conductance map obtained from the inversion into blocks with a half of their initial size, that is, approximately 17.5 by 17.5 km. The new interpolated conductance is obtained by employing a simple 5 by 5 block 2D moving average filter. The resulting model (see Figure 6(b)) has a more significant meaning as it offers the possibility of correlation with the geological features in the area. In general, our result shows that the middle and lower crust between 20 and 30 km depth are mostly conductive. We can also observe relatively narrow conductive zones (between 1000 and 1500 Siemens/m) which separates more resistant blocks (from 100 to 1000 Siemens/m).
The result confirms the existence of the well-known conductivity anomalies, that is, the Bohemian Massif anomaly (BMA) and the West Carpathian anomaly (WCA). More detailed examination of the result shows that there is a lateral shift between the surface curved features of the West Carpathians and the conductive anomaly from the model. This anomaly is interpreted rather as the effect of the root of the dipping contact between geological entities forming the West Carpathians at depth (more than 10 km). This may also reflect the fact that we determine the integrated conductivity of the thin plate, which may include the contribution of anomalies located at different depths between 20 and 30 km.
Our model also highlights conductive zones other than previously known anomalies. By considering also the qualitative interpretation of the induction vector maps, other less obvious anomalies found coincide with the tectonic features present in this area such as the Moravian-Silesian lineament (MSL) and the Labe lineament (LL) [18]. Furthermore, we can observe a significant conductivity anomaly at an area to the south of the West Carpathians. Although the inverse model in that area is not sufficiently constrained by the data, the existence of such anomaly would correspond to significant conductivity heterogeneities in the northern part of the Pannonian Basin [19]. Nevertheless, we could not interpret rigorously the significance of the conductive zones at northern and eastern parts of the study area and the resistive one to the north-east. These variations of conductance are lacking of geological evidences and are not well constrained by the observed data.
The resistivity corresponding to the conductance of the anomaly of a thin plate with thickness of 10 km is less than 10 Ohm·m. Other estimates indicate the resistivity lower than 5 Ohm·m resistivity. Several hypotheses have been proposed to explain the presence of such a conducting medium in a very resistive host (between 500 and 1000 Ohm·m) more than 20 km deep. One that is most plausible is to assume that the conductive body is sedimentary rocks saturated with saline water. The presence of saline water at high temperature in porous rocks greatly increases the conductivity of the latter. This hypothesis favoured in Jankowski et al. [20] is supported in particular by recent studies on the origin of the conductivity anomaly generally observed in the lower crust [21].
6. Conclusions
The thin-sheet modeling employing MCMC inversion algorithm has been applied to GDS data from the area covering the Bohemian Massif and the West Carpathians with satisfactory result. The conductance map obtained from such inversion modeling correlates well with the regional geology of the study area. The result allows more quantitative analysis of the conductivity anomalies present in the area rather than only qualitative interpretation based on the pattern of the induction vectors or contoured map of the magnetic transfer function values.
The MCMC algorithm applied in this study is relatively generic and would be applicable to most geophysical inverse problems with appropriate parameterization. However, the MCMC algorithm remains a computer intensive method since it has to resolve the forward modeling a substantial number of times to explore the model space. Due to limited number of iterations that can be performed with our current computation resources, the posterior quantities leading to the model uncertainty is still underexploited. With the advances in computation processing power, in the near future, it would be possible to perform more thorough analysis of the inverse model uncertainty from the posterior probability function. It would be also possible to employ a full 3D EM modeling to avoid limitation of the thin-sheet modeling in representing complex subsurface heterogeneities.
Acknowledgments
The authors greatly acknowledge Václav Červ and his colleagues from the Institute of Geophysics, Prague, Czech Republic, for providing the GDS data for this study.