+ +
+ + + + + + + + + + +
![](\"Figures/MTDiamonds-02.png\")
"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "#### Animation of the physical phenomenon"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "HTML(\"Figures/Magnetotelluric_Movie_ThibautAstic.html\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "Details of the physics at:\n",
+ "http://em.geosci.xyz/content/maxwell3_fdem/MT/MT_N_layered_Earth.html"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Physics (quasi-static approximation)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "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 \n",
+ "\n",
+ "$$\n",
+ "\\begin{align}\n",
+ "\\nabla \\times \\textbf{E} + i \\omega \\textbf{B} &= 0 \\\\\n",
+ "\\nabla \\times \\frac{1}{\\mu} \\textbf{B} - \\sigma \\textbf{E} &= \\textbf{s}_E \n",
+ "\\end{align}\n",
+ "$$\n",
+ "\n",
+ "The system can be reduced to a single equation \n",
+ "\n",
+ "$$\n",
+ "\\begin{align}\n",
+ "(\\nabla \\times \\frac{1}{\\mu} \\nabla \\times + i \\omega \\sigma) \\textbf{E}_s &= - i \\omega \\Delta \\sigma \\textbf{E}_p \n",
+ "\\end{align}\n",
+ "$$\n",
+ "\n",
+ "which is then solved numerically using finite volume methods, in SimPEG."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Load the discrezited geological model"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "#### Original isosurfaces model\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We start with a geological model which is based on a surfaces derived from diamond exploration. \n",
+ "\n",
+ "![](\"Figures/GeoModel.png\")
\n",
+ "\n",
+ "We need to discritize the Earth and correlate the geological units with physical property (conductivity) need to generate the data. \n",
+ "\n",
+ "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."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "# Load the geological discretized model\n",
+ "mesh, modelDict = simpeg.Mesh.TensorMesh.readVTK('./datafiles/nsmesh_model.vtr')\n",
+ "sigma = modelDict['S/m']"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "# Print model information\n",
+ "print mesh.nC,\" cells\"\n",
+ "print mesh"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "# Define the area of interest in UTM (meters)\n",
+ "bw, be = 557100, 557580\n",
+ "bs, bn = 7133340, 7133960\n",
+ "bb, bt = 0,480"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "# View the model\n",
+ "slice=20\n",
+ "fig,ax=plt.subplots(figsize=(10,8))\n",
+ "modelplot = mesh.plotSlice(np.log10(sigma),normal='Y',ind=slice,ax=ax,grid=True, pcolorOpts={\"cmap\":\"viridis\"})\n",
+ "ax.set_xlim([bw,be])\n",
+ "ax.set_ylim([0,bt])\n",
+ "ax.set_aspect('equal')\n",
+ "plt.colorbar(modelplot[0])\n",
+ "ax.set_title(\"Discretized Model\",fontsize=24)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "#### Paraview view"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "![](\"Figures/NSEMorignalModelPV.png\")
"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Load stations locations and frequency range\n",
+ "\n",
+ "Below with load in the information needed for solving the problem"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "# Load stations locations and frequency range\n",
+ "locs = np.load('./datafiles/MTlocations.npy')\n",
+ "freqList = np.load('./datafiles/MTfrequencies.npy')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "# View a scatter plot of the locations at the surface\n",
+ "plt.scatter(locs[:,0],locs[:,1])"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "# List of the frequencies used for the problem\n",
+ "print freqList"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Types of data\n",
+ "\n",
+ "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. "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "#### Impedance\n",
+ "\n",
+ "Is the ratio of the horizontal electrical over magnetic fields. It is a complex tensor, defined as"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\\begin{align}\n",
+ "\\textbf{Z}(\\omega) =\n",
+ "\\begin{bmatrix}\n",
+ "Z_{xx}(\\omega) & Z_{xy}(\\omega)\\\\\n",
+ "Z_{yx}(\\omega) & Z_{yy}(\\omega)\n",
+ "\\end{bmatrix}\n",
+ "\\end{align} \n",
+ "\\qquad where \\qquad\n",
+ "\\begin{bmatrix}\n",
+ "E_{x}(\\omega)\\\\\n",
+ "E_{y}(\\omega) \n",
+ "\\end{bmatrix} = \n",
+ "\\textbf{Z}(\\omega) \n",
+ "\\begin{bmatrix}\n",
+ "H_{x}(\\omega)\\\\\n",
+ "H_{y}(\\omega) \n",
+ "\\end{bmatrix}\n",
+ "$$\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "#### Tipper\n",
+ "\n",
+ "Is the ratio of the vertical magnetic field over horizontal. It is a complex vector."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\\begin{align}\n",
+ "\\textbf{T}(\\omega) =\n",
+ "\\begin{bmatrix}\n",
+ "T_{zx}(\\omega) \\\\\n",
+ "T_{zy}(\\omega)\n",
+ "\\end{bmatrix} \n",
+ "\\qquad where \\qquad\n",
+ "H_{z}(\\omega) = \n",
+ "\\textbf{T}(\\omega ) \\begin{bmatrix}\n",
+ "H_{x}(\\omega) \\\\\n",
+ "H_{y}(\\omega)\n",
+ "\\end{bmatrix}\n",
+ "\\end{align}$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Run Forward modeling script on a cluster\n",
+ "\n",
+ "Running the numerical codes to solve the proble are resource heavy. The file below is the script we used to run the forward modelling on a cluster."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "[Forward Modeling](datafiles/findDiam_MTforward_HKPK1.py)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Data Visualization\n",
+ "\n",
+ "The data are visulized as sections. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false,
+ "scrolled": true
+ },
+ "outputs": [],
+ "source": [
+ "# Load the data - stored as numpy.recArray\n",
+ "mtData = np.load('./datafiles/MTdata.npy')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false,
+ "scrolled": false
+ },
+ "outputs": [],
+ "source": [
+ "# Plot data\n",
+ "fig, axes, csList = MT_poster_utils.pseudoSect_OffDiagTip_RealImag(mtData,{'y':7133627.5},colBarMode='each')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "# Make the plot\n",
+ "fig, axes, csList = MT_poster_utils.pseudoSect_FullImpTip_RealImag(mtData,{'y':7133627.5},colBarMode='each')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Run the inversions on a cluster\n",
+ "\n",
+ "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. \n",
+ "\n",
+ "We present 2 inversion results, for the impedance and the tipper. "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Links to the files used to run the inversions.\n",
+ "\n",
+ "[Off-diagonal impedance inversion](datafiles/findDiam_inversion_offdiagonal.py)\n",
+ "\n",
+ "[Tipper inversion](datafiles/findDiam_inversion_tipper.py)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Load Inversion results"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "#Load Model Off-diagonal\n",
+ "mesh,inv= simpeg.Mesh.TensorMesh.readVTK('./datafiles/recoveredMod_off_it18.vtr')\n",
+ "siginvoff=inv['S/m']"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "#Load Model Tipper\n",
+ "mesh,invtip= simpeg.Mesh.TensorMesh.readVTK('./datafiles/recoveredMod_tip_it36.vtr')\n",
+ "siginvtip=invtip['S/m']"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "MT_poster_utils.CompareInversion(mesh,sigma,siginvoff,siginvtip,slice_ver=20,slice_hor=40)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 2",
+ "language": "python",
+ "name": "python2"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 2
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython2",
+ "version": "2.7.11"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}
diff --git a/notebooks/NSEM/Figures/Diamond.png b/notebooks/NSEM/Figures/Diamond.png
new file mode 100644
index 0000000..e604e02
Binary files /dev/null and b/notebooks/NSEM/Figures/Diamond.png differ
diff --git a/notebooks/NSEM/Figures/GeoModel.png b/notebooks/NSEM/Figures/GeoModel.png
new file mode 100644
index 0000000..15dba69
Binary files /dev/null and b/notebooks/NSEM/Figures/GeoModel.png differ
diff --git a/notebooks/NSEM/Figures/MTDiamonds-02.png b/notebooks/NSEM/Figures/MTDiamonds-02.png
new file mode 100644
index 0000000..6d4a300
Binary files /dev/null and b/notebooks/NSEM/Figures/MTDiamonds-02.png differ
diff --git a/notebooks/NSEM/Figures/MT_Data_Locations.png b/notebooks/NSEM/Figures/MT_Data_Locations.png
new file mode 100644
index 0000000..630aae7
Binary files /dev/null and b/notebooks/NSEM/Figures/MT_Data_Locations.png differ
diff --git a/notebooks/NSEM/Figures/MT_conceptsketch.png b/notebooks/NSEM/Figures/MT_conceptsketch.png
new file mode 100644
index 0000000..63c9c52
Binary files /dev/null and b/notebooks/NSEM/Figures/MT_conceptsketch.png differ
diff --git a/notebooks/NSEM/Figures/MT_offtip_results.png b/notebooks/NSEM/Figures/MT_offtip_results.png
new file mode 100644
index 0000000..6dc74d8
Binary files /dev/null and b/notebooks/NSEM/Figures/MT_offtip_results.png differ
diff --git a/notebooks/NSEM/Figures/Magnetotelluric_Movie_ThibautAstic.html b/notebooks/NSEM/Figures/Magnetotelluric_Movie_ThibautAstic.html
new file mode 100644
index 0000000..0f06756
--- /dev/null
+++ b/notebooks/NSEM/Figures/Magnetotelluric_Movie_ThibautAstic.html
@@ -0,0 +1,278 @@
+
+
+
+