notebooks/dcip/DC_SurveyDataInversion.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/syckdog/opt/anaconda3/envs/e350_notebook/lib/python3.9/site-packages/discretize/utils/code_utils.py:247: FutureWarning: ExtractCoreMesh has been deprecated, please use extract_core_mesh. It will be removed in version 1.0.0 of discretize.\n",
      "  warnings.warn(\n"
     ]
    }
   ],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "from geoscilabs.dcip.DC_cylinder import cylinder_app\n",
    "from geoscilabs.dcip.DCLayers import plot_layer_potentials_app\n",
    "from geoscilabs.dcip.DC_Pseudosections import (\n",
    "    DC2DPseudoWidget, MidpointPseudoSectionWidget, DC2DfwdWidget\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Understanding currents, fields, charges and potentials"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cylinder app\n",
    "\n",
    " - **survey**: Type of survey\n",
    " - **A**: (+) Current electrode  location\n",
    " - **B**: (-) Current electrode  location\n",
    " - **M**: (+) Potential electrode  location\n",
    " - **N**: (-) Potential electrode  location\n",
    " - **r**: radius of cylinder\n",
    " - **xc**: x location of cylinder center\n",
    " - **zc**: z location of cylinder center\n",
    " - **$\\rho_1$**: Resistivity of the halfspace\n",
    " - **$\\rho_2$**: Resistivity of the cylinder\n",
    " - **Field**: Field to visualize\n",
    " - **Type**: which part of the field\n",
    " - **Scale**: Linear or Log Scale visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "04569e8a845b4a32bc77d043b2950dd3",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "MyApp(children=(ToggleButtons(description='survey', options=('Dipole-Dipole', 'Dipole-Pole', 'Pole-Dipole', 'P…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "cylinder_app()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Potential differences and Apparent Resistivities\n",
    "\n",
    "Using the widgets contained in this notebook you will develop a better understand of what values are actually measured in a DC resistivity survey and how these measurements can be processed, plotted, inverted, and interpreted.\n",
    "\n",
    "## Computing Apparent Resistivity\n",
    "\n",
    "In practice we cannot measure the potentials everywhere, we are limited to those locations where we place electrodes. For each source (current electrode pair) many potential differences are measured between M and N electrode pairs to characterize the overall distribution of potentials. The widget below allows you to visualize the potentials, electric fields, and current densities from a dipole source in a simple model with 2 layers. For different electrode configurations you can measure the potential differences and see the calculated apparent resistivities. \n",
    "\n",
    "In a uniform halfspace the potential differences can be computed by summing up the potentials at each measurement point from the different current sources based on the following equations:\n",
    "\n",
    "\\begin{align}\n",
    "    V_M = \\frac{\\rho I}{2 \\pi} \\left[ \\frac{1}{AM} - \\frac{1}{MB} \\right] \\\\\n",
    "    V_N = \\frac{\\rho I}{2 \\pi} \\left[ \\frac{1}{AN} - \\frac{1}{NB} \\right] \n",
    "\\end{align} \n",
    "where $AM$, $MB$, $AN$, and $NB$ are the distances between the corresponding electrodes. \n",
    "\n",
    "The potential difference $\\Delta V_{MN}$ in a dipole-dipole survey can therefore be expressed as follows,\n",
    "\n",
    "\\begin{equation}\n",
    "    \\Delta V_{MN} = V_M - V_N = \\rho I \\underbrace{\\frac{1}{2 \\pi} \\left[ \\frac{1}{AM} - \\frac{1}{MB} - \\frac{1}{AN} + \\frac{1}{NB} \\right]}_{G}\n",
    "\\end{equation}\n",
    "\n",
    "and the resistivity of the halfspace $\\rho$ is equal to,\n",
    "\n",
    "$$\n",
    "    \\rho = \\frac{\\Delta V_{MN}}{IG}\n",
    "$$\n",
    "\n",
    "In this equation $G$ is often referred to as the geometric factor. \n",
    "\n",
    "In the case where we are not in a uniform halfspace the above equation is used to compute the apparent resistivity ($\\rho_a$) which is the resistivity of the uniform halfspace which best reproduces the measured potential difference.\n",
    "\n",
    "In the top plot the location of the A electrode is marked by the red +, the B electrode is marked by the blue -, and the M/N potential electrodes are marked by the black dots. The $V_M$ and $V_N$ potentials are printed just above and to the right of the black dots. The calculted apparent resistivity is shown in the grey box to the right. The bottom plot can show the resistivity model, the electric fields (e), potentials, or current densities (j) depending on which toggle button is selected. Some patience may be required for the plots to update after parameters have been changed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Two layer app"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- **A**: (+) Current electrode  location\n",
    "- **B**: (-) Current electrode  location\n",
    "- **M**: (+) Potential electrode  location\n",
    "- **N**: (-) Potential electrode  location\n",
    "- **$\\rho_1$**: Resistivity of the top layer\n",
    "- **$\\rho_2$**: Resistivity of the bottom layer\n",
    "- **h**: thickness of the first layer\n",
    "- **Plot**: Field to visualize\n",
    "- **Type**: which part of the field"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d1384e6efce843ca873e1a2f8571cc74",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "MyApp(children=(FloatSlider(value=-30.0, continuous_update=False, description='A', max=40.0, min=-40.0, step=1…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_layer_potentials_app()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Building Pseudosections \n",
    "\n",
    "2D profiles are often plotted as pseudo-sections by extending $45^{\\circ}$ lines downwards from the A-B and M-N midpoints and plotting the corresponding $\\Delta V_{MN}$, $\\rho_a$, or misfit value at the intersection of these lines as shown below. For pole-dipole or dipole-pole surveys the $45^{\\circ}$ line is simply extended from the location of the pole. By using this method of plotting, the long offset electrodes plot deeper than those with short offsets. This provides a rough idea of the region sampled by each data point, but the vertical axis of a pseudo-section is not a true depth.\n",
    "\n",
    "In the widget below the red dot marks the midpoint of the current dipole or the location of the A electrode location in a pole-dipole array while the green dots mark the midpoints of the potential dipoles or M electrode locations in a dipole-pole array. The blue dots then mark the location in the pseudo-section where the lines from Tx and Rx midpoints intersect and the data is plotted. By stepping through the Tx (current electrode pairs) using the slider you can see how the pseudo section is built up.\n",
    "\n",
    "The figures shown below show how the points in a pseudo-section are plotted for pole-dipole, dipole-pole, and dipole-dipole arrays. The color coding of the dots match those shown in the widget.\n",
    "<br />\n",
    "<br />\n",
    "<img style=\"float: center; width: 60%; height: 60%\" src=\"https://github.com/geoscixyz/geosci-labs/blob/main/images/dc/PoleDipole.png?raw=true\">\n",
    "<center>Basic skematic for a uniformly spaced pole-dipole array.\n",
    "<br />\n",
    "<br />\n",
    "<br />\n",
    "<img style=\"float: center; width: 60%; height: 60%\" src=\"https://github.com/geoscixyz/geosci-labs/blob/main/images/dc/DipolePole.png?raw=true\">\n",
    "<center>Basic skematic for a uniformly spaced dipole-pole array. \n",
    "<br />\n",
    "<br />\n",
    "<br />\n",
    "<img style=\"float: center; width: 60%; height: 60%\" src=\"https://github.com/geoscixyz/geosci-labs/blob/main/images/dc/DipoleDipole.png?raw=true\">\n",
    "<center>Basic skematic for a uniformly spaced dipole-dipole array.\n",
    "<br />\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pseudo-section app"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d645ed44e78b485a93927aa419e818e8",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "MyApp(children=(IntSlider(value=0, description='i', max=17), Output()), layout=Layout(align_items='stretch', d…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "MidpointPseudoSectionWidget()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## DC pseudo-section app"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " - **$\\rho_1$**: Resistivity of the first layer (thickness of the first layer is 5m)\n",
    " - **$\\rho_2$**: Resistivity of the cylinder\n",
    " - resistivity of the second layer is 1000 $\\Omega$m\n",
    " - **xc**: x location of cylinder center\n",
    " - **zc**: z location of cylinder center\n",
    " - **r**: radius of cylinder\n",
    " - **surveyType**: Type of survey"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6a25181275ec473a99eb1e64715b50d2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(BoundedFloatText(value=1000.0, description='$\\\\rho_1$', max=1000.0, min=10.0), BoundedFl…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "DC2DPseudoWidget()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Parametric Inversion\n",
    "\n",
    "In this final widget you are able to forward model the apparent resistivity of a cylinder embedded in a two layered earth. Pseudo-sections of the apparent resistivity can be generated using dipole-dipole, pole-dipole, or dipole-pole arrays to see how survey geometry can distort the size, shape, and location of conductive bodies in a pseudo-section.  Due to distortion and artifacts present in pseudo-sections trying to interpret them directly is typically difficult and dangerous due to the risk of misinterpretation. Inverting the data to find a model which fits the observed data and is geologically reasonable should be standard practice.   \n",
    "\n",
    "By systematically varying the model parameters and comparing the plots of observed vs. predicted apparent resistivity a parametric inversion can be preformed by hand to find the \"best\" fitting model. Normalized data misfits, which provide a numerical measure of the difference between the observed and predicted data, are useful for quantifying how well and inversion model fits the observed data. The manual inversion process can be difficult and time consuming even with small examples sure as the one presented here. Therefore, numerical optimization algorithms are typically utilized to minimized the data misfit and a model objective function, which provides information about the model structure and complexity, in order to find an optimal solution.\n",
    "\n",
    "## Parametric DC inversion app\n",
    "\n",
    "Definition of variables:\n",
    "- **$\\rho_1$**: Resistivity of the first layer \n",
    "- **$\\rho_2$**: Resistivity of the cylinder\n",
    "- **xc**: x location of cylinder center\n",
    "- **zc**: z location of cylinder center\n",
    "- **r**: radius of cylinder\n",
    "- **predmis**: toggle which allows you to switch the bottom pannel from predicted apparent resistivity to normalized data misfit\n",
    "- **suveyType**: toggle which allows you to switch between survey types.\n",
    "\n",
    "Knonw information\n",
    "- resistivity of the second layer is 1000 $\\Omega$m\n",
    "- thickness of the first layer is known: 5m\n",
    "\n",
    "Unknowns are: $\\rho_1$, $\\rho_2$, xc, zc, and r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "79e1c6477a08484db56eb9cd47a1916a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "MyApp(children=(BoundedFloatText(value=1000.0, description='$\\\\rho_1$', max=1000.0, min=10.0), BoundedFloatTex…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "DC2DfwdWidget()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.13"
  },
  "widgets": {
   "state": {
    "0202a4dd21c24302a603e710b4d9ff2c": {
     "views": [
      {
       "cell_index": 4
      }
     ]
    },
    "1363fe755ac4491f9397e3641ef9f882": {
     "views": [
      {
       "cell_index": 7
      }
     ]
    }
   },
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}