.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/fdomain/magnetotelluric.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_fdomain_magnetotelluric.py: Magnetotelluric data ==================== The magnetotelluric (MT) method is a passive method using as a source variations in Earth's magnetic field, which create telluric currents in the Earth. The variation of Earth's magnetic field has many origins, e.g., lightning or the interaction between the Earth's magnetic field and solar wind. Ubiquitous in MT is the plane-wave approximation, hence, that the source signal is a plane wave hitting the Earth's surface. Having a 1D source (vertical plane wave) for a layered, 1D model reduces the computation of the impedances and from there the apparent resistivity and apparent phase a simple recursion algorithm. As such it does not make sense to use a full EM wavefield algorithm with three-dimensional sources such as empymod to compute MT responses. However, it is still interesting to see if we can compute MT impedances with a three-dimensional point source. As background theory we reproduce here Equations (11) to (17) from Pedersen and Hermance (1986), with surrounding text. For a more in-depth read we refer to Chave and Jones (2012). -------- If we define the impedance as .. math:: :label: ph-11 Z = E_x / H_y \, , [...] we can develop a recursive relationship for the impedance at the top of the j-th layer looking down .. math:: :label: ph-12 Z_j = z_{oj} \frac{1-R_j \exp(-2\gamma_j t_j)}{1+R_j \exp(-2\gamma_j t_j)} \, , \quad j = N-1, \dots, 1 \, , where .. math:: :label: ph-13 z_{oj} \equiv \text{intrinsic impedance} \equiv \sqrt{\rm{i}\omega \mu \rho_j} \, , .. math:: :label: ph-14 R_j \equiv \text{reflection coefficient} \equiv \frac{z_{oj} - Z_{j+1}}{z_{oj} + Z_{j+1}} \, , .. math:: t_j = \text{thickness of layer} j \, , and the impedance at the surface of the deepest layer (:math:j=N) is given by .. math:: :label: ph-15 Z_N = z_{oN}\, . The surface impedance, :math:Z_j, is found by applying Equation :eq:ph-12 recursively from the top of the bottom half-space, :math:j = N, and propagating upwards. From the surface impedance, :math:Z_1, we can then calculate the apparent resistivity, :math:\rho_a, and phase, :math:\theta_a, as .. math:: :label: ph-16 \rho_a = \frac{|Z_1|^2}{\omega \mu} \, , .. math:: :label: ph-17 \theta_a = \tan^{-1}\frac{\Im(Z1)}{\Re(Z_1)} \, . This calculation is repeated for a range of periods and is used to model the magnetotelluric response of the layered structure. -------- **Reference**: - Chave, A., and Jones, A. (Eds.). (2012). The Magnetotelluric Method: Theory and Practice. Cambridge: Cambridge University Press; https://doi.org/10.1017/CBO9781139020138. - Pedersen, J., Hermance, J.F. Least squares inversion of one-dimensional magnetotelluric data: An assessment of procedures employed by Brown University. Surv Geophys 8, 187–231 (1986); https://doi.org/10.1007/BF01902413. - This example was strongly motivated by Andrew Pethicks blog post https://www.digitalearthlab.com/tutorial/tutorial-1d-mt-forward. .. GENERATED FROM PYTHON SOURCE LINES 103-109 .. code-block:: default import empymod import numpy as np import matplotlib.pyplot as plt from scipy.constants import mu_0 plt.style.use('ggplot') .. GENERATED FROM PYTHON SOURCE LINES 110-112 Define model parameter and frequencies -------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 112-118 .. code-block:: default resistivities = np.array([2e14, 300, 2500, 0.8, 3000, 2500]) depths = np.array([0, 200, 600, 640, 1140]) frequencies = np.logspace(-4, 5, 101) omega = 2 * np.pi * frequencies .. GENERATED FROM PYTHON SOURCE LINES 119-123 1D-MT recursion following Pedersen & Hermance ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Using the variable names as in the paper. .. GENERATED FROM PYTHON SOURCE LINES 123-150 .. code-block:: default # Initiate recursive formula with impedance of the deepest layer. Z_j = np.sqrt(2j * np.pi * frequencies * mu_0 * resistivities[-1]) # Move up the stack of layers till the top (without air). for j in range(depths.size-1, 0, -1): # Thickness t_j = depths[j] - depths[j-1] # Intrinsic impedance z_oj = np.sqrt(1j * omega * mu_0 * resistivities[j]) # Reflection coefficient R_j = (z_oj - Z_j) / (z_oj + Z_j) # Exponential factor gamma_j = np.sqrt(1j * omega * mu_0 / resistivities[j]) exp_j = np.exp(-2 * gamma_j * t_j) # Impedance at this layer Z_j = z_oj * (1 - R_j * exp_j) / (1 + R_j * exp_j) # Step 3. Compute apparent resistivity last impedance apres_mt1d = abs(Z_j)**2/(omega * mu_0) phase_mt1d = np.arctan2(Z_j.imag, Z_j.real) .. GENERATED FROM PYTHON SOURCE LINES 151-161 1D MT using empymod ------------------- The above derivation and code assume plane waves. We can "simulate" plane waves by putting the source _very_ far away. In this example, we set it one million km away in all directions. As the impedance is the ratio of the Ex and Hy fields the source type (electric or magnetic) and the source orientation do not matter; hence, it can be an arbitrarily rotated electric or magnetic source, the apparent resistivity and apparent phase will always be the same. .. GENERATED FROM PYTHON SOURCE LINES 161-183 .. code-block:: default dist = 1_000_000_000 # 1 million kilometer (!) inp = { 'src': (-dist, -dist, -dist), 'rec': (0, 0, 0.1), 'res': resistivities, 'depth': depths, 'freqtime': frequencies, 'verb': 1, } # Get Ex, Hy. ex = empymod.dipole(ab=11, **inp) hy = empymod.dipole(ab=51, **inp) # Impedance. Z = ex/hy # Apparent resistivity and apparent phase. apres_empy = abs(Z)**2 / (omega * mu_0) phase_empy = np.arctan2(Z.imag, Z.real) .. GENERATED FROM PYTHON SOURCE LINES 184-187 Plot results ------------ .. GENERATED FROM PYTHON SOURCE LINES 187-211 .. code-block:: default fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 4)) ax1.set_title('Apparent resistivity') ax1.plot(frequencies, apres_mt1d, label='MT-1D') ax1.plot(frequencies, apres_empy, '--', label='empymod') ax1.set_xscale('log') ax1.set_yscale('log') ax1.set_xlabel('Frequency (Hz)') ax1.set_ylabel(r'Apparent resistivity ($\Omega\,$m)') ax1.legend() ax2.set_title('Phase') ax2.plot(frequencies, phase_mt1d*180/np.pi) ax2.plot(frequencies, phase_empy*180/np.pi, '--') ax2.set_xscale('log') ax2.yaxis.tick_right() ax2.set_xlabel('Frequency (Hz)') ax2.set_ylabel('Phase (degree)') ax2.yaxis.set_label_position("right") fig.tight_layout() fig.show() .. image-sg:: /gallery/fdomain/images/sphx_glr_magnetotelluric_001.png :alt: Apparent resistivity, Phase :srcset: /gallery/fdomain/images/sphx_glr_magnetotelluric_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 212-213 .. code-block:: default empymod.Report() .. raw:: html
 Wed May 31 09:36:18 2023 UTC OS Linux CPU(s) 2 Machine x86_64 Architecture 64bit RAM 7.6 GiB Environment Python File system ext4 Python 3.8.6 (default, May 10 2023, 17:44:10) [GCC 7.5.0] numpy 1.24.3 scipy 1.10.1 numba 0.57.0 empymod 2.2.2 IPython 8.12.2 matplotlib 3.7.1

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