I constantly ask myself where the hell are the geophysics tutorials? Coming from a programming background I find the extreme lack of geophysics tutorials concerning. If you want to solve the most trivial issues in computing the internet has it covered. When you look for trivial geophysics solutions they are nowhere to be found. As geoscientists, we live in a small industry and I expect less free online resources but come on, this is ridiculous.

As a start I include the following tutorial on 1D MT forward modelling. Hopefully it is easy to follow.

*...and oh yeah.... the source code is downloadable and is licensed under WTFPL. If you don't know what WTFPL is, ill give you a hint "do what the _ you want to". So in other words, distribute, copy, modify, sell, integrate, delete, make erroneous.... I really don't care what you do with it.*

# Introduction

The magnetotelluric (MT) method is an electromagnetic geophysical exploration technique. The technique utilizes natural sources of electromagnetic energy for very deep soundings of the earth for crustal, global, geothermal, mineral and hydrocarbon exploration.

We understand that magnetotelluric and telluric fields are large-scale low frequency electromagnetic fields generated by natural source electromagnetic source. The terms "magnetotelluric" and "telluric" are used to describe the fields and currents respectively. These MT sources generally occur outside of the earth including solar emission resulting from the diurnal variation of the earth's magnetic field (i.e., Arora's).

The point of this article is not to go in depth into MT interpretation or background but the formulation of the 1D solution. The derivation and background has been taken from several sources including Niwas et al. (2005) and Pedersen and Hermance (1986)

For a uniform geo-electrical earth the measured surface impedance is related to the true resistivity of the earth.

(1)

where,

- is the surface impedance

- is the angular frequency in radians

- Electric field x-component in V

- Magnetic Field y-component in T

- Magnetic permeability of free space in Henry/m

- Layer resistivity in

There are several derivations to compute the apparent resistivity (). Cagniard (1953) relates the apparent resistivity to the geo-electrical impedance () (). Weidelt (1972) worked with response functions where the corresponding definition of apparent resistivity being . For the following case we are using Max Meju's solution, one of the more simpler solutions to follow.

# Main Formulation

Most, if not all, 1D MT formulations are an iterative process. For the benefit of this formulation consider Figure 2. Figure 2 shows the labeling conventions employed. The formulation starts from the basement (i.e., layer n) and iterates through the above layers. In other words the solution is computed and passed to the above layer (from layer n to layer 1).

The process begins by solving the homogeneous infinite half space impedance. The impedance can be thought of as the proportionality between and ,

(F1)

So starting at the basement (i.e., layer n) layer becomes , so we can see that in F1,

(F2)

is a complex number and can be thought of as the theoretical impedance of the basement. So if there were no above layers, this would be the measured impedance. So to solve the above layers we need to know the energy that is reflected and transmitted at each layer boundary. This is done by computing the reflection coeficient . Prior to computing the reflection coeficient the induction parameter, , exponential factor, and intrinsic impedance must be computed.

Induction Parameter (F3)

Exponential Factor (F4)

Intrinsic Impedance (F5)

where,

-the jth layer

-Layer conductivity

The reflection coeficient is then computed by,

Reflection Coeficient= (F6)

The impedance of the layer is then computed from the reflection coeficient,

Impedance (F7)

Once all impedances are computed (i.e., layers j=n ... 1) take the top layer impedance (layer j=1) and compute apparent resistivity and phase,

Apparent Resistivity (F8)

Phase (F9)

# Coding the Solution: Matlab

Let's start this by creating a function called **modelMT**. The function modelMT will compute the apparent resistivity given a 1D geoelectrical model and MT frequency. The method will take in three arguments, the layer **resistivities**, **thicknesses** and a single frequency. The output of the program is an **apparent resistivity** in Ohm.m and **phase** in radians.

1 | function [apparentResistivity, phase] = modelMT(resistivities, thicknesses,frequency) |

The next step is to pre-compute constants which we will use a number of times in the following computation.This includes the **magnetic permeability**, **angular frequency** and **number of layers**.

2 3 4 | mu = 4*pi*1E-7; %Magnetic Permeability (H/m) w = 2 * pi * frequency; %Angular Frequency (Radians); n=length(resistivities); %Number of Layers |

As we have seen in the formulation, each layer has its own impedance. So to begin we will initialise an array called **impedances**.

5 | impedances = zeros(n,1); |

The first step is to compute the halfspace impedance **(see F2)** and assign it to the final value of the **impedances** array corresponding to the halfspace,

6 7 | Zn = sqrt(sqrt(-1)*w*mu*resistivities(n)); impedances(n) = Zn; |

So to solve the 1D MT problem you need to iterate through layers starting from layer j=n-1 (i.e. the layer above the basement) and finish in the top layer (j=1)

8 9 10 | for j = n-1:-1:1 resistivity = resistivities(j); thickness = thicknesses(j); |

within this loop you need to compute equations F3 to F5. These include the **induction parameter**, **exponential factor** and **intrinsic impedance**. Ensure that when you compute the conductivity a float and not an integer is used. So use **1.0/resistivity** rather than **1/resistivity**.

11 12 13 | dj = sqrt(sqrt(-1)* (w * mu * (1.0/resistivity))); wj = dj * resistivity; ej = exp(-2*thickness*dj); |

The next step is to calculate the **reflection coeficient (F6)** and **impedance (F7)** using the current layer intrinsic impedance and the prior computer layer impedance **j+1**.

14 15 16 17 18 19 | belowImpedance = impedances(j + 1); rj = (wj - belowImpedance)/(wj + belowImpedance); re = rj*ej; Zj = wj * ((1 - re)/(1 + re)); impedances(j) = Zj; end |

Finally you can compute the **apparent resistivity F8** and **phase F9** and assign them to the output matrix!

20 21 22 23 | Z = impedances(1); absZ = abs(Z); apparentResistivity = (absZ * absZ)/(mu * w); phase = atan2(imag(Z),real(Z)); |

**Scared of Programming?** Download the Matlab Code

**1D Magnetotelluric Forward Modelling (Matlab) 2.42 KB**

# Coding the Solution: Python

Import the required libraries **cmath** and **math**. The **math** library contains real number mathematical functions and contains the same functions found in the ANSI c library. The **cmath** library contains mathematical functions for complex number calculation.

1 2 | import math import cmath |

Initialize constants, geo-electrical and recording parameters. **mu** is the Magnetic Permeability (H/m) and is only computed once at the start of the program. The **resistivities** array contains each layer resistivity in Ohm.m, starting with the top layer in index 0 (left) and the half-space resistivity in index 3 (right). The **thicknesses **array contains the thickness in m for each of the layers, starting with the top layer thickness in index 0 and ending with the thickness of the layer above the half-space. This program will accept any number of layers, but you must ensure that if N is the number of resistivities, the thicknesses array must have N-1 elements (See Figure 1 for resulting geo-electrical model). The nlayers variable computes the number of layers (including the halfspace) this is so it does not have to be computed for each iteration. Finally **frequencies** contains all the frequencies, in Hz, that the apparent resistivity will be calculated for.

*Note: The program will work with any number of layers and frequencies.*

3 4 5 6 7 | mu = 4*math.pi*1E-7; #Magnetic Permeability (H/m) resistivities = [1, 2, 3]; thicknesses = [50, 300]; n = len(resistivities); frequencies = [0.0001,0.005,0.01,0.05,0.1,0.5,1,5,10,50,100,500,10000]; |

Using a for loop, iterate over each frequency found in the **frequencies** array and compute the apparent resistivity. The method requires the angular frequency (**w**) and the **impedance** array is reset each iteration for robustness.

8 9 10 | for frequency in frequencies: w = 2*math.pi*frequency; impedances = list(range(n)); |

The first step is to compute the halfspace impedance **(see F2)** and assign it to the final value of the **impedances** array corresponding to the halfspace. The basement impedance is computed using the cmath.sqrt operator.

*Note: 1j is the equivalent of the complex number 0 + 1i*

10 | impedances[nlayers-1] = cmath.sqrt(w*mu*resistivities[nlayers-1]*1j); |

So to solve the 1D MT problem you need to iterate through layers starting from layer j=n-1 (i.e. the layer above the basement) and finish in the top layer (j=1)

Within this loop you need to compute equations F3 to F5. These include the **induction parameter**, **exponential factor** and **intrinsic impedance**.

11 12 13 14 15 16 | for j in range(n-2,-1,-1): resistivity = resistivities[j]; thickness = thicknesses[j]; dj = cmath.sqrt((w * mu * (1/resistivity))*1j); wj = dj * resistivity; ej = cmath.exp(-2*thickness*dj); |

The next step is to calculate the **reflection coeficient (F6)** and **impedance (F7)** using the current layer intrinsic impedance and the prior computer layer impedance **j+1**.

17 18 19 20 21 | belowImpedance = impedances[j + 1]; rj = (wj - belowImpedance)/(wj + belowImpedance); re = rj*ej; Zj = wj * ((1 - re)/(1 + re)); impedances[j] = Zj; |

Finally you can compute the **apparent resistivity F8** and **phase F9** and print the resulting data!

22 23 24 25 26 | Z = impedances[0]; absZ = abs(Z); apparentResistivity = (absZ * absZ)/(mu * w); phase = math.atan2(Z.imag, Z.real); print(frequency, '\t', apparentResistivity, '\t', phase); |

**Scared of Programming?** Download the Python Code

**1D Magnetotelluric Forward Modelling (Python) 2.06 KB**

# References

Basokur A T 1994 Definitions of apparent resistivity for the presentation of magnetotelluric sounding data; Geophysical Prospecting 42 141–149.

Cagniard L 1953 Basic theory of magnetotelluric method of geophysical prospecting; Geophysics 18 605–635.

Niwas S, Gupta P K and Gaur V K 2005 Normalized impedance function and the straightforward inversion scheme for magnetotelluric data; J. Earth Syst. Sci. 114 5 523-531.

Pederesen J F and Hermance 1986 Least-square inversion of one-dimensional magnetotelluric data: An assessment of procedures employed by Brown University; Surv. Geophys.

8 187–231.

Weidelt P 1972 The inverse problem of geomagnetic induction; Z. fur. Geophys. 38 257–289.

## Leave A Comment