{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "da5c2724-57ac-4684-a1a3-ed7dc8a9ce11",
   "metadata": {},
   "source": [
    "## Benchmark loading efficiency\n",
    "### Comparison to results of Bakker (2016)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b68a70e-2986-456c-9beb-8a006e79726a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "import timflow.transient as tft"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c237a98-5b44-4910-b962-9387a464d128",
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 20  # m/d\n",
    "H = 10  # m\n",
    "Ss = 1e-4  # 1/m\n",
    "Sp = 0.1  # -\n",
    "c = 500  # d\n",
    "tsandhstar = [(0, 1)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35f29c29-917b-42ca-9092-6bd16f8af518",
   "metadata": {},
   "source": [
    "### Figure 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e080bceb-3ae4-47db-9f5e-f792ec451190",
   "metadata": {},
   "outputs": [],
   "source": [
    "def model(zland, leffaq=0, cland=1e12, Sland=Ss, land_top=\"conf\"):\n",
    "    ml = tft.ModelXsection(naq=1, tmin=1e-4, tmax=1e2)\n",
    "\n",
    "    tft.XsectionMaq(\n",
    "        model=ml,\n",
    "        x1=-np.inf,  # river extends to infinitiy\n",
    "        x2=0,\n",
    "        z=[1, 0, -H],\n",
    "        kaq=k,\n",
    "        Saq=Ss,\n",
    "        c=c,\n",
    "        leffaq=leffaq,\n",
    "        topboundary=\"semi\",\n",
    "        tsandhstar=tsandhstar,\n",
    "        name=\"river\",\n",
    "    )\n",
    "\n",
    "    tft.XsectionMaq(\n",
    "        model=ml,\n",
    "        x1=0,\n",
    "        x2=np.inf,  # hinterland extends to infinity\n",
    "        kaq=k,\n",
    "        z=zland,\n",
    "        Saq=Sland,\n",
    "        c=cland,\n",
    "        topboundary=land_top,\n",
    "        name=\"land\",\n",
    "    )\n",
    "\n",
    "    ml.solve()\n",
    "    return ml"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5332c9bd-8dd1-4f61-8b0d-ed13079f6a19",
   "metadata": {},
   "outputs": [],
   "source": [
    "ml_leff = []\n",
    "ml_noleff = []\n",
    "# land is confined\n",
    "ml_leff.append(model([0, -H], leffaq=1, Sland=Ss, land_top=\"conf\"))\n",
    "ml_noleff.append(model([0, -H], leffaq=0, Sland=Ss, land_top=\"conf\"))\n",
    "# land is semi-confined\n",
    "ml_leff.append(model([1, 0, -H], leffaq=1, cland=c, Sland=Ss, land_top=\"semi\"))\n",
    "ml_noleff.append(model([1, 0, -H], leffaq=0, cland=c, Sland=Ss, land_top=\"semi\"))\n",
    "# land is phreatic\n",
    "ml_leff.append(model([0, -H], leffaq=1, Sland=Sp, land_top=\"phre\"))\n",
    "ml_noleff.append(model([0, -H], leffaq=0, Sland=Sp, land_top=\"phre\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8053886-ec7c-4a6a-8009-f02bb3a82e1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(-300, 300, 100)\n",
    "t = [1 / 24, 1, 24]\n",
    "plt.figure(figsize=(6, 6))\n",
    "title = [\"confined\", \"semi-confined\", \"unconfined\"]\n",
    "for i in range(len(ml_leff)):\n",
    "    ha = ml_leff[i].headalongline(x, 0, t).squeeze()\n",
    "    hb = ml_noleff[i].headalongline(x, 0, t).squeeze()\n",
    "    ax = plt.subplot(3, 1, i + 1)\n",
    "    legendtext = [\"1 hr\", \"1 d\", \"24 d\"]\n",
    "    for it in range(len(t)):\n",
    "        plt.plot(x, ha[it], color=\"C\" + str(it), label=legendtext[it])\n",
    "        plt.plot(x, hb[it], \"--\", color=\"C\" + str(it))\n",
    "    plt.axvline(0, linestyle=\":\", color=\"k\")\n",
    "    plt.xlim(-300, 300)\n",
    "    plt.ylim(0, 1)\n",
    "    plt.legend(loc=\"lower left\", bbox_to_anchor=(1, 0.2), frameon=False)\n",
    "    plt.ylabel(\"head (m)\")\n",
    "    plt.text(0.6, 0.7, title[i], transform=ax.transAxes)\n",
    "    plt.text(1.05, 0.9, \"loading eff = 1 (solid)\", transform=ax.transAxes)\n",
    "    plt.text(1.05, 0.75, \"no loading eff (dashed)\", transform=ax.transAxes)\n",
    "    plt.grid()\n",
    "plt.xlabel(\"x (m)\")\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e5165da5-a2ba-47d4-8b6d-e36adc515b5b",
   "metadata": {},
   "source": [
    "### Figure 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52e8542e-6f9d-4ca9-8afb-99bf46d4c161",
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.logspace(-2, 2, 100)\n",
    "x = [0, 100, 200]\n",
    "plt.figure(figsize=(5.5, 6))\n",
    "for i in range(len(x)):\n",
    "    ax = plt.subplot(3, 1, i + 1)\n",
    "    legendtext = [\"confined\", \"semi-confined\", \"phreatic\"]\n",
    "    for im in range(len(ml_leff)):\n",
    "        ha = ml_leff[im].head(x[i], 0, t).squeeze()\n",
    "        hb = ml_noleff[im].head(x[i], 0, t).squeeze()\n",
    "        plt.semilogx(t, ha, color=\"C\" + str(im), label=legendtext[im])\n",
    "        plt.semilogx(t, hb, \"--\", color=\"C\" + str(im))\n",
    "    plt.axvline(0, linestyle=\":\")\n",
    "    plt.xlim(0.01, 100)\n",
    "    plt.ylim(0, 1)\n",
    "    plt.legend(loc=\"lower left\", bbox_to_anchor=(1, 0.2), frameon=False)\n",
    "    plt.ylabel(\"head (m)\")\n",
    "    plt.text(0.1, 0.85, f\"at x={x[i]} m\", transform=ax.transAxes)\n",
    "    plt.text(1.05, 0.9, \"loading eff = 1 (solid)\", transform=ax.transAxes)\n",
    "    plt.text(1.05, 0.75, \"no loading eff (dashed)\", transform=ax.transAxes)\n",
    "    plt.grid()\n",
    "plt.xlabel(\"time (d)\")\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c9c7dfb-3aeb-42e3-9479-bfcfdba6bd82",
   "metadata": {},
   "source": [
    "### Figure 6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9ca80d9-4a8c-4265-ae34-7f55d19e7955",
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 20  # m/d\n",
    "H = 10  # m\n",
    "Ss = 1e-4  # 1/m\n",
    "Sp = 0.1  # -\n",
    "c = 100  # d\n",
    "tsandhstar = [(0, 1)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "015f738a-8101-44ce-966f-571f892f8a64",
   "metadata": {},
   "outputs": [],
   "source": [
    "def model(zland, leffaq=0, cland=1e12, Sland=Ss, land_top=\"conf\", W=20):\n",
    "    ml = tft.ModelXsection(naq=1, tmin=1e-4, tmax=1e2)\n",
    "\n",
    "    tft.XsectionMaq(\n",
    "        model=ml,\n",
    "        x1=-W / 2,  # river extends to infinitiy\n",
    "        x2=W / 2,\n",
    "        z=[1, 0, -H],\n",
    "        kaq=k,\n",
    "        Saq=Ss,\n",
    "        c=c,\n",
    "        leffaq=leffaq,\n",
    "        topboundary=\"semi\",\n",
    "        tsandhstar=tsandhstar,\n",
    "        name=\"river\",\n",
    "    )\n",
    "\n",
    "    tft.XsectionMaq(\n",
    "        model=ml,\n",
    "        x1=-np.inf,\n",
    "        x2=-W / 2,  # hinterland extends to infinity\n",
    "        kaq=k,\n",
    "        z=zland,\n",
    "        Saq=Sland,\n",
    "        c=cland,\n",
    "        topboundary=land_top,\n",
    "        name=\"left\",\n",
    "    )\n",
    "\n",
    "    tft.XsectionMaq(\n",
    "        model=ml,\n",
    "        x1=W / 2,\n",
    "        x2=np.inf,  # hinterland extends to infinity\n",
    "        kaq=k,\n",
    "        z=zland,\n",
    "        Saq=Sland,\n",
    "        c=cland,\n",
    "        topboundary=land_top,\n",
    "        name=\"right\",\n",
    "    )\n",
    "\n",
    "    ml.solve()\n",
    "    return ml"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3cbfa57-39b6-4b1f-bf23-b3318cb14a3d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_models(W=20):\n",
    "    ml_leff = []\n",
    "    ml_noleff = []\n",
    "    # land is confined\n",
    "    ml_leff.append(model([0, -H], leffaq=1, Sland=Ss, land_top=\"conf\", W=W))\n",
    "    ml_noleff.append(model([0, -H], leffaq=0, Sland=Ss, land_top=\"conf\", W=W))\n",
    "    # land is semi-confined\n",
    "    ml_leff.append(\n",
    "        model(\n",
    "            [1, 0, -H],\n",
    "            leffaq=1,\n",
    "            cland=c,\n",
    "            Sland=Ss,\n",
    "            land_top=\"semi\",\n",
    "            W=W,\n",
    "        )\n",
    "    )\n",
    "    ml_noleff.append(\n",
    "        model(\n",
    "            [1, 0, -H],\n",
    "            leffaq=0,\n",
    "            cland=c,\n",
    "            Sland=Ss,\n",
    "            land_top=\"semi\",\n",
    "            W=W,\n",
    "        )\n",
    "    )\n",
    "    # land is phreatic\n",
    "    ml_leff.append(model([0, -H], leffaq=1, Sland=Sp, land_top=\"phreatic\", W=W))\n",
    "    ml_noleff.append(model([0, -H], leffaq=0, Sland=Sp, land_top=\"phreatic\", W=W))\n",
    "    return ml_leff, ml_noleff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92e362e6-10d6-41aa-aa0e-7eb7e948bad2",
   "metadata": {},
   "outputs": [],
   "source": [
    "W = [100, 50, 20]\n",
    "t = np.logspace(-2, 2, 100)\n",
    "plt.figure(figsize=(5.5, 6))\n",
    "for i in range(len(W)):\n",
    "    ml_leff, ml_noleff = create_models(W[i])\n",
    "    ax = plt.subplot(3, 1, i + 1)\n",
    "    legendtext = [\"confined\", \"semi-confined\", \"phreatic\"]\n",
    "    for im in range(len(ml_leff)):\n",
    "        ha = ml_leff[im].head(W[i] / 2, 0, t).squeeze()\n",
    "        hb = ml_noleff[im].head(W[i] / 2, 0, t).squeeze()\n",
    "        plt.semilogx(t, ha, color=\"C\" + str(im), label=legendtext[im])\n",
    "        plt.semilogx(t, hb, \"--\", color=\"C\" + str(im))\n",
    "    plt.axvline(0, linestyle=\":\")\n",
    "    plt.xlim(0.01, 100)\n",
    "    plt.ylim(0, 1)\n",
    "    plt.legend(loc=\"lower left\", bbox_to_anchor=(1, 0.2), frameon=False)\n",
    "    plt.ylabel(\"head (m)\")\n",
    "    plt.text(0.1, 0.85, f\"W={W[i]} m\", transform=ax.transAxes)\n",
    "    plt.text(1.05, 0.9, \"loading eff = 1 (solid)\", transform=ax.transAxes)\n",
    "    plt.text(1.05, 0.75, \"no loading eff (dashed)\", transform=ax.transAxes)\n",
    "    plt.grid()\n",
    "plt.xlabel(\"time (d)\")\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a5a1cee-1736-428d-b6ca-70a2debb5b79",
   "metadata": {},
   "source": [
    "### Reference\n",
    "Bakker, M., 2016. The effect of loading efficiency on the groundwater response to water level changes in shallow lakes and streams. Water Resources Research, 52(3), pp.1705-1715, [https://doi.org/10.1002/2015WR017977](https://doi.org/10.1002/2015WR017977)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}