{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "# 128.0x: Tidal fluctuation of open water."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "This notebook shows comparisons of tidal wave propagation in a 'Timflow' ModelXsection against known formulas from literature (Bruggeman, 1999: `Analytical solutions of geohydrological problems`).\n",
    "\n",
    "Three situations are shown:\n",
    "1. Bruggeman 128.01: tidal fluctation of open water, in a confined aquifer with open boundary at x=0\n",
    "2. Bruggeman 128.03: tidal fluctation of open water, in a leaky aquifer with open boundary at x=0\n",
    "3. Bruggeman 128.04: tidal fluctation of open water, in a leaky aquifer with entry resistance at x=0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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_128_01, bruggeman_128_03, bruggeman_128_04"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "View the functions. Since the latter equation is the most generic one (with the others being more specific), we will show that one:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bruggeman_128_04"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "View the docstring to get a description of the input parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(bruggeman_128_04)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Set up the comparison in `timflow`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The analytic solution is for a river with varying head that has been varying for ever. \n",
    "In `timflow`, head is simulated for several days to get the model to spin-up. \n",
    "\n",
    "By using a `sine` function, the head starts at zero (as should be the case in Timflow)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def head_river(h, t, tau, tp=0):\n",
    "    return h * np.sin(2 * np.pi * (t - tp) / tau)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below the `sine` (water level on river) is split into steps as shown in the example [How to model a fluctuating head boundary](https://timflow.readthedocs.io/en/latest/transient/00userguide/howtos/howto_fluctuating_head_boundary.html) from the `timflow` documentation. We use a small timestep ($\\Delta t$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters of tidal wave\n",
    "tau = 0.5  # tidal period, d\n",
    "h_tidal = 1  # amplitude of tidal fluctuation, m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmax = 5  # day\n",
    "delt = 0.01  # day\n",
    "t = np.arange(0, tmax, delt)\n",
    "hexact = head_river(h_tidal, t, tau)\n",
    "hmid = 0.5 * (\n",
    "    head_river(h_tidal, t - delt / 2, tau) + head_river(h_tidal, t + delt / 2, tau)\n",
    ")\n",
    "tmid = np.hstack((0, 0.5 * (t[:-1] + t[1:])))\n",
    "# plot\n",
    "plt.figure(figsize=(10, 3))\n",
    "plt.plot(t, hexact, \"k\", label=\"observed\")\n",
    "plt.step(tmid, hmid, where=\"post\", label=\"mid\")\n",
    "plt.xlim(0, 1)\n",
    "plt.xlabel(\"time (d)\")\n",
    "plt.ylabel(\"head (m)\")\n",
    "plt.legend()\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Bruggeman 128.01 - Confined, no entry resistance\n",
    "\n",
    "Create the Timflow model for comparison with Bruggeman 128.01 (confined, no entry resistance).\n",
    "\n",
    "In this example, we use the 1D inhomogeneities because of the nice plotting. The same result could be reached with a 'normal' Timflow model.\n",
    "In all the Timflow models, the left-hand inhomogeneity (`river`) is shown for illustration purposes. That inhomogeneity is only created to satisfy the requirement of ModelXsection to have inhoms for  $-\\infty<x<\\infty$.\n",
    "\n",
    "The hydraulic boundary `river_hls` is modeled as a River1D at x=0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# parameters\n",
    "k = 20.0  # hydraulic conductivity, m/d\n",
    "H = 50.0  # thickness of aquifer, m\n",
    "T = k * H  # transmissivity, m^2/d\n",
    "S = 0.001  # storage coefficient of aquifer, [-]\n",
    "\n",
    "Saq = S / H  # specific storage [1/m]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tft.ModelXsection(naq=1, tmin=1e-4, tmax=1e2)\n",
    "\n",
    "riv = tft.XsectionMaq(\n",
    "    model=ml,\n",
    "    x1=-np.inf,\n",
    "    x2=0,\n",
    "    kaq=k,\n",
    "    z=[0, -H],\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=[0, -H],\n",
    "    Saq=Saq,\n",
    "    topboundary=\"confined\",\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, xls=0 - small, tsandh=list(zip(tmid, hmid, strict=True))\n",
    ")\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",
    "ax.set_xlim(-100, 100)\n",
    "ax.set_ylim(ymax=12);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compare the head at several locations.\n",
    "\n",
    "As can be seen from the lowest graph $x = 500$, the `timflow` model needs some spinup time, but after a few days the graphs are similar."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xlocs = [1, 50, 500, 1000]\n",
    "tp = t\n",
    "\n",
    "f, axes = plt.subplots(\n",
    "    len(xlocs), 1, figsize=(10, 3 * len(xlocs)), sharex=True, sharey=True\n",
    ")\n",
    "\n",
    "for i in range(len(xlocs)):\n",
    "    h = ml.head(xlocs[i], 0, tp)\n",
    "    axes[i].plot(tp, h[0], label=\"Timflow - confined\", linestyle=\"--\")\n",
    "\n",
    "    h_128_01 = bruggeman_128_01(xlocs[i], tp, h_tidal, S, k, H, tau)\n",
    "    axes[i].plot(tp, h_128_01, \"r\", label=\"Bruggeman 128.01\", linestyle=\":\")\n",
    "\n",
    "    axes[i].set_title(f\"x = {xlocs[i]} m\", loc=\"right\")\n",
    "    axes[i].legend(loc=(0, 1), frameon=False, ncol=2)\n",
    "    axes[i].grid()\n",
    "    axes[i].set_xlabel(\"time (d)\")\n",
    "    axes[i].set_ylabel(\"head (m)\")\n",
    "    axes[i].set_xlim([0, 3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Bruggeman 128.03 Leaky, no entry resistance\n",
    "\n",
    "Create a Timflow model for comparison with Bruggeman 128.03 (leaky, no entry resistance)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Additional parameters for this model:\n",
    "c = 50  # leakage factor, d\n",
    "\n",
    "# The model\n",
    "ml2 = tft.ModelXsection(naq=1, tmin=1e-4, tmax=1e2)\n",
    "\n",
    "riv2 = tft.XsectionMaq(\n",
    "    model=ml2,\n",
    "    x1=-np.inf,\n",
    "    x2=0,\n",
    "    kaq=k,\n",
    "    z=[0, -H],\n",
    "    Saq=1e-20,\n",
    "    Sll=1e-20,\n",
    "    topboundary=\"confined\",\n",
    "    name=\"river\",\n",
    ")\n",
    "\n",
    "land2 = tft.XsectionMaq(\n",
    "    model=ml2,\n",
    "    x1=0,\n",
    "    x2=np.inf,\n",
    "    kaq=k,\n",
    "    z=[5, 0, -H],\n",
    "    c=c,\n",
    "    Saq=Saq,\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_hls2 = tft.River1D(\n",
    "    model=ml2, xls=0 - small, tsandh=list(zip(tmid, hmid, strict=True))\n",
    ")\n",
    "\n",
    "\n",
    "ml2.solve()\n",
    "\n",
    "ax = riv2.plot(params=True, names=True, labels=False)\n",
    "land2.plot(ax=ax, params=True, names=True, labels=False)\n",
    "river_hls2.plot(ax=ax)\n",
    "ax.set_xlim(-100, 100)\n",
    "ax.set_ylim(ymax=12);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, compare heads at several locations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xlocs = [1, 50, 500, 1000]\n",
    "tp = t\n",
    "\n",
    "f, axes = plt.subplots(\n",
    "    len(xlocs), 1, figsize=(10, 3 * len(xlocs)), sharex=True, sharey=True\n",
    ")\n",
    "\n",
    "for i in range(len(xlocs)):\n",
    "    h = ml2.head(xlocs[i], 0, tp)\n",
    "    axes[i].plot(tp, h[0], label=\"Timflow - semi-confined\", linestyle=\"--\")\n",
    "\n",
    "    h_128_03 = bruggeman_128_03(xlocs[i], tp, h_tidal, S, k, H, tau, c)\n",
    "    axes[i].plot(tp, h_128_03, \"g\", label=\"Bruggeman 128.03\", linestyle=\":\")\n",
    "\n",
    "    axes[i].set_title(f\"x = {xlocs[i]} m\", loc=\"right\")\n",
    "    axes[i].legend(loc=(0, 1), frameon=False, ncol=2)\n",
    "    axes[i].grid()\n",
    "    axes[i].set_xlabel(\"time (d)\")\n",
    "    axes[i].set_ylabel(\"head (m)\")\n",
    "    axes[i].set_xlim([0, 3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Bruggeman 128.04 Leaky, with entry resistance\n",
    "\n",
    "Create a Timflow model for comparison with Bruggeman 128.04 (leaky, with entry resistance at x=0). \n",
    "\n",
    "The entry resistance is added to the River1D `river_hls` by adding the keyword argument `res=w`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Additional parameters for this model:\n",
    "w = 20  # entry resistance at x=0, d\n",
    "\n",
    "# The model\n",
    "ml3 = tft.ModelXsection(naq=1, tmin=1e-4, tmax=1e2)\n",
    "\n",
    "riv3 = tft.XsectionMaq(\n",
    "    model=ml3,\n",
    "    x1=-np.inf,\n",
    "    x2=0,\n",
    "    kaq=k,\n",
    "    z=[0, -H],\n",
    "    Saq=1e-20,\n",
    "    Sll=1e-20,\n",
    "    topboundary=\"confined\",\n",
    "    name=\"river\",\n",
    ")\n",
    "\n",
    "land3 = tft.XsectionMaq(\n",
    "    model=ml3,\n",
    "    x1=0,\n",
    "    x2=np.inf,\n",
    "    kaq=k,\n",
    "    z=[5, 0, -H],\n",
    "    c=c,\n",
    "    Saq=Saq,\n",
    "    topboundary=\"semi\",\n",
    "    name=\"hinterland\",\n",
    ")\n",
    "\n",
    "\n",
    "# Use a small offset to avoid a singular matrix.\n",
    "small = 1e-5\n",
    "\n",
    "river_hls3 = tft.River1D(\n",
    "    model=ml3,\n",
    "    xls=0 - small,\n",
    "    tsandh=list(zip(tmid, hmid, strict=True)),\n",
    "    res=w,\n",
    ")\n",
    "\n",
    "ml3.solve()\n",
    "\n",
    "ax = riv3.plot(params=True, names=True, labels=False)\n",
    "land3.plot(ax=ax, params=True, names=True, labels=False)\n",
    "river_hls3.plot(ax=ax)\n",
    "\n",
    "ax.set_xlim(-100, 100)\n",
    "ax.set_ylim(ymax=12);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xlocs = [1, 50, 500, 1000]\n",
    "tp = t\n",
    "\n",
    "f, axes = plt.subplots(\n",
    "    len(xlocs), 1, figsize=(10, 3 * len(xlocs)), sharex=True, sharey=True\n",
    ")\n",
    "\n",
    "for i in range(len(xlocs)):\n",
    "    h = ml3.head(xlocs[i], 0, tp)\n",
    "    axes[i].plot(tp, h[0], label=\"Timflow - leaky + entry resistance\", linestyle=\"--\")\n",
    "    h_128_04 = bruggeman_128_04(xlocs[i], tp, h_tidal, S, k, H, tau, c, w)\n",
    "    axes[i].plot(tp, h_128_04, \"k\", label=\"Bruggeman 128.04\", linestyle=\":\")\n",
    "    axes[i].set_title(f\"x = {xlocs[i]} m\", loc=\"right\")\n",
    "    axes[i].legend(loc=(0, 1), frameon=False, ncol=2)\n",
    "    axes[i].grid()\n",
    "    axes[i].set_xlabel(\"time (d)\")\n",
    "    axes[i].set_ylabel(\"head (m)\")\n",
    "    axes[i].set_xlim([0, 3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As can be seen, all three cases are replicated accurately after approximately 2 days spinup time.\n",
    "\n",
    "In the case with entry resistance, the head in the aquifer is already considerably damped compared to the open water fluctuation, even at a small distance (x=1 m)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Comparison of the three equations\n",
    "\n",
    "Below is a comparison shown between the three equations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(0, 1000, 101)\n",
    "y = np.zeros_like(x)\n",
    "tg = (10 + np.array([0.25, 0.5, 0.75, 1])) * tau  # time in days\n",
    "\n",
    "f, axes = plt.subplots(\n",
    "    len(xlocs), 1, figsize=(10, 3 * len(tg)), sharex=True, sharey=True\n",
    ")\n",
    "\n",
    "\n",
    "for i in range(len(tg)):\n",
    "    h_128_01 = bruggeman_128_01(x, tg[i], h_tidal, S, k, H, tau)\n",
    "    axes[i].plot(x, h_128_01, label=\"Bruggeman 128.01\", linestyle=\"-\")\n",
    "    h_128_03 = bruggeman_128_03(x, tg[i], h_tidal, S, k, H, tau, c)\n",
    "    axes[i].plot(x, h_128_03, label=\"Bruggeman 128.03\", linestyle=\"-\")\n",
    "    h_128_04 = bruggeman_128_04(x, tg[i], h_tidal, S, k, H, tau, c, w)\n",
    "    axes[i].plot(x, h_128_04, \"k\", label=\"Bruggeman 128.04\", linestyle=\"-\")\n",
    "\n",
    "    h_riv = head_river(h_tidal, tg[i], tau)\n",
    "    axes[i].plot([-50, 0], [h_riv, h_riv], \"b\")\n",
    "\n",
    "    axes[i].legend(loc=(0, 1), frameon=False, ncol=6, fontsize=\"small\")\n",
    "    axes[i].set_xlabel(\"x [m]\")\n",
    "    axes[i].set_ylabel(\"drawdown [m]\")\n",
    "    axes[i].grid()\n",
    "    axes[i].set_title(f\"t = {tg[i]:.3f} d\", loc=\"right\")\n",
    "    # axes[i].set_ylim(0, h_step * 1.2)\n",
    "    # axes[i].axvline(0, color=\"grey\", linewidth=3)\n",
    "\n",
    "    axes[i].set_xlim(-50, 1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "bruggeman (broken)",
   "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": 4
}