A Bayesian Inference Approach to Reduce Uncertainty in Magnetotelluric Inversion: A Synthetic Case Study

Show more

1. Introduction

The basic principle behind any geophysical inversion is to explore the subsurface properties of geological structure. This involves measuring the geophysical properties at earth’s surface and subsequently inverting the data in order to reveal the subsurface properties. Several geophysical techniques are used in this process, among which the magnetotelluric method has been applied extensively in many fields. The magnetotelluric is a natural source method where both magnetic and electric fields are measured in orthogonal directions giving two impedance tensors and their phase. Later, the measured data is processed and inverted to derive the subsurface resistivity hence revealing the subsurface physical properties of the earth.

In quantifying the geophysical inverse problem, an initial model is assumed by means of forward modeling process. Forward modeling allows us to predict possible future values of the field’s measurable outputs and is calculated from the prior knowledge of what we know about the geological and geophysical structure under investigation. During the inversion process, the assumed initial model result is compared to the real observed measured data of the geological surface whereby the best fit model is chosen in the end. However, due to problems such as random noise, over-parameterization of the subsurface models, incomplete and sparse data collection, all model solutions are subjected to errors, under-determination and uncertainty. Hence, all other inversion techniques in the field of science and engineering suffer from a common problem where multiple models or descriptions can be obtained that fit or interpret equally very well the observed data (Sen & Stoffa, 1996) .

The deterministic geophysical inversion methods are also dominant (Chen, Hoversten, & Nordquist, 2010) when inverting magnetotelluric data and thus the results largely depends on the assumed initial model and usually, only a single representative solution is obtained. Therefore, the results provide limited statistical information about the solution and a non-uniqueness solution arises. A possible solution to this problem is to utilize Bayesian inversion. In this method, the inversion parameters are expressed in the form of probability density functions and are regarded as random variables. In addition, the uncertainty of the earth model parameters can be included in the form of prior distributions and quantified by examining the posterior distribution model. Normally, this task is achieved by utilizing Markov Chain Monte Carlo (MCMC) methods (Conway, Simpson, Didana, Rugari, & Heinson, 2018) .

Therefore, in this paper, we present the use of stochastic method for inverting magnetotelluric data. Based on MCMC method, a one-dimension Bayesian inversion code was developed for a simple layered earth model and was tested using synthetic model data. A number of studies have applied 1D Bayesian inversion of magnetotelluric data. For example, Grandis, Menvielle and Roussignol (1999) solved the Bayesian magnetotelluric inverse problem using MCMC in layered situations by dividing the domain under study into homogeneous layers. In addition, Guo, Dosso, Liu, Dettmer and Tong (2011) applied a 1D magnetotelluric Bayesian approach to examine non-linearity for the inverse problem. More recently, Conway et al. (2018) successfully applied 1D MCMC magnetotelluric inversion using a No-U-Turns Sampler. However, Mandolesi, Ogaya, Campany and Piana (2018) attest that magnetotelluric data are sensitive to the depth-distribution of rock resistivity. Hence, apart from providing the possible resistivity posterior distribution, we have also incorporated the estimation of depth posterior distribution in our code in order to reduce its uncertainty. Furthermore, the developed Bayesian inversion code has low dependency on the initial model eliminate non-uniqueness problem. The sampling of the posterior distribution was achieved using a Metropolis-Hastings algorithm.

2. Method and Principle

2.1. Electromagnetic and Magnetotelluric Theory

The electromagnetic methods are an excellent means of mapping the subsurface by delineating the electrical conductivity and resistivity of earth materials. The fundamental laws of electromagnetic fields governing geophysical applications have been discussed by many authors (see Vozoff, 1991; Ward & Hohmann, 1988; Zhdanov, 2009 ) and they are best described using the Maxwell’s equations. These equations govern all electromagnetic phenomena for any conductive medium, and are describe as:

$\nabla \times E=-\frac{\partial B}{\partial t}$ (Faraday’s Law) (1)

$\nabla \times H=J+\frac{\partial D}{\partial t}$ (Ampere’s Law) (2)

$\nabla \cdot D=q$ (Gauss’ Law for Electricity) (3)

$\nabla \cdot B=0$ (Gauss’ Law for Magnetism) (4)

where the four field parameters are: E is electric field intensity in (V/m), D is electric flux density (C/m^{2}), H is magnetic field intensity (A/m), B is magnetic flux density (T or Wb/m^{2}), J is the electric current density (A/m^{2}), q is electric charge density (C) and t is the time (s). In general, Equation (1) is referred to as Faraday’s law of electromagnetic induction, Equation (2) is the Ampere’s law while Equation (3) and Equation (4) is the Gauss’ law for electricity and Gauss’ law concerning magnetism respectively.

Nevertheless, in order to reduce the number of basic vector field functions for a linear isotropic medium, the Maxwell’s equations can sometimes be coupled by the empirical constitutive relations as:

$B=\mu H$ (5)

$D=\epsilon E$ (6)

$J=\sigma E$ (7)

where $\mu $ is magnetic permeability (Vs/Am), $\epsilon $ is dielectric permittivity (As/Vm) and $\sigma $ is electrical conductivity (S/m). The magnetic permeability and dielectric permittivity are further defined as:

$\mu ={\mu}_{r}{\mu}_{0}$ (8)

$\epsilon ={\epsilon}_{r}{\epsilon}_{0}$ (9)

${\mu}_{r}$ is the relative magnetic permeability and ${\mu}_{0}$ is the magnetic permeability of the free space whose value is $4\text{\pi}\times {10}^{-7}$ H/m. ${\epsilon}_{r}$ is the relative electric permittivity and ${\epsilon}_{0}$ is the electric permittivity of the free space whose value is $8.85\times {10}^{-12}$ F/m. The parameters $\mu $ , $\epsilon $ and $\sigma $ vary with position in the earth, describing the physical properties of the material through which the EM wave is propagating (Rosenkjaer, 2011) . These parameters, including frequency f and resistivity $\rho $ (reciprocal of $\sigma $ ), will also determine the depth of penetration of EM field through Earth materials. The depth of penetration, called the skin depth ( $\delta $ ) is given by:

$\delta =\sqrt{\frac{2}{\omega \sigma \mu}}=503\sqrt{f\sigma}=503\sqrt{\frac{\rho}{f}}$ (10)

where $\omega $ is the angular frequency. Since $\rho $ is proportional to and f is inversely proportional to $\delta $ , it follows that deeper depth of investigation will be achieved with increase in $\rho $ and a decrease in f.

Magnetotelluric is a passive electromagnetic exploration technique for determining the electrical conductivity of earth geological structure by utilizing the natural occurring broad range spectrum of electromagnetic fields. This is achieved by measuring the secondary orthogonal electric and magnetic fields of the earth structure caused by a change in the electromagnetic field from the naturally occurring primary source. The information about the resistivity (or impedance $Z$ ) profile of the earth is now calculated from the ratio of electric and magnetic fields. The orthogonal components of the horizontal $E$ and $H$ fields are related via a complex impedance tensor which has the form:

$Z=\frac{E}{H}$ or $\left(\begin{array}{c}{E}_{x}\\ {E}_{y}\end{array}\right)=\left(\begin{array}{c}{Z}_{xx}{Z}_{xy}\\ {Z}_{yx}{Z}_{yy}\end{array}\right)\left(\begin{array}{c}{H}_{x}\\ {H}_{y}\end{array}\right)$ (11)

It is not easy to visualize impedance tensors and thus, a useful and widely accepted tool is to transform these values into apparent resistivity ( ${\rho}_{a}$ ) and phase ( $\varnothing $ ) respectively:

${\rho}_{a}=\frac{1}{\mu \omega}{\left|Z\right|}^{2}$ (12)

$\varnothing =\mathrm{arctan}\left(\frac{im\left[Z\right]}{re\left(\left[Z\right]\right)}\right)$ (13)

where im and re represents the imaginary and real part of the impedance elements respectively. Equation (12) is referred to as Cagniard resistivity, who laid the foundation of magnetotelluric theory in Cagniard (1953) .

2.2. Bayesian Inversion Method

Bayesian inversion is a probabilistic method for finding model parameters. The Bayes’ rule or Bayes’ theorem is derived from a conditional probability of A given that B has occurred which is expressed as:

$P\left(A|B\right)=\frac{P\left(B|A\right)\times P\left(A\right)}{P\left(B\right)}$ (14)

We can substitute A and B in Equation (14) in terms of model parameters (m) and observed data vectors (d) respectively as:

$P\left(m|d\right)=\frac{P\left(d|m\right)\times P\left(m\right)}{P\left(d\right)}$ (15)

where $P\left(m|d\right)$ is called the posterior distribution function, $P\left(d|m\right)$ is called the model likelihood function, $P\left(m\right)$ is the regularization referred to as the prior distribution and finally, $P\left(d\right)$ is the marginal probability of the measurements sometimes referred to as the evidence. Usually, $P\left(d\right)$ acts as a normalizing constant over all possible models to ensure that the posterior distribution function $P\left(m|d\right)$ integrates to 1 (Ray & Key, 2012) . Therefore, most of the times $P\left(d\right)$ can be ignored which results into the following relationship:

$P\left(m|d\right)\propto P\left(d|m\right)\times P\left(m\right)$ (16)

The likelihood function $P\left(d|m\right)$ , often written as $L\left(m\right)$ , is used to measure the difference between a set of data and a proposed model $m$ . Therefore, $L\left(m\right)$ depends on the statistics of the Gaussian noise distribution, of which its value depends on the model $m$ being sampled and its misfit (Ray, Alumbaugh, Hoversten, & Key, 2013) given by:

$L\left(m\right)\propto \mathrm{exp}\left(\frac{-{\left[d-f\left(m\right)\right]}^{\text{T}}{C}_{d}^{-1}\left[d-f\left(m\right)\right]}{2}\right)$ (17)

where
${C}_{d}^{-1}$ is the data covariance matrix and
${\left[d-f\left(m\right)\right]}^{\text{T}}{C}_{d}^{-1}\left[d-f\left(m\right)\right]$ is the Chi-square (X^{2}) misfit for the evaluated model
$m$ . A small misfit between the modeled response
$f\left(m\right)$ and the real observed data
$d$ results in a higher likelihood for a model
$m$ than for that with a large misfit (Buland & Kolbjørnsen, 2012) . Obtaining a large likelihood value indicates how closer the model response is to the experimental data and thus the goal in Bayesian inversion.

2.3. Markov Chain Monte Carlo: Metropolis-Hastings Algorithm

MCMC are sampling techniques that have been successfully applied in many fields including geophysical inverse problems. The integration schemes in MCMC makes use of pseudo random number generators: the integrand is evaluated at points chosen uniformly at random (Sen & Stoffa, 1996) . One of the most widely used sampling method in the Bayesian inversion approach is the Metropolis-Hastings (M-H) algorithm. The M-H is an MCMC simulation method in which a ratio of posterior distributions is used to move through the solution space (creating a Markov chain) to generate random samples. In the constructed Markov chain, the M-H algorithm helps in sampling the best candidate from the unknown distribution by evaluating the posterior distributions ratio of two states and hence forming a new model. In order to determine if the new model (m_{2}) or the old model (m_{1}) is the best, an acceptance probability
$\alpha $ is calculated from:

$\alpha \left({m}_{2}|{m}_{1}\right)=\mathrm{min}\left[1,\frac{p\left({m}_{2}\right)}{p\left({m}_{1}\right)}\times \frac{p\left(d|{m}_{2}\right)}{p\left(d|{m}_{1}\right)}\times \frac{q\left({m}_{1}|{m}_{2}\right)}{q\left({m}_{2}|{m}_{1}\right)}\right]$ (18)

where $\frac{p\left({m}_{2}\right)}{p\left({m}_{1}\right)}$ is the prior ratio, $\frac{p\left(d|{m}_{2}\right)}{p\left(d|{m}_{1}\right)}$ is the likelihood ratio and $\frac{q\left({m}_{1}|{m}_{2}\right)}{q\left({m}_{2}|{m}_{1}\right)}$ is the proposal ratio. It is important to note that the ratio of the

likelihood function of the models is referenced whenever the new model is examined. For a given set of impedances, the likelihood function is calculated from Chi-square (X^{2}), which for a one dimension case is given by:

${X}^{2}={\displaystyle {\sum}_{i}\frac{\left({Z}_{i}-{\stackrel{\xaf}{Z}}_{i}\right)}{2{\sigma}_{i}^{2}}}$ (19)

The new model is better than the old one when
$\alpha \ge 1$ , and thus is always accepted. However, sometimes we would like to accept the rejected models. This is an important feature of the MCMC methods and ensures that there is a full sampling of the parameter space. For
$\alpha <1$ , a random number is selected between [0, 1] and if that number lies within the interval
$\left[\alpha ,1\right]$ then the new model is accepted, else the old model is retained (Ray et al., 2013) . This process is repeated for a large number of iterations. Over times, the new state moves towards a better solution: one for which the misfit between data and forward models is small. After sufficient warm-up period, the new states are collected together to provide a distribution of the parameter space (Malinverno, 2002) . Since consecutive states will be highly correlated, only every N^{th} state is used in the final result.

2.4. Bayesian Inversion Algorithm Flow

The pseudo code for the Bayesian inversion algorithm flow we developed is shown in Figure 1. We used MATLAB as a platform for this pseudo code. The core of this algorithm is the LOOPS part as is shown in the dotted box. Generally, the Bayesian inversion algorithm can be summarized in four main parts which are: 1) create a new model with Jump; 2) evaluate the MT response with Forward code; 3) calculate the likelihood probability by evaluating chi-square; and finally 4) make a decision using the M-H for which model is best. Thereafter, the result from the M-H is compared with all different kinds of input data and prior information before the final output result is given after all iterations.

3. Results and Discussion

In order to test the effectiveness of the developed 1D Bayesian inversion magnetotelluric pseudo code, we adopted a three layer model to carry out a synthetic

Figure 1. Bayesian inversion algorithm flow.

study. We chose different models namely H, K, A and Q type ones to perform the synthetic test whose resistivity values were taken from the ones available in literature. The generated frequency bandwidth for the magnetotelluric impedance data ranged from 0.01 to 100 Hz. A 5% Gaussian noise was added at each frequency in order to simulate errors to our synthetic results according to (Conway et al., 2018; Grandis et al., 1999; Mandolesi et al., 2018) . We set the number of iterations to 30,000 times for each model. Taking into account the problem of the convergence rate, samples for each model were drawn starting from 15,000 and every 10^{th} sample was used to avoid correlation between consecutive models.

3.1. H-Type Model ( $\rho 1>\rho 2<\rho 3$ )

First, we examined the response of the H-type model to our pseudo code. An H-type model has a bowl shaped curve whereby the middle layer resistivity is less than the top and the third layer resistivities. The profile composition of this model was adopted from Mohammad, Srigutomo, Sutarno and Sumintadiredja (2013) with a 300 m thick top layer at 500 Ω.m, middle layer thickness is 700 m at 5 Ω.m and 50 Ω.m for the third layer while its thickness was set at infinity. The distribution of resistivity values of the model parameters is shown in Figure 2 while Figure 3 shows the distribution values of thickness with the correct synthetic layer resistivity and thickness value displayed with a horizontal red line.

It can be observed from the results that the distribution of the thickness after 30,000 iterations were near to the initial set real model values. This was also the same case with resistivity distribution values. The convergence of the resistivities values to the initial synthetic true value for all layers starting from top to bottom is also shown in Figure 4 and we can clearly observe that the third layer converged

Figure 2. Resistivity distribution of H-type model.

Figure 3. Thickness distribution of H-type model.

Figure 4. Convergence trend of H-type model.

to an average resistivity value of 50 Ω.m after 500 samples and hence the results should improve with more iterations. This H-type model is a typical example of geothermal site with a conductive layer between two resistive layers (Grandis & Sumintadireja, 2017; Mohammad et al., 2013) . The inversion results we have obtained indicated a fairly good compatibility with the initial synthetic model and thus it can be applied to geothermal explorations.

3.2. K-Type Model ( $\rho 1<\rho 2>\rho 3$ )

Second, we looked at the performance of the K-type model with our pseudo code algorithm. A K-type model is the opposite of H-type model whereby its middle layer resistivity is high than the top and the third layer resistivity’s hence its curve is bell shaped. We set the top layer thickness at 1000 m and 500 Ω.m, middle layer thickness is 100 m at 3000 Ω.m while the bottom layer resistivity is 100 Ω.m according to Mohammad et al. (2013) . The resistivity and thickness distribution values of the model parameters are shown in Figure 5 and Figure 6 respectively, again with the correct synthetic values displayed with a horizontal red line.

The results for both average resistivity and thickness value after 30,000 iterations showed that they were closer to the synthetic true model values for all layers apart from the second layer resistivity. This is also depicted in Figure 7 where it shows the convergence of the resistivities values to the set synthetic true value for all layers with the second layer averaging around 2990 Ω.m. According to Mohammad et al. (2013) , this K-type layered earth model is a typical example of a hydrocarbon site whereby a very thin resistive layer (100 m) is embedded between relatively conductive layers. Similar to the H-type result, the general observation from the inversion results is that it can be applied to hydrocarbon explorations as the results obtained indicated good compatibility with the initial synthetic model data.

3.3. A-Type Model ( $\rho 1<\rho 2<\rho 3$ )

Third, we examined the response of A-type model which has a continuous increasing resistivities with depth. We set the synthetic model according to Shireesha & Harinarayana (2013) in which the top layer thickness is 500 m at 10 Ω.m, middle layer thickness is 1000 m at 20 Ω.m while the bottom layer resistivity is 1000 Ω.m. The inversion results for the resistivity and thickness distribution values of the model parameters are given in Figure 8 and Figure 9 respectively.

Figure 5. Resistivity distribution of K-type model.

Figure 6. Thickness distribution of K-type model.

Figure 7. Convergence trend of K-type model.

Figure 8. Resistivity distribution of A-type model.

Figure 9. Thickness distribution of A-type model.

It can observed (Figure 8 and Figure 9) that the distribution peaks were centered near the true synthetic value for both resistivity and thickness layered earth. This is also depicted in Figure 10, which shows the convergence of the resistivities values to the synthetic true value for all layers starting from first to third layer. Overall, an A-type curve is usually obtained in hard rock with thin conductive top soils and our synthetic model analysis indicates good compatibility results.

3.4. Q-Type Model ( $\rho 1>\rho 2>\rho 3$ )

Finally, we examined the performance response of Q-type model according to Shireesha & Harinarayana (2013) example. Q-type model has a continuous decreasing resistivities with depth. The top layer thickness of the synthetic model was 500 m at 1000 Ω.m, middle layer thickness is 1000 m at 100 Ω.m while the bottom layer resistivity is 10 Ω.m. The resistivity distribution values and thickness distribution values of the model parameters are shown in Figure 11 and Figure 12 respectively, again with the correct synthetic layer value displayed in horizontal red line.

The distribution peaks for both resistivity and thickness layered earth models were observed to be centered near the true synthetic values. This is also depicted in Figure 13, which shows uniform convergence of the resistivities values to the initial synthetic true value for all layers starting from first to third layer. Q-type curves are commonly obtained in coastal region where a compact top layer is underlain by a thick clay layer or saline water aquifer and is also called conductive basement.

4. Conclusion

In this paper, we developed a pseudo code algorithm to invert magnetotelluric data based on Bayesian inference framework. As a test, we applied the inversion

Figure 10. Convergence Trend of H-type model.

Figure 11. Resistivity distribution of Q-type model.

Figure 12. Thickness distribution of Q-type model.

Figure 13. Convergence trend of Q-type model.

algorithm to synthetic model data obtained from available literature which was based on a three layer model (K, H, A and Q-type models). The results were given as a model distribution of possible resistivity values for each and every layer rather than a linear model. This helps in solving the non-uniqueness and eliminating uncertainty problems in geophysical inversion. The developed algorithm has been successfully applied to all types of models and results obtained have demonstrated a good compatibility with the initial synthetic model data. Thus, the developed algorithm has shown promising results for future work in magnetotelluric Bayesian inversion as it can be used to accurately estimate the resistivity uncertainty of subsurface geological structures. However, addition work should be performed to assess how the developed algorithm will perform using observed real field data.

References

[1] Buland, A., & Kolbjørnsen, O. (2012). Bayesian Inversion of CSEM and Magnetotelluric Data. Geophysics, 77, E33-E42.

https://doi.org/10.1190/geo2010-0298.1

[2] Cagniard, L. (1953). Basic Theory of the Magnetotelluruc Method of Geophysical Prospecting. Geophysics, 18, 605-635.

https://doi.org/10.1190/1.1437915

[3] Chen, J., Hoversten, G. M., & Nordquist, G. (2010). Stochastic Inversion of 2D Magnetotelluric Data Using Sharp Boundary Parameterization. Denver: SEG 2010 Annual Meeting.

[4] Conway, D., Simpson, J., Didana, Y., Rugari, J., & Heinson, G. (2018). Probabilistic Magnetotelluric Inversion with Adaptive Regularisation Using the No-U-Turns Sampler. Pure and Applied Geophysics.

https://doi.org/10.1007/s00024-018-1870-5

[5] Grandis, H., & Sumintadireja, P. (2017). Improved Pseudo-Section Representation for CSAMT Data in Geothermal Exploration. IOP Conference Series: Earth and Environmental Science.

https://doi.org/10.1088/1755-1315/62/1/012035

[6] Grandis, H., Menvielle, M., & Roussignol, M. (1999). Bayesian Inversion with Markov Chains-I. The Magnetotelluric One-Dimensional Case. Geophysical Journal International, 138, 757-768.

https://doi.org/10.1046/j.1365-246x.1999.00904.x

[7] Guo, R., Dosso, S. E., Liu, J., Dettmer, J., & Tong, X. (2011). Non-Linearity in Bayesian 1-D Magnetotelluric Inversion. Geophysical Journal International, 185, 663-675.

https://doi.org/10.1111/j.1365-246X.2011.04996.x

[8] Malinverno, A. (2002). Parsimonious Bayesian Markov Chain Monte Carlo Inversion in a Nonlinear Geophysical Problem. Geophysical Journal International, 151, 675-688.

https://doi.org/10.1046/j.1365-246X.2002.01847.x

[9] Mandolesi, E., Ogaya, X., Campany, J., & Piana, N. (2018). A Reversible-Jump Markov Chain Monte Carlo Algorithm for 1D Inversion of Magnetotelluric Data. Computers and Geosciences, 113, 94-105.

https://doi.org/10.1016/j.cageo.2018.01.011

[10] Mohammad, I. H., Srigutomo, W., Sutarno, D., & Sumintadiredja, P. (2013). Interpretation of 1D Vector Controlled-Source Audio-Magnetotelluric (CSAMT) Data Using Full Solution Modeling. Journal of Mathematical and Fundamental Sciences, 45, 172-188.

https://doi.org/10.5614/j.math.fund.sci.2013.45.2.7

[11] Ray, A., & Key, K. (2012). Bayesian Inversion of Marine CSEM Data with a Trans-Dimensional Self Parametrizing Algorithm. Geophysical Journal International, 191, 1135-1151.

https://doi.org/10.1111/j.1365-246X.2012.05677.x

[12] Ray, A., Alumbaugh, D. L., Hoversten, G. M., & Key, K. (2013). Robust and Accelerated Bayesian Inversion of Marine Controlled-Source Electromagnetic Data Using Parallel Tempering. Geophysics, 78.

https://doi.org/10.1190/geo2013-0128.1

[13] Rosenkjaer, G. K. (2011). Electromagnetic Methods in Geothermal Exploration. 1D and 3D Inversion of TEM and MT Data from a Synthetic Geothermal Area and the Hengill Geothermal Area, SW Iceland. University of Iceland.

[14] Sen, M. K., & Stoffa, P. L. (1996). Bayesian Inference, Gibbs a€TM Sampler and Uncertainty Estimation in Geophysical Inversion. Geophysical Prospecting, 44, 313-350.

https://doi.org/10.1111/j.1365-2478.1996.tb00152.x

[15] Shireesha, M., & Harinarayana, T. (2013). Acta Geophysica Increased Resolution of Subsurface Parameters from 1D Magnetotelluric Modeling. Acta Geophysica, 61, 569-582.

https://doi.org/10.2478/s11600-012-0099-4

[16] Vozoff, K. (1991). The Magnetotelluric Method. In M. N. Nabighian (Ed.), Electromagnetic Methods in Applied Geophysics (Vol. 22, pp. 1943-1961).

https://doi.org/10.1190/1.9781560802686.ch8

[17] Ward, S. H., & Hohmann, G. W. (1988). Electromagnetic Theory for Geophysical Applications. Society of Exploration Geophysicists.

[18] Zhdanov, M. S. (2009). Methods in Geochemistry and Geophysics: Geophysical Electromagnetic Theory and Methods. Netherlands: Elsevier.