{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d84c9e3a",
   "metadata": {},
   "source": [
    "# Modeling infiltration in a cross-section\n",
    "\n",
    "This notebook shows how to model infiltration in a transient cross-section model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee675c08-ba03-4d14-9d77-73c8dc790d8f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "import timflow.transient as tft\n",
    "\n",
    "plt.rcParams[\"figure.figsize\"] = (10, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "edcb6ae1-626a-447f-baec-36e5bafd02cb",
   "metadata": {},
   "source": [
    "## Single aquifer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04bafb55-0a73-497e-a254-ce5709c170a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "k = [20.0]\n",
    "H = 10.0\n",
    "Saq = [0.1]\n",
    "L = 100.0\n",
    "N = 1e-3\n",
    "z = [0, -H]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3fc16ea-80fc-489b-b9a3-fefd76bdc96d",
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tft.ModelXsection(naq=1, tmin=0.1, tmax=1e3)\n",
    "\n",
    "left = tft.XsectionMaq(\n",
    "    model=ml,\n",
    "    x1=-np.inf,\n",
    "    x2=-L / 2,\n",
    "    kaq=k,\n",
    "    z=z,\n",
    "    Saq=Saq,\n",
    "    topboundary=\"conf\",\n",
    "    phreatictop=True,\n",
    ")\n",
    "rech = tft.XsectionMaq(\n",
    "    model=ml,\n",
    "    x1=-L / 2,\n",
    "    x2=L / 2,\n",
    "    kaq=k,\n",
    "    z=z,\n",
    "    Saq=Saq,\n",
    "    topboundary=\"conf\",\n",
    "    phreatictop=True,\n",
    "    tsandN=[(0.0, N)],\n",
    ")\n",
    "right = tft.XsectionMaq(\n",
    "    model=ml,\n",
    "    x1=L / 2,\n",
    "    x2=np.inf,\n",
    "    kaq=k,\n",
    "    z=z,\n",
    "    Saq=Saq,\n",
    "    topboundary=\"conf\",\n",
    "    phreatictop=True,\n",
    ")\n",
    "\n",
    "ml.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bcd575db-4c42-40a5-bacb-564dfbd4e391",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(-100, 100, 100)\n",
    "y = np.zeros(100)\n",
    "plt.axvspan(-50, 50, color=3 * [0.9])\n",
    "for t in np.logspace(-1, 1, 10):\n",
    "    h = ml.headalongline(x, y, t)\n",
    "    plt.plot(x, h[0, 0], label=f\"t={t:.2f}\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "206b66c5-0ab4-4ce1-b446-654e7979529e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# check solution inside\n",
    "x = 25\n",
    "y = 0\n",
    "t = 7\n",
    "d = 1e-2\n",
    "\n",
    "# discharge vector\n",
    "Qxtim, _ = ml.disvec(x, y, t)\n",
    "Qxnum = rech.T * (ml.head(x - d, y, t) - ml.head(x + d, y, t)) / (2 * d)\n",
    "print(\"Qxtim, Qxnum: \", Qxtim, Qxnum)\n",
    "\n",
    "# deq\n",
    "d2hdx2 = (\n",
    "    ml.head(x - d, y, t)\n",
    "    + ml.head(x + d, y, t)\n",
    "    + ml.head(x, y - d, t)\n",
    "    + ml.head(x, y + d, t)\n",
    "    - 4 * ml.head(x, y, t)\n",
    ") / d**2\n",
    "dhdt = (ml.head(x, y, t + d) - ml.head(x, y, t - d)) / (2 * d)\n",
    "print(\"lhs: \", d2hdx2)\n",
    "print(\"rhs: \", rech.Saq / rech.T * dhdt - N / rech.T)\n",
    "\n",
    "# check solution outside\n",
    "x = 75\n",
    "y = 0\n",
    "t = 7\n",
    "d = 1e-2\n",
    "\n",
    "# discharge vector\n",
    "Qxtim, _ = ml.disvec(x, y, t)\n",
    "Qxnum = rech.T * (ml.head(x - d, y, t) - ml.head(x + d, y, t)) / (2 * d)\n",
    "print(\"Qxtim, Qxnum: \", Qxtim, Qxnum)\n",
    "\n",
    "# deq\n",
    "d2hdx2 = (\n",
    "    ml.head(x - d, y, t)\n",
    "    + ml.head(x + d, y, t)\n",
    "    + ml.head(x, y - d, t)\n",
    "    + ml.head(x, y + d, t)\n",
    "    - 4 * ml.head(x, y, t)\n",
    ") / d**2\n",
    "dhdt = (ml.head(x, y, t + d) - ml.head(x, y, t - d)) / (2 * d)\n",
    "print(\"lhs: \", d2hdx2)\n",
    "print(\"rhs: \", rech.Saq / rech.T * dhdt)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42eb3f80-aaef-4e87-8905-9655b5a68baa",
   "metadata": {},
   "source": [
    "## 2 aquifers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a400239-1f34-4d1b-b69b-b4deb028334c",
   "metadata": {},
   "outputs": [],
   "source": [
    "k = [10.0, 20.0]\n",
    "z = [0, -10, -12, -20]\n",
    "c = [500]\n",
    "Saq = [0.1, 1e-4]\n",
    "L = 100.0\n",
    "N = 1e-3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a547698-5068-4c15-b37c-cf69bd44fba6",
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tft.ModelXsection(naq=2, tmin=0.1, tmax=1e3)\n",
    "\n",
    "left = tft.XsectionMaq(\n",
    "    model=ml,\n",
    "    x1=-np.inf,\n",
    "    x2=-L / 2,\n",
    "    kaq=k,\n",
    "    z=z,\n",
    "    Saq=Saq,\n",
    "    c=c,\n",
    "    topboundary=\"conf\",\n",
    "    phreatictop=True,\n",
    ")\n",
    "\n",
    "inf = tft.XsectionMaq(\n",
    "    model=ml,\n",
    "    x1=-L / 2,\n",
    "    x2=L / 2,\n",
    "    kaq=k,\n",
    "    z=z,\n",
    "    Saq=Saq,\n",
    "    c=c,\n",
    "    topboundary=\"conf\",\n",
    "    phreatictop=True,\n",
    "    tsandN=[(0.0, N)],\n",
    ")\n",
    "\n",
    "right = tft.XsectionMaq(\n",
    "    model=ml,\n",
    "    x1=L / 2,\n",
    "    x2=np.inf,\n",
    "    kaq=k,\n",
    "    z=z,\n",
    "    Saq=Saq,\n",
    "    c=c,\n",
    "    topboundary=\"conf\",\n",
    "    phreatictop=True,\n",
    ")\n",
    "\n",
    "ml.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0abfd9d0-1bcb-43d0-a540-e19fc2aff301",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(-100, 100, 100)\n",
    "y = np.zeros(100)\n",
    "plt.axvspan(-50, 50, color=3 * [0.9])\n",
    "t = np.logspace(-1, 1, 3)\n",
    "h = ml.headalongline(x, y, t)\n",
    "for i in range(len(t)):\n",
    "    plt.plot(x, h[0, i], color=\"C0\")\n",
    "    plt.plot(x, h[1, i], color=\"C1\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5198e9b-9126-412b-acbe-d937b4bcf7dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"check solution inside\")\n",
    "x = 25\n",
    "y = 0\n",
    "t = 0.7\n",
    "d = 1e-3\n",
    "dt = 0.001\n",
    "\n",
    "# discharge vector\n",
    "Qxtim, _ = ml.disvec(x, y, t)\n",
    "Qxnum = inf.Tcol * (ml.head(x - d, y, t) - ml.head(x + d, y, t)) / (2 * d)\n",
    "print(\"Qxtim: \", Qxtim[:, 0])\n",
    "print(\"Qxnum: \", Qxnum[:, 0])\n",
    "\n",
    "# deq\n",
    "d2hdx2 = (\n",
    "    ml.head(x - d, y, t)\n",
    "    + ml.head(x + d, y, t)\n",
    "    + ml.head(x, y - d, t)\n",
    "    + ml.head(x, y + d, t)\n",
    "    - 4 * ml.head(x, y, t)\n",
    ") / d**2\n",
    "dhdt = (ml.head(x, y, t + dt) - ml.head(x, y, t - dt)) / (2 * dt)\n",
    "h = ml.head(x, y, t)\n",
    "rhs1 = (\n",
    "    inf.Saq[0] / inf.T[0] * dhdt[0, 0]\n",
    "    + (h[0, 0] - h[1, 0]) / (inf.T[0] * inf.c[1])\n",
    "    - N / inf.T[0]\n",
    ")\n",
    "rhs2 = inf.Saq[1] * inf.Haq[1] / inf.T[1] * dhdt[1, 0] - (h[0, 0] - h[1, 0]) / (\n",
    "    inf.T[1] * inf.c[1]\n",
    ")\n",
    "print(\"lhs: \", d2hdx2[:, 0])\n",
    "print(\"rhs: \", rhs1, rhs2)\n",
    "\n",
    "print(\"check solution outside\")\n",
    "x = 75\n",
    "y = 0\n",
    "t = 0.7\n",
    "d = 0.01\n",
    "dt = 0.01\n",
    "\n",
    "# discharge vector\n",
    "Qxtim, _ = ml.disvec(x, y, t)\n",
    "Qxnum = inf.Tcol * (ml.head(x - d, y, t) - ml.head(x + d, y, t)) / (2 * d)\n",
    "print(\"Qxtim: \", Qxtim[:, 0])\n",
    "print(\"Qxnum: \", Qxnum[:, 0])\n",
    "\n",
    "# deq\n",
    "d2hdx2 = (\n",
    "    ml.head(x - d, y, t)\n",
    "    + ml.head(x + d, y, t)\n",
    "    + ml.head(x, y - d, t)\n",
    "    + ml.head(x, y + d, t)\n",
    "    - 4 * ml.head(x, y, t)\n",
    ") / d**2\n",
    "dhdt = (ml.head(x, y, t + dt) - ml.head(x, y, t - dt)) / (2 * dt)\n",
    "h = ml.head(x, y, t)\n",
    "rhs1 = inf.Saq[0] / inf.T[0] * dhdt[0, 0] + (h[0, 0] - h[1, 0]) / (inf.T[0] * inf.c[1])\n",
    "rhs2 = inf.Saq[1] * inf.Haq[1] / inf.T[1] * dhdt[1, 0] - (h[0, 0] - h[1, 0]) / (\n",
    "    inf.T[1] * inf.c[1]\n",
    ")\n",
    "print(\"lhs: \", d2hdx2[:, 0])\n",
    "print(\"rhs: \", rhs1, rhs2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83d50fa9-5112-4d89-9398-9f0031e2ed63",
   "metadata": {},
   "source": [
    "## Example river, 1 aquifer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15eb1fa0-a6bc-4e31-a4e5-da72dc46ef38",
   "metadata": {},
   "outputs": [],
   "source": [
    "k = [20.0]\n",
    "z = [11, 10, 0]\n",
    "Saq = [1e-3]\n",
    "hstar = 2\n",
    "L = 100\n",
    "\n",
    "ml = tft.ModelXsection(naq=1, tmin=1e-4, tmax=1e3)\n",
    "\n",
    "left = tft.XsectionMaq(\n",
    "    model=ml,\n",
    "    x1=-np.inf,\n",
    "    x2=-L / 2,\n",
    "    kaq=k,\n",
    "    z=z,\n",
    "    Saq=Saq,\n",
    "    c=[100],\n",
    "    topboundary=\"semi\",\n",
    ")\n",
    "\n",
    "riv = tft.XsectionMaq(\n",
    "    model=ml,\n",
    "    x1=-L / 2,\n",
    "    x2=L / 2,\n",
    "    kaq=k,\n",
    "    z=z,\n",
    "    Saq=Saq,\n",
    "    c=[100],\n",
    "    topboundary=\"semi\",\n",
    "    tsandhstar=[(0.0, hstar), (5, 2 * hstar)],\n",
    ")\n",
    "\n",
    "right = tft.XsectionMaq(\n",
    "    model=ml,\n",
    "    x1=L / 2,\n",
    "    x2=np.inf,\n",
    "    kaq=k,\n",
    "    z=z,\n",
    "    Saq=Saq,\n",
    "    c=[100],\n",
    "    topboundary=\"semi\",\n",
    ")\n",
    "\n",
    "ml.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83f1e071-5eda-4833-a3a4-29733ce1bf7b",
   "metadata": {},
   "outputs": [],
   "source": [
    "lab = np.sqrt(riv.T[0] * riv.c[0])\n",
    "print(f\"leakage factor: {lab} m\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8316af2c-a2c5-4f91-b95c-f7ed46cecc98",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(-100, 100, 100)\n",
    "y = np.zeros(100)\n",
    "plt.axvspan(-50, 50, color=3 * [0.9])\n",
    "for t in np.logspace(-1, 1, 10):\n",
    "    h = ml.headalongline(x, y, t)\n",
    "    plt.plot(x, h[0, 0], label=f\"t={t:.2f}\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01260162-6e08-467f-88b2-052ebeae172e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# mf6 time step 0.1 day, cell size 10 m\n",
    "hmf6 = np.loadtxt(\"./data/mf6_ttim_xsec_riv1.txt\")  # load mf6 output to compare"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9563f8db-6d70-4653-bd61-825f9551eca2",
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.linspace(0.01, 10, 100)\n",
    "h = ml.head(25, 0, t)\n",
    "plt.plot(t, h[0], label=\"timflow\")\n",
    "plt.plot(hmf6[0], hmf6[1], \"--\", label=\"mf6\")\n",
    "plt.xticks(np.linspace(0, 10, 5))\n",
    "plt.legend()\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1913c367-e966-482a-b127-ebba9502d959",
   "metadata": {},
   "outputs": [],
   "source": [
    "# check solution inside\n",
    "x = 25\n",
    "y = 0\n",
    "t = 1\n",
    "d = 0.01\n",
    "dt = 1e-3\n",
    "\n",
    "# discharge vector\n",
    "Qxtim, _ = ml.disvec(x, y, t)\n",
    "Qxnum = riv.T * (ml.head(x - d, y, t) - ml.head(x + d, y, t)) / (2 * d)\n",
    "print(\"Qxtim, Qxnum: \", Qxtim, Qxnum)\n",
    "\n",
    "# deq\n",
    "h = ml.head(x, y, t)\n",
    "d2hdx2 = (\n",
    "    ml.head(x - d, y, t)\n",
    "    + ml.head(x + d, y, t)\n",
    "    + ml.head(x, y - d, t)\n",
    "    + ml.head(x, y + d, t)\n",
    "    - 4 * ml.head(x, y, t)\n",
    ") / d**2\n",
    "dhdt = (ml.head(x, y, t + dt) - ml.head(x, y, t - dt)) / (2 * dt)\n",
    "print(\"lhs: \", d2hdx2)\n",
    "print(\n",
    "    \"rhs: \",\n",
    "    (riv.Saq[0] * riv.Haq[0]) / riv.T * dhdt\n",
    "    + h / (riv.c[0] * riv.T[0])\n",
    "    - hstar / (riv.c[0] * riv.T[0]),\n",
    ")\n",
    "# print(\n",
    "#     \"rhs: \",\n",
    "#     (ml.aq.Saq[0] * ml.aq.Haq[0]) / ml.aq.T * dhdt\n",
    "#     + (h - hstar) / (ml.aq.c[0] * ml.aq.T[0]),\n",
    "# )  # - hstar / (ml.aq.c[0] * ml.aq.T[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f444c47-a827-4e59-acda-49850d277666",
   "metadata": {},
   "source": [
    "## Very wide river"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d8c6380-c88d-4d46-a3f0-f4211a7396a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "k = [20.0]\n",
    "z = [11, 10, 0]\n",
    "c = 100\n",
    "Saq = [1e-3]\n",
    "hstar = 2\n",
    "L = 2000\n",
    "\n",
    "ml = tft.ModelXsection(naq=1, tmin=1e-2, tmax=1e3)\n",
    "\n",
    "left = tft.XsectionMaq(\n",
    "    model=ml,\n",
    "    x1=-np.inf,\n",
    "    x2=-L / 2,\n",
    "    kaq=k,\n",
    "    z=z,\n",
    "    Saq=Saq,\n",
    "    c=[100],\n",
    "    topboundary=\"semi\",\n",
    ")\n",
    "\n",
    "riv = tft.XsectionMaq(\n",
    "    model=ml,\n",
    "    x1=-L / 2,\n",
    "    x2=L / 2,\n",
    "    kaq=k,\n",
    "    z=z,\n",
    "    Saq=Saq,\n",
    "    c=[100],\n",
    "    topboundary=\"semi\",\n",
    "    tsandhstar=[(0.0, hstar)],\n",
    ")\n",
    "\n",
    "right = tft.XsectionMaq(\n",
    "    model=ml,\n",
    "    x1=L / 2,\n",
    "    x2=np.inf,\n",
    "    kaq=k,\n",
    "    z=z,\n",
    "    Saq=Saq,\n",
    "    c=[100],\n",
    "    topboundary=\"semi\",\n",
    ")\n",
    "\n",
    "ml.solve()\n",
    "\n",
    "x = np.linspace(-L, L, 100)\n",
    "y = np.zeros(100)\n",
    "plt.axvspan(-L / 2, L / 2, color=3 * [0.9])\n",
    "for t in np.logspace(-1, 2, 10):\n",
    "    h = ml.headalongline(x, y, t)\n",
    "    plt.plot(x, h[0, 0], label=f\"t={t:.2f}\")\n",
    "plt.legend()\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b415f60b-8685-4293-92ff-58251582992b",
   "metadata": {},
   "source": [
    "## Two aquifer "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16d3a7eb-e74d-471f-9d09-9b568413cbcc",
   "metadata": {},
   "outputs": [],
   "source": [
    "k = [20.0, 40]\n",
    "z = [11, 10, 0, -2, -12]\n",
    "c = 100\n",
    "Saq = [1e-3, 1e-3]\n",
    "hstar = 2\n",
    "L = 100\n",
    "\n",
    "ml = tft.ModelXsection(naq=2, tmin=1e-4, tmax=1e3)\n",
    "\n",
    "left = tft.XsectionMaq(\n",
    "    model=ml,\n",
    "    x1=-np.inf,\n",
    "    x2=-L / 2,\n",
    "    kaq=k,\n",
    "    z=z,\n",
    "    Saq=Saq,\n",
    "    c=c,\n",
    "    topboundary=\"semi\",\n",
    ")\n",
    "\n",
    "riv = tft.XsectionMaq(\n",
    "    model=ml,\n",
    "    x1=-L / 2,\n",
    "    x2=L / 2,\n",
    "    kaq=k,\n",
    "    z=z,\n",
    "    Saq=Saq,\n",
    "    c=[100],\n",
    "    topboundary=\"semi\",\n",
    "    tsandhstar=[(0.0, hstar), (5, 2 * hstar)],\n",
    ")\n",
    "\n",
    "right = tft.XsectionMaq(\n",
    "    model=ml,\n",
    "    x1=L / 2,\n",
    "    x2=np.inf,\n",
    "    kaq=k,\n",
    "    z=z,\n",
    "    Saq=Saq,\n",
    "    c=[100],\n",
    "    topboundary=\"semi\",\n",
    ")\n",
    "\n",
    "ml.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "958771a5-cd77-47a4-8969-ac6c86ffa433",
   "metadata": {},
   "outputs": [],
   "source": [
    "# check solution inside\n",
    "x = 25\n",
    "y = 0\n",
    "t = 1\n",
    "d = 0.01\n",
    "dt = 1e-3\n",
    "\n",
    "# discharge vector\n",
    "Qxtim, _ = ml.disvec(x, y, t)\n",
    "Qxnum = riv.Tcol * (ml.head(x - d, y, t) - ml.head(x + d, y, t)) / (2 * d)\n",
    "print(\"Qxtim, Qxnum: \")\n",
    "print(Qxtim)\n",
    "print(Qxnum)\n",
    "\n",
    "# deq\n",
    "h = ml.head(x, y, t)\n",
    "d2hdx2 = (\n",
    "    ml.head(x - d, y, t)\n",
    "    + ml.head(x + d, y, t)\n",
    "    + ml.head(x, y - d, t)\n",
    "    + ml.head(x, y + d, t)\n",
    "    - 4 * ml.head(x, y, t)\n",
    ") / d**2\n",
    "dhdt = (ml.head(x, y, t + dt) - ml.head(x, y, t - dt)) / (2 * dt)\n",
    "rhs0 = (\n",
    "    (riv.Saq[0] * riv.Haq[0]) / riv.T[0] * dhdt[0, 0]\n",
    "    + h[0, 0] / (riv.c[0] * riv.T[0])\n",
    "    - hstar / (riv.c[0] * riv.T[0])\n",
    "    + (h[0, 0] - h[1, 0]) / (riv.c[1] * riv.T[0])\n",
    ")\n",
    "rhs1 = (riv.Saq[1] * riv.Haq[1]) / riv.T[1] * dhdt[1, 0] + (h[1, 0] - h[0, 0]) / (\n",
    "    riv.c[1] * riv.T[1]\n",
    ")\n",
    "print(\"lhs: \", d2hdx2[:, 0])\n",
    "print(\"rhs: \", rhs0, rhs1)\n",
    "# print(\n",
    "#     \"rhs: \",\n",
    "#     (riv.Saq[0] * riv.Haq[0]) / riv.T * dhdt + (h - hstar) / (riv.c[0] * riv.T[0]),\n",
    "# )  # - hstar / (riv.c[0] * riv.T[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a12cd826-850d-49b5-9a1c-8315be649438",
   "metadata": {},
   "outputs": [],
   "source": [
    "# mf6 time step 0.1 day, cell size 10 m\n",
    "hmf6 = np.loadtxt(\"./data/mf6_ttim_xsec_riv2.txt\")  # load mf6 output to compare"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3be3b98-86b4-4cc2-b6c9-6598f08987ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.linspace(0.01, 10, 100)\n",
    "h = ml.head(25, 0, t)\n",
    "plt.plot(t, h[0], label=\"timflow lay 0\")\n",
    "plt.plot(t, h[1], label=\"timflow lay 1\")\n",
    "plt.plot(hmf6[0], hmf6[1], \"--\", label=\"mf6 lay 0\")\n",
    "plt.plot(hmf6[0], hmf6[2], \"--\", label=\"mf6 lay 1\")\n",
    "plt.xticks(np.linspace(0, 10, 5))\n",
    "plt.legend()\n",
    "plt.grid()"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}