{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "edd8871e-c5c8-41c0-b21a-8cebc129b30e",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "# HowTo: Make a cross-section model with `timflow`\n",
    " \n",
    "This HowTo describes how to make a transient cross-section model with `timflow`.\n",
    "\n",
    "We start by importing the required packages:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64c3b56a-ec12-4da3-ad88-64ae7fd4193d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "import timflow.transient as tft\n",
    "\n",
    "#\n",
    "plt.rcParams[\"figure.figsize\"] = (4, 3)  # set default figure size\n",
    "plt.rcParams[\"font.size\"] = 8  # set default font size"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8abce62-4971-441d-b45f-5a3dfb9486fa",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "We will create an xsection model for a two-aquifer system with a river with a leaky river bottom (see Figure).\n",
    "The model is divided in three sections: the left side of the river, the river, and the right side of the river.\n",
    "The resistance of the top semi-confining layer is different than the resistance of the leaky river bottom. \n",
    "The aquifer parameters are defined in the table and code cell below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0278d562-ffdf-4495-a842-1eb9f9815329",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.colors as mcolors\n",
    "\n",
    "\n",
    "def gradient_fill(x, y, fill_color=\"C0\", ax=None, **kwargs):\n",
    "    if ax is None:\n",
    "        ax = plt.gca()\n",
    "\n",
    "    z = np.empty((100, 1, 4), dtype=float)\n",
    "    rgb = mcolors.colorConverter.to_rgb(fill_color)\n",
    "    z[:, :, :3] = rgb\n",
    "    z[:, :, -1] = np.linspace(0.1, 0.5, 100)[:, None]\n",
    "\n",
    "    xmin, xmax, ymin, ymax = x[0], x[1], y[0], y[1]\n",
    "    ax.imshow(\n",
    "        z,\n",
    "        aspect=\"auto\",\n",
    "        extent=[xmin, xmax, ymin, ymax],\n",
    "        origin=\"lower\",\n",
    "    )\n",
    "\n",
    "\n",
    "def arrow(xystart, xyend, text=\"\", arrow=\"<-\", color=\"k\", **kwds):\n",
    "    plt.annotate(\n",
    "        text,\n",
    "        xy=xystart,\n",
    "        xytext=xyend,\n",
    "        arrowprops={\"arrowstyle\": arrow, \"shrinkA\": 0, \"shrinkB\": 0, \"color\": color},\n",
    "        color=color,\n",
    "        **kwds,\n",
    "    )\n",
    "\n",
    "\n",
    "def solution12():\n",
    "    plt.figure(figsize=(8, 3))\n",
    "\n",
    "    plt.fill(\n",
    "        [-0.1, 0.3, 0.3, -0.1],\n",
    "        [0.8, 0.8, 1.2, 1.2],\n",
    "        color=\"grey\",\n",
    "        fill=False,\n",
    "        hatch=\"//\",\n",
    "    )\n",
    "    plt.fill(\n",
    "        [0.3, 0.7, 0.7, 0.3], [0.8, 0.8, 1.0, 1.0], color=\"grey\", fill=False, hatch=\"//\"\n",
    "    )\n",
    "    plt.fill(\n",
    "        [0.7, 1.1, 1.1, 0.7], [0.8, 0.8, 1.2, 1.2], color=\"grey\", fill=False, hatch=\"//\"\n",
    "    )\n",
    "\n",
    "    plt.fill(\n",
    "        [-0.1, 1.1, 1.1, -0.1],\n",
    "        [0, 0, -0.2, -0.2],\n",
    "        color=\"grey\",\n",
    "        fill=False,\n",
    "        hatch=\"///\",\n",
    "    )\n",
    "    plt.fill(\n",
    "        [-0.1, 1.1, 1.1, -0.1],\n",
    "        [-1, -1, -0.8, -0.8],\n",
    "        color=\"grey\",\n",
    "        fill=False,\n",
    "        hatch=\"xxx\",\n",
    "    )\n",
    "    plt.plot([-0.1, 0.3], [1.1, 1.1], \"C0\")\n",
    "    plt.plot([0.3, 0.7], [1.4, 1.4], \"C0\")\n",
    "    plt.plot([0.7, 1.1], [1.1, 1.1], \"C0\")\n",
    "    gradient_fill([0.3, 0.7], [1, 1.4])\n",
    "    plt.plot([0.3, 0.3], [0.8, 1.5], \"k\")\n",
    "    plt.plot([0.7, 0.7], [0.8, 1.5], \"k\")\n",
    "    plt.text(0.1, 1.35, \"$h_0$\", ha=\"center\")\n",
    "    plt.text(0.5, 1.55, \"$h_r(t)$\", ha=\"center\")\n",
    "    plt.text(0.9, 1.35, \"$h_0$\", ha=\"center\")\n",
    "    plt.plot([0.5, 0.5], [-1.3, -1.1], \"k\")\n",
    "    arrow((0.5, -1.2), (0.7, -1.2), \"$x$\", va=\"center\")\n",
    "    plt.text(-0.1, 0.4, \"left\", va=\"center\")\n",
    "    plt.text(0.5, 0.25, \"river\", va=\"center\", ha=\"center\")\n",
    "    plt.text(1.1, 0.4, \"right\", ha=\"right\", va=\"center\")\n",
    "    arrow((0.3, 0.45), (0.7, 0.45), arrow=\"<->\")\n",
    "    plt.text(0.5, 0.55, \"$2L$\", ha=\"center\")\n",
    "    plt.ylim(-1.3, 1.5)\n",
    "    plt.axis(\"off\")\n",
    "\n",
    "\n",
    "solution12()  # plt.savefig('riverxsection.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9c645d3-a790-42ea-be9b-509ac94aec48",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "| Parameter            | left/right | river |\n",
    "| :------------: | :------: | ----: |\n",
    "| $c_0$ (d)      |   500    | 200 |\n",
    "| $k_0$ (m/d)    |   10     | 10 |\n",
    "| $S_{s0}$ (m$^{-1}$)    |   1e-4     | 1e-4 |\n",
    "| $c_1$ (d)      |  1000    | 1000 |\n",
    "| $k_1$ (m/d)    |   20     | 20|\n",
    "| $S_{s1}$ (m$^{-1}$)    |   1e-4     | 1e-4 |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "675937ac-6786-4ffb-9c7d-bba705009055",
   "metadata": {},
   "outputs": [],
   "source": [
    "zside = [2, 0, -10, -12, -22]  # elevation of top and bottoms of layers, m\n",
    "zriver = [1, 0, -10, -12, -22]  # elevation of top and bottoms of layers, m\n",
    "cside0 = 500  # resistance of top semi-confining layer, d\n",
    "criver0 = 200  # resistance of leaky river bottom, d\n",
    "k0 = 10  # hydraulic conductivity of first aquifer, m/d\n",
    "Ss0 = 1e-4  # specific storage of first aquifer, m^(-1)\n",
    "c1 = 1000  # resistance of leaky layer between aquifers, d\n",
    "k1 = 20  # hydraulic conductivity of second aquifer, m/d\n",
    "Ss1 = 1e-4  # specific storage of second aquifer, m^(-1)\n",
    "L = 50  # half width of river, m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86d97286-0870-4cea-9b47-fae0be9284ab",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "An cross-sectional model is created with two aquifer layer, \n",
    "and is stored in the variable `ml`. Next, `XsectionMaq` sections are defined for the three sections of the model. Make sure that the entire cross-section is covered by an `XsectionMaq`. The head above the semi-confining top is fixed on the left and right side, which means that the head change equals zero. The head in the river varies with time and is specified with the `tsandhstar` keyword. `tsandhstar` must be a list of 2-tuples, where each tuple consists of `(time, head in river starting at that time)`.  Alternatively, `tsandhstar` can be an array with $N$ rows and two columns, where $N$ is the number of times that a different water level in the river is specified. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88304f3d-5614-474f-a7b5-98cf6600393b",
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tft.ModelXsection(naq=2, tmin=1e-3, tmax=1e2)\n",
    "\n",
    "tft.XsectionMaq(\n",
    "    model=ml,\n",
    "    x1=-np.inf,\n",
    "    x2=-L,\n",
    "    kaq=[k0, k1],\n",
    "    z=zside,\n",
    "    Saq=[Ss0, Ss1],\n",
    "    c=[cside0, c1],\n",
    "    topboundary=\"semi\",\n",
    "    name=\"left\",\n",
    ")\n",
    "\n",
    "tft.XsectionMaq(\n",
    "    model=ml,\n",
    "    x1=-L,\n",
    "    x2=L,\n",
    "    kaq=[k0, k1],\n",
    "    z=zriver,\n",
    "    Saq=[Ss0, Ss1],\n",
    "    c=[criver0, c1],\n",
    "    topboundary=\"semi\",\n",
    "    tsandhstar=[(0, 1)],\n",
    "    name=\"river\",\n",
    ")\n",
    "\n",
    "tft.XsectionMaq(\n",
    "    model=ml,\n",
    "    x1=L,\n",
    "    x2=np.inf,\n",
    "    kaq=[k0, k1],\n",
    "    z=zside,\n",
    "    Saq=[Ss0, Ss1],\n",
    "    c=[cside0, c1],\n",
    "    topboundary=\"semi\",\n",
    "    name=\"right\",\n",
    ")\n",
    "\n",
    "ml.solve()  # solve the model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4290951-4131-4f52-8f80-32fa703c3707",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Compute the head at 100 points in space between $-200$ and $+200$ for 5 times"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bca9d78a-d7cb-45e2-8c1f-782f16c84c41",
   "metadata": {},
   "outputs": [],
   "source": [
    "ng = 100\n",
    "nt = 5\n",
    "x = np.linspace(-200, 200, ng)\n",
    "y = np.zeros(ng)\n",
    "t = np.logspace(-2, 1, nt)\n",
    "h = ml.headalongline(x, y, t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "062663d5-7b27-4a2e-96f6-81b5e78673f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8, 4))\n",
    "for i, itime in enumerate(range(5)):\n",
    "    plt.plot(x, h[0, itime], \"C\" + str(i), label=f\"{t[itime]:.3f}\")\n",
    "    plt.plot(x, h[1, itime], \"C\" + str(i), ls=\"--\")\n",
    "plt.legend(title=\"time (d)\")\n",
    "plt.axhline(0, color=\"k\", lw=1)\n",
    "plt.axvspan(-L, L, color=[0.9, 0.9, 0.9]);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ce042e0-e249-4ff3-aaff-15081642f163",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Compute and plot the head at $x=50$ for time varying from 0.001 to 10 with 99 steps, equall spaced in log-space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2fe572d-b7c4-43f0-8a71-7a369cd31bb8",
   "metadata": {},
   "outputs": [],
   "source": [
    "nt = 100\n",
    "t = np.logspace(-3, 1, nt)\n",
    "h = ml.head(x=50, y=0, t=t)\n",
    "plt.semilogx(t, h[0])\n",
    "plt.semilogx(t, h[1])\n",
    "plt.xlabel(\"time (d)\")\n",
    "plt.ylabel(\"head (m)\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4837450-b195-4d95-bb1e-cd87a8984f9f",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "### Sythetic calibration example\n",
    "A synthetic example is presented \n",
    "The head in the river is varied 20 times (daily) over a period of 20 days. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0d89768-445c-49ab-adb0-bddb4b2c27a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "ts = np.arange(0, 20)\n",
    "rng = np.random.default_rng(seed=22)\n",
    "hstar = rng.integers(low=0, high=10, size=20) * 0.1\n",
    "tsandh = list(zip(ts, hstar, strict=True))\n",
    "plt.figure(figsize=(8, 3))\n",
    "plt.step(ts, hstar, where=\"post\")\n",
    "plt.xlabel(\"time (d)\")\n",
    "plt.ylabel(\"river level (m)\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "980e6691-f4f9-40b6-953d-eccc14d63bfe",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "The head is computed at $x=50$ in the top aquifer from day 10 till day 20 with intervals of 0.1 day. This is the synthetic truth. A normally distributed error is added to the synthetic truth and this is called the (synthetic) observed head `hsyn`, which is plotted. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bdc9457b-f263-4edc-a7a5-c9a08eb43c59",
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tft.ModelXsection(naq=2, tmin=1e-3, tmax=1e2)\n",
    "\n",
    "tft.XsectionMaq(\n",
    "    model=ml,\n",
    "    x1=-np.inf,\n",
    "    x2=-L,\n",
    "    kaq=[k0, k1],\n",
    "    z=zside,\n",
    "    Saq=[Ss0, Ss1],\n",
    "    c=[cside0, c1],\n",
    "    topboundary=\"semi\",\n",
    "    name=\"left\",\n",
    ")\n",
    "\n",
    "tft.XsectionMaq(\n",
    "    model=ml,\n",
    "    x1=-L,\n",
    "    x2=L,\n",
    "    kaq=[k0, k1],\n",
    "    z=zriver,\n",
    "    Saq=[Ss0, Ss1],\n",
    "    c=[criver0, c1],\n",
    "    topboundary=\"semi\",\n",
    "    tsandhstar=tsandh,\n",
    "    name=\"river\",\n",
    ")\n",
    "\n",
    "tft.XsectionMaq(\n",
    "    model=ml,\n",
    "    x1=L,\n",
    "    x2=np.inf,\n",
    "    kaq=[k0, k1],\n",
    "    z=zside,\n",
    "    Saq=[Ss0, Ss1],\n",
    "    c=[cside0, c1],\n",
    "    topboundary=\"semi\",\n",
    "    name=\"right\",\n",
    ")\n",
    "\n",
    "ml.solve()\n",
    "\n",
    "tsyn = np.linspace(10, 20, 101)\n",
    "xsyn = 50\n",
    "hexact = ml.head(xsyn, 0, tsyn, layers=0)\n",
    "hsyn = hexact[0] + 0.01 * rng.standard_normal(len(tsyn))\n",
    "#\n",
    "plt.figure(figsize=(8, 3))\n",
    "plt.plot(tsyn, hsyn, \"k.\", label=\"layer 0\")\n",
    "plt.xlabel(\"time (d)\")\n",
    "plt.ylabel(\"head (m)\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29b4b6b5-97dd-420d-8e21-0a773e3780df",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Next, we pretend that we don't quite know the resistance of the semi-confining layer on\n",
    "the left and right sides nor do we know the resistance of the leaky river bed. We will\n",
    "create a `timflow` model and use the observed head to calibrate the model and estimate\n",
    "the resistance of the semi-confining layer the leakystream bed. As as first guess, we\n",
    "use $c=100$ d for both."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79930651-2658-4fe7-97bb-000d8ff6e205",
   "metadata": {},
   "outputs": [],
   "source": [
    "crivguess = 100.0\n",
    "csideguess = 100.0\n",
    "\n",
    "ml = tft.ModelXsection(naq=2, tmin=1e-3, tmax=1e2)\n",
    "\n",
    "tft.XsectionMaq(\n",
    "    model=ml,\n",
    "    x1=-np.inf,\n",
    "    x2=-L,\n",
    "    kaq=[k0, k1],\n",
    "    z=zside,\n",
    "    Saq=[Ss0, Ss1],\n",
    "    c=[csideguess, c1],\n",
    "    topboundary=\"semi\",\n",
    "    name=\"left\",\n",
    ")\n",
    "\n",
    "tft.XsectionMaq(\n",
    "    model=ml,\n",
    "    x1=-L,\n",
    "    x2=L,\n",
    "    kaq=[k0, k1],\n",
    "    z=zriver,\n",
    "    Saq=[Ss0, Ss1],\n",
    "    c=[crivguess, c1],\n",
    "    topboundary=\"semi\",\n",
    "    tsandhstar=tsandh,\n",
    "    name=\"river\",\n",
    ")\n",
    "\n",
    "tft.XsectionMaq(\n",
    "    model=ml,\n",
    "    x1=L,\n",
    "    x2=np.inf,\n",
    "    kaq=[k0, k1],\n",
    "    z=zside,\n",
    "    Saq=[Ss0, Ss1],\n",
    "    c=[csideguess, c1],\n",
    "    topboundary=\"semi\",\n",
    "    name=\"right\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "564f41b8-a18e-4a9b-8387-a94bc5283459",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Next, we create a `Calibrate` object for the model. We add the series to use for calibration. Next, we add two parameters to be calibrated. The first one on the resistance of the leaky stream bed. The second one is the resistance of the semi-confining layer on the left and right sides; the same value is used for both sides by specifying `inhoms=(\"left\", \"right\")`. Note that the inhoms are specified by the name (string) used to create them. Finally, the model is calibrated using the `fit` function. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42fe85a8-528a-40e8-9d0d-d92726df22d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "cal = tft.Calibrate(ml)\n",
    "cal.series(\n",
    "    name=\"synthetic\",\n",
    "    x=xsyn,\n",
    "    y=0,\n",
    "    layer=0,\n",
    "    t=tsyn,\n",
    "    h=hsyn,\n",
    ")\n",
    "\n",
    "cal.set_parameter(\n",
    "    name=\"c\",\n",
    "    layers=[0],\n",
    "    initial=12,\n",
    "    pmin=1.0,\n",
    "    pmax=1000.0,\n",
    "    inhoms=(\"river\"),\n",
    ")\n",
    "\n",
    "cal.set_parameter(\n",
    "    \"c\",\n",
    "    layers=[0],\n",
    "    initial=csideguess,\n",
    "    pmin=1.0,\n",
    "    pmax=1000.0,\n",
    "    inhoms=(\"left\", \"right\"),\n",
    ")\n",
    "\n",
    "cal.fit(report=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96a8d512-c5e7-4248-97ff-66de137ed2f4",
   "metadata": {},
   "source": [
    "All inhomogeneities are stored in a dictionary `inhomdict`, which is stored in `ml.aq`. For example, the resistance of the leaky layer on the left and right side of the river may be retrieved as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10d9dcc9-a07c-43c9-b01f-3c2ba47322d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"resistance values on left side:  \", ml.aq.inhomdict[\"left\"].c)\n",
    "print(\"resistance falues on right side: \", ml.aq.inhomdict[\"right\"].c)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e04a20bf-115a-47a8-be05-c62f927b749d",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Finally, the fitted head is plotted together with the 'observed' heads and the 'truth'. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e1a203f-6590-4cbd-b26c-ae7da65f6f8e",
   "metadata": {},
   "outputs": [],
   "source": [
    "hmodel = ml.head(xsyn, 0, tsyn, 0)\n",
    "plt.figure(figsize=(8, 3))\n",
    "plt.plot(tsyn, hsyn, \"k.\", label=\"synthetic\")\n",
    "plt.plot(tsyn, hmodel[0], label=\"fit\")\n",
    "plt.plot(tsyn, hexact[0], \"--\", label=\"truth\")\n",
    "plt.xlabel(\"time (d)\")\n",
    "plt.ylabel(\"head (m)\")\n",
    "plt.legend()\n",
    "plt.grid()"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}