Refined models of the conductivity distribution at the transition from the Bohemian Massif to the West Carpathians using stochastic MCMC thin sheet inversion of the geomagnetic induction data

Václav Červ, Michel Menvielle, Svetlana Kováčiková, Josef Pek

► To cite this version:

Václav Červ, Michel Menvielle, Svetlana Kováčiková, Josef Pek. Refined models of the conductivity distribution at the transition from the Bohemian Massif to the West Carpathians using stochastic MCMC thin sheet inversion of the geomagnetic induction data. Geophysical Journal International, 2019, 218 (3), pp.1983-2000. 10.1093/gji/ggz265. insu-02157175

HAL Id: insu-02157175

https://insu.hal.science/insu-02157175

Submitted on 14 Oct 2021

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Refined models of the conductivity distribution at the transition from the Bohemian Massif to the West Carpathians using stochastic MCMC thin sheet inversion of the geomagnetic induction data

Václav Červ<sup>D, 1</sup> Michele Menvielle,<sup>2</sup> Svetlana Kováčiková<sup>D 1</sup> and Josef Pek<sup>1</sup>

Accepted 2019 June 9. Received 2019 March 13; in original form 2018 August 29

## SUMMARY

Although volume 3-D modelling solutions has become widespread in recent time, thin sheet approximation of Earth's conductivity distribution can still serve as a useful tool when quasi-3-D conductivity structures in the heterogeneous subsurface are investigated and the available database of observations is limited to long-period electromagnetic induction data from large-scale deep sounding arrays. We present results of stochastic Monte Carlo Markov Chains (MCMC) inversion of long-period induction arrows based on the Bayesian statistics strategy.

We concentrated on the different methodological aspects of MCMC for Gibbs sampling and for adaptive Metropolis algorithm together with convergence of these methods. The results are presented on a case study from the transition zone between the Bohemian Massif and the West Carpathians where a phantom effect caused by superposition of the prominent SW-NE trending Carpathian Conductivity Anomaly and NW-SE trending anomalous structure related to the fault system at the eastern boundary of the Bohemian Massif appears.

**Key words:** Europe; Geomagnetic induction; Magnetotellurics; Inverse theory; Numerical modelling.

## 1 INTRODUCTION

In some regions, information on crustal electrical structures on a regional scale is available mainly from long-period geomagnetic transfer functions, or induction arrows. In a number of areas, targeted geomagnetic studies based on data from specially designed simultaneously recording magnetometer arrays were carried out (see reviews by Lilley 1975; Alabi 1983).

Over several past decades, the region of Central and Eastern Europe has been covered more or less uniformly by a large collection of long-period induction arrows which helped in tracing boundaries of large-scale tectonic units and blocks over the area (e.g. Praus & Pěčová 1995).

As a result of electromagnetic investigations at the eastern margin of the Variscan Bohemian Massif and its contact with the Alpine-Carpathian orogenic system, long-period induction arrows were collected extensively for more than 20 yr since 1960s, and maps of these arrows were used to trace boundaries of regional tectonic units and blocks (Praus & Pěčová 1995). Induction arrows over broad frequency ranges have still provided important interpretation parameters in a number of recent electromagnetic induction experiments in the region. Deep electrical conductivity structure at the contact zone between the Variscan Bohemian Massif and the Alpine Western Carpathians is mainly known from long-period

electromagnetic induction investigations. In the past, magnetotelluric (MT) experiments and geomagnetic depth soundings (GDS) in the area of interest were limited by technical parameters of analogous instrumentation and now, they are limited due to ubiquitous high-level industrial noise all over the area. Only recently, new MT data have been collected in the area, partly with commercial motivation (Vozár 2005; Bezák et al. 2014). Long-period geomagnetic induction data at periods of the order of thousands of seconds, with penetration depths starting at several tens of kilometres, cannot recover conductivity structure beneath the studied region in detail. Another difficulty is their low sensitivity with respect to the vertical conductivity distribution of the layered background. However, they make it possible to model regional lateral conductivity anomalies in the lower crust and upper mantle. In such conditions, a quasi-3-D thin sheet approximation of crustal conductivity structures (e.g. Vasseur & Weidelt 1977) can be a helpful modelling technique which can help to smartly avoid difficulties of time and memory-consuming full 3-D solutions. At present time, two methods of inversion of geomagnetic induction data for integrated conductivity (conductance) in a thin sheet are known. One is based on the non-linear conjugate gradient method (Wang & Lilley 2002) and the other on a Bayesian strategy with Monte Carlo Markov Chains (MCMC) used as a tool to carry out a stochastic sampling (Grandis et al. 2002, 2013; Grandis & Menvielle 2015).

<footnote n="1"><sup>1</sup>Institute of Geophysics, Czech Acad. Sci., Prague, Czech Republic. E-mail: vcv@ig.cas.cz</footnote>

<footnote n="2"><sup>2</sup>Université Versailles St-Quentin; CNRS/INSU, LATMOS-IPSL, Guyancourt, France</footnote>
In this work, we present results of a case study of the conductivity distribution at the transition zone between the Bohemian Massif and the West Carpathians using the MCMC method.

## 2 GEOLOGICAL SETTING

The transition from the Bohemian Massif to the West Carpathians is located in the junction zone of complexes originated from various tectonogeneses. The Cadomian units formed by mid-neoproterozoic to early Cambrian tectonic events (about 750–550 Ma) at the edge of the ancient Gondwana continent (Linnemann et al. 2007) and exposed only in fragments and the Variscides (Hercynides) originating from collision of the Gondwana and the Laurussia, which culminated through Late Devonian to Permian (380–280 Ma, Franke 2000) are overprinted by the Alpine orogenic phase—the product of Mesozoic (145–65 Ma) and Cenozoic (to-recent) collision between Eurasia with northward migrating African, Arabian and Indian plates). The main geological units of the former Czechoslovakia covered by the experimental electromagnetic (EM) data are displayed in a schematic map in Fig. 1. The Bohemian Massif (BM) represents the easternmost fragment of the European Hercynian (Variscan) realm occupying Bohemia and the western part of Moravia. It is built of several roughly SW-NE trending independent regional units separated by sutures detected by the deep seismic profiling results, gravity, geomagnetic and heat flow data (Suk et al. 1984). This fault pattern is intersected by the second-order NW-SE trending fault system creating a complicated blocky structure. The West Carpathians (WCP) belong to the Carpathian orogenic belt which was formed by the young Alpine folding processes and represent the northeast branch of the European Alpine orogeny. At the southwest, their boundary with the Eastern Alps is formed by the Danube Valley. Their northern boundary with the East European Platform is formed by an accretionary wedge of the Outer Carpathians flysch nappes. Apart from the Alps with their typical collisional structure, the movement in the Carpathians had an oblique character, with invasion of the terranes known as Tisza-Dacia and ALCAPA to the backarc space, docking and rotation, the fan-shaped collisional front and features characteristic for subduction observed in the Eastern Carpathians (Vrancea zone – Alasonati Tasarova et al. 2009).

The West Carpathians consist of two belts—the Inner and the Outer (Flysch) Carpathians divided by the Pieniny Klippen Belt (PKB). This very narrow subvertical stripe of intensely deformed rocks represents the deep-seated contact of the West Carpathians with the consolidated European Platform units (e.g. Stránik et al. 1993 or Hrušecký et al. 2006). The basement of both the Bohemian Massif and the West Carpathians is formed by the Cadomian complexes which appear as an independent unit known as the Brunovistulicum (BV) between them (Dudek 1980; Suk 1995) and is assumed to be a promontory of the African terrane.

North of the Carpathian Chain the transition from the Palaeozoic western terranes and the Precambrian East European Platform is formed by the Trans-European Suture Zone (TESZ) assumed to represent a Caledonian deformation front. A significant part of TESZ is represented by a system of faults of the Teisseyre-Tornquist Zone (TTZ, e.g. Grad et al. 1999). This significant boundary stretches from the Baltic Sea in the northwest to the Black Sea in the southeast hidden beneath the Eastern Carpathians and is characterized by the Moho depth increase from 30 km beneath the Variscan terranes to 45 km beneath the East European Platform, as well as by the heat

flow difference between relatively cold Precambrian Craton and the Western Europe with relatively high heat flow.

## 3 GEOMAGNETIC INDUCTION DATA

### 3.1 Experimental data

Experimental magnetic field variations data at 150 temporary stations has been recorded with analogous equipment in cooperation with Slovak and Polish colleagues in 1970s and 1980s (e.g. Praus & Pěčová 1991). The study polygon covers an area of about 500 km x 250 km in the east and north directions respectively between 15° and 23° east longitudes and 48° and 50.5° north latitudes. As a result of data processing and analysis single-station geomagnetic transfer functions have been estimated within the interval of periods of 1200–5840 s. They are components of the geomagnetic induction arrows $T_x$ and $T_y$ relating the horizontal ($H_x$, $H_y$) and vertical ($H_z$) magnetic field at each single site as functions of period and location (Wiese 1962)

$$H_z(\mathbf{r}, T) = T_x(\mathbf{r}, T) H_x(\mathbf{r}, T) + T_y(\mathbf{r}, T) H_y(\mathbf{r}, T),$$

The vertical magnetic field component at the Earth's surface is purely secondary if a uniform plane wave is assumed a source field. The vertical magnetic field thus indicates existence of lateral inhomogeneities in the crustal electrical conductivity distribution. Above a 2-D geoelectrical structure the induction arrows are perpendicular to the anomalous structure axis and according to the Wiese (1962) definition they point out from conductive zones towards a resistive medium. In this way, reversals of the induction arrows can help to map zones of anomalous induction which can be related somehow to shear zones, fault systems, important tectonic boundaries or zones of recent tectonic activation. As an example, the real and imaginary induction arrows at the period of $T = 3840$ s are presented in Fig. 2.

### 3.2 Two major zones of anomalous induction

Analysis of the induction arrows distribution made it possible to distinguish two major geomagnetic induction anomalies in the transition from the Bohemian Massif to the West Carpathians:

(1) *The West Carpathian Anomaly, WCA* – has been a well-known phenomenon since late 1950s (Wiese 1962) characterized by clear 180° reversals of the induction arrows above the central part of the conductivity anomaly. The anomaly axis runs along the Carpathians chain at the external side of the Inner and Outer Carpathians boundary (corresponding to PKB, Figs 1 and 2). The induction arrows are perpendicular to the anomaly strike and vanish above the axis what suggests a 2-D character of the WCA. The amplitudes of the arrows located to the south of the zero line are larger (by 25–50 per cent) than those located to the north of the anomaly. The presence of highly conductive structures beneath the orogeny generating the anomaly can have various explanations. Hot mineral waters in fracture systems or fluid migration in shear zones related to important tectonic boundaries can play their role in the anomaly genesis (Ádám & Pospíšil 1984; Jankowski et al. 1984, 1991; Hvoždara & Vozár, 2004, 2008). Association with metamorphosed carbon-bearing layers or graphitization zones is discussed in Glover & Vine (1992, 1995) or Žytko (1997).

(2) *The Eastern Bohemian Massif Anomaly, EBMA* – the anomalous zone on the eastern margin of the BM has a different character compared to the WCA. The induction arrows orientation
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Simplified geology of the Central Europe in the vicinity of Czech, Slovak republic and Poland.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Real and imaginary experimental induction arrows in studied area for period 3840 s; WCA, the West Carpathian Anomaly axis; EBMA, the Eastern Bohemian Massif Anomaly axis.</figcaption></figure>

is nearly or fully parallel to the SW-NE strike of the anomaly and their reversals are weakly manifested (Petr et al. 1987). The mentioned features clearly indicate a three-dimensionality of the anomaly and corresponding conductivity distribution beneath the zone. The anomaly source seems to be connected with the

processes related to the collision and submerging of the BM beneath the Carpathians although, discussions about induction processes, physical and geological mechanisms generating specific 3-D character of the anomaly still remain speculative (e.g. Kováčiková et al. 2005).
## 3.3 Previous modelling results

First estimates of the anomalies source depths using the line current approach suggested the depth of the WCA source at 18–27 km (Pěčová 1982) and similarly, the EBMA source at 17–26 km (Pěčova & Praus 1996). Jankowski et al. (1985) estimated the WCA source depth using the method of separation of fields into external and internal components and obtained the depths reaching 16 km. The analysis of equivalent current systems made it possible to locate the top of the anomaly source at the depth of about 15 km in the WCA area and at 10 km in the case of the EBMA (Červ et al. 1997).

In the 2-D joint forward modelling applied by Jankowski et al. (1991) for interpretation of both MT apparent resistivities and phases and induction arrows from the Polish part of the Carpathians, the WCA source was situated at depths less than 10 km and assumed to be related to sediments filling the Carpathian Foredeep. Another regional electrical model along the N–S profile across the High Tatra Mts in the WCP was suggested by Ernst et al. (1997) who also localized the WCA body with the resistivity of the order of 1 Ωm to depths between 8 and 16 km.

A 2-D inversion of the Carpathian data from the Polish–Slovak section of the WCA was first used by Červ & Pek (1981), and it basically confirmed the results obtained previously by fitting the data via forward simulations along N–S profiles crossing the northern segment of the WCA. The model of the conductivity distribution along the Deep Seismic Profile (DSS) No VI, traversing the BM, BV and WCP from NW to SE, based on MT results and on induction arrows, was presented in (Červ et al. 1984), and refined later in (Červ et al. 2001) by inverting the data using the 2-D REBOCC inversion procedure (Siripunvaraporn & Egbert 2000). The first rough attempt to apply the MCMC procedure to the magnetovariational data from the BM-to-WCP transition zone was introduced in Grandis et al. (2013). Our models were calculated with substantially longer Markov chains and more detailed grid of cells resulting in better convergence, less scattered mean conductances, and consequently in more precise resolution of the regional geology. We compare the latest 2-D models across the West Carpathians with our thin sheet model in Section 7 (Fig. 9c).

## 4 THIN SHEET APPROXIMATION

The subsequent technique of inversion of the experimental induction arrows for a laterally non-uniform conductance is based on a thin sheet approximation. The sheet model is a quasi-3-D approximation of the electromagnetic induction problem for sufficiently long periods of the field. This model is a suitable approximation if the conductivity anomalies may be assumed to be concentrated in a layer of a finite thickness which is much less than the penetration depth of the electromagnetic field. In such a case, the modelling of the data is weakly sensitive to the vertical changes in the conductivity distribution, and primarily aims at revealing large-scale lateral conductivity contrasts rather than at providing a complete 3-D conductivity distribution in the subsurface.

Within the context of regional electromagnetic induction studies, comprehensive theories of the 3-D thin sheet model were published in classical papers by Vasseur & Weidelt (1977), Dawson & Weaver (1979), McKirdy et al. (1985), Singer & Fainberg (1999) and others. Models consist of a laterally non-uniform thin sheet situated at the surface or buried in a 1-D normal layered model. The real 3-D conductivity σ(x, y, z) in the sheet is reduced to a 2-D distribution of the integrated conductivity, or conductance, S(x, y) = ∫°(x, y, z)dz, where the integration goes over the thickness of the sheet.

Building up on the thin sheet theory by Price (1949), Vasseur & Weidelt (1977) showed that the horizontal electric fields inside the anomalous domain of the thin sheet can be obtained by solving linear integral equations with kernels given by a product of the electric Green’s tensor of the normal model and the anomalous sheet conductance. With the anomalous electric fields known throughout the anomalous sheet domain, the electric and magnetic fields can be computed anywhere in the model by simple integrations.

Multisheet models have been also studied, for example by Robertson (1988), Fainberg et al. (1993) and Engels et al. 2002.

## 5 BAYESIAN MONTE CARLO SOLUTION OF THE INVERSE PROBLEM

As modelling of the geomagnetic induction data available in a narrow period band only can face substantial uncertainties, quantifying the conductance estimates as well as their uncertainties in a crustal sheet model may be of primary interest to studies on lateral conductance variations. In what follows, we use a Bayesian probabilistic approach to assessing the anomalous conductance distribution in a regional thin sheet model, essentially in the form presented by Grandis et al. (2002), but with more attention paid to sampling alternatives, effects of a prior selection and the convergence. Although computationally demanding for this type of a problem, the Bayesian approach shows some advantages over common deterministic minimization procedures for non-linear objective functions. In particular, estimation of parameters and their uncertainties is solved as a single problem of determining the posterior probability of the conductance parameters corresponding to the experimental magnetic response data. MCMC technique represents a computationally realizable way to simulate samples from the posterior probability distribution of model parameters conditioned on the experimental observations, which represent a most general solution to the inverse problem of the geomagnetic induction.

### 5.1 Bayesian approach in the inverse problem

The Bayesian formulation of the inverse problem allows for a joint treatment of the experimental data together with the estimation of the model parameters and their uncertainties, and determination of the posterior probability of the conductance model conditioned on the recorded data. It is governed by the Bayes theorem (see, e.g. Gelman et al. 2004, or, in the induction problem context, Grandis et al. 2002),

$$ \text{Prob}(\mathbf{S} | \mathbf{T}^{\text{obs}}, \mathcal{S}) = \frac{\text{Prob}(\mathbf{T}^{\text{obs}} | \mathbf{S}, \mathcal{S}) \times \text{Prob}(\mathbf{S} | \mathcal{S})}{\text{Prob}(\mathbf{T}^{\text{obs}}, \mathcal{S})}, \quad (1) $$

where S are the model parameters, specifically logarithms of the sheet conductances, S = (s₁, ..., sₙ) = (logS₁, ..., logSₙ), T^obs is a vector of components of the observed induction arrows, T^obs = (T₁^obs, ..., Tₘ^obs), and S symbolizes the class of models considered for the approximation of the conductivity structure. The integers N and M are, respectively, the number of model parameters (number of cells in the thin sheet model) and the total number of the data items (number of all components of the experimental induction arrows over all observation sites and all periods considered).

The posterior probability density function, Prob(S | T^obs, S), in (1) represents a solution to the inverse problem. From this function, point estimates for the model parameters can be evaluated, and their credibility intervals assessed. Common point estimators are the
posterior mean, posterior mode (maximum a posteriori, or MAP estimate), or posterior median.

## 5.2 Likelihood and prior

The first fundamental term on the right-hand side of (1), the likelihood, Prob($T^{\text{obs}} | S, S$), represents the probability of obtaining the observed induction arrows giving a particular set of values for the cell conductances. If the true conductance distribution in the sheet is $S^{\text{true}}$, and it generates data $T^{\text{true}}$, then the observations, $T^{\text{obs}}$, will be realizations of the true data disturbed by random noise. If we further assume that the noise is a Gaussian process with zero mean and with standard deviations of $\delta T_j^{\text{obs}}$, $j \in \{1, \dots, M\}$ for the individual observations, the probability density function for the observation vector $T^{\text{obs}}$ is given by

$$ \text{Prob}(T^{\text{obs}} | S^{\text{true}}, S) = \frac{1}{(2\pi)^{M/2} \prod_{j=1}^{M} \delta T_j^{\text{obs}}} \times \exp \left[ -\frac{1}{2} \sum_{j=1}^{M} \frac{(T_j^{\text{obs}} - T_j(S^{\text{true}}))^2}{(\delta T_j^{\text{obs}})^2} \right] \propto \exp \left[ -\frac{1}{2} \Phi(S^{\text{true}}) \right], \quad (2) $$

where

$$ \Phi(S^{\text{true}}) = \sum_{j=1}^{M} \frac{(T_j^{\text{obs}} - T_j(S^{\text{true}}))^2}{(\delta T_j^{\text{obs}})^2} $$

is the common expression for the data versus model $\chi^2$-misfit with residuals normalized by the data variances. If an arbitrary conductance array $S$, rather than $S^{\text{true}}$, is used in the likelihood formula (2), the likelihood function makes it possible to rank different conductance models according to their probability to produce the particular set of observations, $T^{\text{obs}}$. Other than the Gaussian noise assumption can be made in (2), for example a Laplacian biexponential noise model could be used which better captures observations with an excessive amount of outliers.

The prior probability Prob($S | S$) in the general eq. (1) represents the available, or assumed, information about the sheet cell conductances prior to the observed data, where $S$ corresponds to the assumptions made on the a priori conductivity model. In our specific case, $S$ corresponds to a class of thin sheet models embedded in a particular normal layered structure specified in Section 4.

At first, we did not consider any specific a priori information about the cell conductances, and used identical constant (flat) priors for each log-conductance. The constant prior all over the theoretically infinite range of the cell log-conductances, $s_i = \log S_i \in (-\infty, +\infty)$, $i = 1, \dots, N$, is improper since it does not integrate to 1 for any finite prior probability value. Nevertheless, the posterior density is a proper distribution if likelihood (2) is considered, and the stochastic inverse problem is then regular. Moreover, the cell-conductances are often in practice constrained to lie within some reasonable physical bounds, $S_{\min} \le s_i \le S_{\max}$, $i = 1, \dots, N$, which then implies the prior distribution to be of the following form,

$$ \text{Prob}(S | S) = \frac{1}{(s_{\max} - s_{\min})^N} \prod_{i=1}^{N} \Theta(s_i | s_{\min}, s_{\max}) \propto \prod_{i=1}^{N} \Theta(s_i | s_{\min}, s_{\max}), \quad (3) $$

where $\Theta(x | a, b)$ is a box function which is equal to one for $a < x < b$ and zero otherwise. This only slightly informative prior assumption is the same as that used by Grandis et al. (2002, 2013), and is proper and integrates to 1 over the finite conductance support. In practical computations, the lower and upper conductance bounds were mostly chosen $S_{\min} = 1 \div 10$ siemens and $S_{\max} = 8 \div 30$ ksiemens in various runs.

One of the advantages of the Bayesian strategy is that it introduces priors bearing more information about the model parameters into formula (1) if such a priori structural information is at disposal. Regularizing priors are commonly used in Bayesian inverse problems. We can, for example, favour models that are closer to the chosen normal conductance, $S^{\text{norm}}$, by assuming a priori that the cell log-conductances follow Gaussian distribution with the mean $s_i^{\text{norm}} = \log S_i^{\text{norm}}$, $i = 1, \dots, N$, and a standard deviation of $1/\sqrt{\lambda}$, where $\lambda$ characterizes the prior belief in how close, or distant, the cell conductances are from the normal conductance value adopted. Then, the prior will read

$$ \text{Prob}(S | S) = \prod_{i=1}^{N} \mathcal{N}_i (s_i^{\text{norm}}, \lambda^{-1/2}) = \left(\frac{\lambda}{2\pi}\right)^{N/2} \exp \left[ -\frac{\lambda}{2} \sum_{i=1}^{N} (s_i - s_i^{\text{norm}})^2 \right], \quad (4) $$

and may be further truncated, similarly to (3), to the physical conductance bounds, $s^{\min}$ and $s^{\max}$.

If (4) is combined with the likelihood (2), induced by the Gaussian noise model, the posterior will be

$$ \text{Prob}(S | T^{\text{obs}}, S) \propto \lambda^{N/2} \exp \left[ -\frac{1}{2} (\Phi(S) + \lambda \|S - S^{\text{norm}}\|^2) \right], \quad (5) $$

Looking at the exponent in (5), we can see that, for a fixed $\lambda$, maximizing the posterior probability is equivalent to minimizing the target $\Phi^{\text{reg}}(S) = \Phi(S) + \lambda \|S - S^{\text{norm}}\|^2$, which corresponds to the standard formulation of the Tikhonov regularized inverse problem for a minimum norm solution with respect to $S^{\text{norm}}$. Instead of (4), alternative priors can be also defined and used in the Bayesian analysis, which favour different structural features of the conductance model, for example smoothness, minimum total variation, minimum support or minimum gradient support (e.g. Portniaguine & Zhdanov 1999) of the anomalies in the inhomogeneous domain of the thin sheet.

As our prior knowledge as to the degree of the structural penalization in (5) is poor in most cases, the regularization weight, $\lambda$ in (5), is often included into the Bayesian analysis as an additional variable parameter with a vague prior (hyperparameter). An easier, but largely less general way is to adopt a fixed value of the hyperparameter $\lambda$ in (5) a priori from other considerations.

## 5.3 Sampling from the posterior by MCMC

Analytical solutions to the Bayesian inverse problem (1) are rarely possible, and are mostly limited to the simplest statistical models. Application of the Bayesian strategy to non-linear multidimensional models can be based on various approximate approaches, for example, the target probability distributions can be replaced by simpler standard probability densities (Gaussians or their mixtures) or representative samples from the posterior probability can be generated numerically (by Monte Carlo simulation procedures).

In this work, simulation of samples from the posterior probability of thin sheet conductance models based on the experimental
induction arrows is carried out using the Monte Carlo method with Markov chains. Details of the MCMC method can be found, for example, in Gelfand et al. (2004) or, practically applied in the electromagnetic induction studies, in Grandis et al. (1999, 2002). In the MCMC method, an ergodic Markov chain with the limit probability distribution equal to the target posterior probability (1) is constructed and then a particular sample from that Markov chain is obtained and processed. If stabilized and converged, it is approximately distributed according to eq. (1) and makes it possible to carry out probabilistic inference on the model parameters and their uncertainties.

### 5.3.1 Gibbs sampling

MCMC procedure serves as a frame for the generation and construction of the Markov chains. In this study, we first applied a technique suggested by Grandis et al. (2002) who used the classical Gibbs sampling procedure (see, e.g. Geman & Geman 1984; Grandis et al. 1999; Gelman et al. 2004, for a more theoretical background). The Gibbs sampler generates a Markov chain iteratively on a component-by-component basis in the following way. The procedure starts from the latest kth element of the chain ($S^{(k)}$). It successively scans the individual components of the vector $S$. For each component $s_i$, the value is updated from the following univariate conditional probability density

$$s_i^{(k+1)} \sim \text{Prob}(s_i | \underbrace{s_1^{(k+1)}, \dots, s_{i-1}^{(k+1)}}_{\text{updated values}}, \underbrace{s_{i+1}^{(k)}, \dots, s_N^{(k)}}_{\text{not updated values}}, S), \quad i = 1, \dots, N. \quad (6)$$

The whole iteration step terminates if all N components of the parameter vector have been updated. Subsequently, the Gibbs sampler transfers to its new state $S^{(k+1)}$. After the posterior probability density function is approximated by a sufficient number of chain states, marginal probability distributions of the individual model parameters can be visualized by parameter histograms. Summary statistics can be easily obtained by approximately calculating Bayesian integrals from the posterior sample, and both the model parameters and their uncertainties can be assessed by evaluating sample means, variances, quantile estimates, credible intervals, etc.

Since the conditional probabilities in (6) are not elementary distributions which could be sampled directly from, Grandis et al. (1999, 2002) treated the problem as a discrete probability problem by postulating that the cell conductances achieve values from only a set of discrete conductance levels within the admissible conductance range. This is essentially identical with the application of the Griddy-Gibbs sampler suggested by Ritter & Tanner (1992) who evaluate the intractable distribution on the right-hand side of eq. (6) at a finite series of points, form an approximate cumulative distribution function (CDF) and, finally, draw a realization of the parameter by using the inverse CDF method. The choice of the approximation grid is fully in user's hands, but it should capture the shape of the true probability distribution sufficiently well to lead to correct results. The likelihood function (2) is evaluated for each parameter at each of the grid points, what can be extremely time consuming if difficult forward solutions are involved and/or vast domains of the parameter space with low likelihood are to be sampled.

### 5.3.2 Sampling by single component adaptive metropolis (SCAM)

As an alternative to the Gibbs sampler, we have also applied the simplified component-wise adaptive Metropolis algorithm for the

multidimensional MCMC. This sampling method was suggested by Haario et al. (2003, 2005) who applied it for upper atmosphere studies. Contrary to the Gibbs sampling, individual components $s_i$ in SCAM are updated using an adaptive Metropolis rule. For this, first a proposal draw is made from a normal distribution of the parameter centred at its current value, with a data adaptive variance, specifically

$$\tilde{s}_i^{(k+1)} \sim \mathcal{N}(s_i | s_i^{(k)}, \mu (\text{Var}_i^{(k)} + \varepsilon)), \quad (7)$$

where $\text{Var}_i^{(k)}$ is the variance of the $i$th log-conductance, $s_i$, estimated from the previous steps of the sampler, the experimentally selected multiplicative constant $\mu$ optimizes the rejection-to-acceptance ratio of the Metropolis algorithm, and $\varepsilon$ is a small regularizing constant.

The Metropolis acceptance decision consists of accepting appropriate $s_i^{(k+1)} = \tilde{s}_i^{(k+1)}$, with the probability

$$\pi^{(\text{accept})} = \min \left\{ 1, \frac{\text{Prob}(\tilde{s}_i^{(k+1)} | s_1^{(k+1)}, \dots, s_{i-1}^{(k+1)}, s_{i+1}^{(k)}, \dots, s_N^{(k)}, S)}{\text{Prob}(s_i^{(k)} | s_1^{(k+1)}, \dots, s_{i-1}^{(k+1)}, s_{i+1}^{(k)}, \dots, s_N^{(k)}, S)} \right\} \quad (8)$$

while rejection of the proposal keeps the old parameter value $s_i^{(k+1)} = s_i^{(k)}$ with the probability $(1 - \pi^{\text{accept}})$. SCAM operates with longer history of the chain for adapting the variances in (7). The chain thus does not depend solely on its latest state, and evidently has not Markovian character. Nevertheless, according to Haario et al. (2003), the algorithm converges to the target posterior (1) under very general assumptions. The adaptive variance in (7) controls a quasi-optimality of the acceptance-to-rejection ratio for updating the model parameters in the chain evolution, and regulates the convergence of the chain. In our runs, we have mostly selected $\mu = 2.4$ as was suggested by Haario et al. (2003).

The adaptive Metropolis procedure needs only one forward solution per component and per iteration, and is thus cheaper per iteration than the Gibbs sampler. We have observed, however, that the autocorrelation time of the model log-conductances in the SCAM chains is generally longer than that for chains generated by the Gibbs sampler under analogous conditions. Hence, obtaining mutually independent samples from the posterior distribution by SCAM costs more iterations than by the Gibbs sampler.

# 6 STOCHASTIC INVERSION FOR THE THIN SHEET CONDUCTANCE FROM THE BM/WCP TRANSITION ZONE

## 6.1 Parameters of the thin sheet

The thin sheet approach is limited by the specific modelling conditions related to the studied periods, penetration depths of the EM field, the depth and scale of the anomalous structure and embedding medium, and, the applicability of the thin sheet theory should be checked first. Penetration depth for period 1000 s for typical continental crustal section with the resistivity $\rho = 10 \, \Omega m$ should not be less than about 50 km ($z_p = \sqrt{10\rho T}/2\pi \approx 50$ km) what is relevant to studies of anomalous structures localized in the Earth's crust. As a multisheet version of the MCMC technique is not available at present, a single sheet placed at the Earth's surface integrates all anomalous structures at crustal depths in the presented model.

The size of square cells $d$, which discretize the anomalous subdomain of the sheet, should satisfy the requirement $d\sigma\mu S_{\text{max}} \ll 1$ (Bruton 1994; Schmucker 1995). For the shortest applied period
in our study $T_{min}$ = 1200 s and $S_{max}$ = 5000 S, the size of the cells should be less than 30 km, Therefore, typical size of cells of 20–25 km in our models fully fits to the requirements mentioned above.

The embedding rocks beneath the thin sheet differ significantly in the BM region and the WCP, therefore, the simplified normal 1D section is introduced. The top 100 km occupied by a weak conductor with the resistivity varying between several hundreds and 3000 Ωm are underlain by more conductive asthenosphere with resistivities varying between 100 and 500 Ωm in our models. According to 3-D modelling results incorporating seismic data (Praus et al. 1990), a possible effect of the asthenospheric layer topography is insignificant for the period range of thousands of seconds and the resistivity exceeding 10 Ωm.

## 6.2 Inverse problem setting for the reconstruction of the thin sheet conductance

For the numerical experiments with the MCMC inversion of the long-period geomagnetic response data covering the BM-BV-WCP area, we have chosen the following particular parameters of the thin sheet and of experimental data models. They are, in principle, in accord with the earlier thin sheet model discussion in Section 4.

(i) The conductive sheet is situated at the surface of the model. The anomalous domain of the sheet contains 21(Northward) × 30 (Eastward) square tiles with uniform size 25 × 25 km², covering thus a region of about 500 (N) × 750 (E) km². This gives in total 630 variable parameters (log-conductances in the sheet cells) for the inversion. The normal conductance of the sheet outside the anomalous domain was chosen to be 400 siemens.

(ii) The normal background model consists of two layers covered by a uniform thin sheet. Parameters of the layered medium are: 100 km and 3000 Ωm for the top layer, and 300 Ωm for the resistivity of the uniform basement. The basement is used to simulate the asthenospheric layer. The topography of the asthenosphere was disregarded. The layered model is overlain by a uniform infinite thin sheet with the conductance of 400 siemens. The influence of the normal conductance on the inversion was tested for several different values, specifically 300, 400, 800 and 1500 siemens. The final average models showed very similar configurations of conductive/resistive domains throughout the region in all cases. However, the fit to the experimental data was improving if the normal conductance was increasing, especially at sites with extreme induction arrows in the Carpathian region. Here, we present results for the normal conductance of 400 siemens which was considered realistic as an average for central Europe (see 2-D models in Fig. 9). The effect of the distance of the sheet domain covered by experimental data to the external boundary of the anomalous sheet domain was also tested in the inversion experiments. In most calculations we considered the distance from the data domain to the external boundary to be about 150 km. If we extended the anomalous region by additional 150 km in all directions, the effect on the induction arrows in the data domain was practically negligible, but the computation time increased dramatically.

(iii) Before the inversion, components of the 150 experimental induction arrows were interpolated via a kriging procedure (Davis 2002) into the centres of the sheet cells. In this way, the data were also smoothed slightly and effects of inhomogeneities with higher spatial frequencies (i.e. those due to distorting structures of a subcell size) were partly reduced. As reliable data errors were not available for the original experimental transfer functions, we assumed that the

data are contaminated with Gaussian errors with zero mean and the standard deviation of 0.03. This is a reasonable choice considering the relative errors of apparent resistivities and phases as indicated in Gabás & Marcuello (2003). The resulting inversion model may a bit depend in some parts on the selection of the interpolation type. Dependence on the method of interpolation will decrease with the decreasing size of the grid.

(iv) In our long period case ($T > 1200$ s), the spatial pattern of the experimental induction arrows is very similar across the periods up to 5760 s (e.g. Červ et al. 2001). It means that the induction effect is largely due to the same conductivity distribution, and different periods sense largely the same structures at depth. We mostly used the data for the period $T = 3840$ s for various MCMC test runs to keep the computation times in reasonable limits, and only at the end we tested the final model by inverting all the three periods available, 1920, 3840 and 5760 s, simultaneously.

## 6.3 Numerical MCMC experiments for the thin sheet conductance

In our numerical experiments, we applied the Gibbs and SCAM simulation procedures to the thin sheet geomagnetic induction inverse problem, and have focused on several issues that are only marginally addressed in the pioneering studies by Grandis et al. (2002, 2013). Specifically, we try to present a comparison of the two different sampling approaches, study the convergence properties of the MCMC procedures for the thin sheet conductance problem, analyse the effect of introducing various regularizing priors on the final outputs, show the sampling-based uncertainty characteristics of the inverse solution.

### 6.3.1 Gibbs and SCAM sampling experiments

In the first experiment, we ran the Gibbs and SCAM samplers without any regularizing prior used, and then with the minimum norm regularizing prior (6), with $S^{norm} = 400$ siemens, or $s^{norm} = 2.602$. For the Griddy–Gibbs sampling, the conductance of each of the 630 variable cells could assume one of 18 predefined values from the interval of 10–8000 siemens, while the interval of the conductances allowed in the SCAM routine was slightly broader, between 1 and $10^4$ siemens.

The truncation of the normal prior (4) introduces a certain bias to the final conductance distribution. In particular, large cell conductances are eliminated, but this prevents thin sheet models with physically unfeasible parameters to appear in the solution. In our tests, we did not observe any damaging impact of the prior truncation on the MCMC convergence, which has been reported earlier, especially when the truncation region has a small probability with respect to the unconstrained posterior distribution (e.g. Kotecha & Djuric 1999).

In experiments with the regularizing prior, we did not include the regularization weight $\lambda$ into the set of variable parameters, but fixed it with the value of 1.0, which is close to the maximum of the probability (7) found from several test runs with different $\lambda \in (0.1100)$. To make the convergence monitoring meaningful, we ran the chains for several thousands of iteration cycles for both the Gibbs and SCAM algorithms. It was an extremely time-consuming procedure, as one iteration of the Gibbs sampler, with the above parameters and for one single period, needed, on average, 7–10 min to compute on a PC workstation (Intel Core i7 950 CPU, 3.07 GHz, 6 GB RAM, Intel Visual Fortran Compiler XE 12.1 on Windows
7 64 bit OS). Repeated MCMC runs were carried out to simulate
parallel chains and check the convergence of the MCMC process.

Visually, the MCMC convergence to the modus of the target pos-
terior (5) is very fast for our data set. Independently of the starting
model generated from the prior probability, the Markov chain seems
to forget its starting position and achieves values close to the max-
imum of (5) after a few iteration cycles only, typically 10–15 for
the Gibbs sampler and 30–40 for the SCAM. Because of interde-
pendencies of the model parameters (sheet cell log-conductances),
the convergence of the individual parameters is less straightforward,
and shows a more complex pattern which largely depends on the
position of the particular parameters with respect to the domain
covered by the experimental data.

In Fig. 3, we show typical traces of the log-conductances from
Markov chains, for both the Gibbs and SCAM samplers, in three
sheet cells of the model. The chains sampled from the posterior
probability density (5) with the regularizing minimum norm prior,
with $S^{\text{norm}} = 400$ siemens and $\lambda = 1$. Only single period data
were inverted at $T = 3840$ s.

The first two panels, Figs 3(a) and (b), show relatively good con-
vergence of the parameter subchains. The cell (N12, E24) in Fig. 3a
is located over a high-conductivity Carpathian anomaly and is well
surrounded by the experimental data sites (see Fig. 4 for the cell
position). The conductancce histograms are sharp at that cell, giving
a mean conductance and lower and upper bounds of the 95 per cent
(equal tailed) credible interval of 2047, 650 and 8000 siemens, re-
spectively, for the Gibbs sampler. The respective values from the
SCAM sampler are 1634, 528 and 6133 siemens. The autocorrela-
tion plots of the log-conductance in this cell, shown in the bottom
subgraphs in Fig. 3, decrease steadily, but the subsequent samples
in the chain are far from independent. The decorrelation of the sam-
ples takes a relatively long time, specifically about 30 steps for the
Gibbs sampler and 33 steps for the SCAM, measured by the sin-
gle component integrated autocorrelation time, $\tau_a = 1 + \sum_{k=1}^{\infty} \rho_k$,
where $\rho_k$ is the autocorrelation at lag $k$ (Roberts & Rosenthal 2001).

The cell (N16, E27) in Fig. 3(b) is situated far from the data
sites, in the NE sector of the model, and its conductances along
the chain show fast walks all over the conductance range allowed.
The conductance value is not well constrained by the data, the
histograms are broad and their shape is primarily formed by the prior
conductance distribution (4). The mean conductance and the lower
and upper limits of the 95 per cent credible intervals are, respectively,
153, 10 and 6500 siemens from the Gibbs sample, and 102, 3.5 and
3470 siemens from the SCAM sample. The autocorrelation function
drops quickly to zero for both samplers and then oscillates randomly
around the zero level. Estimates of the autocorrelation times make
2 and 4 steps for the Gibbs and SCAM samples, respectively.

The last example, Fig. 3(c), shows the partial chain for the cell
(N11, E17) which is situated at the margin of the data domain
and is immediately surrounded by a single data point only. The
trace plots show long-period wave motions across the whole con-
ductance range, suggesting that both low and high conductances
in that particular sheet cell are allowed with practically the same
probability to explain the data. This is due to relatively high cross-
correlations, up to ±0.4 when estimated from the sample, between
the conductances of several cells in a close vicinity of that particu-
lar cell. The histograms are broad, with the mean, lower and upper
limits of the 95 per cent credible interval being 283, 20 and 3000
siemens from the Gibbs sample, and 190, 6.4 and 4302 siemens from
the SCAM sample. The SCAM histogram even indicates some bi-
modality of the conductance distribution, with the main peak at

about 200 siemens and a side peak at about 1000 siemens. The
autocorrelations decrease only slowly, and they are well above the
zero level even at the lag of 100, especially for the SCAM sample.
The sample estimates of the autocorrelation times are 29 and 88
steps for the Gibbs and SCAM samples, respectively.

Large autocorrelations indicate dependence of the samples in
the chain, reduce the amount of information the sample provides
with regard to the parameters, and decrease thus the accuracy of
the empirical estimates evaluated from the samples in the MCMC
chain. An effective sample size (ESS), $N_{\text{eff}} = N_{\text{chain}} / \tau_a$, is used to
estimate the effective number of independent samples in a chain of
the total length $N_{\text{chain}}$. By using the function effectiveSize from the
R-package 'coda' (Plummer et al. 2006), we calculated the ESS
estimates for the individual sheet log-conductances and present
them in Fig. 4 for both the Gibbs and SCAM generated sample
chains. The ESS estimates are normalized with respect to the total
chain length here, and the plots thus show the spatial distribution
of the inverse autocorrelation times, $\tau_a^{-1}$, throughout the anoma-
lous sheet domain $\Omega_a$. As the autocorrelations are generally larger
for the SCAM sample, the scale of the respective plot (Fig. 4,
right) is reduced by a factor of four. It is clearly seen that the
ESS is largely reduced, to even less than 5 per cent of the total
chain length, in and close to domains covered by the experimen-
tal data. The extreme cases in our MCMC runs, which eventually
control the efficiency of the chain as a whole, are $\tau_a^{-1} = 0.018$ ($N_{\text{eff}} = 72$ out of $N_{\text{chain}} = 4000$) from the Gibbs generated chain,
and $\tau_a^{-1} = 0.008$ ($N_{\text{eff}} = 128$ out of $N_{\text{chain}} = 16 000$) from the
SCAM chain.

6.3.2 Formal convergence diagnostics

Monitoring the parameter traceplots and autocorrelation functions gives a basic visual impression about the MCMC chain behaviour and allows for detecting obvious non-convergence of the MCMC procedure. More formal quantitative procedures to monitor the MCMC convergence have been proposed (see, e.g. Cowles & Carlin 1996, for a critical overview), though none of them is able to give an ultimate answer to this question. The convergence of the MCMC procedure to the target probability is theoretically guaranteed for Markov chains of infinite length only. For finite chain sections, only a violation of the stability of the statistical parameters along or across the chains can be detected, that is only the non-convergence of the chain(s) can be diagnosed, in fact.

In our tests, we have used two well-known convergence diagnos-
tics criteria, specifically those by Geweke (1992) and by Gelman
& Rubin (1992), as they are easy to summarize visually even for
multivariate chains. The Geweke diagnostics tests whether a single
chain has attained an equilibrium state by comparing the means of
the sampled values from two sections of the chain, one from the
beginning (usually the first 10 per cent of the chain length) and the
other from the end (usually the last 50 per cent of the chain length)
of the chain. It produces, for each of the model parameters, a Z
statistics calculated as the difference between the two means di-
vided by the asymptotic standard error of their difference. As the
number of iterations approaches infinity, the Z statistics approaches
the normal distribution N(0, 1) if the chain has converged. For prac-
tical computations, we have used the Geweke diagnostic routines
implemented in the R-packages 'coda' (Plummer et al. 2006) and
'boa' (Smith 2007).

The Gelman–Rubin diagnostics rely on parallel chains to test
whether they all converge to the same posterior distribution.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Conductance traces (top), histograms (middle - grey histograms for normal prior) and autocorrelations (bottom) from partial Markov chains for three sheet cells of the thin sheet model. See Fig. 4 for the cell positions.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Effective sample size (ESS) relative to the total length of the chain.</figcaption></figure>

Gelman–Rubin evaluates the variance within several chains on the one hand and the variance between those chains on the other hand. By taking the ratio of those two values, it provides the so called potential scale reduction factor (PSRF). PSRF that is significantly greater than 1.0, typically greater than 1.1 or 1.2, indicates a notable difference between the chains, and suggests that longer simulation is

needed. The PSRF is monitored for each model parameter individually, but a multivariate extension (MPSRF) has been also proposed (Brooks & Gelman 1998).

We show plots of selected summary parameters along with the convergence statistics for four MCMC experiments in Fig. 5. Specifically, the columns of plots in Fig. 5 show, from top to bottom, the
<figure><img src="image_7.png" /><figcaption>Figure 5. Summary models.</figcaption></figure>

MAP and average conductance models, square roots of the sample variances of the cell log-conductances, the Geweke Z-scores and Gelman–Rubin PSRF factors obtained after 4050 iteration steps of the Gibbs sampler (50 initial iterations discarded for burn-in) and after 16 200 steps of the SCAM sampler (200 initial iterations discarded). The results are shown for sampling without any regularizing prior (except for the truncation of the conductance range, columns 1 and 3) as well as for the minimum norm prior (4) with $S_{\text{norm}} = 400$ siemens and $\lambda = 1$ (columns 2 and 4).

Contrary to standard non-stochastic minimization methods, non-regularized models from stochastic sampling do not suffer excessively from overfitting the data, since the noise effect is largely cancelled due to the marginalization of the model variables and averaging over a whole collection of individual models throughout the stochastic sample. The plots in Fig. 3, columns 1 and 3, show, however, that the convergence of the non-regularized models is poor for both the Gibbs and SCAM procedures. The Geweke’s Z-scores are outside the range (−2, 2) in more than 20 per cent of cells which is by far more than 5 per cent that would be acceptable if $Z \sim N(0, 1)$. For the SCAM sampler in particular, the non-convergent cells tend to form quasi-continuous large-scale clusters of poor convergence.

The Gelman–Rubin PSRF’s, evaluated from three parallel chains starting from largely different conductance models, are unacceptably high, indicating that the chains do not share statistically similar sampling paths in the model space even after thousands of iteration steps. The unstabilized chains produce largely different average and MAP models and, especially the latter ones, show excessive granularity.

### 6.3.3 Effect of regularizing priors

The MCMC convergence improves substantially if the regularizing prior is used, as shown in Fig. 5 in columns 2 and 4. Geweke’s Z-scores fall outside the range (−2, 2) in about 9 per cent of cases, and Gelman-Rubin PSRF from three chains exceeds the 1.2 limit in only one cell for the Gibbs sample, though the multivariate version (MPSRF) of 1.27 for the Gibbs samples and of 1.40 for the SCAM samples indicate that longer simulation runs would be still needed. MAP and average models show largely the same conductance features, except in the model domains situated close to the margin of the data region where large parameter correlations are
expected, specifically in the NW and E parts of the sheet. Large variances of the log-conductances in the non-regularized models suggest that these models are largely overparametrized, especially in the sheet areas that are far from the domain covered by the data. The regularizing prior causes the parameter estimates to be biased to some degree, but it reduces the variances evaluated from the chain samples considerably.

Essentially the same conclusions can be obtained if other than the minimum norm regularizing priors are used in (5). In Fig. 6, we present MCMC results from the SCAM sampler for three alternative priors, specifically the minimum flatness prior (column 2), minimum total variation prior (column 3) and minimum support for anomalous conductance prior (column 4). Formally, application of these priors is carried out by changing the sum in the exponent of (6) to, respectively,

$$ \sum_{i=1}^{N} \sum_{j \in V_i} (s_i - s_j)^2, \quad \sum_{i=1}^{N} \sum_{j \in V_i} |s_i - s_j|, \quad \sum_{i=1}^{N} \frac{(s_i - s^{\text{norm}})^2}{(s_i - s^{\text{norm}})^2 + \beta^2}, $$

where $V_i$ is a vicinity of the ith sheet cell consisting of its four immediate neighbours, and $\beta$ is a small regularizing factor.

All the models presented in Figs 5 and 6 show the same basic conductance features, specifically a continuous high conductance anomaly along the West Carpathians Mts., a high conductance anomaly all over the Sudetes tectonic zone in SW and S Poland, and a conductance anomaly at the eastern margin of the Bohemian Massif, which is well-developed in its southern part, close to the border with Austria, but seems to be disrupted towards the NE. The particular geometrical forms of the anomalies are less well defined, especially in areas with insufficient data coverage where the shapes of the anomalies are influenced by the prior densities to a large degree. The most diffuse anomalies are produced by the flatness prior (Fig. 6, column 2), while the minimum support prior results in a high contrast segmentation of the conductance image (Fig. 6, column 4). A few features in the models are arguable, and may be potentially deemed artifacts due to the regularizing priors or due to the particular spatial distribution of the experimental data. Specifically, this may apply to the conductive bridge connecting the Sudetes and Silesian anomalous zone with the Carpathian anomaly in S Poland, as well as to the conductive connection between the Carpathian anomaly and the anomaly at the SE margin of the Bohemian Massif in S Moravia which lacks any geological justification as a continuous feature (Červ et al. 2001).

### 6.3.4 Model uncertainty

Based on the estimated posterior probability distribution coming from the MCMC samples, the most probable parameter values as well as the corresponding uncertainty bounds can be determined. In the Bayesian analysis, the posterior mean, or MAP or median, and the posterior (co)variance, or posterior standard deviation, are commonly used to summarize the distribution of the model parameters. For models obtained from our MCMC experiments with the Gibbs and SCAM samplers and with different priors, these summary parameters are evaluated directly from the sample chains, and are presented in the maps in Figs 5 and 6.

With samples from the MCMC approximation of the posterior distribution available, Bayesian interval estimates of the model parameters are straightforward to obtain. To find an interval $(s_{\beta}^{-}, s_{\beta}^{+})$ which contains the log-conductance $s_i$ of the ith cell, $i = 1, \dots, N$, with probability $\beta \in (0, 1)$, we have to find two points within the

admissible conductabce range such that

$$ \int_{s_{\beta}^{-}}^{s_{\beta}^{+}} \mathrm{Prob}(s_i | \mathbf{T}^{\mathrm{obs}}, S) d s_i = \beta. \quad (9) $$

The Bayesian interval with probability $\beta$ is called a $\beta$-credible interval, and gives the degree of belief that the true value of $s_i$ lies within that interval. The credible interval (9) is not unique, and an additional condition has to be imposed to constrain it. For example, equi-tailed (central, symmetric) intervals are constrained by requiring that the integrals over each tail of the probability function, that is, $(s_{\min}, s_{\bar{\beta}})$ and $(s_{\beta}^{+}, s_{\max})$, are the same and equal to $(1 - \beta)/2$. The highest posterior density (HPD) intervals are also frequently reported, which are the shortest intervals, containing the most probable content, that meet condition (9).

Fig. 7(b) presents the mean model MG-MN-400 (Gibbs, average model, min. norm prior, normal conductance 400S, see Fig. 5), along with the lower and upper limits of symmetric 95 per cent-credible intervals for the cell conductances in Figs 7(a) and (c), respectively. The credible intervals are large. If induction arrows for only one period are considered ($T = 3840$ s), we can constrain the conductance to a range narrower than one order of magnitude in only 7 per cent of the sheet cells with 95 per cent credibility (Fig. 7d). The same applies to 18 per cent of sheet cells if data for three periods are inverted simultaneously, $T = 1920, 3840, 5760$ s (SCAM sampling example, Fig. 7e). Clearly, the conductances are better constrained mainly in the high conductance zones of the Carpathian anomaly and partly also in Eastern Sudetes, at the northern margin of the data domain. Inversion of data for three periods simultaneously provides a sharper conductance image, indicating a significant conductance increase also at the SE margin of the Bohemian Massif and in the north of the Moravo-Silesian zone. These results are less reliable, however, because of a fairly poor convergence of the MCMC procedure with the three-periods' data set.

Bayesian interval estimates can also help in detecting zones with significant deviations of the conductance from that of a reference model. Fig. 7(f) indicates those regions of the sheet model which differ, with high credibility of 95 per cent, from the normal conductance $S^{\text{norm}} = 400$ S, introduced as a priori conductance value through (4). It shows again a highly anomalous conductance along the Carpathians

### 6.3.5 Data fit

For an individual log-conductance model S from the Markov chain, its fitness to the experimental data is measured by the deviance parameter defined by

$$ D(S) = -2 \log \left[ \mathrm{Prob}(\mathbf{T}^{\mathrm{obs}} | \mathbf{S}, S) \right] \\ = \sum_{i=1}^{M} \left( \frac{T_i(S) - T_i^{\mathrm{obs}}}{\delta T_i^{\mathrm{obs}}} \right)^2, \quad (10) $$

which has a form of a standard $\chi^2$-misfit. A Bayesian measure of model fit is defined as the posterior expectation of the deviance, and is estimated as a mean deviance of models from the converged section of the Markov chain, $\bar{D} = K^{-1} \sum_{i=1}^{K} D(S_i)$, where K is the number of models considered from the chain. Partial deviances evaluated for individual data items or for selected data subsets can be also monitored.

For the model from Fig. 5, column 2 (model with a minimum norm prior, $S_n = 400$ S, $\lambda = 1.0$, 4000 iteration steps of the Gibbs
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Effect of various regularizing priors on the models.</figcaption></figure>

sampler), we show partial mean deviances for the individual real and imaginary induction arrows in, respectively, Figs 8(a) and (b), along with the experimental arrows and those produced by the mean conductance model. The mean deviances for the individual observation sites are relatively large especially in the eastern and central southern parts of the anomalous sheet domain where the model failed in fitting both the large induction arrows and their directional variations simultaneously. Generally, the fit is slightly better for real than for imaginary induction arrows (aggregate $\tilde{D}_{Re} = 441.0$ for all real data, while $\tilde{D}_{lm} = 487.1$ for all imaginary data; for both cases, number of the data items involved is $M_{Re} = M_{lm} = M/2 = 184$). The deviance for the mean conductance model, estimated by averaging the model parameters across the model samples in the whole Markov chain, is $D(\tilde{S}) = 867.0$ (corresponding to a standard RMS = $\sqrt{867.0/368} = 1.53$).

## 7 MEAN CONDUCTANCE MODEL AND ITS RELATION TO THE REGIONAL GEOLOGY

Interpretation of MCMC samples from narrow-band long-period geomagnetic data provides only a large-scale conductance model that can generate the recorded transfer functions. The conductance distribution calculated from the induction arrows observed at the eastern margin of the Bohemian Massif (BM) and the transition

zone between the BM and the West Carpathians (WCP) allows appropriate geological interpretation (Fig. 9).

First, the model provides detailed plausible image of the Carpathian Conductivity Anomaly (CCA) along the whole Carpathians Mts chain. The CCA is assumed to mark a collision zone of two continental plates. Extensive evidence has been brought together suggesting the enhanced conductivity is likely due to deep-seated porous sediments saturated with strongly mineralized hot waters (Jankowski *et al.* 1985, 2008 Hvoždara & Vozár 2004). Post-Oligocene migration and formation of graphite films in porous rocks and fault zones in the crystalline basement may be another reason for the anomalously high conductivity (Žytko 1997), possibly acting in parallel with the fluid mechanism (Hvoždara & Vozár 2004).

As regards the predicted properties of the anomaly in our model, the CCA reaches its maximum width and conductance in the east. This is consistent with an earlier report by Logvinov (2015) who summarized the parameters of the CCA from MT interpretations along a series of profiles across the Carpathians in eastern Slovakia and western Ukraine. He reports an increase in the conductance of the CCA from (4 – 5.5) × 10³ S to (8.5 – 15.8) × 10³ S from the W to the E between longitudes 20°E and 26°E. The conductance rise is accompanied with the CCA broadening from about 15–20 km to 30–45 km, and also with its coming up closer to the surface, from about 18 km in the W to less than 8 km in the E. In Fig. 9, we show the sketch of the principal fault systems across the studied area (Mísař *et al.* 1983) together with a comparison of the sheet
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Top row: Mean model MG-MN-400 (b), and lower (a) and upper (c) limits of the symmetric 95 per cent-credible intervals for the log-conductances in the sheet cells. Bottom row: (d) Cells with conductances constrained to a range narrower than one decade with 95 per cent credibility are indicated by full colours, other cells are shown by tints of the original scale colours. Results from the inversion of induction arrows at one period, *T* = 3840 s. (e) The same as (d), but from the inversion of induction arrows at 3 periods, *T* = 1920, 3840, 5760 s. (f) Conductances that differ with high probability of 95 per cent from the normal conductance *S*<sup>norm</sup> = 400 S. Results from the inversion of induction arrows at one period, *T* = 3840 s.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Maps of partial mean deviances for the individual real (panel a, top) and imaginary (panel b, bottom) induction arrows for the model MG-MN-400.</figcaption></figure>

conductance from our model with the integrated conductivities, down to 30 km, from 2D MT models below the PREPAN95 profile, running roughly N–S along the 22°E meridian (Ádám *et al.* 1997), along the international deep seismic profile DSSVI (Červ *et al.* 2001), profile MT-15 (Bezák *et al.* 2014) as well as below the Tatra Mts profile, deployed along the 19.7°E meridian (Ernst *et al.* 1997). The correspondence between those model results are satisfactory,

both as regards the position of the CCA and the estimates of its conductance.

Two conspicuous conducting northward protrusions in the eastern section of the CCA appear systematically with high confidence in all our models, at approximately 22°E and 21.5°E latitudes. Especially in the models with smoothing priors (Fig. 6), the latter protrusion forms an arc-like conducting zone with a SE–NW trend in southern Poland. These conductors most likely reflect the approaching TTZ from the NW, as also indicated by a distribution of the largest singular values of the reconstructed horizontal magnetic tensor across Poland by Neska *et al.* (2016). Continuation of these conductive features to the south seems to be associated with the Hornád fault zone representing the eastern boundary of the Central Western Carpathians and associated with transtension and oblique horizontal slip during recent past. Towards the west, weakening or even a local interruption in the conductivity anomaly course is observed approximately along 19.2°E (Fig. 7 right-hand panel, Figs 9a, b) at the intersection of the anomaly with the N–S trending Central Slovakia fault zone characterized by discontinuous features in other geophysical fields as well and presumably associated with strike-slip processes (Kováč & Hók 1993). Active strike-slip faults are often accompanied by enhanced conductivity due to fluid migration, although, impermeable resistive segments can exist in such fault systems, controlling fluid hydrodynamics (reference – e.g. Gleeson & Ingebritsen 2017). Further to the W/SW, the CCA closely follows the inner boundary of the Outer Carpathians along the Pieniny Klippen Belt as far as Little Carpathians Mts, west of Bratislava in SW Slovakia.

A SW-ward continuation of the CCA is not confirmed, though attempts were made in the past to trace the CCA beneath the Vienna Basin in NE Austria (Verö *et al.* 1996). Our data set does not cover
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Comparison of conductances from our thin sheet model MG-MN-400 with those obtained along four MT MV⁻¹ profiles situated within the data domain; (a) map of principal fault systems across the studied area and graphite deposits (red stars): EL, Elba Lineament; OL, Odra Lineament; CSFS, Central Slovakia Fault System; HF, Hornád Fault; PKB, Pieniny Klippen Belt; TTZ, Tornquist-Teisseyre tectonic zone; BM, Bohemian Massif; MO, Moldanubicum; SU, Sudetic zone; VB, Viena Basin; BV, Brunovistulicum with the Moravian zone, WCP, The West Carpathians; red stars – graphite deposits according to mindat.org data; (b) conductances from our thin sheet model MG-MN-400; (c) conductances from 2-D profile models within a layer 30 km thick shown by grey filled plots while the thick black lines show conductances derived from our thin sheet model. In the central map (subfigure b), geographical positions of the selected profiles are given.</figcaption></figure>
that region, and no other long period induction data are available here. Our models suggest a moderately conducting zone connecting the SW termination of the CCA and a strong conductor at the SE margin of the Moldanubian zone of the Bohemian Massif, but there is no reason for these two distinct anomalous zones to be physically interconnected. We rather suppose that this conductive ‘bridge’ may be an artifact due to the lack of data to the south of this region. Induction arrows in this area may be also affected by a conducting filling of the Vienna basin, with up to 6 km of Neogene sediments, or effects of graphitized zones which were interpreted in the upper crust at the eastern termination of the Eastern Alps in NE Austria/NW Hungary by Ádám *et al.* (2013), with model conductances of up to several thousands of Siemens.

The induction anomaly running SW-NE along the eastern margin of the Bohemian Massif (Eastern Bohemian Massif Anomaly, EBMA), and indicated primarily by rapid rotations and E-W asymmetry of the induction arrows rather than by simple arrows reversals, does not seem to be produced by a single quasi-linear conductor beneath the anomalous zone. The anomaly is clearly segmented and its segments are separated by and offset along relatively low conductance NW-SE trending strips parallel to the Elbe Zone in Moravia. A possible disruption and SE-NW dextral offset of this anomaly was already hypothesized by Petr *et al.* (1987) purely from tracing the line of close-to-zero real induction arrows. It is likely that graphite smeared within the shear zones of the Moravian thrust zone, a large scale shear zone separating the Moldanubian and Moravian zones in the east of the Bohemian Massif, essentially contributes to the enhanced conductance within this anomalous induction zone, especially in its SW part. A coincidence of the zone of the enhanced conductance with numerous graphite occurrences, including graphite deposits of industrial significance (available at http://www.mindat.org - the database run by the Hudson Institute of Mineralogy, or in Duťa & Rejl 1997, Fig. 9a), lends support to this hypothesis.

Further to the north, a strong anomalous conductivity feature in the NE part of the Bohemian Massif may represent the main source of relatively large, SW pointing induction arrows observed over the whole area to the west of the EBMA. The anomaly is located in SW Poland, and cannot be sufficiently constrained spatially due to lacking induction data in the N of that region. Electromagnetic induction measurements have been avoiding this area, mainly because of a high level of cultural noise due to heavy industry in the Silesia region. The nearest long period induction arrows are available only at latitudes greater than 51°N, more than 100 km to the N of the Czech-Polish border, and they are generally shorter than about 0.2 (at $T = 1800$ s, see, e.g. Jóźwiak 2012), indicating a general decrease of the induction effect of the TESZ towards the south. This general trend in the induction arrows also suggests that the S-SW component of the arrows in the east of the Bohemian Massif, which are up to 0.3 in magnitude, is not likely to be related to the regional induction field of the TESZ. It is more likely that the anomalous zone is related to some tectonic features of the Sudetes, which form the NE margin of the Bohemian Massif of enormous geological complexity. However, a more detailed identification of the anomalous induction sources is not possible because of a considerable spatial ambiguity in the interpretation of the position of the high conductance zone, and also because of a peripheral position of our induction data with respect to the Sudetes region.

The Sudetic region extends between two first-order WNW-ESE trending fault zones, the Middle Odra Fault Zone in the NE and the parallel Elbe Fault Zone in the SW (Aleksandrowski & Mazur 2002). The offset of the segments of the EBMA in Moravia and

sequences of low and high resistivity zones at the eastern margin of the Bohemian Massif in the models in Fig. 9 follow roughly the same direction. Though some influence of the profile arrangement of the data in this area cannot be excluded, the large-scale quasi-anisotropic structure seems to be required by the induction arrows. Such domain would allow to explain the induction processes across the whole transition zone between the BM and WCP and to generate the relatively large induction arrows over the studied region. Similar hypothesis was already mentioned earlier by Kováčiková *et al.* (2005) from results based on a unimodal induction thin sheet modelling. The physical mechanism is similar to that of a segmented mid-crustal conductor suggested by Eisel & Haak (1999) to account for persistent large long-period induction arrows in a regional vicinity of the German Continental Deep Drilling site (KTB). If the induction sources suggested above could be verified on a more detailed scale, the regional induction processes at the eastern margin of the BM may be controlled to a large degree by the NW-SE to W-E orientation of tectonic structures parallel to the Elbe Zone and further influenced by the SW-NE structures corresponding to the tectonic trend of the transition zone to the Western Carpathians.

## 8 CONCLUSION

We studied a stochastic MCMC simulation approach for a regional conductance distribution from long-period geomagnetic induction data (induction arrows). For induction data with periods longer than 1000 s the mathematical model for the crustal conductivity distribution was simplified to a single thin sheet, in the sense of Vasseur & Weidelt (1977), consisting of square tiles of variable conductance (depth integrated conductivity). The study was thus primarily focused on revealing possible lateral configurations of large-scale conducting zones that comply with the observed induction data. As a study polygon, the transition zone from the Bohemian Massif and the Polish Palaeozoic Platform to the West Carpathians in the Czech Republic, Poland and Slovakia was selected, where units of two main European orogenic cycles make contact, the Late Palaeozoic Variscan cycle and the younger, Palaeogene and Neogene, Alpine orogeny. The whole area was systematically covered by long-period geomagnetic induction measurements in the 1970s and 1980s and a collection of 150 induction arrows within the period range of 1200–5769 s was used as input data for the inversion.

MCMC procedure is used as a technique for the numerical solution of the inverse problem formulated in the Bayesian probabilistic sense. An exact probabilistic model for the experimental data distribution is practically impossible to obtain for very old data sets, and, therefore, a Gaussian model of the data error distributions with a conservative choice of the standard deviations was adopted for the likelihood function. As the inverse problem for the sheet conductances is ill-posed, the models were tested with a number of structural prior constraints. Introducing the priors has a regularizing effect on the conductance models and improves the MCMC convergence considerably, but, in some cases, structural features become set off in the regularized models, especially along the periphery of the experimental data domain, which may be, in fact, artifacts due to lacking data.

The MCMC sampling was carried out by using two standard algorithms, the Griddy-Gibbs sampler griddy92 and the SCAM algorithm by Haario *et al.* (2005). The former algorithm is generally slower than the SCAM, as it requires several forward solutions to be evaluated at each internal sweep of the sampler. However, the Gibbs procedure generates samples with shorter autocorrelation times,
providing thus less dependent samples in a shorter time. Otherwise, the sampling results are almost identical for both procedures. Both the above algorithms operate componentwise on the conductances of the sheet tiles during the internal sweeps of the sampling procedures. This allows us to speed up the internal forward solutions in the MCMC procedure by reducing the number of iterations to only a few steps in the solution of the thin sheet integral equations. It is because the current magnetic field in the sheet is always a good starting guess for the solution expected from a model with the conductance changed in a single sheet cell only.

Although the MCMC can appear time-consuming and memory-
demanding, it provides the advantage of a probabilistic solution for
the inverse problem. The method allows the assessment of model
parameters with respect to their values and uncertainty ranges, for
instance in terms of parameter histograms, credibility intervals or
sample estimates of the parameter variance–covariance matrix. A
series of stochastic simulations show that the MCMC inverse tech-
nique allows calculation of medium-size thin sheet models con-
sisting of tens of cells along each horizontal axis with standard
computer equipment, though it becomes less so if full convergence,
in the sense of meeting formal MCMC convergence criteria, is
attempted. It was one of the objectives of this study to compare
long-chain MCMC results with those of short-chain, but fast, so-
lutions as presented earlier by Grandis *et al.* (2013). Short-chain,
clearly non-converged, solutions give models with generally unac-
ceptably large scatter of mean conductances in the individual sheet
cells, though patches of well constrained conductances arise in sheet
domains that are covered sufficiently well by the experimental data.
The largest uncertainty in the mean conductance estimates appears
at the periphery of the data domain where parameter correlations
lead to longer sampling times to explore the parameter space thor-
oughly enough. We thus conclude that a short-chain inference may
be a useful tool for an approximate tectonic zonation of the area cov-
ered by the thin sheet model according to conductance variations
on a large scale.

The conductance model calculated using the MCMC method pro-
vides a plausible geological interpretation of the induction arrows
recorded at the eastern margin of the Bohemian Massif and its tran-
sition to the West Carpathians. Above the Carpathian conductivity
anomaly, the model fits in detail the recorded data features. The
anomaly observed above the eastern slopes of the Bohemian Mas-
sif suggests an alternative interpretation explained by the NW–SE
to W–E trending conductance distribution, which could reflect the
fault systems trend of the area (Odra, Sudetic and Elbe fault zones),
modulated by structures of the SW–NE tectonic style at the contact
zone between the Bohemian Massif and West Carpathians. Never-
theless, a detailed identification of the tectonic structures causative
for the induction anomalies in the studied region is, hardly feasible
by relying only on the long-period geomagnetic responses.

The final model compares well with recent results of Jóźwiak (2012) who suggested a large-scale crustal conductivity pattern in Central Europe based on a formal transformation of induction arrows, from a much larger European data collection than we have used here, into horizontal magnetic tensors. Along with the TESZ, the CCA is clearly one of the dominant conductance features in the East-Central European region. Conductors above the eastern slopes of the Bohemian Massif and in SW Poland are clearly identified in Jóźwiak's maps of the horizontal magnetic tensor invariants, though they are expressed less intensely than in our conductance model in Fig. 9. This may be partly due to a heavy smoothing of the magnetic invariant maps over areas with only coarse data coverage in Jóźwiak (2012).

ACKNOWLEDGEMENTS

We greatly thank Dr Rath and an unknown reviewer for their con-
structive comments and suggestions.

REFERENCES

Ádám, A. & Pospíšil, L., 1984. Crustal conductivity anomalies in the Carpathian region, *Acta Geodaet. Geophys. Montanist. Hung.*, **19**, 19–34.

Ádám, A. et al., 1997. Electromagnetic induction profile (PREPAN95) from the East European Platform (EEP) to the Pannonian Basin, *Acta Geodaet. Geophys. Hung.*, **32**, 203–223.

Ádám, A., Novák, A., Prácser, E., Szarka, L. & Wesztergom, V., 2013. Indication of meta-anthracite by magnetotellurics in the Koszeg-Rechnitz Penninic window: a test area, *Acta Geod. Geophys.*, **48**, 281–292.

Alabi, A.D., 1983. Magnetometer array studies, *Geophys. Surv.*, **6**, 153–172.

Alasonati Tašárová, Z., Afonso, J.C., Bielik, M., Götze, H.J. & Hók, J., 2009. The lithospheric structure of the Western Carpathian–Pannonian Basin region based on the CELEBRATION 2000 seismic experiment and gravity modelling, *Tectonophysics*, **475**(2009), 454–469.

Aleksandrowski, P. & Mazur, S., 2002. Collage tectonics in the north-
easternmost part of the Variscan Belt: the Sudetes, Bohemian Massif, in
*Palaeozoic Amalgamation of Central Europe*, pp. 237–277, eds Winch-
ester, J.A., Pharaoh, T.C. & Verniers, J., Geol. Soc. Spec. Publ., **201**.

Bezák, V., Pek, J., Vozár, J., Bielik, M. & Vozár, J., 2014. Geoelectrical and geological structure of the crust in Western Slovakia. *Stud. Geophys. Geod.*, **58**, 473–488.

Brooks, S.P. & Gelman, A., 1998. General methods for monitoring conver-
gence of iterative simulations, *J. Comput. Graph. Statist.*, **7**, 434–455.

Bruton, P., 1994. Testing Vasseur and Weidelt's thin sheet algorithm, in *Proceedings of the 15th Colloquium 'Electromagnetic Depth Research'*, pp. 54–59, eds Bahr, K. & Junge, A., DGG.

Cowles, M.K. & Carlin, B.P., 1996. Markov chain monte carlo convergence
diagnostics: a comparative review, *J. Am. Stat. Assoc.*, **91**, 883–904.

Červ, V. & Pek, J., 1981. Numerical solution of the two-dimensional inverse geomagnetic induction problem, *Stud. Geoph. Geod.*, **25**, 69–80.

Červ, V., Pek, J. & Praus, O., 1984. Model of geoelectrical anomalies in Czechoslovakia, *J. Geophys.*, **55**, 161–168.

Červ, V., Kováčiková, S., Pek, J., Pěčová, J. & Praus, O., 1997. Model of electrical conductivity distribution across Central Europe, *J. Geomag. Geoelectr.*, **49**, 1585–1600.

Červ, V., Kováčiková, S., Pek, J., Pěčová, J. & Praus, O., 2001. Geoelectrical structure across the Bohemian Massif and the transition zone to the West Carpathians, *Tectonophysics*, **332**, 201–210.

Davis, J.C., 2002. *Statistics and Data Analysis in Geology*, John Willey & Sons, 656pp.

Dawson, T.W. & Weaver, J.T., 1979. *H*-polarization induction in two thin
halfsheets, *Geophys. J. R. astr. Soc.*, **56**, 419–438.

Duda, R. & Rejl, L., 1997. *Minerals (The Great Guide)*, Aventinum (in Czech).

Dudek, A., 1980. The crystalline basement block of the Outer Carpathians in Moravia: Brunovistulicum, *Rozpr. Čs. Akad. Věd, R. mat. a př. věd*, Prague, **90**, 1–85.

Eisel, M. & Haak, V., 1999. Macro-anisotropy of the electrical conductivity of the crust: a magnetotelluric study from the German Continental Deep Drilling site (KTB), *Geophys. J. Int.*, **136**, 109–122.

Engels, M. & Korja, T., the BEAR Working Group, 2002. Multisheet mod-
elling of the electrical conductivity structure in the Fennoscandian Shield,
*Earth Planets Space*, **54**, 559–573.

Ernst, T. et al., 1997. Electromagnetic soundings across the Tatra Mountains, *Acta Geophys. Polon.*, **45**, 33–44.

Fainberg, E.B., Pankratov, O.V. & Singer, B. Sh., 1993. Thin sheet modelling of subsurface and deep inhomogeneities, *Geophys. J. Int.*, **113**, 144–154.

Franke, W., 2000. *The Mid-European Segment of the Variscides: Tectonost-
tratigraphic Units, Terrane Boundaries and Plate Tectonic Evolution.*
Geological Society, London, Special Publications, **179**.
Gabás, A. & Marcuello, A., 2003. The relative influence of different types of magnetotelluric data on joint inversions. *Earth Planets Space, **55***, 243–248.

Gelfand, A.E., Schmidt, A., Banerjee, S. & Sirmans, C.F., 2004. Nonstationary multivariate process modelling through spatially varying coregionalization. *Test **13***, 263–312.

Gelman, A. & Rubin, D.B., 1992. Inference from iterative simulation using multiple sequences, *Stat. Sci.* **7**, 457–511.

Gelman, A., Carlin, J.B., Stern, H.S. & Rubin, D.B., 2004. *Bayesian Data Analysis*, 2nd edn, Chapman & Hall/CRC, 696pp.

Geman, S. & Geman, D., 1984. Stochastic relaxation, gibbs distributions, and the bayesian restoration of images, *IEEE Trans. Pattern Anal. Mach. Intell.*, **PAMI**, **6**, 721–741.

Geweke, J., 1992. Evaluating the accuracy of sampling-based approaches to the calculation of posterior moments, in *Bayesian Statistics, Vol. 4*, pp. 169–193, eds Bernardo, J.M., Berger, J.O., Dawid, A.P. & Smith, A.F.M., Clarendon Press.

Gleeson, T. & Ingebritsen, S.E.eds, 2017. *Crustal Permeability, AGU-Willey*.

Glover, P.W.J. & Vine, F.J., 1992. Electrical conductivity of carbon-bearing granulite at raised temperatures and pressures, *Nature*, **360**, 723–726.

Glover, P.W.J. & Vine, F.J., 1995. Beyond KTB-Electrical conductivity of the deep continental crust, *Surv. Geophys.*, **16**, 5–36.

Grad, M. et al., 1999. Crustal structure of the Mid-Polish Trough beneath Teisseyre-Tornquist Zone seismic profile. *Tectonophysics*, **314**, 145–160.

Grandis, H., Menvielle, M. & Roussignol, M., 1999. Bayesian inversion with Markov chains—I. The magnetotelluric one-dimensional case. *Geophys. J. Int.*, **138**, 757–768.

Grandis, H., Menvielle, M. & Roussignol, M., 2002. Thin-sheet electromagnetic inversion modeling using Monte Carlo Markov Chain (MCMC) algorithm, *Earth Planets Space*, **54**, 511–521.

Grandis, H., Menvielle, M. & Roussignol, M., 2013. Thin-sheet inversion modeling of geomagnetic deep sounding data using MCMC algorithm, *Int. J. Geophys.*, 2013, doi:10.1155/2013/531473.

Grandis, H. & Menvielle, M., 2015. Thin-sheet electromagnetic modeling of magnetovariational data for a regional-scale study, *Earth, Planets Space*, **67**, 67–121.

Haario, H., Saksman, E. & Tamminen, J., 2003. Componentwise adaptation for MCMC, Reports of the Department of Mathematics, Univ. Helsinki, Preprint 342, 20pp.

Haario, H., Saksman, E. & Tamminen, J., 2005. Componentwise adaptation for high dimensional MCMC. *Comput. Stat.*, **20**, 265–273.

Hrušecký, I., Plašienka, D. & Pospíšil, L., 2006. Identification of the North European Platform below the Eastern Part of the Western Carpathian Flysch Belt. *AASP Memoir*, **84**, 717–727.

Hvoždara, M. & Vozár, J., 2004. Laboratory and geophysical implications for explanation of the nature of the Carpathian Conductivity Anomaly, *Acta Geophys. Polon.*, **52**, 497–508.

Jankowski, J., Petr, V., Pěčová, J. & Praus, O., 1984. Geoelectric anomaly in the Czechoslovak-Polish section of the Carpathians on the basis of geomagnetic and magnetotelluric soundings, *Acta Geodaet. Geophys. Montanist. Hung.*, **19**, 81–91.

Jankowski, J., Tarlowski, Z., Praus, O., Pěčová, J. & Petr, V., 1985. The results of deep geomagnetic soundings in the West Carpathians, *Geophys. J. R. astr. Soc.*, **80**, 561–574.

Jankowski, J., Pawliszyn, J., Jóžwiak, W. & Ernst, T., 1991. Synthesis of electric conductivity surveys performed on the Polish part of the Carpathians with geomagnetic and magnetotelluric sounding methods, *Publs. Inst. Geophys. Pol. Acad. Sci.*, A-19(236), 183–214.

Jankowski, J., Jóžwiak, W. & Vozár, J., 2008. Arguments for ionic nature of the Carpathian electric conductivity anomaly, *Acta Geophys.*, **56**, 455–465.

Jóžwiak, W., 2012. Large-scale crustal conductivity pattern in Central Europe and its correlation to deep tectonic structures, *Pure appl. Geophys.*, **169**, 1737–1747.

Kováč, P. & Hók, J., 1993. The Central Slovak Fault System—the field evidence of a strike slip, *Geol. Carpath.*, **44**, 153–159.

Kotecha, H.J. & Djuric, P., 1999. Gibbs sampling approach for generation of truncated multivariate Gaussian random variables. *ICASSP, IEEE International Conference on Acoustics, Speech and Signal Processing - Proceedings*, **3**, 1757–1760.

Kováčiková, S., Červ, V. & Praus, O., 2005. Modelling of the conductance distribution at the eastern margin of the European Hercynides, *Studia geoph. geod.*, **49**, 403–421.

Lilley, F.E.M., 1975. Magnetometer array studies: a review of the interpretation of observed fields, *Phys. Earth planet. Inter.*, **10**, 231–240.

Linnemann, U., Gerdes, A., Drost, K. & Buschmann, B., 2007. The continuum between Cadomian orogenesis and opening of the Rheic Ocean: constraints from LA-ICP-MS U-Pb zircon dating and analysis of plate-tectonic setting (Saxo-Thuringian zone, northeastern Bohemian Massif, Germany). *Geol. Soc. Am. Spec. Paper*, **423**, 61–96.

Logvinov, I.M., 2015. Deep geoelectrical structure of the Central and Western Ukraine, *Acta Geophys.*, **63**, 1216–1230.

McKirdy, D.McA., Weaver, J.T. & Dawson, T.W., 1985. Induction in a thin sheet of variable conductance at the surface of stratified Earth-II. Three-dimensional theory, *Geophys. J. R. astr. Soc.*, **80**, 177–194.

Mísář, Z., Dudek, A., Havlena, V & Weiss, J. 1983. Geologie ČSSR, Český Masiv, SPN Praha: 226 – 236 (in Czech).

Neska, A., Jóžwiak, W. & Ladanivskyy, B., 2016. Crustal structure of the southern Polish Basin imaged by magnetotelluric surveys, in *Proceedings of the 26th Schmucker-Weidelt Colloquium Electromagnetic Depth Research*, pp. 52–62, eds Miensopust, M. & Becken, M., DGG.

Pěčová, J., 1982. Determination of the depth of a conductivitz anomalz bz separating the geomagnetic variation field into its external and internal part, *Travaux Geophys.*, **575**, 179–197.

Pěčová, J. & Praus, O., 1996. Anomalous induction zones in the Czech Republic in relation to large scale European anomalies, *Stud. Geoph. Geod.*, **40**, 50–76.

Petr, V., Pěčová, J., Praus, O. & Pěč, K., 1987. Anomalous induction zone near the eastern margin of the Bohemian Massif, *Phys. Earth planet. Inter.*, **45**, 161–169.

Plummer, M., Best, N., Cowles, K. & Vines, K., 2006. CODA: convergence diagnosis and output analysis for MCMC, *R News*, **6**, 7–11.

Portniaguine, O. & Zhdanov, M.S., 1999. Focusing geophysical inversion images, *Geophysics*, **64**, 874–887.

Praus, O., Pěčová, J., Petr, V., Babuška, V. & Plomerová, J., 1990. Magnetotelluric and seismological determination of the lithosphere-asthenosphere transition in Central Europe, *Phys. Earth planet. Inter.*, **60**, 212–228.

Praus, O. & Pěčová, J., 1991. Anomalous geomagnetic fields of internal origin in Czechoslovakia, *Stud. Geoph. Geod.*, **35**, 81–89.

Praus, O. & Pěčová, J., 1995. Anomalous induction across Europe, *Stud. Geoph. Geod.*, **39**, 286–288.

Price, A.T., 1949. The induction of electric currents in non-uniform thin sheets and shells, *Q. J. Mech. appl. Math.*, **2**, 283–310.

Ritter, Ch. & Tanner, M.A., 1992. Facilitating the Gibbs sampler: the Gibbs stopper and the Griddy-Gibbs sampler, *J. Am. Stat. Assoc.*, **87**, 861–868.

Roberts, G.O. & Rosenthal, J.S., 2001. Optimal scaling for various Metropolis-Hastings algorithms, *Statist. Sci.*, **16**, 351–367.

Robertson, R.C., 1988. The electromagnetic response of a heterogeneous layer modeled by two thin sheets in a uniformly conducting half-space, *IEEE Trans. Geosci.*, **26**, 2–10.

Schmucker, U., 1995. Electromagnetic induction in thin sheets: integral equations and model studies in two dimensions, *Geophys. J. Int.*, **121**, 173–190.

Singer, B.Sh. & Fainberg, E.B., 1999. Modeling electromagnetic fields in thin heterogeneous layers with application to field generation by volcanos, I—Theory and example, *Geophys. J. Int.*, **138**, 125–145.

Siripunvaraporn, W. & Egbert, G., 2000. An efficient data-subspace inversion method for 2-D magnetotelluric data, *Geophysics*, **65**, 791–803.

Smith, G.J., 2007. boa: An R package for MCMC output convergence assessment and posterior inference, *J. Stat. Soft.*, **21**, 1–37.
Stráník, Z., Dvořák, J., Krejčí, O., Muller, P., Přichystal, A., Suk, M. & Tomek, Č., 1993. The Contact of the North European Epivariscan Platform with the West Carpathians. *J. Czech. Geol. Soc.*, **38**(1-2), 21–29.

Suk, M. et al., 1984. *Geological History of the Czech Socialist Republic*, Academia, Geological Survey, Prague, 387p.

Suk, M., 1995. Regional geological division of the Bohemian Massif, *Explor. Geophys.*, *Rem. Sens. Environ.*, **II**(2), 25–30.

Vasseur, G. & Weidelt, P., 1977. Bimodal EM induction in non-uniform thin sheets with an application to the Northern Pyrenean Anomaly, *Geophys. J. R. astr. Soc.*, **51**, 669–690.

Verö, L., Madarasi, A., Seiberl, W. & Varga, G., 1996. Magnetotelluric tracing of the Carpathian conductivity anomaly in the Vienna Basin, *Proceedings of the 58th EAGE Conference and Exhibition Session: Electromagnetics II*, doi:10.3997/2214-4609.201408667.

Vozár, J., 2005. Magnetotelluric Soundings on the Territory of Slovakia, *PhD thesis*, GFÚ SAV Bratislava, 108 pp. (in Slovak).

Wang, L.J. & Lilley, F.E.M., 2002. Inversion of magnetometer array data by thin-sheet modelling, *Geophys. J. Int.*, **137**, 128–138.

Wiese, H., 1962. Geomagnetische Tiefentellurik, Teil II: Die Streichrichtung der Untergrundstrukturen des elektrischen Widerstandes, erschlossen aus geomagnetischen Variationen [Strike direction of underground structures of electric resistivity, inferred from geomagnetic variations], *Pure appl. Geophys.*, **52**, 83–103 (in German).

Žytko, K., 1997, Electrical conductivity anomaly of the Northern Carpathians and the deep structure of the orogen, *Ann. Soc. Geol. Polon.*, **67**, 25–43.