{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Examples reservoir models\n", "*Mark Bakker, Delft University of Technology*\n", "\n", "**Warning:** \n", "This is a proof-of-concept of the Reservoir model. The Reservoir model is still under development and not yet meant for general application. No noise model has been used in the models presented in this notebook. This might lead to wrong estimates of the parameter uncertainties." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy.integrate import solve_ivp\n", "import pandas as pd\n", "import pastas as ps\n", "ps.show_versions()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Single reservoir model as compared to exponential response function" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "jupyter": { "source_hidden": true } }, "outputs": [], "source": [ "def arrow(xystart, xyend, text=\"\", arrow=\"<-\", color='k', **kwds):\n", " plt.annotate(text,\n", " xy=xystart, \n", " xytext=xyend, \n", " arrowprops=dict(arrowstyle=arrow, shrinkA=0, shrinkB=0, color=color),\n", " color=color,\n", " **kwds)\n", "\n", "plt.figure(figsize=(4, 5))\n", "plt.subplot(111, aspect=1)\n", "plt.plot([0.2, 0.2, 1, 1, 1.1, 1.1], [0.8, 0, 0, 0.2, 0.2, 0.4], 'k')\n", "plt.plot([1.06, 1.06, 1, 1], [0.4, 0.24, 0.24, 0.8], 'k' )\n", "plt.plot([0.2, 1], [0.6, 0.6], 'C0', lw=2)\n", "arrow((1.08, 0.35), (1.08, 0.48), arrow='<->', color='C0')\n", "plt.text(0.41, 0.3, '$h(t)$')\n", "plt.text(1.06, 0.5, '$Q(t)=(h-d)/c$')\n", "arrow((1.15, 0), (1.15, 0.4), arrow='<->')\n", "plt.text(1.16, 0.2, '$d$')\n", "arrow((0.6, 0.75), (0.6, 0.65), color='C0')\n", "plt.text(0.6, 0.8, '$R(t)=P(t)-fE(t)$', ha='center')\n", "arrow((0.4, 0), (0.4, 0.6), arrow='<->')\n", "plt.xlim(0, 1.2)\n", "plt.axis('off');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\\begin{equation}\n", "\\frac{\\text{d}h}{\\text{d}t} = \\frac{R}{S} - \\frac{h-d}{cS}\n", "\\end{equation}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Implicit solution:\n", "\\begin{equation}\n", "h_t = \\frac{h_{t-\\Delta t} + R\\Delta t / S + \\Delta t / (cS) d}{1 + \\Delta t / (cS)}\n", "\\end{equation}" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "hobs = pd.read_csv('data/head_nb1.csv', parse_dates=['date'],\n", " index_col='date').squeeze()\n", "prec = pd.read_csv('data/rain_nb1.csv', parse_dates=['date'],\n", " index_col='date').squeeze()\n", "evap = pd.read_csv('data/evap_nb1.csv', parse_dates=['date'],\n", " index_col='date').squeeze()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ml = ps.Model(hobs, name=\"reservoirtest\", constant=False, noisemodel=False)\n", "rsv = ps.ReservoirModel([prec, evap], ps.reservoir.Reservoir1, 'reservoir', ml.oseries.series.mean())\n", "ml.add_stressmodel(rsv)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#%timeit ml.solve(noise=False, solver=ps.LmfitSolve, tmin='1995', tmax='2005', report=False)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ml.solve(noise=False, solver=ps.LmfitSolve, tmin='1995', tmax='2005')\n", "hreservoir = ml.simulate()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Comparison with exponential response function\n", "\\begin{equation}\n", "\\theta_{step}=A(1-\\text{e}^{-t/a})\n", "\\end{equation}\n", "with\n", "\\begin{equation}\n", "A = c\n", "\\end{equation}\n", "\\begin{equation}\n", "a = cS\n", "\\end{equation}" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "mlbase = ps.Model(hobs)\n", "sm = ps.RechargeModel(prec, evap, ps.Exponential, name='rech')\n", "mlbase.add_stressmodel(sm)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#%timeit mlbase.solve(solver=ps.LmfitSolve, noise=False, tmin='1995', tmax='2005', report=False)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "mlbase.solve(solver=ps.LmfitSolve, noise=False, tmin='1995', tmax='2005')\n", "hbase = mlbase.simulate()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "hreservoir.plot(label='reservoir')\n", "hbase.plot(ls='--', label='exponential')\n", "plt.legend();" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "(hreservoir - hbase).plot();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Reservoir model with overflow as compared to Tarso" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "jupyter": { "source_hidden": true } }, "outputs": [], "source": [ "plt.figure(figsize=(4, 5))\n", "plt.subplot(111, aspect=1)\n", "plt.plot([0.2, 0.2, 1, 1, 1.1, 1.1], [0.8, 0, 0, 0.2, 0.2, 0.4], 'k')\n", "plt.plot([1.06, 1.06, 1, 1, 1.05], [0.4, 0.24, 0.24, 0.7, 0.7], 'k' )\n", "plt.plot([1.05, 1, 1], [0.74, 0.74, 0.8], 'k')\n", "plt.plot([0.2, 1], [0.6, 0.6], 'C0', lw=2)\n", "arrow((1.08, 0.35), (1.08, 0.48), arrow='<->', color='C0')\n", "plt.text(0.41, 0.3, '$h(t)$')\n", "plt.text(1.06, 0.5, '$Q(t)=(h-d)/c$')\n", "arrow((1.15, 0), (1.15, 0.4), arrow='<->')\n", "plt.text(1.16, 0.2, '$d$')\n", "arrow((0.6, 0.75), (0.6, 0.65), color='C0')\n", "plt.text(0.6, 0.8, '$R(t)=P(t)-fE(t)$', ha='center')\n", "arrow((0.4, 0), (0.4, 0.6), arrow='<->')\n", "arrow((1, 0.72), (1.1, 0.72), color='C0')\n", "plt.text(1.11, 0.72, '$Q_2(t)=(h-d_2)/c_2$ \\n if $h>d_2$', va='center')\n", "plt.xlim(0, 1.2)\n", "plt.axis('off');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\\begin{equation}\n", "\\frac{\\text{d}h}{\\text{d}t} = \\frac{R}{S} - \\frac{h-d}{cS} \\qquad h