{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib notebook" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# SEG-EAGE 3D Salt Model\n\n[Muld07]_ presented electromagnetic responses for a resistivity model which he\nderived from the seismic velocities of the SEG/EAGE salt model [AmBK97]_.\n\nHere we reproduce and store this resistivity model, and compute electromagnetic\nresponses for it.\n\n## Velocity-to-resistivity transform\n\nQuoting here the description of the velocity-to-resistivity transform used by\n[Muld07]_:\n\n \u00abThe SEG/EAGE salt model (Aminzadeh et al. 1997), originally designed for\n seismic simulations, served as a template for a realistic subsurface model.\n Its dimensions are 13500 by 13480 by 4680 m. The seismic velocities of the\n model were replaced by resistivity values. The water velocity of 1.5 km/s\n was replaced by a resistivity of 0.3 Ohm m. Velocities above 4 km/s,\n indicative of salt, were replaced by 30 Ohm m. Basement, beyond 3660 m\n depth, was set to 0.002 Ohm m. The resistivity of the sediments was\n determined by $(v/1700)^{3.88}$ Ohm m, with the velocity v in m/s\n (Meju et al. 2003). For air, the resistivity was set to $10^8$ Ohm\n m.\u00bb\n\nEquation 1 of [MeGM03]_, is given by\n\n\\begin{align}:label: meju\n\n \\log_{10}\\rho = m \\log_{10}V_P + c \\ ,\\end{align}\n\nwhere $\\rho$ is resistivity, $V_P$ is P-wave velocity, and for\n$m$ and $c$ 3.88 and -11 were used, respectively.\n\nThe velocity-to-resistivity transform uses therefore a Faust model ([Faus53]_)\nwith some additional constraints for water, salt, and basement.\n\n

#### Note

The SEG/EAGE Salt Model is licensed under the CC-BY-4.0\n _.

\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import os\nimport pooch\nimport emg3d\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom matplotlib.colors import LogNorm\nplt.style.use('bmh')\n\n\n# Adjust this path to a folder of your choice.\ndata_path = os.path.join('..', 'download', '')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Fetch the model\n\nRetrieve and load the pre-computed resistivity model.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fname = \"SEG-EAGE-Salt-Model.h5\"\npooch.retrieve(\n 'https://raw.github.com/emsig/data/2021-05-21/emg3d/models/'+fname,\n '6ee10663de588d445332ba7cc1c0dc3d6f9c50d1965f797425cebc64f9c71de6',\n fname=fname,\n path=data_path,\n)\nfmodel = emg3d.load(data_path + fname)['model']\nfgrid = fmodel.grid" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## QC resistivity model\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Limit colour-range to [0.3, 50] Ohm.m\n# (affects only the basement and air, improves resolution in the sediments).\nvmin, vmax = 0.3, 50\n\nfgrid.plot_3d_slicer(\n fmodel.property_x,\n zslice=-2000,\n pcolor_opts={'norm': LogNorm(vmin=vmin, vmax=vmax)}\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Compute some example CSEM data with it\n\n### Survey parameters\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Create a source instance\nsource = emg3d.TxElectricDipole(\n coordinates=[6400, 6600, 6500, 6500, -50, -50],\n strength=1/200. # Normalize for length.\n)\n\n# Frequency (Hz)\nfrequency = 1.0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Initialize computation mesh\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "grid = emg3d.construct_mesh(\n frequency=frequency,\n properties=[0.3, 1, 2, 15],\n center=(6500, 6500, -50),\n seasurface=0,\n domain=([0, 13500], [0, 13500], None),\n vector=(None, None, np.array([-100, -80, -60, -40, -20, 0])),\n min_width_limits=([5, 100], [5, 100], [5, 20]),\n min_width_pps=5,\n stretching=[1.03, 1.05],\n lambda_from_center=True,\n)\ngrid" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Put the salt model onto the modelling mesh\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Interpolate full model from full grid to grid\nmodel = fmodel.interpolate_to_grid(grid)\n\ngrid.plot_3d_slicer(\n model.property_x,\n zslice=-2000,\n zlim=(-4180, 500),\n pcolor_opts={'norm': LogNorm(vmin=vmin, vmax=vmax)}\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solve the system\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "efield = emg3d.solve_source(\n model, source, frequency,\n semicoarsening=False,\n linerelaxation=False,\n verb=1,\n)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "grid.plot_3d_slicer(\n efield.fx.ravel('F'),\n zslice=-2000,\n zlim=(-4180, 500),\n view='abs',\n v_type='Ex',\n pcolor_opts={'norm': LogNorm(vmin=1e-16, vmax=1e-9)}\n)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Interpolate for a \"more detailed\" image\nx = grid.cell_centers_x\ny = grid.cell_centers_y\nrx = np.repeat([x, ], np.size(x), axis=0)\nry = rx.transpose()\nrz = -2000\ndata = efield.get_receiver((rx, ry, rz, 0, 0))\n\n# Colour limits\nvmin, vmax = -16, -10.5\n\n# Create a figure\nfig, axs = plt.subplots(figsize=(8, 5), nrows=1, ncols=2)\naxs = axs.ravel()\nplt.subplots_adjust(hspace=0.3, wspace=0.3)\n\ntitles = [r'|Real|', r'|Imaginary|']\ndat = [np.log10(np.abs(data.real)), np.log10(np.abs(data.imag))]\n\nfor i in range(2):\n plt.sca(axs[i])\n axs[i].set_title(titles[i])\n axs[i].set_xlim(min(x)/1000, max(x)/1000)\n axs[i].set_ylim(min(x)/1000, max(x)/1000)\n axs[i].axis('equal')\n cs = axs[i].pcolormesh(x/1000, x/1000, dat[i], vmin=vmin, vmax=vmax,\n linewidth=0, rasterized=True, shading='nearest')\n plt.xlabel('Inline Offset (km)')\n plt.ylabel('Crossline Offset (km)')\n\n# Colorbar\n# fig.colorbar(cf0, ax=axs, label=r'$\\log_{10}$ Amplitude (V/m)')\n\n# Plot colorbar\ncax, kw = plt.matplotlib.colorbar.make_axes(\n axs, location='bottom', fraction=.05, pad=0.2, aspect=30)\ncb = plt.colorbar(cs, cax=cax, label=r\"$\\log_{10}$ Amplitude (V/m)\", **kw)\n\n# Title\nfig.suptitle(f\"SEG/EAGE Salt Model, depth = {rz/1e3} km.\", y=1, fontsize=16)\n\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## QC resistivity model with PyVista\n\n

#### Note

The following cell is about how to plot the model in 3D using PyVista,\n for which you have to install pyvista.\n\n The code example was created on 2021-05-21 with pyvista=0.30.1.

\n\n.. code-block:: python\n\n import pyvista\n\n dataset = fgrid.to_vtk({'res': np.log10(fmodel.property_x.ravel('F'))})\n\n # Create the rendering scene and add a grid axes\n p = pyvista.Plotter(notebook=True)\n p.show_grid(location='outer')\n\n dparams = {'rng': np.log10([vmin, vmax]), 'show_edges': False}\n # Add spatially referenced data to the scene\n xyz = (5000, 6000, -3200)\n p.add_mesh(dataset.slice('x', xyz), name='x-slice', **dparams)\n p.add_mesh(dataset.slice('y', xyz), name='y-slice', **dparams)\n p.add_mesh(dataset.slice('z', xyz), name='z-slice', **dparams)\n\n # Get the salt body\n p.add_mesh(dataset.threshold([1.47, 1.48]), name='vol', **dparams)\n\n # Show the scene!\n p.camera_position = [\n (27000, 37000, 5800), (6600, 6600, -3300), (0, 0, 1)\n ]\n p.show()\n\n\n.. figure:: ../../_static/images/SEG-EAGE_3D_salt_model.png\n :scale: 66 %\n :align: center\n :alt: SEG-EAGE 3D salt model with PyVista\n :name: salt_model\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Reproduce the resistivity model\n\n

#### Note

The last cell as about how to reproduce the resistivity model. For this\n you have to download the SEG/EAGE salt model, as explained further down.\n\n The code example and the SEG-EAGE-Salt-Model.h5-file used in the\n gallery were created on 2021-05-21.

\n\nTo reduce runtime and dependencies of the gallery build we use a pre-computed\nresistivity model, which was generated with the code provided below.\n\nIn order to reproduce it yourself you have to download the data from the\nSEG-website\n_ or via this\ndirect link\n_.\nThe zip-file is 513.1 MB big. Unzip the archive, and place the file\nSalt_Model_3D/3-D_Salt_Model/VEL_GRIDS/SALTF.ZIP (20.0 MB) in the same\ndirectory as the notebook.\n\nThe following cell loads takes this SALTF.ZIP, carries out the\nvelocity-to-resistivity transform, and stores the resistivity model for later\nuse.\n\n.. code-block:: python\n\n import emg3d\n import zipfile\n import numpy as np\n\n # Dimension of seismic velocities\n nx, ny, nz = 676, 676, 210\n\n # Create a discretize-mesh of the correct dimension\n # (nz: +1, for air)\n fgrid = emg3d.TensorMesh(\n [np.ones(nx)*20., np.ones(ny)*20., np.ones(nz+1)*20.],\n origin=(0, 0, -210*20))\n res = np.zeros(fgrid.shape_cells, order='F')\n\n # Load data\n zipfile.ZipFile('SALTF.ZIP', 'r').extract('Saltf@@')\n with open('Saltf@@', 'r') as file:\n v = np.fromfile(file, dtype=np.dtype('float32').newbyteorder('>'))\n res[:, :, 1:] = v.reshape((nx, ny, nz), order='F')\n\n # Velocity to resistivity transform for whole cube\n res = (res/1700)**3.88 # Sediment resistivity = 1\n\n # Overwrite basement resistivity from 3660 m onwards\n res[:, :, np.arange(fgrid.shape_cells)*20 > 3680] = 500.\n\n # Set sea-water to 0.3\n res[:, :, :16][res[:, :, :16] <= 1500] = 0.3\n # Ensure at least top layer is water\n res[:, :, 1] = 0.3\n\n # Fix salt resistivity\n res[res == 4482] = 30.\n\n # Set air resistivity\n res[:, :, 0] = 1e8\n\n # THE SEG/EAGE salt-model uses positive z downwards; discretize positive\n # upwards. Hence for res, we use np.flip(res, 2) to flip the z-direction\n res = np.flip(res, 2)\n\n # Create the resistivity model\n model = emg3d.Model(fgrid, property_x=res)\n\n # Store the resistivity model\n emg3d.save(\"SEG-EAGE-Salt-Model.h5\", model=model)\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "emg3d.Report()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.12" } }, "nbformat": 4, "nbformat_minor": 0 }