<jupyter_start><jupyter_text>Forward simulation the 1D Magnetotelluric (MT) problemIn the [previous notebook](./MT1D_Simulation.ipynb), we walked through how to discretize and solve the 1D Magnetotelluric (MT) problem using a finite difference approach. In this notebook, we will use the numerical simulation to simulate MT data and explore concepts including- **Mesh Design**: we will compare our numerical result to an analytic and see what happens when cell sizes are too large, or the domain doesn't extend far enough- **Non-uniqueness**: prior to inverting geophysical data, it is important to set expectations on what we hope to recover from those data. Forward modelling is a powerful tool for getting a handle on this.<jupyter_code>import numpy as np
import scipy.sparse as sp
import matplotlib.pyplot as plt
from SimPEG import Mesh, Utils, Solver
from scipy.constants import mu_0, epsilon_0
import matplotlib
matplotlib.rcParams["font.size"] = 13
%matplotlib inline<jupyter_output><empty_output><jupyter_text>In the [previous notebook](./MT1D_Simulation), we generated the function `simulateMT` in the file `MTforward.py`, we will import that and use it in this notebook.<jupyter_code>from MTforward import simulateMT<jupyter_output><empty_output><jupyter_text>Set up Model ParametersTo define an MT simulation, we will set up a conductivity model a frequency range over which we wish to simulate. We will use:- a background resistivity of $\rho = 100 \Omega m$ ($10^{-2} Sm$)- 25 frequencies between $10^{-3}$ Hz and 100 Hz<jupyter_code>rho_half = 100.  # Resistivity of the halfspace in Ohm-m
sigma_half = 1./rho_half  # Conductivity is the inverse of conductivity
frequency = np.logspace(-3, 2, 25)  # frequencies at which to simulate the MT problem<jupyter_output><empty_output><jupyter_text>Skin depth and Mesh DesignWhen setting up a mesh, we want to make sure our fine cells are fine enough to capture the behaviour at the highest frequencies and that the domain extends far enough so that the fields have sufficiently decayed by the time they reach the boundary. To gauge this, we will examine the [skin depth](http://em.geosci.xyz/content/maxwell1_fundamentals/plane_waves_in_homogeneous_media/frequency/analytic_solution.htmlattenuation-and-skin-depth) at the highest and lowest frequencies. Skin depth ($\delta$) is the distance at which the amplitude of an EM wave propagating through a homogeneous medium will have decayed by a factor of $1/e$**Skin Depth**$$\delta = \frac{500}{\sqrt{\sigma f}}$$<jupyter_code>def skin_depth(sigma, f):
    """
    Depth at which the fields propagating through a homogeneous 
    medium have decayed by a factor of 1/e for a given 
    frequency, f and conductivity, sigma
    """
    return 500./np.sqrt(sigma * f)
skin_depth_min = skin_depth(sigma_half, frequency.max())
skin_depth_max = skin_depth(sigma_half, frequency.min())

print("The minimum skin depth is {:1.2f}m".format(skin_depth_min))
print("The maximum skin depth is {:1.2e}m".format(skin_depth_max))<jupyter_output>The minimum skin depth is 500.00m
The maximum skin depth is 1.58e+05m<jupyter_text>Mesh parametersWe start by choosing a mesh with parameters as follows: - smallest cell size: $$\Delta z_{\rm core} = \delta_{\rm min} / 4$$- domain extent: $$z_{\rm max} = 2 \delta_{\rm max} $$For this example, we will be exploring model variations within the top 5km, so we will extend the core region (the region of the having uniform cells with width $\Delta z_{\rm core}$) of the mesh to a depth of 5km.<jupyter_code>cs = skin_depth_min / 4.
core_extent = 5000. 
domain_extent = 2 * skin_depth_max

print("The smallest cell size is {:1.2f}m".format(cs))
print("The core region of the mesh extends {:1.2e}m".format(core_extent))
print("The mesh should extend at least {:1.2e}m".format(domain_extent))<jupyter_output>The smallest cell size is 125.00m
The core region of the mesh extends 5.00e+03m
The mesh should extend at least 3.16e+05m<jupyter_text>PaddingWe will use a tensor mesh, which means we can use non-uniform cells for the padding, that is, expanding the width of the cells with depth. We can get away with this because EM fields and fluxes are diffusive and high frequencies will be attenuated as they move through the conductive earth (eg. in seismic, you wouldn't want to do this as the seismic response is  dominated by wave propagaion). We expand by a factor of 1.3 until we are beyond the desired domain extent. Here, we write a small `while` loop to figure out how many padding cells we should use.<jupyter_code>npad = 1  # start with 1 cell
padding_fact = 1.3  # the amount by which we will expand each cell of the padding

def padding_extent(npad):
    """
    given a number of padding cells, this computes how far 
    the padding extends
    """
    padding_widths = cs*padding_fact**(np.arange(npad) + 1)
    return padding_widths.sum()

# keep adding padding until we are beyond the desired extent
padding_z = padding_extent(npad)
while padding_z < domain_extent:
    npad+=1
    padding_z = padding_extent(npad)
    
print(
    "{:1.0f} padding cells extends {:1.2e}m > {:1.2e}m "
    "(2 skin depths)".format(
        npad, padding_extent(npad), domain_extent
    )
)<jupyter_output>25 padding cells extends 3.82e+05m > 3.16e+05m (2 skin depths)<jupyter_text>Construct the Tensor meshNow that we have defined all of the mesh parameters, we use the [Mesh Class in SimPEG](http://docs.simpeg.xyz) to construct a mesh. This will define all of the geometries, provide the necessary differential operators, and some handy functions like plotting.<jupyter_code>ncz = np.ceil(core_extent / cs)  # number of cells in the core domain
hz = [(cs, npad, -1.3), (cs, ncz)]  # define how to construct the cell widths
mesh = Mesh.TensorMesh([hz], x0='N')  # construct a 1D Tensor Mesh

print(
    "There are {:1.0f} cells in the mesh. The mest extends {:1.2e}m".format(
        ncz, mesh.hx.sum()
    )
) 
# plot the mesh
fig, ax = plt.subplots(1,1, figsize=(8, 3))
mesh.plotGrid(centers=True, faces=True, ax=ax)
ax.legend(["centers", "faces"])
ax.grid(which="both", linewidth=0.5)
ax.invert_xaxis()  # so that the surface is on our left hand side
ax.set_xlabel('z (m)')<jupyter_output><empty_output><jupyter_text>Set up a modelIn the [previous notebook](./MT1D_Simulation.ipynb), we showed that for a half-space, we expect the apparent resistivity computed from impedance data to be the same as the half-space resistivity, and the phase to by $45^\circ$ across the entire frequency range. What happens when we include a conductive target in our model?Here we will use the following model parameters- Resistivity of the background: $\rho_{\text{halfspace}} = 100 \Omega m$ (which is the same as $\sigma_{\text{halfspace}} = 10^{-2} S/m$)- Resistivity of the target: $\rho_{\text{target}} = 10 \Omega m$ (which is the same as $\sigma_{\text{halfspace}} = 10^{-1} S/m$)- Thickness of the target layer: 1000m - Depth of the target layer: 2000m<jupyter_code>rho_target = 10.  # resistivity in Ohm-m
depth = 2000.  # depth to the top of the target in m 
thickness = 1000.  # thickness of the target in m
# put the model on the mesh
sigma = 1./rho_half * np.ones(mesh.nC)

# find the indices of the layer
layer_inds = (
    (mesh.vectorCCx<=-depth) & 
    (mesh.vectorCCx>-(depth+thickness))
)
sigma[layer_inds] = 1./rho_target
# plot the model
fig, ax = plt.subplots(1, 1, figsize=(8, 3))

# trickery to plot from node to node rather than at cell centers
z = np.repeat(mesh.vectorNx[1:-1], 2, axis=0)
z = np.r_[mesh.vectorNx[0], z, mesh.vectorNx[-1]]
sigma_plt = np.repeat(sigma, 2, axis=0)

ax.semilogy(z, sigma_plt,"C0", lw=2)
ax.grid(which="both", linewidth=0.5)
ax.set_xlim([-5000., 0.])
ax.set_ylim([5e-3, 1])
ax.invert_xaxis() # plot the surface on the left

ax.set_xlabel("Elevation (m)", fontsize=14)
ax.set_ylabel("Conductivity (S/m)", fontsize=14)<jupyter_output><empty_output><jupyter_text>Compare with an analyticPrior to drawing conclusions from our simulations, we want to first make sure they are correct! So here, we compare the numerical results computed with our method `simulateMT` with an analytic for the MT response over a layered earth. A more complete description, including the derivation, is available on [EM GeoSci](http://em.geosci.xyz/content/maxwell3_fdem/natural_sources/MT_N_layered_Earth.html)<jupyter_code>from SimPEG.EM.Analytics import MT_LayeredEarth
# the analytic takes the frequencies, layer thicknesses and layer conductivities
sigma_layers = np.r_[
    1./rho_half, 
    1./rho_target, 
    1./rho_half]
h = np.r_[depth, thickness]  

app_res_ana, app_phase_ana = MT_LayeredEarth(
    frequency, h, sigma_layers, 'Res-Phase'
)
# # Uncomment to read the documentation on MT_LayeredEarth
# MT_LayeredEarth??
# numerically compute the response
app_res, app_phase = simulateMT(mesh, sigma, frequency)
def plot_with_analytic(frequency, app_res, app_phase, app_res_ana, app_phase_ana):
    # Plot and compare the results
    fig, ax = plt.subplots(2, 1, figsize=(8, 3*2))

    # apparent resistivity
    ax[0].loglog(frequency, app_res, label='Numeric')
    ax[0].loglog(frequency, app_res_ana, 'k.', label='Analytic')
    ax[0].set_ylabel("$\\rho_a \ (\Omega m)$", fontsize=14)
    ax[0].set_ylim([2e1, 3e2])

    # phase
    ax[1].semilogx(frequency, app_phase, label='Numeric')
    ax[1].semilogx(frequency, app_phase_ana, 'k.', label='Analytic')
    ax[1].set_ylabel("$\phi \ (^{\circ})$", fontsize=13)
    ax[1].set_ylim([0., 90.])

    for a in ax:
        a.grid(True, which='both', linewidth=0.3)
        a.set_xlim(frequency.max(), frequency.min())
        a.set_xlabel("Frequency (Hz)", fontsize=14)
        a.legend(fontsize=11)

    plt.tight_layout()
    plt.show()
plot_with_analytic(frequency, app_res, app_phase, app_res_ana, app_phase_ana)<jupyter_output><empty_output><jupyter_text>Across the entire frequency band, the apparent resistivity and phase agree with the analytic. What happens if1. we use a coarse cell size?2. we reduce how far the padding extends?To answer those questions, we have a small widget.<jupyter_code>import ipywidgets
def mesh_design_app(domain_extent, cs):
    ncz = np.ceil(core_extent / cs)  # number of cells in the core domain
    
    # keep adding padding until we are beyond the desired extent
    npad=0
    padding_z = padding_extent(npad)
    while padding_z < domain_extent*1e3:
        npad+=1
        padding_z = padding_extent(npad)
    
    hz = [(cs, npad, -1.3), (cs, ncz)]  # define how to construct the cell widths
    mesh = Mesh.TensorMesh([hz], x0='N')  # construct a 1D Tensor Mesh

    # put the model on the mesh
    sigma = 1./rho_half * np.ones(mesh.nC)

    # find the indices of the layer
    layer_inds = (
        (mesh.vectorCCx<=-depth) & 
        (mesh.vectorCCx>-(depth+thickness))
    )
    sigma[layer_inds] = 1./rho_target
    
    app_res, app_phase = simulateMT(mesh, sigma, frequency)
    
    plot_with_analytic(frequency, app_res, app_phase, app_res_ana, app_phase_ana)
print(
    "It is helpful to remember that \n"
    "the minimum skin depth is {:1.0f}m \n"
    "and the maximum skin depth is {:1.0f}km".format(
        skin_depth_min, skin_depth_max*1e-3
    )
)<jupyter_output>It is helpful to remember that 
the minimum skin depth is 500m 
and the maximum skin depth is 158km<jupyter_text>- What happens if you make the domain smaller? Which frequencies are impacted first?- What happens if you use a larger cell size? (Pay attention to the phase at high frequencies)    - notice that there will be two regions impacted: high frequencies, and the frequencies impacted by the conductive layer<jupyter_code>skin_depth_min
ipywidgets.interact(
    mesh_design_app,
    domain_extent=ipywidgets.FloatSlider(min=0, max=400, step=25, value=300, description="domain extent (km)"),
    cs=ipywidgets.FloatSlider(min=25, max=1000, step=25, value=125, description="min cell size (m)")
)<jupyter_output><empty_output><jupyter_text>Considerations for the inversionPrior to considering an inversion, first making sure that a suitable mesh is designed is an important step - otherwise, you can introduce artifacts and structures simply because the mesh doesn't capture the physics! Comparing against analytics that model a simplified, but similar, scenario can be a powerful tool for mesh design.   Non-uniquenessThe second concept we will explore with forward modelling is non-uniqueness: when models are different, but the data are similar.In the Magnetotelluric problem, a classic example that demonstrates non-uniqueness of the data is the non-uniquness associated with the conductivity-thickness product of a thin layer. If we start with a layer that has a conductivity of $\sigma$, halve its thickness and double its conductivity, the resulting data will be similar... (but don't just take our word for it!)<jupyter_code># compute models with equivalent conductivity-thickness product
conductivity_thickness = thickness * 1./rho_target
layer_thicknesses = (np.arange(5)+1.)*300.
print("we consider layers having thicknesses (m) {}".format(layer_thicknesses))
sig_targets = conductivity_thickness / layer_thicknesses
print(
    "To preserve the conductivity thickness product,\n"
    "these layer have conductivities (S/m) of {}".format(sig_targets)
)
print("This is the same as resistivities (Ohm m) of {}".format(1./sig_targets))
# construct conductivity models on the mesh
sigma_list = []

# Compute the midpoint and make sure the layers we will compare with are all centered at the same depth 
depth_mid = depth + thickness / 2.

for i in range(len(layer_thicknesses)):
    layer_inds = (
        (mesh.vectorCCx<-depth_mid+layer_thicknesses[i]*0.5) & 
        (mesh.vectorCCx>-depth_mid-layer_thicknesses[i]*0.5)
    )
    sigma_temp = np.ones(mesh.nC) * 1./rho_half
    sigma_temp[layer_inds] = sig_targets[i]
    sigma_list.append(sigma_temp)
# compute the apparent resistivity and phase for each model
app_res_list = []
phase_list = []
for i in range(len(layer_thicknesses)):
    app_res, phase = simulateMT(mesh, sigma_list[i], frequency)
    app_res_list.append(app_res)
    phase_list.append(phase)
def plot_apparent_resistivities_and_phases(
    frequency, sigma_list, layer_thicknesses, app_res_list, phase_list, app_res_uncert=0., phase_uncert=0.
):
    """
    Plot the apparent resistivity and phase for multiple models
    
    frequency: numpy array of frequencies
    sigma_list: list of conductivity models to plot
    layer_thicknesses: thickness of the ~equivalent layers
    app_res_list: list of apparent resistivities computed from the conductivity models
    phase_list: list of phases computed fromt the conductivity models 
    app_res_uncert: uncertainty (in percentage) of the apparent resistivity
    phase_uncert: uncertainty (in degrees) of the phase  
    """
    fig, ax = plt.subplots(3, 1, figsize=(8, 3*3), dpi=400)

    for i in range(len(sigma_list)):
        # trickery to plot from node to node rather than at cell centers
        z = np.repeat(mesh.vectorNx[1:-1], 2, axis=0)
        z = np.r_[mesh.vectorNx[0], z, mesh.vectorNx[-1]]
        sigma_plt = np.repeat(sigma_list[i], 2, axis=0)

        # noise models, apparent resistivity is a percent, phase in units of degrees
        app_res_noise = (app_res_uncert/100.)*np.abs(app_res_list[i])*np.random.randn(len(app_res_list[i]))
        phase_noise = phase_uncert*np.random.randn(len(app_res_list[i]))

        ax[0].plot(z, sigma_plt, 'C{}'.format(i), lw=2, label="{:1.0f} m".format(layer_thicknesses[i]))
        ax[1].loglog(frequency, app_res_list[i]+app_res_noise, 'C{}-'.format(i), lw=1, label="{:1.0f} m".format(layer_thicknesses[i]))
        ax[2].semilogx(frequency, phase_list[i]+phase_noise, 'C{}-'.format(i), lw=1, label="{:1.0f} m".format(layer_thicknesses[i]))   
    # conductivity plot
    ax[0].set_ylabel("Conductivity (S/m)", fontsize=13)
    ax[0].set_xlabel("Depth (m)", fontsize=14)
    ax[0].set_xlim([-5000., 0.])
    ax[0].set_ylim([0, 0.4])
    ax[0].legend(fontsize=10, bbox_to_anchor=(1.18, 1.01))

    # apparent resistivity plot
    ax[1].set_ylabel("Apparent Resistivity $(\Omega m)$", fontsize=13)
    ax[1].set_ylim([2e1, 2e2])

    # phase plot
    ax[2].set_ylabel("Phase $(^{\circ})$", fontsize=13) 
    ax[2].set_ylim([0., 90.])

    # apparent resistivity and phase plots
    for a in ax[1:]:
        a.set_xlabel("Frequency (Hz)", fontsize=13)
        a.set_xlim([frequency.min(), frequency.max()])

    for a in ax:
        a.grid(True, which='both', linewidth=0.3)
        a.invert_xaxis()
        

    plt.tight_layout()
    plt.show()
plot_apparent_resistivities_and_phases(frequency, sigma_list, layer_thicknesses, app_res_list, phase_list)<jupyter_output><empty_output><jupyter_text>What happens if we add noise??<jupyter_code>ipywidgets.interact(
    lambda app_res_uncert, phase_uncert: plot_apparent_resistivities_and_phases(
        frequency, sigma_list, layer_thicknesses, 
        app_res_list, phase_list, app_res_uncert, phase_uncert
    ),
    app_res_uncert = ipywidgets.FloatSlider(
        min=0., max=20., step=1, value=0., description="$\delta \\rho_a (\%)$"
    ),
    phase_uncert = ipywidgets.FloatSlider(
        min=0., max=10., step=0.5, value=0., description="$\delta \phi (^\circ)$" 
    )
)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1063044356001615, 'fraction_numerical': 0.021226279056353464, 'mean_word_length': 4.4113607990012484, 'pattern_counts': {'":': 0, '<': 54, '<?xml version=': 0, '>': 54, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '35220', 'license': 'MIT', 'n_tokens_mistral': 5923, 'n_tokens_neox': 5507, 'n_words': 1856, 'path': 'notebooks/MT_tutorial_2_MT1D_ForwardModellingAndNonuniqueness.ipynb', 'repo_name': 'simpeg/em-notebooks', 'size': '17337'}