avo_class_4_investigation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/aadm/geophysical_notes/blob/master/avo_class_4_investigation.ipynb)\n",
    "\n",
    "# Reference elastic values for AVO class 4\n",
    "\n",
    "This is a little investigation on the reliability of the values given in two papers written by John Castagna et al. (1997, 1998 -- see references below).\n",
    "\n",
    "In these papers I have found Vp, Vs and density for an AVO class 4 interface; I have always put them together with those found for the other 3 AVO classes in Hilterman, 2001 (see below for references) but then I realized there was no porosity assigned to the class 4 sand.\n",
    "\n",
    "I started to look into it and look for the porosity that could be assigned to this sand. But then I realized that depending on the method there is no unique solutions, and that suggests that one of the inputs (either velocities or density) are wrong.\n",
    "\n",
    "My suggestion it this: use the class 3 sand velocities and density in place of class 4 sand. With the same class 4 overburden shale we still get a class 4 response, and none of the absurdities that derive from trying to make sense of the parameters used by Castagna.\n",
    "\n",
    "**REFERENCES**\n",
    "\n",
    "* [Seismic Amplitude Interpretation]((http://library.seg.org/doi/book/10.1190/1.9781560801993)\" (Hilterman, 2001, Distinguished Instructor Short Course, SEG, section 4.13)\n",
    "* \"Principles of AVO crossplotting\" (Castagna et al., 1997, The Leading Edge)\n",
    "* \"Framework for AVO Gradient and Intercept Interpretation\" (Castagna et al., 1998, Geophysics)\n",
    "\n",
    "**THINGS THAT CAN BE DONE**\n",
    "\n",
    "- Use Batzle-Wang to validate fluid moduli considering overburden depths assigned to the various classes\n",
    "- use rock physics models to validate all AVO classes sand properties"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### coding details\n",
    "\n",
    "The next cell block does the standard imports and defines a few other functions. This makes it a rather lenghty block which should be kept hidden.\n",
    "\n",
    "The additional functions are:\n",
    "\n",
    "* `gassmann` to implement Gassmann's equation\n",
    "* `gassmann_phi` to get porosity from elastic parameters using Gassmann's equation\n",
    "* `bulk` to get bulk modulus from Vp, Vs and density\n",
    "* `shuey` to calculate the reflectivity variation with angle\n",
    "* `vels` to calculate velocities from density, bulk and shear moduli\n",
    "* `classref` to plot AVO classes\n",
    "* `hertzmindlin` and `softsand` to implement the Soft Sand rock physics model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "mpl.rcParams['figure.dpi'] = 100\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.patches as patches\n",
    "import ipywidgets as ipw\n",
    "from IPython.display import display\n",
    "\n",
    "def percdiff(start, end):\n",
    "    return (end-start)/start\n",
    "\n",
    "\n",
    "def find_nearest(a, a0):\n",
    "    idx = np.abs(a - a0).argmin()\n",
    "    return idx, a[idx]\n",
    "\n",
    "\n",
    "def bulk(vp, vs, rho):\n",
    "    # converts density to SI (kg/m3)\n",
    "    D = rho*1e3\n",
    "    K = D*vp**2 - 4/3*D*vs**2\n",
    "    return K/1e9\n",
    "\n",
    "\n",
    "def gassmann(vp1, vs1, rho1, rhof1, kfluid1, rhof2, kfluid2, kmin, phi):\n",
    "    # convert density to kg/m3 and elastic moduli to Pa\n",
    "    d1 = rho1*1e3\n",
    "    df1 = rhof1*1e3\n",
    "    df2 = rhof2*1e3\n",
    "    k0 = kmin*1e9\n",
    "    kf1 = kfluid1*1e9\n",
    "    kf2 = kfluid2*1e9\n",
    "    d2 = d1 - phi * df1 + phi * df2\n",
    "    mu1 = d1 * vs1**2\n",
    "    k1 = d1 * vp1**2 - (4/3) * mu1\n",
    "    kd = (k1 * (phi*k0 / kf1 + 1-phi) - k0) / (phi*k0 / kf1 + k1 / k0 - 1-phi)\n",
    "    mu2 = mu1\n",
    "    with np.errstate(divide='ignore', invalid='ignore'):\n",
    "        k2 = kd + (1 - kd/k0 )**2 / (phi/kf2 + (1-phi)/k0 - kd/k0**2)\n",
    "        vp2 = np.sqrt((k2 + 4/3*mu2) / d2)\n",
    "        vs2 = np.sqrt(mu2 / d2)\n",
    "        return vp2, vs2, d2/1e3, k2/1e9, kd/1e9\n",
    "\n",
    "\n",
    "def gassmann_phi(ksat1, ksat2, kf1, kf2, kmin):\n",
    "    a = (kmin-ksat1)*(kmin-ksat2)*(kf1-kf2) \n",
    "    b = (kmin-kf1)*(kmin-kf2)*(ksat1-ksat2)\n",
    "    return a / b\n",
    "\n",
    "\n",
    "def hertzmindlin(K0, G0, sigma, phi_c=0.4, Cn=8.6, f=1):\n",
    "    sigma0 = sigma / 1e3  # converts pressure in same units as solid moduli (GPa)\n",
    "    pr0 = (3*K0-2*G0) / (6*K0+2*G0)  # poisson's ratio of mineral mixture\n",
    "    Khm = (sigma0*(Cn**2*(1 - phi_c)**2*G0**2) / (\n",
    "           18*np.pi**2 * (1 - pr0)**2))**(1/3)\n",
    "    Ghm = ((2+3*f-pr0*(1+3*f)) / (5*(2-pr0))) * (\n",
    "          (sigma0 * (3 * Cn**2 * (1 - phi_c)**2 * G0**2) / (\n",
    "           2 * np.pi**2 * (1 - pr0)**2)))**(1/3)\n",
    "    return Khm, Ghm\n",
    "\n",
    "\n",
    "def softsand(K0, G0, phi, sigma, phi_c=0.4, Cn=8.6, f=1):\n",
    "    Khm, Ghm = hertzmindlin(K0, G0, sigma, phi_c, Cn, f)\n",
    "    Kdry = -4/3 * Ghm + (((phi / phi_c) / (Khm + 4/3 * Ghm)) + (\n",
    "           (1 - phi / phi_c) / (K0 + 4/3 * Ghm)))**-1\n",
    "    gxx = Ghm / 6 * ((9 * Khm + 8 * Ghm) / (Khm + 2 * Ghm))\n",
    "    Gdry = -gxx + ((phi / phi_c) / (Ghm + gxx) + (\n",
    "           (1 - phi / phi_c) / (G0 + gxx)))**-1\n",
    "    return Kdry, Gdry\n",
    "\n",
    "\n",
    "def vels(kdry, gdry, kmin, rho0, kfluid, rhof, phi):\n",
    "    # convert density to kg/m3 and elastic moduli to Pa\n",
    "    d0 = rho0*1e3\n",
    "    df = rhof*1e3\n",
    "    kd = kdry*1e9\n",
    "    gd = gdry*1e9\n",
    "    kf = kfluid*1e9\n",
    "    k0 = kmin*1e9\n",
    "    rho = d0 * (1 - phi) + df * phi\n",
    "    with np.errstate(divide='ignore', invalid='ignore'):\n",
    "        ksat = kd + (1 - kd/k0)**2 / (phi/kf + (1-phi)/k0 - kd/k0**2)\n",
    "        vp = np.sqrt((ksat+4/3*gd)/rho)\n",
    "        vs = np.sqrt(gd/rho)\n",
    "    return vp, vs, rho/1e3, ksat/1e9\n",
    "\n",
    "\n",
    "def shuey(vp1, vs1, rho1, vp2, vs2, rho2, theta, approx=True, terms=False):\n",
    "    a = np.radians(theta)\n",
    "    dvp = vp2-vp1\n",
    "    dvs = vs2-vs1\n",
    "    drho = rho2-rho1\n",
    "    vp  = np.mean([vp1, vp2], axis=0)\n",
    "    vs  = np.mean([vs1, vs2], axis=0)\n",
    "    rho = np.mean([rho1, rho2], axis=0)\n",
    "    R0 = 0.5*(dvp/vp + drho/rho)\n",
    "    G  = 0.5*(dvp/vp) - 2*(vs**2/vp**2)*(drho/rho+2*(dvs/vs))\n",
    "    F =  0.5*(dvp/vp)\n",
    "    # if angles is an array\n",
    "    if a.size>1:\n",
    "        R0 = R0.reshape(-1,1)\n",
    "        G = G.reshape(-1,1)\n",
    "        F = F.reshape(-1,1)\n",
    "    if approx:\n",
    "        R = R0 + G*np.sin(a)**2\n",
    "    else:\n",
    "        R = R0 + G*np.sin(a)**2 + F*(np.tan(a)**2-np.sin(a)**2)\n",
    "    if terms:\n",
    "        return R, R0, G\n",
    "    else:\n",
    "        return R\n",
    "\n",
    "\n",
    "def classref(near=5, far=30, mx=.6, plot_brine=False, plot_colorzones=True):\n",
    "    tmp_shl = np.array([[3094, 1515, 2.40, 0],\n",
    "                        [2643, 1167, 2.29, 0],\n",
    "                        [2192, 818, 2.16, 0],\n",
    "                        [3240, 1620, 2.34, 0]])\n",
    "    tmp_ssg = np.array([[4050, 2526, 2.21, .2],\n",
    "                        [2781, 1665, 2.08, .25],\n",
    "                        [1542, 901, 1.88, .33],\n",
    "                        [1650, 1090, 2.07, .156]])\n",
    "    tmp_ssb = np.array([[4115, 2453, 2.32, .2],\n",
    "                        [3048, 1595, 2.23, .25],\n",
    "                        [2134, 860, 2.11, .33],\n",
    "                        [2590, 1060, 2.21, .156]])\n",
    "    avocl = ['CLASS1', 'CLASS2', 'CLASS3', 'CLASS4']\n",
    "    logs = ['VP', 'VS', 'RHO', 'PHI']\n",
    "    shl = pd.DataFrame(tmp_shl, columns=logs, index=avocl)\n",
    "    ssg = pd.DataFrame(tmp_ssg, columns=logs, index=avocl)\n",
    "    ssb = pd.DataFrame(tmp_ssb, columns=logs, index=avocl)\n",
    "\n",
    "    opttxt = dict(weight='bold', ha='left', va='center')\n",
    "    mrkg = {'ms': 10, 'mew': 2, 'ls': 'none'}\n",
    "    mrkb = {'ms': 10, 'mew': 2, 'ls': 'none', 'mfc':'none'}\n",
    "    mrk_sel = {'marker': '*', 'mec': 'k', 'mfc': 'white', 'ms': 16, 'ls': 'none', 'mew': 2}\n",
    "\n",
    "    angs = np.array([near, far])\n",
    "    tmp = ['C0', 'C1', 'C2', 'C3']\n",
    "    cc = dict(zip(avocl, tmp))\n",
    "    tmp = ['s', 'P', 'v', '^']\n",
    "    mm  = dict(zip(avocl, tmp))\n",
    "\n",
    "    f, ax = plt.subplots(constrained_layout=True)\n",
    "    ax.axhline(0, color='k', lw=3)\n",
    "    ax.axvline(0, color='k', lw=3)\n",
    "    for i, sh in shl.iterrows():\n",
    "        vpsh, vssh, dsh = sh['VP'], sh['VS'], sh['RHO']\n",
    "        vpb, vsb, db = ssb.loc[i, 'VP'], ssb.loc[i, 'VS'], ssb.loc[i, 'RHO']\n",
    "        vpg, vsg, dg = ssg.loc[i, 'VP'], ssg.loc[i, 'VS'], ssg.loc[i, 'RHO']\n",
    "        Ab, Ib, Gb = shuey(vpsh, vssh, dsh, vpb, vsb, db, angs, terms=True)\n",
    "        Ag, Ig, Gg = shuey(vpsh, vssh, dsh, vpg, vsg, dg, angs, terms=True)\n",
    "        ax.plot(Ig, Gg, fillstyle='full', label=sh.name, marker=mm[i], mfc=cc[i], mec=cc[i], **mrkg)\n",
    "        if plot_brine:\n",
    "            ax.plot(Ib, Gb, fillstyle='none', label=None, marker=mm[i], mec=cc[i], **mrkb)\n",
    "    ax.set_xlabel('Intercept')\n",
    "    ax.set_ylabel('Gradient')\n",
    "    ax.legend()\n",
    "    ax.set_xlim(-mx, mx)\n",
    "    ax.set_ylim(-mx, mx)\n",
    "    ax.set_aspect('equal', 'box')\n",
    "    ax.grid()\n",
    "    if plot_colorzones:\n",
    "        opt1 = dict(edgecolor='None', alpha=0.2)\n",
    "        cl1_area = patches.Rectangle((0.02, -1), .98, 1, facecolor=cc['CLASS1'], **opt1)\n",
    "        cl2_area = patches.Rectangle((-0.02, -1), .04, 2, facecolor=cc['CLASS2'], **opt1)\n",
    "        cl3_area = patches.Rectangle((-1, -1), .98, 1, facecolor=cc['CLASS3'], **opt1)\n",
    "        cl4_area = patches.Rectangle((-1, 0), .98, 1, facecolor=cc['CLASS4'], **opt1)\n",
    "        background = patches.Polygon([[-1, 1], [1, -1], [1, 1]], facecolor='w')\n",
    "        ax.add_patch(cl1_area)\n",
    "        ax.add_patch(cl2_area)\n",
    "        ax.add_patch(cl3_area)\n",
    "        ax.add_patch(cl4_area)\n",
    "        ax.add_patch(background)\n",
    "        opt2 = dict(ha='center', va='center', weight='bold', size='large')\n",
    "        ax.text(.15, -.3, 'Class 1', color=cc['CLASS1'], **opt2)\n",
    "        ax.text(0, -.25, 'Class 2/2p', color=cc['CLASS2'], **opt2)\n",
    "        ax.text(-.35, -.3, 'Class 3', color=cc['CLASS3'], **opt2)\n",
    "        ax.text(-.35, .15, 'Class 4', color=cc['CLASS4'], **opt2)\n",
    "    if plot_brine:\n",
    "        ax.set_title('Filled markers: gas, empty=brine')\n",
    "    return ax"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The table below shows the elastic properties for the various AVO classes taken from the references above (Hilterman 2001, Castagna et al., 1997 and 1998):\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    | ?\n",
    "\n",
    "Porosity is unknown for class 4.\n",
    "\n",
    "First we build 3 tables detailing the elastic properties $V_p$, $V_s$ and $\\rho$ plus porosity for shales, gas sands and brine sands for each AVO class:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "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,
   "metadata": {},
   "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>3094.0</td>\n",
       "      <td>1515.0</td>\n",
       "      <td>2.40</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CLASS2</th>\n",
       "      <td>2643.0</td>\n",
       "      <td>1167.0</td>\n",
       "      <td>2.29</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CLASS3</th>\n",
       "      <td>2192.0</td>\n",
       "      <td>818.0</td>\n",
       "      <td>2.16</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CLASS4</th>\n",
       "      <td>3240.0</td>\n",
       "      <td>1620.0</td>\n",
       "      <td>2.34</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            VP      VS   RHO  PHI\n",
       "CLASS1  3094.0  1515.0  2.40  0.0\n",
       "CLASS2  2643.0  1167.0  2.29  0.0\n",
       "CLASS3  2192.0   818.0  2.16  0.0\n",
       "CLASS4  3240.0  1620.0  2.34  0.0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "shale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "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": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sandg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "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>4115.0</td>\n",
       "      <td>2453.0</td>\n",
       "      <td>2.32</td>\n",
       "      <td>0.20</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CLASS2</th>\n",
       "      <td>3048.0</td>\n",
       "      <td>1595.0</td>\n",
       "      <td>2.23</td>\n",
       "      <td>0.25</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CLASS3</th>\n",
       "      <td>2134.0</td>\n",
       "      <td>860.0</td>\n",
       "      <td>2.11</td>\n",
       "      <td>0.33</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CLASS4</th>\n",
       "      <td>2590.0</td>\n",
       "      <td>1060.0</td>\n",
       "      <td>2.21</td>\n",
       "      <td>0.00</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            VP      VS   RHO   PHI\n",
       "CLASS1  4115.0  2453.0  2.32  0.20\n",
       "CLASS2  3048.0  1595.0  2.23  0.25\n",
       "CLASS3  2134.0   860.0  2.11  0.33\n",
       "CLASS4  2590.0  1060.0  2.21  0.00"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sandb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Porosity from densities\n",
    "\n",
    "We can derive porosity from given densities using the simple formula:\n",
    "\n",
    "$$\n",
    "\\phi = \\frac{\\rho - \\rho_m} {\\rho_f- \\rho_m}\n",
    "$$\n",
    "\n",
    "The returned values when compared to the given porosities are a bit off; I need to tweak brine and gas density to get porosities for classes 1-3 \"almost\" correct.\n",
    "\n",
    "For example, starting from brine sands I initially set $\\rho_{B} = 1$ for brine density but I ended up with $1.0135$ to match the given values.\n",
    "\n",
    "Starting from gas sands instead I have to use an unusually high gas density ($\\rho_{G} = 0.35$, usually this is around $0.1$) to get closer to the given values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "STARTING POINT: BRINE SAND DENSITIES\n",
      "CLASS1: calculated=0.2016, given=0.20\n",
      "CLASS2: calculated=0.2566, given=0.25\n",
      "CLASS3: calculated=0.3300, given=0.33\n",
      "CLASS4: calculated=0.2689, given=0.00\n",
      "\n",
      "STARTING POINT: GAS SAND DENSITIES\n",
      "CLASS1: calculated=0.1913, given=0.20\n",
      "CLASS2: calculated=0.2478, given=0.25\n",
      "CLASS3: calculated=0.3348, given=0.33\n",
      "CLASS4: calculated=0.2522, given=0.00\n"
     ]
    }
   ],
   "source": [
    "rho_m = 2.65  # mineral density\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={:.4f}, given={:.2f}'.format(*itms))\n",
    "\n",
    "print('\\nSTARTING POINT: GAS SAND DENSITIES')\n",
    "for i, ss in sandg.iterrows():\n",
    "    rho_fluid = 0.35\n",
    "    phi = (ss['RHO'] - rho_m)/(rho_fluid - rho_m)\n",
    "    itms = (ss.name, phi, ss.PHI)\n",
    "    print('{:s}: calculated={:.4f}, given={:.2f}'.format(*itms))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead of using fixed values for the fluid densities and assuming correct given porosities (at least for classes 1-3) and constant mineral density $\\rho_{m} = 2.65 g/cm^2$, we can scan through realistic values for brine and gas density and select those values that return the given porosity. See next example for class 3:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "i = 'CLASS3'\n",
    "phi = sandb.loc[i, 'PHI']\n",
    "rhoB = sandb.loc[i, 'RHO']\n",
    "rhoG = sandg.loc[i, 'RHO']\n",
    "\n",
    "rho_b = np.linspace(0.9, 1.2)\n",
    "rho_g = np.linspace(0.05, 0.6)\n",
    "\n",
    "phi_from_bri = (rhoB - rho_m)/(rho_b - rho_m)\n",
    "phi_from_gas = (rhoG - rho_m)/(rho_g - rho_m)\n",
    "\n",
    "rho_brine_exact = (rhoB-rho_m*(1-phi)) / phi\n",
    "rho_gas_exact = (rhoG-rho_m*(1-phi)) / phi\n",
    "\n",
    "mrkg = dict(marker='o', color='r', ms=12, mfc='none', ls='none')\n",
    "mrkb = dict(marker='o', color='b', ms=12, mfc='none', ls='none')\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(phi_from_gas, rho_g, '-r', label='from gas')\n",
    "plt.plot(phi_from_bri, rho_b, '-b', label='from brine')\n",
    "if i != 'CLASS4':\n",
    "    plt.axvline(sandb.loc[i, 'PHI'], ls='--', color='k')\n",
    "plt.plot(phi, rho_brine_exact, label='{:.3f}'.format(rho_brine_exact), **mrkb)\n",
    "plt.plot(phi, rho_gas_exact, label='{:.3f}'.format(rho_gas_exact), **mrkg)\n",
    "plt.legend()\n",
    "plt.xlabel('phi')\n",
    "plt.ylabel('rho fluids')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **interactive** validate fluid densities for correct porosity prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def phi_from_dens(avoclass, rho_m=2.65):\n",
    "    phi = sandb.loc[avoclass, 'PHI']\n",
    "    rhoB = sandb.loc[avoclass, 'RHO']\n",
    "    rhoG = sandg.loc[avoclass, 'RHO']\n",
    "\n",
    "    rho_b = np.linspace(0.9, 1.2)\n",
    "    rho_g = np.linspace(0.05, 0.6)\n",
    "\n",
    "    phi_from_bri = (rhoB - rho_m)/(rho_b - rho_m)\n",
    "    phi_from_gas = (rhoG - rho_m)/(rho_g - rho_m)\n",
    "\n",
    "    rho_brine_exact = (rhoB-rho_m*(1-phi)) / phi\n",
    "    rho_gas_exact = (rhoG-rho_m*(1-phi)) / phi\n",
    "\n",
    "    mrkg = dict(marker='o', color='r', ms=12, mfc='none', ls='none')\n",
    "    mrkb = dict(marker='o', color='b', ms=12, mfc='none', ls='none')\n",
    "\n",
    "    f, ax = plt.subplots(constrained_layout=True)\n",
    "    ax.plot(phi_from_gas, rho_g, '-r', label='from gas')\n",
    "    ax.plot(phi_from_bri, rho_b, '-b', label='from brine')\n",
    "    ax.axvline(sandb.loc[avoclass, 'PHI'], ls='--', color='k')\n",
    "    ax.plot(phi, rho_brine_exact, label='{:.3f} g/cc'.format(rho_brine_exact), **mrkb)\n",
    "    ax.plot(phi, rho_gas_exact, label='{:.3f} g/cc'.format(rho_gas_exact), **mrkg)\n",
    "    ax.legend(loc='lower right')\n",
    "    ax.set_xlabel('Porosity')\n",
    "    ax.set_ylabel('Fluid density (g/cc)')\n",
    "    ax.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w_avoc = ipw.ToggleButtons(description='AVO Class', options=['CLASS1', 'CLASS2', 'CLASS3'])\n",
    "w_rhom = ipw.FloatSlider(value=2.65, description='Matrix density', max=2.8, min=2.5, step=0.02)\n",
    "ui0 = ipw.VBox([w_avoc, w_rhom])\n",
    "tool0 = ipw.interactive_output(phi_from_dens, {'avoclass': w_avoc, 'rho_m': w_rhom})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display(ui0, tool0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "The fluid densities that I have found can be stored into dictionaries for later use:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rho_b = {}\n",
    "rho_g = {}\n",
    "for i, ss in sandb.iterrows():\n",
    "    if i != 'CLASS4':\n",
    "        rhoB = sandb.loc[i, 'RHO']\n",
    "        rhoG = sandg.loc[i, 'RHO']\n",
    "        rho_b[i] = (rhoB-rho_m*(1-sandb.loc[i, 'PHI']))/sandb.loc[i, 'PHI']\n",
    "        rho_g[i] = (rhoG-rho_m*(1-sandb.loc[i, 'PHI']))/sandb.loc[i, 'PHI']\n",
    "        print('{:s}: rho brine={:.3f}, rho gas={:.3f}'.format(i, rho_b[i], rho_g[i]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And copy the values for class 3 to class 4 since they should occur at similar depths (but this is unclear from Castagna's papers, so we are hunting in the dark right now):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rho_b['CLASS4'] = rho_b['CLASS3']\n",
    "rho_g['CLASS4'] = rho_g['CLASS3']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Obviously when using these class-specific fluid densities starting from either brine or gas values I will now get the same porosities. The following code is redundant but since this is not a paper I will put it here as well:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i, ss in sandb.iterrows():\n",
    "    phi_from_bri = (ss['RHO'] - rho_m)/(rho_b[i] - rho_m)\n",
    "    phi_from_gas = (sandg.loc[i, 'RHO'] - rho_m)/(rho_g[i] - rho_m)\n",
    "    print('{:s}: given phi={:.2f}, phi from brine={:.2f},  phi from gas={:.2f}'.format(i, ss.PHI, phi_from_bri, phi_from_gas))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Porosity from Gassmann's equation\n",
    "\n",
    "I now try a different approach using Gassmann's formula to derive porosity:\n",
    "\n",
    "$$ \\frac{k_{satB}} {k_m - k_{satG}} = \\frac {k_d} {k_m - k_d} + \\frac{k_B}{\\phi \\cdot (k_m - k_B)} $$\n",
    "\n",
    "and\n",
    "\n",
    "$$ \\frac{k_{satG}} {k_m - k_{satG}} = \\frac {k_d} {k_m - k_d} + \\frac{k_G}{\\phi \\cdot (k_m - k_G)} $$\n",
    "\n",
    "These two combined become:\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})}$$\n",
    "\n",
    "An example of how to code the above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k_b = 2.2\n",
    "k_g = 0.02\n",
    "k_m = 37\n",
    "\n",
    "i = 'CLASS3'\n",
    "q, w, e, phi_given = sandg.loc[i].values\n",
    "ksat_g = bulk(q, w, e)\n",
    "\n",
    "q, w, e, _ = sandb.loc[i].values\n",
    "ksat_b = bulk(q, w, e)\n",
    "\n",
    "phi = gassmann_phi(ksat_b, ksat_g, 2.1, k_g, k_m)\n",
    "print('Porosity: calculated={:.4f}, given={:.4f}'.format(phi, phi_given))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I also need to take into account that the fluids' bulk moduli are depth-dependant (and densities too): I should derive them systematically using Batzle-Wang equations but for now with simple trial and error I have derived the following brine bulk moduli:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Applying the same formula above to class 4:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "i = 'CLASS4'\n",
    "q, w, e, phi_given = sandg.loc[i].values\n",
    "ksat_g = bulk(q, w, e)\n",
    "\n",
    "q, w, e, _ = sandb.loc[i].values\n",
    "ksat_b = bulk(q, w, e)\n",
    "\n",
    "phi = gassmann_phi(ksat_b, ksat_g, k_b, k_g, k_m)\n",
    "print('Porosity: calculated={:.4f}, given={:.4f}'.format(phi, phi_given))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **interactive** graphical display of using Gassmann's to get porosity\n",
    "\n",
    "The idea is to compute dry-rock bulk modulus for a range of porosities starting from brine and sand properties. Where they intersect, that's where we expect to find the 'correct' porosity. This works for classes 1-3 as expected, while for class 4 it obviously returns the same porosity I get by inverting Gassmann (0.163) which is quite clearly wrong.\n",
    "\n",
    "Dry-rock bulk moduli are obviously much more prone to error if we use a denser fluid such as brine, that's why the blue lines below also reach negative (clearly incorrect) values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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,
   "metadata": {},
   "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",
    "ui1 = ipw.VBox([w_avoc, w_km, lrbox, reset_button])\n",
    "\n",
    "tool1 = 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,
   "metadata": {},
   "outputs": [],
   "source": [
    "display(ui1, tool1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The same concept above, shown in a single figure:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "phi = np.linspace(.1, .4)\n",
    "\n",
    "f, ax = plt.subplots(ncols=4, figsize=(8,3), sharey=True,constrained_layout=True)\n",
    "for n in range(4):\n",
    "    i = sandb.index[n]\n",
    "    vpG, vsG, rhoG, _ = sandg.loc[i].values\n",
    "    vp_0, vs_0, rho_0, ksat_0, kdry_0 = gassmann(vpG, vsG, rhoG, rho_g[i], k_g, rho_b[i], k_b, k_m, phi)\n",
    "    vpB, vsB, rhoB, _ = sandb.loc[i].values\n",
    "    vp_1, vs_1, rho_1, ksat_1, kdry_1 = gassmann(vpB, vsB, rhoB, rho_b[i], k_b, rho_g[i], k_g, k_m, phi)\n",
    "\n",
    "    ax[n].plot(phi, kdry_0, '-r', label='from gas')\n",
    "    ax[n].plot(phi, kdry_1, '-b', label='from brine')\n",
    "    ax[n].axhline(0, ls='-', color='k', lw=2)\n",
    "    if i != 'CLASS4':\n",
    "        ax[n].axvline(sandb.loc[i, 'PHI'], ls='--', color='k')\n",
    "    ax[n].set_xlabel('phi')\n",
    "    ax[n].set_ylabel('kdry')\n",
    "    ax[n].grid()\n",
    "    ax[n].set_title(i)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "So using Gassmann I get 16% porosity while using densities I get 26% porosity:\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?\n",
    "* consider that a 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": "code",
   "execution_count": null,
   "metadata": {},
   "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,
   "metadata": {},
   "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",
   "metadata": {},
   "source": [
    "Below is a little snippet of code showing fluid replacement on class 4 starting from both brine and gas sand  with the two different porosities; we observe *greater* fluid effect with 16% porosity which is a bit stange to say the least:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "i = 'CLASS4'\n",
    "vpB, vsB, rhoB, _ = sandb.loc[i].values\n",
    "vpG, vsG, rhoG, _ = sandg.loc[i].values\n",
    "\n",
    "print('\\nBRINE --> GAS')\n",
    "print('initial:             vp={:5.0f}, vs={:5.0f}, rho={:5.2f}'.format(vpB, vsB, rhoB))\n",
    "print('expected:            vp={:5.0f}, vs={:5.0f}, rho={:5.2f}'.format(vpG, vsG, rhoG))\n",
    "\n",
    "phi = 0.26\n",
    "vp2, vs2, rho2, ksat2, kdry2 = gassmann(vpB, vsB, rhoB, rho_b[i], k_b, rho_g[i], k_g, k_m, phi)\n",
    "print('gassmann (phi={:.2f}): vp={:5.0f}, vs={:5.0f}, rho={:5.2f}, ksat={:5.2f}, kdry={:4.2f}'.format(phi, vp2, vs2, rho2, ksat2, kdry2))\n",
    "phi = 0.16\n",
    "vp2, vs2, rho2, ksat2, kdry2 = gassmann(vpB, vsB, rhoB, rho_b[i], k_b, rho_g[i], k_g, k_m, phi)\n",
    "print('gassmann (phi={:.2f}): vp={:5.0f}, vs={:5.0f}, rho={:5.2f}, ksat={:5.2f}, kdry={:4.2f}'.format(phi, vp2, vs2, rho2, ksat2, kdry2))\n",
    "\n",
    "print('\\nGAS --> BRINE')\n",
    "print('initial:             vp={:5.0f}, vs={:5.0f}, rho={:5.2f}'.format(vpG, vsG, rhoG))\n",
    "print('expected:            vp={:5.0f}, vs={:5.0f}, rho={:5.2f}'.format(vpB, vsB, rhoB))\n",
    "\n",
    "phi = 0.26\n",
    "vp2, vs2, rho2, ksat2, kdry = gassmann(vpG, vsG, rhoG, rho_g[i], k_g, rho_b[i], k_b, k_m, phi)\n",
    "print('gassmann (phi={:.2f}): vp={:5.0f}, vs={:5.0f}, rho={:5.2f}, ksat={:5.2f}, kdry={:4.2f}'.format(phi, vp2, vs2, rho2, ksat2, kdry2))\n",
    "phi = 0.16\n",
    "vp2, vs2, rho2, ksat2, kdry = gassmann(vpG, vsG, rhoG, rho_g[i], k_g, rho_b[i], k_b, k_m, phi)\n",
    "print('gassmann (phi={:.2f}): vp={:5.0f}, vs={:5.0f}, rho={:5.2f}, ksat={:5.2f}, kdry={:4.2f}'.format(phi, vp2, vs2, rho2, ksat2, kdry2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **interactive** Gassmann and porosity\n",
    "\n",
    "This thingie shows the numerical results of Gassmann's equation on the various AVO class sands depending on porosity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k_m = 37\n",
    "k_w, rho_w = 2.2, 1.0\n",
    "k_g, rho_g = 0.02, 0.32\n",
    "\n",
    "def frm_check_pors(avoclass, phi):\n",
    "    vpG, vsG, rhoG, phi_given = sandg.loc[avoclass]\n",
    "    vpB, vsB, rhoB, _ = sandb.loc[avoclass]\n",
    "\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)\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)\n",
    "\n",
    "    print('Given porosity:    {:.2f}'.format(phi_given))\n",
    "    print('Assigned porosity: {:.2f}'.format(phi))\n",
    "\n",
    "    print('\\nBRINE --> GAS')\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(vp_1, vs_1, rho_1))   \n",
    "    print('gassmann qc: ksat={:5.2f}, kdry={:4.2f}'.format(ksat_1, kdry_1))\n",
    "\n",
    "    itms = (percdiff(vpB, vp_1), percdiff(vsB, vs_1), percdiff(rhoB, rho_1))\n",
    "    print('Difference:        Vp={:+4.0%}, Vs={:+4.0%}, rho={:+4.0%}'.format(*itms))\n",
    "    \n",
    "    \n",
    "    print('\\nGAS --> BRINE')\n",
    "    print('Input (gas):         Vp={:.0f}, Vs={:.0f}, rho={:.2f}'.format(vpG, vsG, rhoG))\n",
    "    print('Reference (brine):   Vp={:.0f}, Vs={:.0f}, rho={:.2f}'.format(vpB, vsB, rhoB))\n",
    "    print('Output FRM brine:    Vp={:.0f}, Vs={:.0f}, rho={:.2f}'.format(vp_0, vs_0, rho_0))\n",
    "    print('gassmann qc: ksat={:5.2f}, kdry={:4.2f}'.format(ksat_0, kdry_0))\n",
    "    \n",
    "    itms = (percdiff(vpG, vp_0), percdiff(vsG, vs_0), percdiff(rhoG, rho_0))\n",
    "    print('Difference:        Vp={:+4.0%}, Vs={:+4.0%}, rho={:+4.0%}'.format(*itms))\n",
    "\n",
    "ipw.interact(frm_check_pors,\n",
    "             avoclass=ipw.ToggleButtons(options=['CLASS1', 'CLASS2', 'CLASS3', 'CLASS4'], value='CLASS4'),\n",
    "             phi=ipw.FloatSlider(min=0.1, max=0.40, step=0.01, value=0.26))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "What is the conclusion from all of this:\n",
    "\n",
    "* let's not use the values for class 4 given by Castagna in both his 1997 and 1998 papers (brine: Vp=2590 m/s, Vs=1060 m/s, rho=2.21 g/cc; gas: Vp=1650 m/s, Vs=1090 m/s, rho=2.07 g/cc) they are wrong.\n",
    "* class 4 is a further subdivision of class 3 so we have to assume similar petrophysical characteristics (and therefore, similar sedimentological setting perhaps?); what is peculiar of this class is what's *above*.\n",
    "* to prove the point above, if we simply model the AVO response of the same class 3 sand overlain by the class 4 shale, we also get something that plots in the class 4 quadrant!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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)\n",
    "\n",
    "ax = classref(plot_brine=False)\n",
    "ax.plot(i_newcl4, g_newcl4, '^', mec='r', mfc='k', ms=14)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**In the plot above, Castagna's class 4 is the red triangle, the \"new\" class 4 is the black triangle outlined in red.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## find a justification for Class 4 numbers with Rock Physics Modeling\n",
    "\n",
    "\n",
    "_Work in progress. The idea is to set up a Monte Carlo simulation to capture the possibile variations of a real sand using rock physics modeling but I only arrived to setting up an initial porosity distribution. See also the [lightning talk notebook](https://github.com/aadm/geophysical_notes/blob/master/avo_class_4_lightning_talk.ipynb)_\n",
    "\n",
    "Make 500 porosity samples following a normal distribution with:\n",
    "\n",
    "$$ \\mu = 0.2 $$\n",
    "$$ \\sigma = 0.03 $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import norm\n",
    "\n",
    "ns = 500\n",
    "phi_dist = norm.rvs(loc=0.2, scale=0.03, size=ns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.Series(phi_dist).plot(kind='density', figsize=(5,3), xlabel='Porosity',\n",
    "                         xlim=(phi_dist.min(), phi_dist.max()),\n",
    "                         title='Porosities',\n",
    "                         rot=45, color='k', lw=4);"
   ]
  }
 ],
 "metadata": {
  "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.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}