In [ ]:
import SimPEG as simpeg
from SimPEG import NSEM
import MT_poster_utils
from IPython.display import HTML
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
https://drive.google.com/open?id=0B4oB1gnNZnWNcmdncHh3WVZqa1U
Short answer: The magnetotelluric method is a geophysical technique that relays on Natrual Eletromagnetic (EM) waves as a source of energy to image the subsurface. The method is widely used in geophysical exploration, specially imaging deep structures such as in geothermal, kimberlite or techtonic exploration.
Long answer: Keep on reading...
In this notebook with are sharing the work we did for our SciPy 2016 poster, where we used NSEM code that are built with the SimPEG package (http://simpeg.xyz) to make data from a geological model (forward model) and to image the model again (invert the data).
This poster is a part of a community effort, using a diamond exploration example as motivation to use the large spectrum of geophysical methods within SimPEG for the same geological setting. We want to demonstrate the advantages of having a shared language between the governing physics of all these problems, the written codes, and other utilitices used in the work. We will show
The northern lights are generated by interaction of charged particles and the Earth’s magnetic field, and along with lightning strikes around the world, generate electromagnetic energy that induces electrical currents in the ground. These currents are affected by electrical conductivity contrasts and by measuring electromagnetic fields at the surface, we can image the subsurface structures using geophysical numerical simulation and optimization.
In [ ]:
HTML("Figures/Magnetotelluric_Movie_ThibautAstic.html")
Details of the physics at: http://em.geosci.xyz/content/maxwell3_fdem/MT/MT_N_layered_Earth.html
The Maxwell's equations govern the propagation of EM fields in the Earth. In frequency domain (and using the quasi-static approximation and a $e^{i\omega t}$ time relationship), they are given as
$$ \begin{align} \nabla \times \textbf{E} + i \omega \textbf{B} &= 0 \\ \nabla \times \frac{1}{\mu} \textbf{B} - \sigma \textbf{E} &= \textbf{s}_E \end{align} $$The system can be reduced to a single equation
$$ \begin{align} (\nabla \times \frac{1}{\mu} \nabla \times + i \omega \sigma) \textbf{E}_s &= - i \omega \Delta \sigma \textbf{E}_p \end{align} $$which is then solved numerically using finite volume methods, in SimPEG.
We start with a geological model which is based on a surfaces derived from diamond exploration.
We need to discritize the Earth and correlate the geological units with physical property (conductivity) need to generate the data.
Below we use built in functions in SimPEG to read VTK (http://www.vtk.org) files of the physical property model in to mesh and model of the conductivity.
In [ ]:
# Load the geological discretized model
mesh, modelDict = simpeg.Mesh.TensorMesh.readVTK('./datafiles/nsmesh_model.vtr')
sigma = modelDict['S/m']
In [ ]:
# Print model information
print mesh.nC," cells"
print mesh
In [ ]:
# Define the area of interest in UTM (meters)
bw, be = 557100, 557580
bs, bn = 7133340, 7133960
bb, bt = 0,480
In [ ]:
# View the model
slice=20
fig,ax=plt.subplots(figsize=(10,8))
modelplot = mesh.plotSlice(np.log10(sigma),normal='Y',ind=slice,ax=ax,grid=True, pcolorOpts={"cmap":"viridis"})
ax.set_xlim([bw,be])
ax.set_ylim([0,bt])
ax.set_aspect('equal')
plt.colorbar(modelplot[0])
ax.set_title("Discretized Model",fontsize=24)
In [ ]:
# Load stations locations and frequency range
locs = np.load('./datafiles/MTlocations.npy')
freqList = np.load('./datafiles/MTfrequencies.npy')
In [ ]:
# View a scatter plot of the locations at the surface
plt.scatter(locs[:,0],locs[:,1])
In [ ]:
# List of the frequencies used for the problem
print freqList
The natural source of the MT method is unknown and cannot be seporated from the induced signal due to the Earth's conductivity or explicitly modelled numercially. To address this the ratio's of the meaured fields are used, and assuming that the source is the same, the source cancels out, leaving only the secondary signal in the data.
In [ ]:
# Load the data - stored as numpy.recArray
mtData = np.load('./datafiles/MTdata.npy')
In [ ]:
# Plot data
fig, axes, csList = MT_poster_utils.pseudoSect_OffDiagTip_RealImag(mtData,{'y':7133627.5},colBarMode='each')
In [ ]:
# Make the plot
fig, axes, csList = MT_poster_utils.pseudoSect_FullImpTip_RealImag(mtData,{'y':7133627.5},colBarMode='each')
To image the Earth, based on the data we forward modelled above, we use inversion. As with the forward modelling, the inversion is a computationally intensive procedure, and as before we ran yours on a cluster. A single iteration takes at the order of 4 hours to solve using ~ 20Gb of ram.
We present 2 inversion results, for the impedance and the tipper.
In [ ]:
#Load Model Off-diagonal
mesh,inv= simpeg.Mesh.TensorMesh.readVTK('./datafiles/recoveredMod_off_it18.vtr')
siginvoff=inv['S/m']
In [ ]:
#Load Model Tipper
mesh,invtip= simpeg.Mesh.TensorMesh.readVTK('./datafiles/recoveredMod_tip_it36.vtr')
siginvtip=invtip['S/m']
In [ ]:
MT_poster_utils.CompareInversion(mesh,sigma,siginvoff,siginvtip,slice_ver=20,slice_hor=40)
In [ ]: