{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "# Test tidal response with analytic solutions\n",
    "\n",
    "This notebook shows some validation of tidal wave propagation in a 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",
    "\n",
    "import timflow.transient as tft"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we define some solutions from the book of Bruggeman. \n",
    "\n",
    "_(Note 1: Note how easily `numpy` deals with the complex algebra (splitting into an imaginary and a real part)._\n",
    "\n",
    "_(Note 2: We might move the Bruggeman code to a new repository: https://github.com/dbrakenhoff/bruggeman)._"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bruggeman_128_01(x, t, h, S, k, D, tau):\n",
    "    \"\"\"Tidal fluctuation open water, confined aquifer with open boundary (x = 0).\n",
    "\n",
    "    From Bruggeman 128.01\n",
    "\n",
    "    h = amplitude of tidal fluctuation, [m]\n",
    "    k = hydraulic conductivity [m/d]\n",
    "    D = aquifer thickness [m]\n",
    "    S = storage coefficient [-]\n",
    "    tau = tidal period [d]\n",
    "    \"\"\"\n",
    "    beta = np.sqrt(S / (k * D))\n",
    "    omega = 2 * np.pi / tau\n",
    "    omega_accent = beta * np.sqrt(omega / 2)\n",
    "\n",
    "    return h * np.exp(-omega_accent * x) * np.sin(omega * t - omega_accent * x)\n",
    "\n",
    "\n",
    "def bruggeman_128_03(x, t, h, S, k, D, tau, c):\n",
    "    \"\"\"Tidal fluctuation open water, leaky aquifer with open boundary (x = 0).\n",
    "\n",
    "    From Bruggeman 128.03\n",
    "\n",
    "    h = amplitude of tidal fluctuation, [m]\n",
    "    k = hydraulic conductivity [m/d]\n",
    "    D = aquifer thickness [m]\n",
    "    S = storage coefficient [-]\n",
    "    tau = tidal period [d]\n",
    "    c = leakance [d]\n",
    "    \"\"\"\n",
    "    beta = np.sqrt(S / (k * D))\n",
    "    eta = 1 / (c * S)\n",
    "    omega = 2 * np.pi / tau\n",
    "\n",
    "    a = np.real(np.sqrt(eta + 1j * omega))\n",
    "    b = np.imag(np.sqrt(eta + 1j * omega))\n",
    "\n",
    "    return h * np.exp(-beta * a * x) * np.sin(omega * t - beta * b * x)\n",
    "\n",
    "\n",
    "def bruggeman_128_04(x, t, h, S, k, D, tau, c, w):\n",
    "    \"\"\"Tidal fluctuation open water, leaky aquifer with entrance resistance (x = 0).\n",
    "\n",
    "    From Bruggeman 128.04\n",
    "\n",
    "    h = amplitude of tidal fluctuation, [m]\n",
    "    k = hydraulic conductivity [m/d]\n",
    "    D = aquifer thickness [m]\n",
    "    S = storage coefficient [-]\n",
    "    tau = tidal period [d]\n",
    "    c = leakance [d]\n",
    "    w = entry resistance at x=0 [d]\n",
    "    \"\"\"\n",
    "    beta = np.sqrt(S / (k * D))\n",
    "    eta = 1 / (c * S)\n",
    "    omega = 2 * np.pi / tau\n",
    "    theta = 1 / (np.power(beta, 2) * np.power(k, 2) * np.power(w, 2))\n",
    "\n",
    "    a = np.real(np.sqrt(eta + 1j * omega))\n",
    "    b = np.imag(np.sqrt(eta + 1j * omega))\n",
    "\n",
    "    return (\n",
    "        h\n",
    "        * np.sqrt(theta)\n",
    "        * np.exp(-beta * a * x)\n",
    "        * np.sin(omega * t - beta * b * x - np.arctan(b / (a + np.sqrt(theta))))\n",
    "        / (np.sqrt(np.power((a + np.sqrt(theta)), 2) + np.power(b, 2)))\n",
    "    )"
   ]
  },
  {
   "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://ttim.readthedocs.io/en/latest/00tutorials/howto_fluctuating_head_boundary.html). 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",
    "    phreatictop=False,\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",
    "    phreatictop=False,\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 = 1000  # 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",
    "    phreatictop=False,\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",
    "    phreatictop=False,\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\n",
    "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",
    "    phreatictop=False,\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",
    "    phreatictop=False,\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)."
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}