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

The deterministic geophysical inversion methods are dominant when inverting magnetotelluric data whereby its results largely depends on the assumed initial model and only a single representative solution is obtained. A common problem to this approach is that all inversion techniques suffer from non-uniqueness since all model solutions are subjected to errors, under-determination and uncertainty. A statistical approach in nature is a possible solution to this problem as it can provide extensive information about unknown parameters. In this paper, we developed a 1D Bayesian inversion code based Metropolis-Hastings algorithm whereby the uncertainty of the earth model parameters were quantified by examining the posterior model distribution. As a test, we applied the inversion algorithm to synthetic model data obtained from available literature based on a three layer model (K, H, A and Q). The frequency for the magnetotelluric impedance data was generated from 0.01 to 100 Hz. A 5% Gaussian noise was added at each frequency in order to simulate errors to the synthetic results. 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.


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.Journal of Geoscience and Environment Protection 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.

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: 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: where µ is magnetic permeability (Vs/Am), ε is dielectric permittivity (As/Vm) and σ is electrical conductivity (S/m).The magnetic permeability and dielec-Journal of Geoscience and Environment Protection tric permittivity are further defined as: r µ is the relative magnetic permeability and 0 µ is the magnetic permeabili- ty of the free space whose value is 7 4π 10 − × H/m.r ε is the relative electric permittivity and 0 ε is the electric permittivity of the free space whose value is 12 8.85 10 − × F/m.The parameters µ , ε and σ 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 ρ (reciprocal of σ ), will also determine the depth of penetra- tion of EM field through Earth materials.The depth of penetration, called the skin depth ( δ ) is given by: where ω is the angular frequency.Since ρ is proportional to and f is inversely proportional to δ , it follows that deeper depth of investigation will be achieved with increase in ρ 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: It is not easy to visualize impedance tensors and thus, a useful and widely accepted tool is to transform these values into apparent resistivity ( a ρ ) and phase ( ∅ ) respectively: 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).

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: We can substitute A and B in Equation ( 14) in terms of model parameters (m) and observed data vectors (d) respectively as: The likelihood function ( ) L m , is used to measure the difference between a set of data and a proposed model m .Therefore,

( )
L m 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: where ) misfit for the evaluated model m .A small misfit between the modeled response ( ) f m 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.

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 α is calculated from:  ( ) q m m 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: ( ) The new model is better than the old one when 1 α ≥ , and thus is always ac- cepted.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 1 α < , a random number is selected be- tween [0, 1] and if that number lies within the interval [ ] ,1 α 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.

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.

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 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.

H-Type Model
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   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.

K-Type Model
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.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.

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

Figure 2 .
Figure 2. Resistivity distribution of H-type model.

Figure 3 .
Figure 3. Thickness distribution of H-type model.

Figure 4 .
Figure 4. Convergence trend of H-type model.

Figure 5 .
Figure 5. Resistivity distribution of K-type model.

Figure 6 .
Figure 6.Thickness distribution of K-type model.

Figure 7 .
Figure 7. Convergence trend of K-type model.

Figure 8 .
Figure 8. Resistivity distribution of A-type model.

Figure 9 .
Figure 9. Thickness distribution of A-type model.
Figure 12 respectively, again with the correct synthetic layer value displayed in horizontal red line.

Figure 10 .
Figure 10.Convergence Trend of H-type model.

Figure 11 .
Figure 11.Resistivity distribution of Q-type model.

Figure 12 .
Figure 12.Thickness distribution of Q-type model.

Figure 13 .
Figure 13.Convergence trend of Q-type model.