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.

$Z(\omega)=\frac{E_x(\omega)}{H_y(\omega)}=\sqrt{i\omega\mu_0\rho}$           (1)

where,

$Z$ - is the surface impedance
$\omega$ - is the angular frequency in radians $\omega=2\pi f$
$E_x$ - Electric field x-component in V
$H_y$ - Magnetic Field y-component in T
$\mu_0$ - Magnetic permeability of free space in Henry/m $4\pi \times 10^{-7}$
$\rho$ - Layer resistivity in $\Omega \cdot m$

There are several derivations to compute the apparent resistivity ($\rho_a$). Cagniard (1953) relates the apparent resistivity to the geo-electrical impedance ($Z$) ($\rho_a|Z|=\frac{1}{\omega\mu_0}|Z|^2$). Weidelt (1972) worked with response functions $C(\omega)=\frac{E_x(\omega)}{dE_x{\omega}/dx}$ where the corresponding definition of apparent resistivity being $\rho_a|C|=\omega \mu_0 |C(\omega)|^2=\rho_a|Z|$. 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).

Figure 2 - This Magnetotelluric forward modeling solution is conducted over a homogeneous 1D layered earth with n layers. Each layer has its own resistivity and thickness. As the naturally occurring electromagnetic waves propagate through each layer they must obey a number of boundary conditions. Each layer boundary can be considered to transmit and reflect waves.

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

$Z_{xy}=\frac{E_x}{H_y}=\sqrt{i\omega\mu_0\rho}$         (F1)

So starting at the basement (i.e., layer n) layer $Z_{xy}$ becomes  $Z_n$, so we can see that in F1,

$Z_n=\sqrt{i\omega\mu\rho_n}$          (F2)

$Z_n$ 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 $R_j$. Prior to computing the reflection coeficient the induction parameter, $\gamma_j$, exponential factor, $E_j$ and intrinsic impedance $w_j$ must be computed.

Induction Parameter$=\gamma_j=\sqrt{i\omega\mu_0\sigma_j}$          (F3)

Exponential Factor$=E_j=\exp{-2\gamma_j h_j}$          (F4)

Intrinsic Impedance$=w_j=\gamma_j\rho_j$          (F5)

where,

$j$-the jth layer
$\sigma$-Layer conductivity $S/m$

The reflection coeficient is then computed by,

Reflection Coeficient=$R_j=\frac{w_j-Z_j+1}{w_j+Z_j+1}$          (F6)

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

Impedance$=Z_j=w_j\frac{1-R_j E_j}{1+R_j E_j}$          (F7)

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

Apparent Resistivity $=\rho_a=\frac{1}{\omega}|Z_1|^2$          (F8)

Phase $=tan^{-1}{\frac{Im(Z_1)}{Re(Z_1)}}$          (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 $Z_n$ (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));

# 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];

Figure 1: A schematic of the resulting geo-electrical model.

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 $Z_n$ (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);