{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Synthetic Pumping Test - Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "### What is `timflow.transient`?\n",
    "\n",
    "`timflow.transient` uses the Laplace Transform Analytic Element Method (AEM) to derive a solution to transient flow in multi-layered systems. `timflow` uses discrete features such as Wells and Line-elements to represent real-world aquifer features of interest.\n",
    "\n",
    "\n",
    "In this notebook we demonstrate:\n",
    "\n",
    "* How to specify a simple conceptual model consisting of one confining layer and one well\n",
    "* Simulate the model\n",
    "* Calibrate aquifer parameters by providing data from observation wells"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To achieve this we will create a Model, sample some observations add noise and try to find the model parameters through fitting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We begin by importing the required libraries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 1: Import the Required Libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt  # Plotting library\n",
    "import numpy as np  # Numpy\n",
    "\n",
    "from timflow import transient as tft\n",
    "\n",
    "plt.rcParams[\"figure.figsize\"] = (6, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2: Set Model Parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We set the model parameters with dimensions in **meters** and time in **days**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H = 7  # aquifer thickness [m]\n",
    "k = 70  # hydraulic conductivity [m/d]\n",
    "S = 1e-4  # specific storage fo the aquifer\n",
    "Q = 788  # constant discharge [m/d]\n",
    "d1 = 30  # distance of observation well 1 to pumpinq well\n",
    "d2 = 90  # distance of observation well 2 to pumpinq well (same as Oude Korendijk)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3: Loading Data:\n",
    "\n",
    "- Since we are modelling a synthetic example. We will just borrow the time interval from another pumping test\n",
    "\n",
    "* Data consistis of two columns:\n",
    "    * First column is time in minutes\n",
    "    * Second column is the piezometer level in meters above mean sea level [amsl]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data1 = np.loadtxt(\"../04pumpingtests/data/piezometer_h30.txt\", skiprows=1)\n",
    "t = data1[:, 0] / 60 / 24  # convert min to days"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As seen above, we have loaded the data as a numpy array. That is the preffered format for loading data into `timflow`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 4: Creating a Conceptual Model\n",
    "\n",
    "We will model our aquifer using ModelMaq, which is the 2d model interface from\n",
    "`timflow.transient`. It assumes horizontal stacking of layers consiting of one aquifer\n",
    "and a leaky aquitard."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tft.ModelMaq(\n",
    "    kaq=k,  # Hydraulic Conductivity of the aquifer\n",
    "    # Top and bottom dimensions of the aquifer layer.\n",
    "    # The leaky aquitard can have 0 thickness:\n",
    "    z=[-18, -25],\n",
    "    Saq=1e-4,  # Specific storage of the aquifers\n",
    "    # the minimum time for which heads can be computed after any change\n",
    "    # in boundary condition:\n",
    "    tmin=1e-5,\n",
    "    tmax=1,  # The maximum time for which heads will be computed\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we add the Well element at position (0,0) with screen radius of 10 cm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = tft.Well(\n",
    "    ml,  # Model where we add the well element\n",
    "    xw=0,  # Position x\n",
    "    yw=0,  # Position y\n",
    "    rw=0.1,  # Well radius,\n",
    "    tsandQ=[(0, 788)],  # Tuple describing starting time and discharge\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now 'solve' the model and compute the heads at the two observation locations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml.solve(silent=\"False\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# To compute the heads at the specified time intervals and location, we use the 'head'\n",
    "# method\n",
    "h1 = ml.head(\n",
    "    x=d1,  # location of observation well 1\n",
    "    y=0,\n",
    "    t=t,  # Time array that the heads will be returned for\n",
    ")\n",
    "h2 = ml.head(\n",
    "    x=d2, y=0, t=t\n",
    ")  # Computing heads at distance d2, and time array t for observation well 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can take a look at the data:\n",
    "h1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The head method output a numpy array with dimensions [number of aquifer layers, number of time data]. In this case we only have one row"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 5: Demonstration of Calibration with `timflow`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To demonstrate the capability of `timflow` for deriving aquifer parameters with drawdown data, we will first add noise to the sampled data. We will test the model performance with two standard devitations: 0.02 and 0.05"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# np.savetxt('data/syn_30_0.0.txt', h1[0])\n",
    "# np.savetxt('data/syn_90_0.0.txt', h2[0])\n",
    "# print(h2[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Creating the arrays with noise for sigma = 0.02"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rnd = np.random.default_rng(5)  # Adding a Random seed\n",
    "he12 = h1[0] - rnd.random(len(t)) * 0.02\n",
    "he22 = h2[0] - rnd.random(len(t)) * 0.02\n",
    "np.savetxt(\"data/syn_p30_0.02.txt\", he12)\n",
    "np.savetxt(\"data/syn_p90_0.02.txt\", he22)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Creating the arrays with noise for sigma = 0.05"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rnd = np.random.default_rng(4)  # Adding a Random seed\n",
    "he15 = h1[0] - rnd.random(len(t)) * 0.05\n",
    "he25 = h2[0] - rnd.random(len(t)) * 0.05\n",
    "np.savetxt(\"data/syn_p30_0.05.txt\", he15)\n",
    "np.savetxt(\"data/syn_p90_0.05.txt\", he25)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting and checking the noise added data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.semilogx(t, he12, \".\", label=\"obs at 30 m with $\\sigma$ =0.02\")\n",
    "plt.semilogx(t, h1[0], label=\"modelled heads at 30 m (obs 1)\")\n",
    "plt.semilogx(t, he22, \".\", label=\"obs at 90 m with $\\sigma$ =0.02\")\n",
    "plt.semilogx(t, h2[0], label=\"modelled heads at 90 m (obs 2)\")\n",
    "plt.legend()\n",
    "plt.xlabel(\"time (d)\")\n",
    "plt.ylabel(\"drawdown (m)\")\n",
    "plt.title(\"Drawdown model and observations with noise: $\\sigma$ = 0.02\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.semilogx(t, he15, \".\", label=\"obs at 30 m with $\\sigma$ =0.05\")\n",
    "plt.semilogx(t, h1[0], label=\"modelled heads at 30 m (obs 1)\")\n",
    "plt.semilogx(t, he25, \".\", label=\"obs at 90 m with $\\sigma$ =0.05\")\n",
    "plt.semilogx(t, h2[0], label=\"modelled heads at 90 m (obs 2)\")\n",
    "plt.legend()\n",
    "plt.xlabel(\"time (d)\")\n",
    "plt.ylabel(\"drawdown (m)\")\n",
    "plt.title(\"Drawdown model and observations with noise: $\\sigma$ = 0.05\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 5.1: Calibration of the model using the noisy data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='sig002obs1'></a>\n",
    "Calibrate the model using the $\\sigma $ = 0.02 noisy data from observation well 1.\n",
    "* We calibrate the model by creating a `Calibrate` object with the initial model as input.\n",
    "* Then we set the parameter initial values using the `set_parameter` method\n",
    "* we add the observation data with the `series` method\n",
    "* And fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ca23 = tft.Calibrate(ml)  # class for model calibration\n",
    "# Setting initial parameters values for calibration\n",
    "ca23.set_parameter(name=\"kaq\", initial=10, layers=0)  # Hydraulic Conductivity\n",
    "ca23.set_parameter(name=\"Saq\", initial=1e-3, layers=0)  # Specific Storage\n",
    "\n",
    "ca23.series(  # Adding the observations for calibration\n",
    "    name=\"obs1\",  # Observation well 1\n",
    "    x=d1,  # Location\n",
    "    y=0,\n",
    "    t=t,  # Time Array\n",
    "    h=he12,  # Drawdown noisy data for well 1\n",
    "    layer=0,  # Aquifer layer where we have the observations from\n",
    ")\n",
    "ca23.fit(report=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can already see the results from the calibration if we set `report = True` in the `fit` method. Besides fit statistics, we also get the Variables, with confidence interval and the correlations between variables during fitting process. Note that these confidence intervals can be estimated reasonably well here because the errors are random and uncorrelated, which is generally not the case for pumping tests."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also check the calibrated parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ca23.parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now compare the model performance for both observation wells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"rmse:\", ca23.rmse())  # Return the RMSE error for the calibration\n",
    "h123 = ml.head(d1, 0, t)  # Compute drawdown of the calibrated model for obs 1\n",
    "h223 = ml.head(d2, 0, t)  # Compute drawdown for obs 2\n",
    "plt.semilogx(t, he12, \".\", label=\"obs at 30 m with $\\sigma$ = 0.02\")\n",
    "plt.semilogx(t, h123[0], label=\"modelled drawdown at 30 m\")\n",
    "plt.semilogx(t, he22, \".\", label=\"obs at 90 m with $\\sigma$ = 0.02\")\n",
    "plt.semilogx(t, h223[0], label=\"modelled drawdown at 90 m\")\n",
    "plt.xlabel(\"time (d)\")\n",
    "plt.ylabel(\"drawdown (m)\")\n",
    "plt.title(\n",
    "    \"timflow analysis with synthetic data, \"\n",
    "    \"Calibrated model with $\\sigma$ = 0.02 errors at 30 m.\"\n",
    ")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we do the same procedure as [before](#sig002obs1) to calibrate the model, but now using the observation data from well 2:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ca29 = tft.Calibrate(ml)\n",
    "ca29.set_parameter(name=\"kaq\", initial=10, layers=0)\n",
    "ca29.set_parameter(name=\"Saq\", initial=1e-3, layers=0)\n",
    "ca29.series(name=\"obs2\", x=d2, y=0, t=t, h=he22, layer=0)\n",
    "ca29.fit(report=True)\n",
    "display(ca29.parameters)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can once again check the model performance with the new calibrated model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"rmse:\", ca29.rmse())\n",
    "h129 = ml.head(d1, 0, t)\n",
    "h229 = ml.head(d2, 0, t)\n",
    "plt.semilogx(t, he15, \".\", label=\"obs at 30 m with sig=0.02\")\n",
    "plt.semilogx(t, h129[0], label=\"timflow at 30 m\")\n",
    "plt.semilogx(t, he25, \".\", label=\"obs at 90 m with sig=0.02\")\n",
    "plt.semilogx(t, h229[0], label=\"timflow at 90 m\")\n",
    "plt.legend()\n",
    "plt.xlabel(\"time (d)\")\n",
    "plt.ylabel(\"drawdown (m)\")\n",
    "plt.title(\"timflow analysis with synthetic data, sig=0.02 errors at 90 m.\")\n",
    "plt.legend(loc=\"best\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 5.2: Calibrate with two datasets simultaneously"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`timflow` can also analyse the pumping tests using drawdown data from more than one borehole."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we calibrate the model using the data without error from both observation wells."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ca0 = tft.Calibrate(ml)\n",
    "ca0.set_parameter(name=\"kaq\", initial=10, layers=0)\n",
    "ca0.set_parameter(name=\"Saq\", initial=1e-3, layers=0)\n",
    "ca0.series(name=\"obs1\", x=d1, y=0, t=t, h=h1[0], layer=0)\n",
    "ca0.series(name=\"obs2\", x=d2, y=0, t=t, h=h2[0], layer=0)\n",
    "ca0.fit(report=True)\n",
    "display(ca0.parameters)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The obtained results match exactly the input parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"rmse:\", ca0.rmse())\n",
    "h1n = ml.head(d1, 0, t)\n",
    "h2n = ml.head(d2, 0, t)\n",
    "plt.semilogx(t, h1[0], \"b.\", label=\"obs at 30 m no errors\")\n",
    "plt.semilogx(t, h1n[0], color=\"r\", label=\"timflow at 30 m\")\n",
    "plt.semilogx(t, h2[0], \"g.\", label=\"obs at 90 m no errors\")\n",
    "plt.semilogx(t, h2n[0], color=\"orange\", label=\"timflow at 90 m\")\n",
    "plt.xlabel(\"time (d)\")\n",
    "plt.ylabel(\"drawdown (m)\")\n",
    "plt.title(\"timflow analysis with synthetic data without errors.\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We do the same, but now with drawdowns with errors with $\\sigma=0.02$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ca2 = tft.Calibrate(ml)\n",
    "ca2.set_parameter(name=\"kaq\", initial=10, layers=0)\n",
    "ca2.set_parameter(name=\"Saq\", initial=1e-3, layers=0)\n",
    "ca2.series(name=\"obs1\", x=d1, y=0, t=t, h=he12, layer=0)\n",
    "ca2.series(name=\"obs2\", x=d2, y=0, t=t, h=he22, layer=0)\n",
    "ca2.fit()\n",
    "display(ca2.parameters)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"rmse:\", ca2.rmse())\n",
    "h12 = ml.head(d1, 0, t)\n",
    "h22 = ml.head(d2, 0, t)\n",
    "plt.semilogx(t, he12, \"b.\", label=\"obs at 30 m, sig=0.02\")\n",
    "plt.semilogx(t, h12[0], color=\"r\", label=\"timflow at 30 m\")\n",
    "plt.semilogx(t, he22, \"g.\", label=\"obs at 90 m, sig=0.02\")\n",
    "plt.semilogx(t, h22[0], color=\"orange\", label=\"timflow at 90 m\")\n",
    "plt.xlabel(\"time (d)\")\n",
    "plt.ylabel(\"drawdown (m)\")\n",
    "plt.title(\"timflow analysis with synthetic data and errors with sig=0.02\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Drawdowns with errors with $\\sigma=0.05$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ca5 = tft.Calibrate(ml)\n",
    "ca5.set_parameter(name=\"kaq\", initial=10, layers=0)\n",
    "ca5.set_parameter(name=\"Saq\", initial=1e-3, layers=0)\n",
    "ca5.series(name=\"obs1\", x=d1, y=0, t=t, h=he15, layer=0)\n",
    "ca5.series(name=\"obs2\", x=d2, y=0, t=t, h=he25, layer=0)\n",
    "ca5.fit()\n",
    "display(ca5.parameters)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"rmse:\", ca5.rmse())\n",
    "h15 = ml.head(d1, 0, t)\n",
    "h25 = ml.head(d2, 0, t)\n",
    "plt.semilogx(t, he15, \"b.\", label=\"obs at 30 m\")\n",
    "plt.semilogx(t, h15[0], color=\"r\", label=\"timflow at 30 m\")\n",
    "plt.semilogx(t, he25, \"g.\", label=\"obs at 90 m\")\n",
    "plt.semilogx(t, h25[0], color=\"orange\", label=\"timflow at 90 m\")\n",
    "plt.xlabel(\"time (d)\")\n",
    "plt.ylabel(\"drawdown (m)\")\n",
    "plt.title(\"timflow analysis with synthetic data and errors with sig=0.05\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see the model also performed well for $\\sigma = 0.05$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Final Remarks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example we have succesfully:\n",
    "\n",
    "* Initiated a timflow model instance using the `Model` class\n",
    "* Sampled observation data from the model and defined the synthetic data by adding noise.\n",
    "* Calibrated the model with the `Calibrate` class using one and two calibration wells\n",
    "* Inspected the calibration performance"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}