{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(_huitrelle_sompuis)=\n",
    "# How to simulate a threshold effect in groundwater level time series?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example is taken from the Gardenia tutorial {cite:p}`2013:thiery_didacticiel`.\n",
    "\n",
    "This example shows how to simulate groundwater level time series with overflow threshold, fractures or permeability increase near the surface, using two groundwater flow components.\n",
    "\n",
    "The Huitrelle River basin is located in the Marne French department between the Troyes and Châlons-en-Champagne cities in the Chalk Champagne. Its drainage area at the gauging station Huitrelle at Lhuitre (H1503510) is equal to 524 km².\n",
    "\n",
    "The Sompuis piézometer (BSS : 02255X0003) is located at about 10 km from the middle of the basin. Time series show a very marked threshold effect on the groundwater levels at an altitude of around 145 m NGF.\n",
    "\n",
    "The following data are available:\n",
    "\n",
    "- Observed mean decadal river flow of the Huitrelle at Lhuitre from July 1997 to January 2008 (m3/s)\n",
    "- Observed groundwater levels at Sompuis (m NGF) from 1969 to 2008\n",
    "- Decadal rainfall and potential evapotranspiration from 1969 to 2008 (mm/decade)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## TOML configuration file\n",
    "\n",
    "```{literalinclude} model.toml\n",
    ":language: toml\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optimization without threshold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import rameau as rm\n",
    "\n",
    "# Load model from a toml file\n",
    "model = rm.Model.from_toml(f\"model.toml\")\n",
    "\n",
    "# Run a simulation with parameters defined in the toml file\n",
    "sim = model.run_optimization()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Optimized parameters are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = sim.tree.watersheds[0]\n",
    "print(\"Base level = \", round(w.groundwater.base_level.value, 4), \" m NGF\")\n",
    "print(\"Storage coefficient = \", round(w.groundwater.storage.coefficient.value, 4), \" %\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Metrics are a priori correct:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print the riverflow metrics\n",
    "scores = sim.get_metrics(\"riverflow\")\n",
    "print(scores)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print the watertable metrics\n",
    "scores = sim.get_metrics(\"watertable\")\n",
    "print(scores)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We plot river flows and groundwater levels time series from this first optimization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the riverflow simulation\n",
    "riv_sim = sim.get_output(\"riverflow\")\n",
    "\n",
    "# Get the riverflow observation\n",
    "riv_obs = model.get_input(\"riverobs\")\n",
    "\n",
    "# Plot river flow\n",
    "df = pd.DataFrame(\n",
    "    {\n",
    "        \"Obs\":riv_obs.iloc[:, 0],\n",
    "        \"Sim\":riv_sim.iloc[:, 0],\n",
    "    }\n",
    ")\n",
    "df.loc[\"1976-01-01\":, :].plot(\n",
    "    grid=True,\n",
    "    title=\"Simulated and observed river flows\",\n",
    "    ylabel=\"m3/s\",\n",
    "    style=[\"-\", \"+\"]\n",
    ")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the watertable simulation\n",
    "wtl_sim = sim.get_output(\"watertable\")\n",
    "\n",
    "# Get the watertable observation\n",
    "wtl_obs = model.get_input(\"groundwaterobs\")\n",
    "\n",
    "# Plot watertable\n",
    "df = pd.DataFrame({\"Sim\":wtl_sim.iloc[:, 0], \"Obs\":wtl_obs.iloc[:, 0]})\n",
    "df.loc[\"1976-01-01\":, :].plot(\n",
    "    grid=True,\n",
    "    title=\"Simulated and observed groundwater levels\",\n",
    "    ylabel=\"m NGF\",\n",
    "    style=['-', '+']\n",
    ")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This plot shows that the simulated groundwater levels exeed the observed maximal level of around 145 m NGF (1983, 1988, 1994, 2001). We propose to perform a second optimization run by letting the possibility for the model to have an overflow threshold corresponding approximatively to this 145 m NGF value."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optimization with groundwater threshold\n",
    "\n",
    "This is exactly the same configuration but now by introducting a threshold\n",
    "in the groundwater reservoir, i.e. a flow that appears only if the groundwater \n",
    "reservoir level exceeds a threshold value. \n",
    "\n",
    "The overflow threshold value corresponding to the 145 m NGF is deduced as\n",
    "follows:\n",
    "\n",
    "{math}`H_o = (145 - H_b)S`\n",
    "\n",
    "{math}`H_o = (145 - 132.0037) \\times 0.000979 = 12.7 \\times 10^{-3} \\ \\mathrm{mNGF} = 12.7 \\ \\mathrm{ mm}`\n",
    "\n",
    "We modify the model to start from this value during the optimization. The initial value for the overflow halflife time is set to 0.2 months. This quick value allows to limit the rising of groundwater levels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run a seoncd optimization with optimization threshold parameters now optimized\n",
    "model.tree.watersheds[0].groundwater.reservoirs[0].overflow.threshold = rm.Parameter(\n",
    "    value = 12.7, lower=0.0, upper=500, opti=True\n",
    ")\n",
    "model.tree.watersheds[0].groundwater.reservoirs[0].overflow.halflife = rm.Parameter(\n",
    "    value = 0.2, lower=0.05, upper=15, opti=True\n",
    ")\n",
    "sim = model.run_optimization()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Optimized parameters are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = sim.tree.watersheds[0]\n",
    "print(\"Base level = \", round(w.groundwater.base_level.value, 4), \" m NGF\")\n",
    "print(\"Storage coefficient = \", round(w.groundwater.storage.coefficient.value, 4), \" %\")\n",
    "print(\"Overflow threshold = \", round(w.groundwater.reservoirs[0].overflow.threshold.value, 4), \" mm\")\n",
    "print(\"Overflow halflife time = \", round(w.groundwater.reservoirs[0].overflow.halflife.value, 4), \" mm\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Metrics for river flow stays correct while metrcis for grondwater levels are improved:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print the riverflow metrics\n",
    "scores = sim.get_metrics(\"riverflow\")\n",
    "print(scores)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print the watertable metrics\n",
    "scores = sim.get_metrics(\"watertable\")\n",
    "print(scores)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We plot river flows and groundwater level time series. As expected the simulated groundwater levels do not exceed the threshold value noted in the observations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the riverflow simulation from optimization 1 \n",
    "riv_sim = sim.get_output(\"riverflow\")\n",
    "\n",
    "# Plot river flow\n",
    "df = pd.DataFrame(\n",
    "    {\n",
    "        \"Obs\":riv_obs.iloc[:, 0],\n",
    "        \"Sim\":riv_sim.iloc[:, 0],\n",
    "    }\n",
    ")\n",
    "df.loc[\"1976-01-01\":, :].plot(\n",
    "    grid=True,\n",
    "    title=\"Simulated and observed river flows\",\n",
    "    ylabel=\"m3/s\",\n",
    "    style=[\"+\", \"-\", \"-\"]\n",
    ")\n",
    "plt.show() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the watertable simulation from optimization 1 \n",
    "wtl_sim = sim.get_output(\"watertable\")\n",
    "\n",
    "# Plot watertable\n",
    "df = pd.DataFrame(\n",
    "    {\n",
    "        \"Obs\":wtl_obs.iloc[:, 0],\n",
    "        \"Sim\":wtl_sim.iloc[:, 0],\n",
    "    }\n",
    ")\n",
    "df.loc[\"1976-01-01\":, :].plot(\n",
    "    grid=True,\n",
    "    title=\"Simulated and observed groundwater levels\",\n",
    "    ylabel=\"m NGF\",\n",
    "    style=[\"+\", \"-\", \"-\"]\n",
    ")\n",
    "plt.show() "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "dev",
   "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}