One-dimension forward modeling of MT(python, numerical)

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
import scipy.linalg as la

def mt1dte(freq, dz, sig):
    mu = 4.0E-7 * math.pi
    ii = cm.sqrt(-1)
    omega = 2.0 * math.pi * freq
    nz = len(sig)
    sig = np.append(sig, sig[nz - 1])
    # dz.append(math.sqrt(2.0 / (sig[nz] * omega * mu)))
    dz = np.append(dz, math.sqrt(2.0 / (sig[nz] * omega * mu)))
    diagA = []
    offdiagA = []

    for ki in range(nz):
        diagA.append(ii*omega*mu*(sig[ki]*dz[ki]+sig[ki+1]*dz[ki+1])-2/dz[ki]-2/dz[ki+1])

    for ki in range(nz-1):
        offdiagA.append(2/dz[ki+1])

    mtxA = np.diag(diagA) + np.diag(offdiagA, k=-1) + np.diag(offdiagA, k=1)
    rhs = np.zeros(nz)
    rhs[0] = - 2.0 / dz[0]
    a = la.inv(mtxA)
    a = np.dot(a, rhs)
    #a = la.solve(mtxA, rhs)
    ex = [1.00]
    ex.extend(a)
    return ex

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 math
import matplotlib.pyplot as plt
import lab2
import lab3

freq = np.logspace(-4, 2, num=20)
h = [100, 200]
sig = [0.01, 1.0, 0.02]
rho1, phs1, zxy1 = lab2.mt1dan(freq, h, sig)

# numerical solution
air = [100000.0, 30000.0, 10000.0, 3000.0, 1000.0, 300.0, 100.0, 30.0, 30.0, 10.0]
ground = [10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0]
layer2 = [20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0]
base = [30.0, 50.0, 80.0, 100.0, 300.0, 300.0, 500.0, 700.0, 800.0, 1000.0, 1000.0, 2000.0, 2000.0, 3000.0, 5000.0, 8000.0, 10000.0, 15000.0, 20000.0]
dz = []
dz.extend(air)
dz.extend(ground)
dz.extend(layer2)
dz.extend(base)

sig = []
for i in range(10):
    sig.append(0)
for i in range(10):
    sig.append(0.01)
for i in range(10):
    sig.append(1.0)
for i in range(19):
    sig.append(0.02)


nza = 10

# forward modeling
mu = 4e-7*math.pi
omega = 2.0*math.pi*freq
zxy2 = []
rho2 = []
phs2 = []
rho_error = []
phase_error = []

for k in range(len(freq)):
    ex = lab3.mt1dte(freq[k], dz, sig)
    exs = ex[nza+1]
    hys = -(ex[nza+2]-ex[nza+1]) / dz[nza+1] / cm.sqrt(-1) / omega[k]/mu
    zxy2.append(exs/hys)

    rho2.append(np.abs(np.power(zxy2[k], 2))/(omega[k]*mu))
    rho_error.append(rho2[k]-rho1[k])
    phs2.append(180-math.atan2(zxy2[k].imag, zxy2[k].real)*180.0/math.pi)
    phase_error.append(phs2[k]-phs1[k])

plt.subplot(2, 1, 1)
plt.title("The analytic and numerical solution(TE) of 1D MT")
plt.semilogx(freq, rho1, 'b-', freq, rho2, 'ro')
plt.semilogx(freq, rho_error, 'y--')
ax = plt.gca()
ax.invert_xaxis()
plt.xlabel("frequency/Hz")
plt.ylabel("apparent resistivity/$\Omega\cdot m$")
plt.legend(['analytic', 'numerical', 'diff'])
plt.grid()

plt.subplot(2, 1, 2)
plt.semilogx(freq, phs1, 'b-', freq, phs2, 'ro')
plt.semilogx(freq, phase_error, 'y--')
ax = plt.gca()
ax.invert_xaxis()
plt.xlabel("frequency/Hz")
plt.ylabel("Phase/degree")
plt.legend(['analytic', 'numerical', 'diff'])
plt.grid()
plt.show()
plt.grid()

And here is the result.