{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib notebook" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\nDC apparent resistivity\n=======================\n\nDC apparent resistivity, dipole-dipole configuration.\n\nThere are various DC sounding layouts, the most common ones being Schlumberger,\nWenner, pole-pole, pole-dipole, and **dipole-dipole**, at which we have a look\nhere.\n\nDipole-dipole layout as shown in figure 8.32 in Kearey et al. (2002):\n\n![](../../_static/figures/Fig_from_8-32.jpg)\n\n\nThe apparent resistivity $\\rho_a$ of the *plotting point* is then\ncomputed with\n\n\\begin{align}\\rho_a = \\frac{V}{I} \\pi na(n+1)(n+2)\\ ,\\end{align}\n\nwhere $V$ is measured Voltage, $I$ is source strength, $a$ is\ndipole length, and $n$ is the factor of source-receiver separation.\n\n**References**\n\n**Kearey, P., M. Brooks, and I. Hill, 2002**, An introduction to geophysical\nexploration, 3rd ed.: Blackwell Scientific Publications, ISBN: 0 632 04929 4.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import empymod\nimport numpy as np\nimport matplotlib.pyplot as plt\nplt.style.use('ggplot')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compute $\\boldsymbol{\\rho_a}$\n-----------------------------------\n\nFirst we define a function to compute apparent resistivity for a given model\nand given source and receivers.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def comp_appres(depth, res, a, n, srcpts=1, recpts=1, verb=1):\n \"\"\"Return apparent resistivity for dipole-dipole DC measurement\n\n rho_a = V/I pi a n (n+1) (n+2).\n\n Returns die apparent resistivity due to:\n - Electric source, inline (y = 0 m).\n - Source of 1 A strength.\n - Source and receiver are located at the air-interface.\n - Source is centered at x = 0 m.\n\n Note: DC response can be obtained by either t->infinity s or f->0 Hz. f = 0\n Hz is much faster, as there is no Fourier transform involved and only\n a single frequency has to be computed. By default, the minimum\n frequency in empymod is 1e-20 Hz. The difference between the signals\n for 1e-20 Hz and 0 Hz is very small.\n\n For more explanation regarding input parameters see `empymod.model`.\n\n Parameters\n ----------\n depth : Absolute depths of layer interfaces, without air-interface.\n res : Resistivities of the layers, one more than depths (lower HS).\n a : Dipole length.\n n : Separation factors.\n srcpts, recpts : If < 3, bipoles are approximated as dipoles.\n verb : Verbosity.\n\n Returns\n -------\n rho_a : Apparent resistivity.\n AB2 : Src/rec-midpoints\n\n \"\"\"\n\n # Get offsets between src-midpoint and rec-midpoint, AB\n AB = (n+1)*a\n\n # Collect model, putting source and receiver slightly (1e-3 m) into the\n # ground.\n model = {\n 'src': [-a/2, a/2, 0, 0, 1e-3, 1e-3],\n 'rec': [AB-a/2, AB+a/2, AB*0, AB*0, 1e-3, 1e-3],\n 'depth': np.r_[0, np.array(depth, ndmin=1)],\n 'freqtime': 1e-20, # Smaller f would be set to 1e-20 be empymod.\n 'verb': verb, # Setting it to 1e-20 avoids warning-message.\n 'res': np.r_[2e14, np.array(res, ndmin=1)],\n 'strength': 1, # So it is NOT normalized to 1 m src/rec.\n 'htarg': {'pts_per_dec': -1},\n }\n\n return np.real(empymod.bipole(**model))*np.pi*a*n*(n+1)*(n+2), AB/2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot-function\n-------------\n\nSecond we create a plot-function, which includes the call to `comp_appres`,\nto use for a couple of different models.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def plotit(depth, a, n, res1, res2, res3, title):\n \"\"\"Call `comp_appres` and plot result.\"\"\"\n\n # Compute the three different models\n rho1, AB2 = comp_appres(depth, res1, a, n)\n rho2, _ = comp_appres(depth, res2, a, n)\n rho3, _ = comp_appres(depth, res3, a, n)\n\n # Create figure\n plt.figure()\n\n # Plot curves\n plt.loglog(AB2, rho1, label='Case 1')\n plt.plot(AB2, rho2, label='Case 2')\n plt.plot(AB2, rho3, label='Case 3')\n\n # Legend, labels\n plt.legend(loc='best')\n plt.title(title)\n plt.xlabel('AB/2 (m)')\n plt.ylabel(r'Apparent resistivity $\\rho_a (\\Omega\\,$m)')\n\n plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Model 1: 2 layers\n~~~~~~~~~~~~~~~~~\n\n+--------+---------------------+---------------------+\n|layer | depth (m) | resistivity (Ohm m) |\n+========+=====================+=====================+\n|air | $-\\infty$ - 0 | 2e14 |\n+--------+---------------------+---------------------+\n|layer 1 | 0 - 50 | 10 |\n+--------+---------------------+---------------------+\n|layer 2 | 50 - $\\infty$ | 100 / 10 / 1 |\n+--------+---------------------+---------------------+\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plotit(\n 50, # Depth\n 20, # a (src- and rec-lengths)\n np.arange(3, 500), # n\n [10, 100], # Case 1\n [10, 10], # Case 2\n [10, 1], # Case 3\n 'Model 1: 2 layers')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Model 2: layer embedded in background\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n\n+--------+----------------------+---------------------+\n|layer | depth (m) | resistivity (Ohm m) |\n+========+======================+=====================+\n|air | $-\\infty$ - 0 | 2e14 |\n+--------+----------------------+---------------------+\n|layer 1 | 0 - 50 | 10 |\n+--------+----------------------+---------------------+\n|layer 2 | 50 - 500 | 100 / 10 / 1 |\n+--------+----------------------+---------------------+\n|layer 3 | 500 - $\\infty$ | 10 |\n+--------+----------------------+---------------------+\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plotit(\n [50, 500], # Depth\n 20, # a (src- and rec-lengths)\n np.arange(3, 500), # n\n [10, 100, 10], # Case 1\n [10, 10, 10], # Case 2\n [10, 1, 10], # Case 3\n 'Model 2: layer embedded in background')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Model 3: 3 layers\n~~~~~~~~~~~~~~~~~\n\n+--------+----------------------+----------------------+\n|layer | depth (m) | resistivity (Ohm m) |\n+========+======================+======================+\n|air | $-\\infty$ - 0 | 2e14 |\n+--------+----------------------+----------------------+\n|layer 1 | 0 - 50 | 10 |\n+--------+----------------------+----------------------+\n|layer 2 | 50 - 500 | 100 / 10 / 1 |\n+--------+----------------------+----------------------+\n|layer 3 | 500 - $\\infty$ | 1000 / 10 / 0.1 |\n+--------+---------------------+-----------------------+\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plotit(\n [50, 500], # Depth\n 20, # a (src- and rec-lengths)\n np.arange(3, 500), # n\n [10, 100, 1000], # Case 1\n [10, 10, 10], # Case 2\n [10, 1, 0.1], # Case 3\n 'Model 3: 3 layers')" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "empymod.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.7.3" } }, "nbformat": 4, "nbformat_minor": 0 }