.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/models/SEG-EAGE_3D_salt_model.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_gallery_models_SEG-EAGE_3D_salt_model.py: SEG-EAGE 3D Salt Model ====================== [Muld07]_ presented electromagnetic responses for a resistivity model which he derived from the seismic velocities of the SEG/EAGE salt model [AmBK97]_. Here we reproduce and store this resistivity model, and compute electromagnetic responses for it. Velocity-to-resistivity transform --------------------------------- Quoting here the description of the velocity-to-resistivity transform used by [Muld07]_: «The SEG/EAGE salt model (Aminzadeh et al. 1997), originally designed for seismic simulations, served as a template for a realistic subsurface model. Its dimensions are 13500 by 13480 by 4680 m. The seismic velocities of the model were replaced by resistivity values. The water velocity of 1.5 km/s was replaced by a resistivity of 0.3 Ohm m. Velocities above 4 km/s, indicative of salt, were replaced by 30 Ohm m. Basement, beyond 3660 m depth, was set to 0.002 Ohm m. The resistivity of the sediments was determined by :math:`(v/1700)^{3.88}` Ohm m, with the velocity v in m/s (Meju et al. 2003). For air, the resistivity was set to :math:`10^8` Ohm m.» Equation 1 of [MeGM03]_, is given by .. math:: :label: meju \log_{10}\rho = m \log_{10}V_P + c \ , where :math:`\rho` is resistivity, :math:`V_P` is P-wave velocity, and for :math:`m` and :math:`c` 3.88 and -11 were used, respectively. The velocity-to-resistivity transform uses therefore a Faust model ([Faus53]_) with some additional constraints for water, salt, and basement. .. note:: The SEG/EAGE Salt Model is licensed under the `CC-BY-4.0 `_. .. GENERATED FROM PYTHON SOURCE LINES 47-60 .. code-block:: Python import os import pooch import emg3d import numpy as np import matplotlib.pyplot as plt from matplotlib.colors import LogNorm plt.style.use('bmh') # Adjust this path to a folder of your choice. data_path = os.path.join('..', 'download', '') .. GENERATED FROM PYTHON SOURCE LINES 61-65 Fetch the model --------------- Retrieve and load the pre-computed resistivity model. .. GENERATED FROM PYTHON SOURCE LINES 65-77 .. code-block:: Python fname = "SEG-EAGE-Salt-Model.h5" pooch.retrieve( 'https://raw.github.com/emsig/data/2021-05-21/emg3d/models/'+fname, '6ee10663de588d445332ba7cc1c0dc3d6f9c50d1965f797425cebc64f9c71de6', fname=fname, path=data_path, ) fmodel = emg3d.load(data_path + fname)['model'] fgrid = fmodel.grid .. rst-class:: sphx-glr-script-out .. code-block:: none Data loaded from «/home/dtr/Codes/emsig/emg3d-gallery/examples/download/SEG-EAGE-Salt-Model.h5» [emg3d v1.0.0rc3.dev5+g0cd9e09 (format 1.0) on 2021-05-21T13:15:50.617756]. .. GENERATED FROM PYTHON SOURCE LINES 78-80 QC resistivity model -------------------- .. GENERATED FROM PYTHON SOURCE LINES 80-92 .. code-block:: Python # Limit colour-range to [0.3, 50] Ohm.m # (affects only the basement and air, improves resolution in the sediments). vmin, vmax = 0.3, 50 fgrid.plot_3d_slicer( fmodel.property_x, zslice=-2000, pcolor_opts={'norm': LogNorm(vmin=vmin, vmax=vmax)} ) .. image-sg:: /gallery/models/images/sphx_glr_SEG-EAGE_3D_salt_model_001.png :alt: SEG EAGE 3D salt model :srcset: /gallery/models/images/sphx_glr_SEG-EAGE_3D_salt_model_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 93-98 Compute some example CSEM data with it -------------------------------------- Survey parameters ````````````````` .. GENERATED FROM PYTHON SOURCE LINES 98-108 .. code-block:: Python # Create a source instance source = emg3d.TxElectricDipole( coordinates=[6400, 6600, 6500, 6500, -50, -50], strength=1/200. # Normalize for length. ) # Frequency (Hz) frequency = 1.0 .. GENERATED FROM PYTHON SOURCE LINES 109-111 Initialize computation mesh ``````````````````````````` .. GENERATED FROM PYTHON SOURCE LINES 111-127 .. code-block:: Python grid = emg3d.construct_mesh( frequency=frequency, properties=[0.3, 1, 2, 15], center=(6500, 6500, -50), seasurface=0, domain=([0, 13500], [0, 13500], None), vector=(None, None, np.array([-100, -80, -60, -40, -20, 0])), min_width_limits=([5, 100], [5, 100], [5, 20]), min_width_pps=5, stretching=[1.03, 1.05], lambda_from_center=True, center_on_edge=True, ) grid .. raw:: html
TensorMesh 2,097,152 cells
MESH EXTENT CELL WIDTH FACTOR
dir nC min max min max max
x 128 -71.48 13,830.69 55.13 195.38 1.02
y 128 -71.48 13,830.69 55.13 195.38 1.02
z 128 -4,579.62 12,968.70 20.00 625.54 1.05


.. GENERATED FROM PYTHON SOURCE LINES 128-130 Put the salt model onto the modelling mesh `````````````````````````````````````````` .. GENERATED FROM PYTHON SOURCE LINES 130-141 .. code-block:: Python # Interpolate full model from full grid to grid model = fmodel.interpolate_to_grid(grid) grid.plot_3d_slicer( model.property_x, zslice=-2000, zlim=(-4180, 500), pcolor_opts={'norm': LogNorm(vmin=vmin, vmax=vmax)} ) .. image-sg:: /gallery/models/images/sphx_glr_SEG-EAGE_3D_salt_model_002.png :alt: SEG EAGE 3D salt model :srcset: /gallery/models/images/sphx_glr_SEG-EAGE_3D_salt_model_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 142-144 Solve the system ```````````````` .. GENERATED FROM PYTHON SOURCE LINES 144-152 .. code-block:: Python efield = emg3d.solve_source( model, source, frequency, semicoarsening=False, linerelaxation=False, verb=1, ) .. rst-class:: sphx-glr-script-out .. code-block:: none :: emg3d :: 4.8e-07; 5(10); 0:00:48; CONVERGED .. GENERATED FROM PYTHON SOURCE LINES 153-163 .. code-block:: Python grid.plot_3d_slicer( efield.fx.ravel('F'), zslice=-2000, zlim=(-4180, 500), view='abs', v_type='Ex', pcolor_opts={'norm': LogNorm(vmin=1e-16, vmax=1e-9)} ) .. image-sg:: /gallery/models/images/sphx_glr_SEG-EAGE_3D_salt_model_003.png :alt: SEG EAGE 3D salt model :srcset: /gallery/models/images/sphx_glr_SEG-EAGE_3D_salt_model_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 164-208 .. code-block:: Python # Interpolate for a "more detailed" image x = grid.cell_centers_x y = grid.cell_centers_y rx = np.repeat([x, ], np.size(x), axis=0) ry = rx.transpose() rz = -2000 data = efield.get_receiver((rx, ry, rz, 0, 0)) # Colour limits vmin, vmax = -16, -10.5 # Create a figure fig, (ax1, ax2) = plt.subplots( 1, 2, figsize=(8, 5), sharex=True, constrained_layout=True) titles = [r'|Real|', r'|Imaginary|'] def log10abs(inp): """Log10 of absolute values, avoiding zero-division.""" tiny = np.finfo(float).tiny return np.log10(np.where(abs(inp) < tiny, tiny, abs(inp))) dat = [log10abs(data.real), log10abs(data.imag)] for i, ax in enumerate([ax1, ax2]): ax.set_title(titles[i]) cs = ax.pcolormesh(x/1000, x/1000, dat[i], vmin=vmin, vmax=vmax, linewidth=0, rasterized=True, shading='nearest') ax.set_xlim([min(x)/1000, max(x)/1000]) ax.axis('equal') ax.set_xlabel('Inline Offset (km)') ax.set_ylabel('Crossline Offset (km)') # Plot colorbar cb = plt.colorbar(cs, ax=[ax1, ax2], orientation='horizontal', aspect=30) cb.set_label("log10 Amplitude (V/m)") # Title fig.suptitle(f"SEG/EAGE Salt Model, depth = {rz/1e3} km.", fontsize=16) .. image-sg:: /gallery/models/images/sphx_glr_SEG-EAGE_3D_salt_model_004.png :alt: SEG/EAGE Salt Model, depth = -2.0 km., |Real|, |Imaginary| :srcset: /gallery/models/images/sphx_glr_SEG-EAGE_3D_salt_model_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 209-251 QC resistivity model with PyVista --------------------------------- .. note:: The following cell is about how to plot the model in 3D using PyVista, for which you have to install ``pyvista``. The code example was created on 2021-05-21 with ``pyvista=0.30.1``. .. code-block:: python import pyvista dataset = fgrid.to_vtk({'res': np.log10(fmodel.property_x.ravel('F'))}) # Create the rendering scene and add a grid axes p = pyvista.Plotter(notebook=True) p.show_grid(location='outer') dparams = {'rng': np.log10([vmin, vmax]), 'show_edges': False} # Add spatially referenced data to the scene xyz = (5000, 6000, -3200) p.add_mesh(dataset.slice('x', xyz), name='x-slice', **dparams) p.add_mesh(dataset.slice('y', xyz), name='y-slice', **dparams) p.add_mesh(dataset.slice('z', xyz), name='z-slice', **dparams) # Get the salt body p.add_mesh(dataset.threshold([1.47, 1.48]), name='vol', **dparams) # Show the scene! p.camera_position = [ (27000, 37000, 5800), (6600, 6600, -3300), (0, 0, 1) ] p.show() .. figure:: ../../_static/images/SEG-EAGE_3D_salt_model.png :scale: 66 % :align: center :alt: SEG-EAGE 3D salt model with PyVista :name: salt_model .. GENERATED FROM PYTHON SOURCE LINES 254-329 Reproduce the resistivity model ------------------------------- .. note:: The last cell as about how to reproduce the resistivity model. For this you have to download the SEG/EAGE salt model, as explained further down. The code example and the ``SEG-EAGE-Salt-Model.h5``-file used in the gallery were created on 2021-05-21. To reduce runtime and dependencies of the gallery build we use a pre-computed resistivity model, which was generated with the code provided below. In order to reproduce it yourself you have to download the data from the `SEG-website `_ or via this `direct link `_. The zip-file is 513.1 MB big. Unzip the archive, and place the file ``Salt_Model_3D/3-D_Salt_Model/VEL_GRIDS/SALTF.ZIP`` (20.0 MB) in the same directory as the notebook. The following cell loads takes this ``SALTF.ZIP``, carries out the velocity-to-resistivity transform, and stores the resistivity model for later use. .. code-block:: python import emg3d import zipfile import numpy as np # Dimension of seismic velocities nx, ny, nz = 676, 676, 210 # Create a discretize-mesh of the correct dimension # (nz: +1, for air) fgrid = emg3d.TensorMesh( [np.ones(nx)*20., np.ones(ny)*20., np.ones(nz+1)*20.], origin=(0, 0, -210*20)) res = np.zeros(fgrid.shape_cells, order='F') # Load data zipfile.ZipFile('SALTF.ZIP', 'r').extract('Saltf@@') with open('Saltf@@', 'r') as file: v = np.fromfile(file, dtype=np.dtype('float32').newbyteorder('>')) res[:, :, 1:] = v.reshape((nx, ny, nz), order='F') # Velocity to resistivity transform for whole cube res = (res/1700)**3.88 # Sediment resistivity = 1 # Overwrite basement resistivity from 3660 m onwards res[:, :, np.arange(fgrid.shape_cells[2])*20 > 3680] = 500. # Set sea-water to 0.3 res[:, :, :16][res[:, :, :16] <= 1500] = 0.3 # Ensure at least top layer is water res[:, :, 1] = 0.3 # Fix salt resistivity res[res == 4482] = 30. # Set air resistivity res[:, :, 0] = 1e8 # THE SEG/EAGE salt-model uses positive z downwards; discretize positive # upwards. Hence for res, we use np.flip(res, 2) to flip the z-direction res = np.flip(res, 2) # Create the resistivity model model = emg3d.Model(fgrid, property_x=res) # Store the resistivity model emg3d.save("SEG-EAGE-Salt-Model.h5", model=model) .. GENERATED FROM PYTHON SOURCE LINES 331-333 .. code-block:: Python emg3d.Report() .. raw:: html
Thu Jul 04 15:42:23 2024 CEST
OS Linux (Ubuntu 22.04) CPU(s) 16 Machine x86_64
Architecture 64bit RAM 31.0 GiB Environment Python
File system ext4
Python 3.12.4 | packaged by conda-forge | (main, Jun 17 2024, 10:23:07) [GCC 12.3.0]
numpy 1.26.4 scipy 1.14.0 numba 0.60.0
emg3d 1.8.3 empymod 2.3.1 xarray 2024.6.1.dev31+g6c2d8c33
discretize 0.10.0 h5py 3.11.0 matplotlib 3.8.4
tqdm 4.66.4
Intel(R) oneAPI Math Kernel Library Version 2023.2-Product Build 20230613 for Intel(R) 64 architecture applications


.. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 53.631 seconds) **Estimated memory usage:** 2873 MB .. _sphx_glr_download_gallery_models_SEG-EAGE_3D_salt_model.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: SEG-EAGE_3D_salt_model.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: SEG-EAGE_3D_salt_model.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_