One-dimension forward modeling of MT(Python, Analytical)

One-dimension forward modeling of MT

Assumptions:

• Layer medium
• Uniform, homogeneous

Here is python implementation:

import cmath as cm
import numpy as np
import math

def mt1dan(freq, dz, sig):

    mu = 4e-7*math.pi
    ii = cm.sqrt(-1)

    nf = len(freq)
    nz = len(sig)

    zxy = []
    rho = []
    phs = []

    for kf in range(nf):
        omega = 2.0*math.pi*freq[kf]
        z = cm.sqrt(-ii*omega*mu/sig[nz-1])
        for m in range(nz-2, -1, -1):
            km = cm.sqrt(-ii*omega*mu*sig[m])
            z0 = -ii*omega*mu/km
            z = cm.exp(-2*km*dz[m])*(z-z0)/(z+z0)
            z = z0*(1+z)/(1-z)
        zxy.append(z)
        rho.append(np.abs(np.power(zxy[kf], 2))/(omega*mu))
        phs.append(math.atan2(zxy[kf].imag, zxy[kf].real)*180/math.pi)
    return rho, phs, zxy

This function enables us to calculate the conductivity, phase and impedance of underground medium while providing the model.

Note that freq is the frequency we use to observe, dz is the depth of layer(Unit:m), and sig is electrical conductance(Unit: S/m).

The formula we use could be easily found in any MT textbook.

$$Z_m=Z_{om}\frac{1+\frac{Z_{m+1}-Z_m}{Z_{m+1}+Z_m}{e}^{-2k_m(z_{m+1}-z_m)}}{1-\frac{Z_{m+1}-Z_m}{Z_{m+1}+Z_m}{e}^{-2k_m(z_{m+1}-z_m)}}$$

where $k_m=\sqrt{i\omega {\mu }_0{\sigma }_m}$, and $Z_{om}=\sqrt{i\omega {\mu }_0/{\sigma }_m}$.

the test file is also provided:

import numpy as np
import lab2
import matplotlib.pyplot as plt

# the model is described as following:
##########
# h1 = 100m,  sig1 = 0.01 S/m
##########
# h2 = 200m,  sig2 = 1.0 S/m
##########
# h3 = inf, sig3 = 0.02 S/m

# define the observation model
freq = np.logspace(-5, 8, num=20)

# comput
h = [100, 200]
sig = [0.01, 1.0, 0.02]

# realization of python, function mt1dan
rho1, phs1, zxy1 = lab2.mt1dan(freq, h, sig)

# plot the data.
plt.title("The analytic solution(TE) of 1D MT")
plt.subplot(2, 1, 1)
plt.semilogx(freq, rho1)
plt.semilogx(freq, rho1, 'o')

plt.xlabel("frequency/Hz")
plt.ylabel("apparent resistivity/$\Omega\cdot m$")

plt.subplot(2, 1, 2)
plt.semilogx(freq, phs1)
plt.semilogx(freq, phs1,'o')
plt.xlabel("frequency/Hz")
plt.ylabel('phase/degree')
plt.ylim(-80, 17)

plt.show()

And here is the result.