3-injectivity.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "6iJGW9Ee1v0UFTmVNXBU",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "***\n",
    "\n",
    "# Geothermal Well Test Analysis with Python\n",
    "### Notebook 3: Determining well capacity as injectivity index\n",
    "#### Irene Wallis and Katie McLean \n",
    "#### Software Underground, Transform 2021\n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "5uEy8MtapeDwl2BFVfq5",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "### Google Colab Setup\n",
    "\n",
    "If you are using Google Colab to run this notebook, we assume you have already followed the Google Colab setup steps outlined [here](https://github.com/ICWallis/T21-Tutorial-WellTestAnalysis).\n",
    "\n",
    "Because we are importing data, we need to \"mount your Google Drive\", which is where we tell this notebook to look for the data files. You will need to mount the Google Drive into each notebook.  \n",
    "\n",
    "1. Run the cell below if you are in Google Colab. If you are not in Google Colab, running the cell below will just return an error that says \"No module named 'google'\". If you get a Google Colab error that says \"Unrecognised runtime 'geothrm'; defaulting to 'python3' Notebook settings\", just ignore it. \n",
    "\n",
    "2. Follow the link generated by running this code. That link will ask you to sign in to your google account (use the one where you have saved these tutorial materials in) and to allow this notebook access to your google drive. \n",
    "\n",
    "3. Completing step 2 above will generate a code. Copy this code, paste below where it says \"Enter your authorization code:\", and press ENTER. \n",
    "\n",
    "Congratulations, this notebook can now import data!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "2yc0A14MocmIUBJeyxJS",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "from google.colab import drive\n",
    "drive.mount('/content/drive')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "yUDqmVG46PXRPDkJOroN",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "***\n",
    "\n",
    "# 7. Injectivity index\n",
    "\n",
    "The injectivity index is a measure of the amount of fluid the well can accept during injection. Specifically, it is the change in mass rate (t/hr) per change in pressure (bar), hence it has the units t/hr/bar. \n",
    "\n",
    "If the well is destined to be used as an injection well, then the injectivity is a direct measure of the future well performance (though it must be corrected for different injectate temperatures). \n",
    "\n",
    "If the well is destined to be used as a production well then the injectivity is used to give an indication of future productivity. Productivity index also has the units t/h/bar though it refers to the flow rate out of the well during production, divided by the pressure drop (pressure drawdown) downhole during production. \n",
    " \n",
    "This notebook gives a workflow to: \n",
    "- Import and check data \n",
    "- Select the stable pressure values to use for each flow rate\n",
    "- Create the flow rate vs stable pressure plot (t/h vs bar)\n",
    "- Use linear regression to find the slope (t/h/bar)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "mKBqf6Z7RNrIFeXo8UvP",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "***\n",
    "\n",
    "# 8. Import, munge and check data\n",
    "\n",
    "## 8.1 Use bespoke functions to import and munge data\n",
    "\n",
    "Install all packages required for this notebook. \n",
    "\n",
    "If you do not already have iapws in your environment, then you will need to pip install it. This will need to be done in Google Colab. If you followed the Anaconda setup instructions to make an environment with the environment.yml, then you will not need to do this. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "bHC0oDQlQoHOzqf3ALqI",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "!pip install iapws"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "80RxZjPjjfjCCP2THrOd",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "import iapws # steam tables\n",
    "import openpyxl\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy import stats\n",
    "from datetime import datetime\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.dates as mdates\n",
    "from IPython.display import Image\n",
    "from ipywidgets import interactive, Layout, FloatSlider"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "sdYoDOiIv6MK0laKujbZ",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "def timedelta_seconds(dataframe_col, test_start):\n",
    "    '''\n",
    "    Make a float in seconds since the start of the test\n",
    "\n",
    "    args:   dataframe_col: dataframe column containing datetime objects\n",
    "            test_start: test start time formatted '2020-12-11 09:00:00'\n",
    "\n",
    "    returns: float in seconds since the start of the test\n",
    "    '''\n",
    "    test_start_datetime = pd.to_datetime(test_start)\n",
    "    list = []\n",
    "    for datetime in dataframe_col:\n",
    "        time_delta = datetime - test_start_datetime\n",
    "        seconds = time_delta.total_seconds()\n",
    "        list.append(seconds)\n",
    "    return list\n",
    "\n",
    "\n",
    "\n",
    "def read_flowrate(filename):\n",
    "    ''' \n",
    "    Read PTS-2-injection-rate.xlsx in as a pandas dataframe and munge for analysis\n",
    "\n",
    "    args: filename is r'PTS-2-injection-rate.xlsx'\n",
    "\n",
    "    returns: pandas dataframe with local NZ datetime and flowrate in t/hr\n",
    "    '''\n",
    "    df = pd.read_excel(filename, header=1) \n",
    "    df.columns = ['raw_datetime','flow_Lpm']\n",
    "\n",
    "    list = []\n",
    "    for date in df['raw_datetime']:\n",
    "        newdate = datetime.fromisoformat(date)\n",
    "        list.append(newdate)\n",
    "    df['ISO_datetime'] = list \n",
    "\n",
    "    list = []\n",
    "    for date in df.ISO_datetime:\n",
    "        newdate = pd.to_datetime(datetime.strftime(date,'%Y-%m-%d %H:%M:%S'))\n",
    "        list.append(newdate)\n",
    "    df['datetime'] = list\n",
    "\n",
    "    df['flow_tph'] = df.flow_Lpm * 0.060\n",
    "\n",
    "    df['timedelta_sec'] = timedelta_seconds(df.datetime, '2020-12-11 09:26:44.448')\n",
    "\n",
    "    df.drop(columns = ['raw_datetime', 'flow_Lpm', 'ISO_datetime'], inplace = True)\n",
    "\n",
    "    return df\n",
    "\n",
    "\n",
    "\n",
    "def read_pts(filename):\n",
    "    '''\n",
    "    Read PTS-2.xlsx in as a Pandas dataframe and munge for analysis\n",
    "\n",
    "    args: filename is r'PTS-2.xlsx'\n",
    "\n",
    "    returns: Pandas dataframe with datetime (local) and key coloumns of PTS data with the correct dtype\n",
    "    '''\n",
    "    df = pd.read_excel(filename)\n",
    "\n",
    "    dict = {\n",
    "        'DEPTH':'depth_m',\n",
    "        'SPEED': 'speed_mps',\n",
    "        'Cable Weight': 'cweight_kg',\n",
    "        'WHP': 'whp_barg',\n",
    "        'Temperature': 'temp_degC',\n",
    "        'Pressure': 'pressure_bara',\n",
    "        'Frequency': 'frequency_hz'\n",
    "    }\n",
    "    df.rename(columns=dict, inplace=True)\n",
    "\n",
    "    df.drop(0, inplace=True)\n",
    "    df.reset_index(drop=True, inplace=True)\n",
    "\n",
    "    list = []\n",
    "    for date in df.Timestamp:\n",
    "        newdate = openpyxl.utils.datetime.from_excel(date)\n",
    "        list.append(newdate)\n",
    "    df['datetime'] = list\n",
    "\n",
    "    df.drop(columns = ['Date', 'Time', 'Timestamp','Reed 0',\n",
    "       'Reed 1', 'Reed 2', 'Reed 3', 'Battery Voltage', \n",
    "       'PRT Ref Voltage','SGS Voltage', 'Internal Temp 1', \n",
    "       'Internal Temp 2', 'Internal Temp 3','Cal Temp', \n",
    "       'Error Code 1', 'Error Code 2', 'Error Code 3',\n",
    "       'Records Saved', 'Bad Pages',], inplace = True)\n",
    "    \n",
    "    df[\n",
    "        ['depth_m', 'speed_mps','cweight_kg','whp_barg','temp_degC','pressure_bara','frequency_hz']\n",
    "    ] = df[\n",
    "        ['depth_m','speed_mps','cweight_kg','whp_barg','temp_degC','pressure_bara','frequency_hz']\n",
    "        ].apply(pd.to_numeric)\n",
    "    \n",
    "    df['timedelta_sec'] = timedelta_seconds(df.datetime, '2020-12-11 09:26:44.448')\n",
    "\n",
    "    return df\n",
    "\n",
    "\n",
    "\n",
    "def append_flowrate_to_pts(flowrate_df, pts_df):\n",
    "    '''\n",
    "    Add surface flowrate to pts data\n",
    "\n",
    "    Note that the flowrate data is recorded at a courser time resolution than the pts data\n",
    "    The function makes a linear interpolation to fill the data gaps\n",
    "    Refer to bonus-combine-data.ipynb to review this method and adapt it for your own data\n",
    "\n",
    "    Args:   flowrate and pts dataframes generated by the read_flowrate and read_pts functions\n",
    "\n",
    "    Returns: pts dataframe with flowrate tph added\n",
    "    \n",
    "    '''\n",
    "    flowrate_df = flowrate_df.set_index('timedelta_sec')\n",
    "    pts_df = pts_df.set_index('timedelta_sec')\n",
    "    combined_df = pts_df.join(flowrate_df, how = 'outer',  lsuffix = '_pts', rsuffix = '_fr')\n",
    "    combined_df.drop(columns = ['datetime_fr'], inplace = True)\n",
    "    combined_df.columns = ['depth_m', 'speed_mps', 'cweight_kg', 'whp_barg', 'temp_degC',\n",
    "       'pressure_bara', 'frequency_hz', 'datetime', 'flow_tph']\n",
    "    combined_df['interpolated_flow_tph'] = combined_df['flow_tph'].interpolate(method='linear')\n",
    "    trimmed_df = combined_df[combined_df['depth_m'].notna()]\n",
    "    trimmed_df.reset_index(inplace=True)\n",
    "    return trimmed_df\n",
    "\n",
    "\n",
    "\n",
    "def find_index(value, df, colname):\n",
    "    '''\n",
    "    Find the dataframe index for the exact matching value or nearest two values\n",
    "\n",
    "    args:   value: (float or int) the search term\n",
    "            df: (obj) the name of the dataframe that is searched\n",
    "            colname: (str) the name of the coloum this is searched\n",
    "\n",
    "    returns:  dataframe index(s) for the matching value or the two adjacent values\n",
    "              rows can be called from a df using df.iloc[[index_number,index_number]]\n",
    "    '''\n",
    "    exactmatch = df[df[colname] == value]\n",
    "    if not exactmatch.empty:\n",
    "        return exactmatch.index\n",
    "    else:\n",
    "        lowerneighbour_index = df[df[colname] < value][colname].idxmax()\n",
    "        upperneighbour_index = df[df[colname] > value][colname].idxmin()\n",
    "        return [lowerneighbour_index, upperneighbour_index] \n",
    "\n",
    "    \n",
    "\n",
    "def overview_fig(pts_df,flowrate_df,title=''):\n",
    "    fig, (ax1, ax2, ax3, ax4, ax5, ax6) = plt.subplots(6, 1,figsize=(10,15),sharex=True)\n",
    "    ax1.set_title(title,y=1.1,fontsize=15)\n",
    "\n",
    "    ax1.plot(flowrate_df.datetime, flowrate_df.flow_tph, label='Surface pump flowrate', \n",
    "        c='k', linewidth=0.8, marker='.')\n",
    "    ax1.set_ylabel('Surface flowrate [t/hr]')\n",
    "    ax1.set_ylim(0,150)\n",
    "    \n",
    "    ax2.plot(pts_df.datetime, pts_df.depth_m, label='PTS tool depth', \n",
    "        c='k', linewidth=0.8)\n",
    "    ax2.set_ylabel('PTS tool depth [m]')\n",
    "    ax2.set_ylim(1000,0)\n",
    "    \n",
    "    ax3.plot(pts_df.datetime, pts_df.pressure_bara, label='PTS pressure', \n",
    "        c='tab:blue', linewidth=0.8)\n",
    "    ax3.set_ylabel('PTS pressure [bara]')\n",
    "    \n",
    "    ax4.plot(pts_df.datetime, pts_df.temp_degC, label='PTS temperature', \n",
    "        c='tab:red', linewidth=0.8)\n",
    "    ax4.set_ylabel('PTS temperature')\n",
    "    \n",
    "    ax5.plot(pts_df.datetime, pts_df.frequency_hz, label='PTS impeller frequency', \n",
    "        c='tab:green', linewidth=0.8)\n",
    "    ax5.set_ylim(-30,30)\n",
    "    ax5.set_ylabel('PTS impeller frequency [hz]')\n",
    "    # 1 hz = 60 rpm\n",
    "\n",
    "    ax6.plot(pts_df.datetime, pts_df.speed_mps, label='PTS tool speed', \n",
    "        c='tab:orange', linewidth=0.8)\n",
    "    ax6.set_ylim(-2,2)\n",
    "    ax6.set_ylabel('PTS tool speed [mps]')\n",
    "    \n",
    "    ax6.set_xlabel('Time [hh:mm]')\n",
    "    plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%H:%M'))\n",
    "    \n",
    "    for ax in [ax1,ax2,ax3,ax4,ax5,ax6]:\n",
    "        ax.grid()\n",
    "    \n",
    "    return plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "Qm0lajQWUB998JCUjf3F",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "The cells below will take a little while to run because it includes all steps required to import and munge the data (i.e., everything we did in notebook 1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "A4W82lnuPbxjRnkih13j",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "# Use this method if you are running this notebook in Google Colab\n",
    "flowrate = read_flowrate(r'/content/drive/My Drive/T21-Tutorial-WellTestAnalysis-main/Data-FlowRate.xlsx')\n",
    "\n",
    "# Use this method if you are running this notebook locally (Anaconda)\n",
    "#flowrate = read_flowrate(r'Data-FlowRate.xlsx')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "2xKu79NCPc3Uc98xcYVG",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "# Use this method if you are running this notebook in Google Colab\n",
    "pts = read_pts(r'/content/drive/My Drive/T21-Tutorial-WellTestAnalysis-main/Data-PTS.xlsx')\n",
    "\n",
    "# Use this method if you are running this notebook locally (Anaconda)\n",
    "#pts = read_pts(r'Data-PTS.xlsx')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "2On4PwPSZAeAnXmcFrKJ",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "We use surface pump flow rate and downhole pressure together for our analysis, so we append flow rate to our pts dataframe to make the method simpler. This is done using the append_flowrate_to_pts function from the collection of helper functions above. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "rb9a3TbH9IZNzofj8cax",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "pts = append_flowrate_to_pts(flowrate, pts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "rXJCeaRg7QyrUjPaOzhS",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "## 8.2 Import stable temperature/pressure data for comparison\n",
    "\n",
    "After drilling and injection testing, geothermal wells heat up over time. \n",
    "So for comparison, we have provided a temperature log acquired when the well had reached its stable condition.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "MO0wNqo4xUJIbFco80QX",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "# Use this method if you are running this notebook in Google Colab\n",
    "heating_37days = pd.read_csv(r'/content/drive/My Drive/T21-Tutorial-WellTestAnalysis-main/Data-Temp-Heating37days.csv') \n",
    "\n",
    "# Use this method if you are running this notebook locally (Anaconda)\n",
    "#heating_37days = pd.read_csv('Data-Temp-Heating37days.csv')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "5pbUb1lrKx07UgA9wFUm",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "# Convert bar gauge to bar atmosphere\n",
    "\n",
    "heating_37days['pressure_bara'] = heating_37days.pres_barg - 1\n",
    "\n",
    "heating_37days.head(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "HfECgu82Qhka9Q1QKiAc",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "## 8.3 Check the data\n",
    "\n",
    "It is good practice to check your data after import. \n",
    "\n",
    "You can use the Pandas methods listed in Section 2.1.1 (1-intro-and-data.ipynb) to check your data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "DBoGLVWSGnjUmCgeF6IN",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "pts.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "0oyWggUcUEwCyKnGzHG7",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "flowrate.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "pGpKRfLajjDeN7BQxMHN",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "heating_37days.head(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "qOruxRo8LVnLUi69LTzi",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "We made the plot below to check that we have imported what we expected to import. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "qA75D5hJTEMpEaEyrpa1",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(1, 2,figsize=(24,8),sharey=True)\n",
    "\n",
    "ax1.scatter(pts.pressure_bara, pts.depth_m, c = pts.timedelta_sec, s = 5, linewidths = 0)\n",
    "\n",
    "ax1.plot(heating_37days.pressure_bara, heating_37days.depth_m, c = 'k', label = 'Heating 37 days')\n",
    "\n",
    "ax1.legend()\n",
    "\n",
    "ax2.scatter(pts.datetime, pts.depth_m, c = pts.timedelta_sec, s = 5, linewidths = 0)\n",
    "\n",
    "ax3 = ax2.twinx()\n",
    "ax3.plot(flowrate.datetime, flowrate.flow_tph, \n",
    "    c='k', linestyle = '-', linewidth = 3, alpha = 0.3, \n",
    "    label='Surface pump flowrate')\n",
    "\n",
    "ax1.set_ylim(1000,0)  # 940,400 \n",
    "ax1.set_xlim(0,55)     # 20,500\n",
    "\n",
    "ax1.set_xlabel('Pressure [bara]')\n",
    "ax1.set_ylabel('Depth [m]')\n",
    "\n",
    "ax2.set_xlabel('Time [hh:mm]')\n",
    "ax2.xaxis.set_major_formatter(mdates.DateFormatter('%H:%M'))\n",
    "\n",
    "ax3.set_ylabel('Flowrate [t/hr]')\n",
    "\n",
    "for ax in [ax1, ax2]:\n",
    "    ax.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "HpdPHqZesnUuvImI5aiX",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "***\n",
    "\n",
    "# 9. Pressure inside geothermal wells\n",
    "\n",
    "Note that the pressure measured inside the well does not equate to the pressure in the reservoir. It varies depending on the liquid level and the density (temperature & phase) of the fluid inside the well. As the well heats up after drilling, pressure profiles will appear to pivot around a single point. That pressure of the pivot point equates to the reservoir pressure at the depth of pivoting.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "wYFX7D7BaS997ED83Wh4",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "***\n",
    "\n",
    "# 10. Select the most stable pressure values for each flow rate \n",
    "\n",
    "## 10.1 Interactive plot with ipywidgets\n",
    "\n",
    "Use the interactive plot to select the times within the completion test program when the pressure down-hole is most likely to be stable for a given flow rate. \n",
    "\n",
    "The most stable pressure is usually just after the pts tool returned to the programmed hanging depth that was for pressure transient after the well passes are complete and before the pump rate is changed. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "xuSK8ON5isw9Wony7y3n",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "min_timestamp = pts.timedelta_sec.iloc[0]\n",
    "max_timestamp = pts.timedelta_sec.iloc[-1]\n",
    "\n",
    "def subselect_plot(first_value, second_value, third_value):\n",
    "    f,ax1 = plt.subplots(1,1, figsize = (20,6))\n",
    "    ax1.plot(pts.timedelta_sec, pts.depth_m, c = 'k', label = 'PTS tool depth')\n",
    "    ax2 = ax1.twinx()\n",
    "    ax2.plot(flowrate.timedelta_sec, flowrate.flow_tph, c='k', linestyle = ':', label='Surface pump flowrate')\n",
    "    ymin = pts.depth_m.min()\n",
    "    ymax = pts.depth_m.max() + 100\n",
    "    ax1.vlines(first_value, ymin, ymax, color='tab:green')\n",
    "    ax1.vlines(second_value, ymin, ymax, color='tab:orange')\n",
    "    ax1.vlines(third_value, ymin, ymax, color='tab:red')\n",
    "    ax1.set_ylim(pts.depth_m.max() + 100, 0)\n",
    "    ax1.set_xlabel('Time elapsed since the test started [sec]')\n",
    "    ax2.set_ylabel('Surface pump flowrate [t/hr]')\n",
    "    ax1.set_ylabel('PTS tool depth [mMD]')\n",
    "\n",
    "result = interactive(subselect_plot,\n",
    "         \n",
    "         first_value = FloatSlider\n",
    "         (\n",
    "             value = (max_timestamp - min_timestamp)/6 + min_timestamp,\n",
    "             description = '1st value',\n",
    "             min = min_timestamp, \n",
    "             max = max_timestamp, \n",
    "             step = 10, \n",
    "             continuous_update=False,\n",
    "             layout = Layout(width='80%'),\n",
    "             ),\n",
    "          \n",
    "          second_value = FloatSlider\n",
    "          (\n",
    "             value = (max_timestamp - min_timestamp)/4 + min_timestamp, \n",
    "             description = '2nd value',\n",
    "             min = min_timestamp, \n",
    "             max = max_timestamp, \n",
    "             step = 10, \n",
    "             continuous_update=False,\n",
    "             layout = Layout(width='80%')\n",
    "             ),\n",
    "\n",
    "          third_value = FloatSlider\n",
    "          (\n",
    "             value = (max_timestamp - min_timestamp)/2 + min_timestamp, \n",
    "             description = '3rd value',\n",
    "             min = min_timestamp, \n",
    "             max = max_timestamp, \n",
    "             step = 10, \n",
    "             continuous_update=False,\n",
    "             layout = Layout(width='80%')\n",
    "             )\n",
    ")\n",
    "\n",
    "display(result);\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "xx06lhQWwgjPhYs7ixo7",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "## 10.2 Call results from the interactive plot\n",
    "\n",
    "After you place the 1st (green), 2nd (orange), and 3rd (red) line locations, run the cell below to call the results. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "95nWlAGj9C9g5kIdjSKl",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "# extract pressure and flow rate at the marked points\n",
    "print(\n",
    "    'first_timestamp =',result.children[0].value, \n",
    "    '\\nsecond_timestamp =', result.children[1].value, \n",
    "    '\\nthird_timestamp =', result.children[2].value, \n",
    "    )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "ayifs1FmLCEbbmaOo6hX",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "## 10.3 Store analysis\n",
    "\n",
    "Because result.children will change each time you move the sliders in the plot above or re-run this Jupyter Notebook, we copy-paste our selection below. This records your choice and will be the values you do the rest of the interpretation with. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "ZulLwvtPRngBKEad39pc",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "#### Our selected data and some metadata\n",
    " \n",
    "The third value was selected before the tool made passes because the pumps were shut off so quickly \n",
    "that it is difficult to reliably pick the value after the tool passes at this plot scale. \n",
    "\n",
    "If we were concerned that the pressure was not stable before the logging passes, \n",
    "we could take the time to adjust the scale of our plot and pick the correct location."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "HGkb9vgz9BX90XXGbBnU",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "# define a set of objects using our elapsed time\n",
    "# we will use these to extract the data\n",
    "\n",
    "first_timestamp = 3620.0 \n",
    "second_timestamp = 12370.0 \n",
    "third_timestamp = 18740.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "uvY08GQfo4IZfjUx4Jkp",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "***\n",
    "\n",
    "# 11. Generate an array of stable pressure and flow rate\n",
    "\n",
    "Recall that we used the append_flowrate_to_pts function from our utilities.py in Section 8.1 above to append surface pump flow rate to our pts dataframe.\n",
    "\n",
    "## 11.1 Make new dataframes containing the pts log values at our timestamps\n",
    "\n",
    "Making our new dataframes is a two-step process:\n",
    "1. Find the index value where the timestamp is either an exact match to the value we pass in or is the nearest neighbour above and below using the helper function at the start of this notebook\n",
    "2. Make new dataframes using the .iloc method and the timestamp returned by step 1 (i.e., the find_index function)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "Pswur9Zh6qQrdCNQVPzz",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "first_pts = pts.iloc[find_index(first_timestamp, pts, 'timedelta_sec')]\n",
    "second_pts = pts.iloc[find_index(second_timestamp, pts, 'timedelta_sec')]\n",
    "third_pts = pts.iloc[find_index(third_timestamp, pts, 'timedelta_sec')]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "EHHyQZ5QjxaabJ2lXVcH",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "first_pts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "1q3rqBwNw0evW28Xgg9j",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "second_pts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "DeBRHC89PNL9y8lWb8x5",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "third_pts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "gBVdIH8NzTVhQElaKtHM",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "## 11.2 Make pressure and flow rate arrays\n",
    "\n",
    "Now we either use the exact match value or the mean of the two neighbouring values to make an array of pressure and flow rate that we can use in our injectivity index analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "ORTWVJoN5HXjvrsTQFnB",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "# make array of mean pressure values\n",
    "\n",
    "list = []\n",
    "for df in [first_pts, second_pts, third_pts]:\n",
    "    mean_pressure = df['pressure_bara'].mean()\n",
    "    list.append(mean_pressure)\n",
    "pressure_array = np.array(list)\n",
    "\n",
    "print('pressure data =', pressure_array)\n",
    "print('object type =', type(pressure_array))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "wuo1cKaHbAt5jUIFQgcM",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "# make array of mean flowrate values\n",
    "\n",
    "list = []\n",
    "for df in [first_pts, second_pts, third_pts]:\n",
    "    mean_flowrate = df['interpolated_flow_tph'].mean()\n",
    "    list.append(mean_flowrate)\n",
    "flowrate_array = np.array(list)\n",
    "\n",
    "print('flowrate data =', flowrate_array)\n",
    "print('object type =', type(flowrate_array))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "lJ2DIRsF7VGCz8VUAn5f",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "***\n",
    "\n",
    "# 12. Generate a linear regression to find injectivity index\n",
    "\n",
    "The injectivity is a measure of total well capacity during injection. It is calculated from the change in pressure (bar) that occurs in response to changing the mass rate (t/hr) that is injected into the well. \n",
    "\n",
    "If there is nothing unusual going on downhole, then we can find this with simple linear regression. There are many cases where a single linear model is not the best approach, such as where changes in pressure or thermal conditions in the well change the permeability. In these cases, a two-slope approach may be better. \n",
    "\n",
    "## 12.1 Make the linear model\n",
    "\n",
    "There are many ways to do a simple linear regression. We selected the stats.linregress method from Scipy because it is relatively simple to use. It's also fast because it is tooled specifically for our two-variable use case. The stats.linregress method returns the slope, intersect and R value that we need for our analysis. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "cJNE7hdVvIE8oYMXf1dR",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "linear_model = stats.linregress(pressure_array, flowrate_array) \n",
    "\n",
    "linear_model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "TVla75NTuyUX1qfVTE9u",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "# Define some sensibly named objects for reuse below\n",
    "\n",
    "slope = linear_model[0]\n",
    "intercept = linear_model[1]\n",
    "rvalue = linear_model[2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "ifTNzsjj2SxcccY9W18M",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "# Print a nicely formatted string describing of our model\n",
    "\n",
    "print(\"The linear model for our data is y = {:.5} + {:.5} * pressure\"\n",
    "      .format(intercept, slope))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "xcnGi9tggYnhfUgXieyc",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "## 12.2 Plot the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "WODF6RAAHji9MWZYXUhN",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "f,ax = plt.subplots(1,1,figsize=(6,4))\n",
    "\n",
    "ax.scatter(pressure_array, flowrate_array, \n",
    "           s=100, c='k', label = 'Data')\n",
    "\n",
    "ax.plot(pressure_array, # x values\n",
    "    slope * pressure_array + intercept, # use the model to generate y values\n",
    "    color='tab:orange', linestyle='-', linewidth=4, alpha=0.5, label='Linear fit')\n",
    "\n",
    "\n",
    "ax.set_title(\"Sample well injectivity index = {:} t/hr.bar\".format(round(slope)))\n",
    "\n",
    "#ax.set_ylim(0,130)\n",
    "#ax.set_xlim(0,40)\n",
    "\n",
    "ax.legend()\n",
    "ax.set_ylabel('Rate [t/hr]')\n",
    "ax.set_xlabel('Pressure [bara]')  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "4NMvebhNSQxhD4xSzwoD",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "***\n",
    "\n",
    "# 13. Interpreting injectivity index\n",
    "\n",
    "### 13.1 Understanding the scale of injectivity index\n",
    "\n",
    "Based on data from New Zealand geothermal fields, pre-eminent reservoir engineer Malcolm Grant determined that injectivity index forms a lognormal distribution with a median value around 20 t/hr.bar (Figure 5)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "1Tl17lWE3qtSjy5Xlpjg",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     },
     "outputId": null
    }
   },
   "outputs": [],
   "source": [
    "Image('https://raw.githubusercontent.com/ICWallis/T21-Tutorial-WellTestAnalysis/main/Figures/Figure5.png',width = 500,)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "RSLHZ5ynZRX9pz8DiQOK",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "_Figure 5: Like many natural phenomena, geothermal well injectivity has a lognormal distribution. Figure adapted from [Grant (2008)](https://pangea.stanford.edu/ERE/pdf/IGAstandard/SGW/2008/grant.pdf)._\n",
    "\n",
    "There is no 1:1 relationship between the injectivity index (II) and productivity index (PI) of a geothermal well. Each reservoir is subject to a unique set of conditions that influence the relationship between these indices, such as depth to the liquid level, and enthalpy. \n",
    "\n",
    "The table below reflects a general rule of thumb to be used for a reservoir where the local relationship between II and PI has not already been established.\n",
    "\n",
    "|Permeability magnitude | II \\[t/hr.bar\\] |\n",
    "| --- | --- |\n",
    "| Very low permeability (near-conductive) | < 1 |\n",
    "| Poor permeability, usable in special cases | 1 - 5 |\n",
    "| Likely productive | 5 - 20 |\n",
    "| Median production well | 20 |\n",
    "| Reliably economic production for well if T > 250$^{\\circ}$C | 20-50 |\n",
    "| High permeability | 50 - 100 |\n",
    "| Very high permeability (unusual) | > 100 |\n",
    "\n",
    "\n",
    "### 13.2 Injectivity index and temperature\n",
    "\n",
    "The injectivity index is only truly valid for the injectate temperature at which it was measured (i.e., ambient temperature for completion testing). If you inject at a higher temperature (i.e., injection of separated geothermal water or condensate in operational injection wells), then the injectivity index will be less. Empirical corrections are available to adjust injectivity index for temperature - for more information refer to [Siega et al. (2014)](https://www.geothermal-energy.org/pdf/IGAstandard/NZGW/2014/109.Siega.pdf)\n",
    "\n",
    "\n",
    "# 14. Injectivity in our case study well\n",
    "\n",
    "With an II of 252 t/hr.bar, our case study well is extremely permeable! Yah! \n",
    "\n",
    "It is worth noting that as II increases, so does the level of uncertainty in the actual value. This is because the pressure changes resulting from the flow changes become so small. This is a good \"problem\" to have. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "A6n7OYBhdxS1DUb0sHgV",
      "project": "mmReuqVTAa9JzPpNr22I",
      "version": 1
     }
    }
   },
   "source": [
    "***\n",
    "\n",
    "### Cited references\n",
    "\n",
    "Grant, M.A. (2008) Decision tree analysis of possible drilling outcomes to optimise drilling decisions: Proceedings, Thrity-Third Workshop of Geothermal Reservoir Engineering. Stanford University, CA. \n",
    "\n",
    "Siega, C., Grant, M.A., Bixley, P. Mannington, W. (2014) Quantifying the effect of temperature on well injectivity: Proceedings 36th New Zealand Geothermal Workshop. Auckland, NZ.\n",
    "\n",
    "Zarrouk, S.J. and McLean, K. (2019): Geothermal well test analysis: fundamentals, applications, and advanced techniques. 1st edition, Elsevier. \n",
    "\n",
    "***\n",
    "\n",
    "© 2021 [Irene Wallis](https://www.cubicearth.nz/) and [Katie McLean](https://www.linkedin.com/in/katie-mclean-25994315/) \n",
    "\n",
    "Licensed under the Apache License, Version 2.0\n",
    "\n",
    "***"
   ]
  }
 ],
 "metadata": {
  "iooxa": {
   "id": {
    "block": "CZKZK8UMST12MbM5OpK9",
    "project": "mmReuqVTAa9JzPpNr22I",
    "version": 1
   }
  },
  "kernelspec": {
   "display_name": "geothrm",
   "language": "python",
   "name": "geothrm"
  },
  "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}