{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "fd242557",
   "metadata": {},
   "source": [
    "# 126.33: Steady state after step head change."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "import-packages",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import timflow.transient as tft\n",
    "\n",
    "from bruggeman.flow1d import bruggeman_126_33"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cfa264a9",
   "metadata": {},
   "source": [
    "View the function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74014f25",
   "metadata": {},
   "outputs": [],
   "source": [
    "bruggeman_126_33"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5df8bc7d",
   "metadata": {},
   "source": [
    "View the docstring to get a description of the input parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d2d0ae3",
   "metadata": {},
   "outputs": [],
   "source": [
    "help(bruggeman_126_33)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f53cf151",
   "metadata": {},
   "source": [
    "Define some aquifer parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "parameters",
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 20.0  # hydraulic conductivity, m/d\n",
    "D = 50.0  # thickness of aquifer, m\n",
    "T = k * D  # transmissivity, m^2/d\n",
    "c = 1000  # leakage factor, d\n",
    "w = 20  # entry resistance at x=0, d\n",
    "h_step = 1.0  # head at x=0, m\n",
    "S = 0.001  # storage coefficient of aquifer, [-]\n",
    "\n",
    "Saq = S / D  # specific storage [1/m]\n",
    "Sll = 0.001  # specific storage of leaky layer, [-]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91a403f1",
   "metadata": {},
   "source": [
    "Set up a `timflow` model. \n",
    "\n",
    "The ModelXsection to be able to plot the sections, but the same result could be reached with the 'normal' ModelMaq model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "model-setup",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The model\n",
    "ml = tft.ModelXsection(naq=1, tmin=1e-4, tmax=1e3)\n",
    "\n",
    "riv = tft.XsectionMaq(\n",
    "    model=ml,\n",
    "    x1=-np.inf,\n",
    "    x2=0,\n",
    "    kaq=k,\n",
    "    z=[0, -D],\n",
    "    Saq=1e-20,\n",
    "    Sll=1e-20,\n",
    "    topboundary=\"confined\",\n",
    "    name=\"river\",\n",
    ")\n",
    "\n",
    "land = tft.XsectionMaq(\n",
    "    model=ml,\n",
    "    x1=0,\n",
    "    x2=np.inf,\n",
    "    kaq=k,\n",
    "    z=[5, 0, -D],\n",
    "    c=c,\n",
    "    Saq=Saq,\n",
    "    Sll=Sll,\n",
    "    topboundary=\"semi\",\n",
    "    name=\"hinterland\",\n",
    ")\n",
    "\n",
    "# Use a small offset to avoid a singular matrix.\n",
    "small = 1e-5\n",
    "\n",
    "river_hls = tft.River1D(\n",
    "    model=ml,\n",
    "    xls=0 - small,\n",
    "    tsandh=[0, h_step],\n",
    "    res=w,\n",
    ")\n",
    "\n",
    "ml.solve()\n",
    "\n",
    "ax = riv.plot(params=True, names=True, labels=False)\n",
    "land.plot(ax=ax, params=True, names=True, labels=False)\n",
    "river_hls.plot(ax=ax)\n",
    "\n",
    "ax.set_xlim(-100, 100)\n",
    "ax.set_ylim(ymax=12)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10f3a609",
   "metadata": {},
   "source": [
    "Compare `timflow` implementation to the analytical solution from Bruggeman."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ae8b0cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(0, 1000, 101)\n",
    "y = np.zeros_like(x)\n",
    "t = np.logspace(-1, 1, 3)\n",
    "\n",
    "plt.figure(figsize=(10, 3))\n",
    "for i in range(len(t)):\n",
    "    h = ml.headalongline(x, y, t[i])\n",
    "    plt.plot(x, h.squeeze(), label=f\"t={t[i]:.2f} d\")\n",
    "    ha = bruggeman_126_33(x, h_step, k, D, c, w)\n",
    "    plt.plot(x, ha, \"k:\")\n",
    "plt.plot([], [], \"k:\", label=\"Bruggeman 126.33\")\n",
    "plt.legend(loc=(0, 1), frameon=False, ncol=6, fontsize=\"small\")\n",
    "plt.xlabel(\"x [m]\")\n",
    "plt.ylabel(\"drawdown [m]\")\n",
    "plt.grid()\n",
    "plt.ylim(0, h_step * 1.2)\n",
    "plt.axvline(0, color=\"grey\", linewidth=3)\n",
    "plt.plot([-50, 0], [h_step, h_step], \"b\")\n",
    "plt.xlim(-50, 1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3414d28",
   "metadata": {},
   "source": [
    "As can be seen from the graph, the `storativity` causes a delay in reaching the steady-state value. When the storativity approaches zero, the steady-state is immediately reached.\n",
    "\n",
    "Furthermore, the head in the aquifer never reaches the head in the river, because of the entry resistance."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "bruggeman (3.13.11)",
   "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.13.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}