avo_class_4_lightning_talk.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f4c5cf76",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# An investigation on AVO class 4\n",
    "\n",
    "A little investigation on the the values given for Vp, Vs, density in \"Principles of AVO crossplotting\" (Castagna et al., 1997, The Leading Edge) and \"Framework for AVO Gradient and Intercept Interpretation\" (Castagna et al., 1998, Geophysics).\n",
    "\n",
    "**Notes**\n",
    "\n",
    "* This is a version of the notebook made for the [Transform 2021](https://softwareunderground.org/events/transform-2021) Lightning Talks. It is intendended to be run in presentation mode using [Rise](https://rise.readthedocs.io/en/stable/). You should download the notebook, intall Rise (`conda install -c conda-forge rise`), open it with jupyter-notebook and launch the presentation with `Alt-R`.\n",
    "* Since this was supposed to be a presentation it has minimal comments. For additional background and information look at the companion notebook: [Avo Class 4 Investigation](https://nbviewer.jupyter.org/github/aadm/geophysical_notes/blob/master/avo_class_4_investigation.ipynb). There may be overlaps and some differences though.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbdf9e11",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "![anticline avo class 3](https://github.com/aadm/geophysical_notes/blob/master/avo_class_3_anticline.png?raw=true)\n",
    "\n",
    "![anticline avo class 4](https://github.com/aadm/geophysical_notes/blob/bb6e028f15ea411d7a4e1ee91b3e9c447c3d8872/avo_class_4_anticline.png?raw=true)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8ba65a5",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "All the AVO classes in one table:\n",
    "\n",
    "facies        | AVO |  Vp   | Vs    | density | porosity\n",
    "--------------|-----|-------|-------|---------|----------\n",
    "shale         | 1   | 3094  | 1515  | 2.40    |\n",
    "gas sand      | 1   | 4050  | 2526  | 2.21    | 0.20\n",
    "brine sand    | 1   | 4115  | 2543  | 2.32    | 0.20\n",
    "shale         | 2   | 2643  | 1167  | 2.29    |\n",
    "gas sand      | 2   | 2781  | 1665  | 2.08    | 0.25\n",
    "brine sand    | 2   | 3048  | 1595  | 2.23    | 0.25\n",
    "shale         | 3   | 2192  | 818   | 2.16    |\n",
    "gas sand      | 3   | 1543  | 901   | 1.88    | 0.33\n",
    "brine sand    | 3   | 2134  | 860   | 2.11    | 0.33\n",
    "shale         | 4   | 3240  | 1620  | 2.34    |\n",
    "gas sand      | 4   | 1650  | 1090  | 2.07    | ?\n",
    "brine sand    | 4   | 2590  | 1060  | 2.21    | ?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ee9ae42",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Taken from:\n",
    "\n",
    "* Seismic Amplitude Interpretation, Hilterman, 2001, Distinguished Instructor Short Course, SEG\n",
    "* Principles of AVO crossplotting, Castagna et al., 1997, The Leading Edge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2c542e54",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from avo_class_4_investigation_funcs import *\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "mpl.rcParams['figure.dpi'] = 100\n",
    "\n",
    "import ipywidgets as ipw\n",
    "from IPython.display import display\n",
    "from IPython.display import Image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "7834d3fa",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "tmp_shale = np.array([[3094,1515,2.40,0], [2643,1167,2.29,0], [2192,818,2.16,0], [3240,1620,2.34,0]])\n",
    "tmp_sandg = np.array([[4050,2526,2.21,.2], [2781,1665,2.08,.25], [1542,901,1.88,.33], [1650,1090,2.07,0]])\n",
    "tmp_sandb = np.array([[4115,2453,2.32,.2], [3048,1595,2.23,.25], [2134,860,2.11,.33], [2590,1060,2.21,0]])\n",
    "avocl = ['CLASS1', 'CLASS2', 'CLASS3', 'CLASS4']\n",
    "logs = ['VP', 'VS', 'RHO', 'PHI']\n",
    "shale = pd.DataFrame(tmp_shale, columns=logs, index=avocl)\n",
    "sandg = pd.DataFrame(tmp_sandg, columns=logs, index=avocl)\n",
    "sandb = pd.DataFrame(tmp_sandb, columns=logs, index=avocl)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "912d2f7a",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>VP</th>\n",
       "      <th>VS</th>\n",
       "      <th>RHO</th>\n",
       "      <th>PHI</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>CLASS1</th>\n",
       "      <td>4050.0</td>\n",
       "      <td>2526.0</td>\n",
       "      <td>2.21</td>\n",
       "      <td>0.20</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CLASS2</th>\n",
       "      <td>2781.0</td>\n",
       "      <td>1665.0</td>\n",
       "      <td>2.08</td>\n",
       "      <td>0.25</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CLASS3</th>\n",
       "      <td>1542.0</td>\n",
       "      <td>901.0</td>\n",
       "      <td>1.88</td>\n",
       "      <td>0.33</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CLASS4</th>\n",
       "      <td>1650.0</td>\n",
       "      <td>1090.0</td>\n",
       "      <td>2.07</td>\n",
       "      <td>0.00</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            VP      VS   RHO   PHI\n",
       "CLASS1  4050.0  2526.0  2.21  0.20\n",
       "CLASS2  2781.0  1665.0  2.08  0.25\n",
       "CLASS3  1542.0   901.0  1.88  0.33\n",
       "CLASS4  1650.0  1090.0  2.07  0.00"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sandg"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25bb9a5d",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Porosity from densities\n",
    "\n",
    "$$\n",
    "\\phi = \\frac{\\rho - \\rho_m} {\\rho_f- \\rho_m}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "e33a0398",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "STARTING POINT: BRINE SAND DENSITIES\n",
      "CLASS1: calculated=0.20, given=0.20\n",
      "CLASS2: calculated=0.26, given=0.25\n",
      "CLASS3: calculated=0.33, given=0.33\n",
      "CLASS4: calculated=0.27, given=0.00\n",
      "\n",
      "STARTING POINT: GAS SAND DENSITIES\n",
      "CLASS1: calculated=0.19, given=0.20\n",
      "CLASS2: calculated=0.24, given=0.25\n",
      "CLASS3: calculated=0.33, given=0.33\n",
      "CLASS4: calculated=0.25, given=0.00\n"
     ]
    }
   ],
   "source": [
    "rho_m = 2.65\n",
    "\n",
    "print('STARTING POINT: BRINE SAND DENSITIES')\n",
    "for i, ss in sandb.iterrows():\n",
    "    rho_fluid = 1.0135\n",
    "    phi = (ss['RHO'] - rho_m)/(rho_fluid - rho_m)\n",
    "    itms = (ss.name, phi, ss.PHI)\n",
    "    print('{:s}: calculated={:.2f}, given={:.2f}'.format(*itms))\n",
    "\n",
    "print('\\nSTARTING POINT: GAS SAND DENSITIES')\n",
    "for i, ss in sandg.iterrows():\n",
    "    rho_fluid = 0.32\n",
    "    phi = (ss['RHO'] - rho_m)/(rho_fluid - rho_m)\n",
    "    itms = (ss.name, phi, ss.PHI)\n",
    "    print('{:s}: calculated={:.2f}, given={:.2f}'.format(*itms))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee9de3b1",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Porosity from densities\n",
    "\n",
    "$$\n",
    "\\phi = \\frac{\\rho - \\rho_m} {\\rho_f- \\rho_m}\n",
    "$$\n",
    "\n",
    "...for AVO Class 4, set $\\phi=0.26$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "39c6e100",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "sandb.loc['CLASS4', 'PHI'] = 0.26\n",
    "sandg.loc['CLASS4', 'PHI'] = 0.26"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "34c1a78f",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>VP</th>\n",
       "      <th>VS</th>\n",
       "      <th>RHO</th>\n",
       "      <th>PHI</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>CLASS1</th>\n",
       "      <td>4050.0</td>\n",
       "      <td>2526.0</td>\n",
       "      <td>2.21</td>\n",
       "      <td>0.20</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CLASS2</th>\n",
       "      <td>2781.0</td>\n",
       "      <td>1665.0</td>\n",
       "      <td>2.08</td>\n",
       "      <td>0.25</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CLASS3</th>\n",
       "      <td>1542.0</td>\n",
       "      <td>901.0</td>\n",
       "      <td>1.88</td>\n",
       "      <td>0.33</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CLASS4</th>\n",
       "      <td>1650.0</td>\n",
       "      <td>1090.0</td>\n",
       "      <td>2.07</td>\n",
       "      <td>0.26</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            VP      VS   RHO   PHI\n",
       "CLASS1  4050.0  2526.0  2.21  0.20\n",
       "CLASS2  2781.0  1665.0  2.08  0.25\n",
       "CLASS3  1542.0   901.0  1.88  0.33\n",
       "CLASS4  1650.0  1090.0  2.07  0.26"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sandg"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b95f960b",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## porosity from Gassmann's equation\n",
    "\n",
    "$$ \\phi = \\frac{(K_m - K_{satB}) \\cdot (K_m - K_{satG}) \\cdot (K_B - K_G)}{(K_m - K_B) \\cdot (K_m - K_G) \\cdot (K_{satB} - K_{satG})}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "6db5f83a",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "k_m = 37\n",
    "\n",
    "k_w = 2.2\n",
    "rho_w = 1.0\n",
    "\n",
    "k_g = 0.02\n",
    "rho_g = 0.32"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "436ca078",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Porosity: calculated=0.3388, given=0.3300\n"
     ]
    }
   ],
   "source": [
    "avoclass = 'CLASS3'\n",
    "\n",
    "vpG, vsG, rhoG, phi_given = sandg.loc[avoclass]\n",
    "ksat_g = bulk(vpG, vsG, rhoG)\n",
    "\n",
    "vpB, vsB, rhoB, _ = sandb.loc[avoclass]\n",
    "ksat_b = bulk(vpB, vsB, rhoB)\n",
    "\n",
    "phi = gassmann_phi(ksat_b, ksat_g, k_w, k_g, k_m)\n",
    "\n",
    "print('Porosity: calculated={:.4f}, given={:.4f}'.format(phi, phi_given))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "746e6697",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "def phi_from_gassmann(avoclass):\n",
    "    vpG, vsG, rhoG, phi_given = sandg.loc[avoclass]\n",
    "    ksat_g = bulk(vpG, vsG, rhoG)\n",
    "\n",
    "    vpB, vsB, rhoB, _ = sandb.loc[avoclass]\n",
    "    ksat_b = bulk(vpB, vsB, rhoB)\n",
    "\n",
    "    phi = gassmann_phi(ksat_b, ksat_g, k_w, k_g, k_m)\n",
    "    \n",
    "    vp2, vs2, rho2, *_ = gassmann(vpB, vsB, rhoB, rho_w, k_w, rho_g, k_g, k_m, phi)\n",
    "    \n",
    "    print('Input (brine):     Vp={:.0f}, Vs={:.0f}, rho={:.2f}'.format(vpB, vsB, rhoB))\n",
    "    print('Reference (gas):   Vp={:.0f}, Vs={:.0f}, rho={:.2f}'.format(vpG, vsG, rhoG))\n",
    "    print('Output FRM gas:    Vp={:.0f}, Vs={:.0f}, rho={:.2f}'.format(vp2, vs2, rho2))   \n",
    "    itms = (percdiff(vpG, vp2), percdiff(vsG, vs2), percdiff(rhoG, rho2))\n",
    "    print('Difference:        Vp={:+4.1%}, Vs={:+4.1%}, rho={:+4.1%}'.format(*itms))\n",
    "    print('\\nCalculated porosity:  {:.2f}'.format(phi))\n",
    "    print('Given porosity:       {:.2f}'.format(phi_given))\n",
    "    \n",
    "    return phi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "a33a3964",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Input (brine):     Vp=2134, Vs=860, rho=2.11\n",
      "Reference (gas):   Vp=1542, Vs=901, rho=1.88\n",
      "Output FRM gas:    Vp=1550, Vs=911, rho=1.88\n",
      "Difference:        Vp=+0.5%, Vs=+1.1%, rho=-0.0%\n",
      "\n",
      "Calculated porosity:  0.34\n",
      "Given porosity:       0.33\n"
     ]
    }
   ],
   "source": [
    "phi_from_gassmann('CLASS3');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "cc7eb7dc",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Input (brine):     Vp=2590, Vs=1060, rho=2.21\n",
      "Reference (gas):   Vp=1650, Vs=1090, rho=2.07\n",
      "Output FRM gas:    Vp=1643, Vs=1088, rho=2.10\n",
      "Difference:        Vp=-0.4%, Vs=-0.2%, rho=+1.4%\n",
      "\n",
      "Calculated porosity:  0.16\n",
      "Given porosity:       0.26\n"
     ]
    }
   ],
   "source": [
    "phi_from_gassmann('CLASS4');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "c332fbe4",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "def phi_finder(avoclass, rho_g, k_g, rho_w, k_w, k_m):\n",
    "    phi_array = np.linspace(.1, .4)\n",
    "   \n",
    "    vpG, vsG, rhoG, phi_given = sandg.loc[avoclass]\n",
    "    vpB, vsB, rhoB, _ = sandb.loc[avoclass]\n",
    "    \n",
    "    vp_0, vs_0, rho_0, ksat_0, kdry_0 = gassmann(vpG, vsG, rhoG, rho_g, k_g, rho_w, k_w, k_m, phi_array)\n",
    "    vp_1, vs_1, rho_1, ksat_1, kdry_1 = gassmann(vpB, vsB, rhoB, rho_w, k_w, rho_g, k_g, k_m, phi_array)\n",
    "\n",
    "    # find intersection and return porosity\n",
    "    idx, val = find_nearest(kdry_0,kdry_1)\n",
    "    phi_calc = phi_array[idx]\n",
    "    \n",
    "    f, ax = plt.subplots(figsize=(7, 5), constrained_layout=True)\n",
    "    ax.plot(phi_array, kdry_0, '-r', lw=4, alpha=0.5, label='K_dry from gas')\n",
    "    ax.plot(phi_array, kdry_1, '-b', lw=4, alpha=0.5, label='K_dry from brine')\n",
    "    ax.axvline(phi_given, ls=':', color='k', label='Given phi={:.3f}'.format(phi_given))\n",
    "    ax.plot(phi_calc, val, marker='o', mec='k', mfc='none', ms=15, ls='none', label='Calculated phi={:.3f}'.format(phi_calc))\n",
    "    ax.set_xlabel('Porosity')\n",
    "    ax.set_ylabel(r'$K_{dry}$ (GPa)')\n",
    "    ax.legend(bbox_to_anchor=(1.01, 1))\n",
    "    ax.set_ylim(-5, 25)\n",
    "    ax.set_title(avoclass)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "541345ce",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "def reset_values(b):\n",
    "    w_avoc.value = 'CLASS3'\n",
    "    w_km.value = 37\n",
    "    w_rhog.value = 0.3\n",
    "    w_kg.value = 0.02\n",
    "    w_rhow.value = 1.0\n",
    "    w_kw.value = 2.2\n",
    "\n",
    "sty = {'description_width': 'initial'}\n",
    "w_avoc = ipw.ToggleButtons(value='CLASS3', description='AVO Class', options=['CLASS1', 'CLASS2', 'CLASS3', 'CLASS4'])\n",
    "w_km = ipw.FloatSlider(value=37, description='K matrix (GPa)', max=40, min=30, step=1, style=sty)\n",
    "w_rhog = ipw.FloatSlider(value=0.3, description='rho gas (g/cc)', max=0.4, min=0.1, step=0.05, style=sty)\n",
    "w_kg = ipw.FloatSlider(value=0.02, description='K gas (GPa)', max=0.2, min=0.01, step=0.01, style=sty)\n",
    "w_rhow = ipw.FloatSlider(value=1.0, description='rho water (g/cc)', max=1.1, min=0.9, step=0.05, style=sty)\n",
    "w_kw = ipw.FloatSlider(value=2.2, description='K water (GPa)', max=3.5, min=2, step=0.05, style=sty)\n",
    "reset_button = ipw.Button(description = \"Reset\")\n",
    "reset_button.on_click(reset_values)\n",
    "\n",
    "lbox = ipw.VBox([w_kg, w_rhog])\n",
    "rbox = ipw.VBox([w_kw, w_rhow])\n",
    "lrbox = ipw.HBox([lbox, rbox])\n",
    "ui0 = ipw.VBox([w_avoc, w_km, lrbox, reset_button])\n",
    "\n",
    "tool0 = ipw.interactive_output(phi_finder, {\n",
    "    'avoclass': w_avoc,\n",
    "    'rho_g': w_rhog,\n",
    "    'k_g': w_kg,\n",
    "    'rho_w': w_rhow,\n",
    "    'k_w': w_kw,\n",
    "    'k_m': w_km})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83d8eae1",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "display(ui0, tool0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a84adb3b",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Gassmann does not work\n",
    "\n",
    "* porosity from density equation: $0.26$\n",
    "* porosity from Gassmann: $0.16$\n",
    "* is it realistic for a $\\phi=0.16$ sand to increase its Vp of 940 m/s (+57%) going from gas to brine?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2620a549",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "check_vpb = sandb.loc['CLASS4', 'VP']\n",
    "check_vpg = sandg.loc['CLASS4', 'VP']\n",
    "print(check_vpb-check_vpg)\n",
    "print(percdiff(check_vpg, check_vpb))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b24c48e",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "check_vpb = sandb.loc['CLASS3', 'VP']\n",
    "check_vpg = sandg.loc['CLASS3', 'VP']\n",
    "print(check_vpb-check_vpg)\n",
    "print(percdiff(check_vpg, check_vpb))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1aeb9c4f",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* Class 3 sand with $\\phi=0.33$ shows a +38% increase!\n",
    "* ...if Gassmann does not work then one of the inputs must be wrong!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95c96a46",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## find a justification for Class 4 numbers with Rock Physics Modeling\n",
    "\n",
    "* Make Vp, Vs, rho using Soft Sand model on a range of porosities\n",
    "* see if I match the input Vp, Vs, rho for AVO class 4 sand\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a0b81b3",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "def match_class4_rpm(avoclass, sigma, phic):\n",
    "    # assign elastic properties for fluid (gas) and rock matrix\n",
    "    Df, Kf = 0.3,  0.02\n",
    "    Dqz, Kqz, Gqz = 2.65, 37, 44\n",
    "    f=1\n",
    "\n",
    "    # relation critical porosity-coordination number by Murphy (1982)\n",
    "    # Mavko et al. (2009), The Rock Physics Handbook, Cambridge UP (p.246)\n",
    "    Cn = 20-34*phic+14*phic**2\n",
    "\n",
    "    # make vp, vs, rho for sand\n",
    "    phi_dist = np.linspace(0.1, 0.4)\n",
    "    ns = phi_dist.size\n",
    "    Kdry, Gdry = softsand(Kqz, Gqz, phi_dist, sigma, phic, Cn, f)\n",
    "    vp, vs, rho, _ = vels(Kdry, Gdry, Kqz, Dqz, Kf, Df, phi_dist)\n",
    "\n",
    "    # get vp, vs, rho for class 4\n",
    "    sand_vp, sand_vs, sand_rho, _ = sandg.loc[avoclass]\n",
    "   \n",
    "    # do the plot\n",
    "    f, ax = plt.subplots(ncols=3, sharex=True, figsize=(8,4), constrained_layout=True)\n",
    "    ax[0].plot(phi_dist, vp, '-k')\n",
    "    ax[1].plot(phi_dist, vs, '-k')\n",
    "    ax[2].plot(phi_dist, rho, '-k')\n",
    "    ax[0].plot(0.26, sand_vp, 'or', alpha=0.5)\n",
    "    ax[1].plot(0.26, sand_vs, 'or', alpha=0.5)\n",
    "    ax[2].plot(0.26, sand_rho, 'or', alpha=0.5)\n",
    "    ax[0].set_xlabel('Porosity')\n",
    "    ax[0].set_ylabel('Vp (m/s)')\n",
    "    ax[1].set_ylabel('Vs (m/s)')\n",
    "    ax[2].set_ylabel('Density (g/cc)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d91ffed",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "w_sigma = ipw.FloatSlider(value=20, description='Effective pressure (MPa)', max=40, min=5, step=5, style=sty)\n",
    "w_phic = ipw.FloatSlider(value=0.3, description='Critical porosity', max=0.6, min=0.2, step=0.05, style=sty)\n",
    "w_fudg = ipw.FloatSlider(value=1.0, description='Fudge factor', max=1.0, min=0., step=0.1, style=sty)\n",
    "\n",
    "box = ipw.HBox([w_sigma, w_phic])\n",
    "ui1 = ipw.VBox([w_avoc, box])\n",
    "tool1 = ipw.interactive_output(match_class4_rpm, {\n",
    "    'avoclass': w_avoc,\n",
    "    'sigma': w_sigma,\n",
    "    'phic': w_phic})\n",
    "#    'f': w_fudg})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e1f9615",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "display(ui1, tool1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0df777e4",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "caprock_vp, caprock_vs, caprock_rho, _ = shale.loc['CLASS4']\n",
    "sand_vp, sand_vs, sand_rho, _ = sandg.loc['CLASS3']\n",
    "_, i_newcl4, g_newcl4 = shuey(caprock_vp, caprock_vs, caprock_rho, sand_vp, sand_vs, sand_rho, np.arange(30), terms=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7748d72",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# in short\n",
    "\n",
    "* somebody would please ask John Castagna where those numbers come from?\n",
    "* class 4 AVO is only peculiar because of what's above so:\n",
    "    * take a class 3 sand\n",
    "    * put class 4 shale on top\n",
    "    * you still get something that plots in the right quadrant"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3277aae3",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "ax = classref(plot_brine=False)\n",
    "ax.plot(i_newcl4, g_newcl4, '^', mec='r', mfc='k', ms=14)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "769dd44b",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# anything else?\n",
    "\n",
    "* some more RPM to validate those Vp, Vs, rho for Class 4 sand\n",
    "* ...perhaps using `scipy.optimize.least_squares` and see if the parameters are ok"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb6ca302",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**notebook available here:** https://github.com/aadm/geophysical_notes\n",
    "\n",
    "`avo_class_4_lightning_talk.ipynb`"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.10.6"
  },
  "rise": {
   "controls": false,
   "enable_chalkboard": true,
   "footer": "<h3>AVO Class 4 Investigation by AAdM, T2021 Lightning Talks</h3>",
   "show_buttons_on_startup": false,
   "theme": "simple",
   "transition": "slide"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}