Head in leaky layer

docs/transient/02examples/head_in_vertical.ipynb

@@ -0,0 +1,214 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Head propagation in leaky layers"
+   ]
+  },
+  {
+   "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\n",
+    "\n",
+    "plt.rcParams[\"font.size\"] = 8.0"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Consider a three-aquifer system. Aquifer properties are given in Table 1. All aquifers have elastic storage. A well is located at $(x,y)=(0,0)$ and is screened in layer 1. The well starts pumping at time $t=0$ with a discharge $Q=1000$ m$^3$/d. The radius of the well is 0.2 m.\n",
+    "\n",
+    "**Table 1 - Aquifer properties for exercise 1**\n",
+    "|Layer          | $k$ (m/d) | $c$ (d) | $S_s$ (m$^{-1}$) | $z_t$ (m) | $z_b$ (m)|\n",
+    "|---------------| ---------:| -------:| -----:| ---------:| --------:|\n",
+    "|Aquifer 0      |      1    |    -    |0.0001 |   25      |        20|\n",
+    "|Leaky layer 1  |      -    |  1000   |0      |   20      |        18|\n",
+    "|Aquifer 1      |     20    |    -    |0.0001 |   18      |        10|\n",
+    "|Leaky layer 2  |      -    |  2000   |0      |   10      |         8|\n",
+    "|Aquifer 2      |      2    |    -    |0.0001 |    8      |         0|\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Confined top"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ml = tft.ModelMaq(\n",
+    "    kaq=[5, 10, 20],\n",
+    "    z=[25, 22, 18, 13, 8, 0],\n",
+    "    c=[1000, 2000],\n",
+    "    Saq=[1e-4, 1e-4, 1e-4],\n",
+    "    Sll=[1e-3, 1e-3],\n",
+    "    phreatictop=False,\n",
+    "    tmin=0.1,\n",
+    "    tmax=1000,\n",
+    ")\n",
+    "w1 = tft.Well(ml, xw=0, yw=0, rw=0.2, tsandQ=[(0, 1000)], layers=1)\n",
+    "w2 = tft.Well(ml, xw=20, yw=0, rw=0.2, tsandQ=[(0, 177)], layers=1)\n",
+    "ml.solve()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ml.head(10, 10, 1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "z = np.linspace(0, 25, 100)\n",
+    "h = ml.headinvertical(10, 10, z, [0.1, 1])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "z = np.linspace(0, 25, 100)\n",
+    "t = [0.1, 1, 10, 100]\n",
+    "h = ml.headinvertical(10, 10, z, t)\n",
+    "for i in range(len(t)):\n",
+    "    plt.plot(h[:, i], z)\n",
+    "for i in range(ml.aq.naq - 1):\n",
+    "    plt.axhspan(ml.aq.zaqtop[i + 1], ml.aq.zaqbot[i], color=[0.9, 0.9, 0.9])\n",
+    "plt.axhspan(25, 30, color=[0.6, 0.6, 0.6])\n",
+    "plt.ylim(z[0], 30);"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Semi-confined top"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ml = tft.ModelMaq(\n",
+    "    kaq=[5, 10, 20],\n",
+    "    z=[30, 25, 22, 18, 13, 8, 0],\n",
+    "    c=[500, 1000, 2000],\n",
+    "    Saq=[1e-4, 1e-4, 1e-4],\n",
+    "    Sll=[1e-3, 1e-3, 1e-3],\n",
+    "    topboundary=\"semi\",\n",
+    "    tmin=0.1,\n",
+    "    tmax=1000,\n",
+    ")\n",
+    "w1 = tft.Well(ml, xw=0, yw=0, rw=0.2, tsandQ=[(0, 1000)], layers=1)\n",
+    "w2 = tft.Well(ml, xw=20, yw=0, rw=0.2, tsandQ=[(0, 177)], layers=1)\n",
+    "ml.solve()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "z = np.linspace(0, 30, 100)\n",
+    "t = [0.1, 1, 10, 100]\n",
+    "h = ml.headinvertical(10, 10, z, t)\n",
+    "for i in range(len(t)):\n",
+    "    plt.plot(h[:, i], z)\n",
+    "for i in range(len(ml.aq.z[::2]) - 1):\n",
+    "    plt.axhspan(ml.aq.z[2 * i + 1], ml.aq.z[2 * i], color=[0.9, 0.9, 0.9])\n",
+    "plt.ylim(z[0], z[-1]);"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Leaky layer top"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ml = tft.ModelMaq(\n",
+    "    kaq=[5, 10, 20],\n",
+    "    z=[30, 25, 22, 18, 13, 8, 0],\n",
+    "    c=[500, 1000, 2000],\n",
+    "    Saq=[1e-4, 1e-4, 1e-4],\n",
+    "    Sll=[1e-3, 1e-3, 1e-3],\n",
+    "    topboundary=\"lea\",\n",
+    "    tmin=0.1,\n",
+    "    tmax=1000,\n",
+    ")\n",
+    "w1 = tft.Well(ml, xw=0, yw=0, rw=0.2, tsandQ=[(0, 1000)], layers=1)\n",
+    "w2 = tft.Well(ml, xw=20, yw=0, rw=0.2, tsandQ=[(0, 177)], layers=1)\n",
+    "ml.solve()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "z = np.linspace(0, 30, 100)\n",
+    "t = [0.1, 1, 10, 100]\n",
+    "h = ml.headinvertical(10, 10, z, t)\n",
+    "for i in range(len(t)):\n",
+    "    plt.plot(h[:, i], z)\n",
+    "for i in range(len(ml.aq.z[::2]) - 1):\n",
+    "    plt.axhspan(ml.aq.z[2 * i + 1], ml.aq.z[2 * i], color=[0.9, 0.9, 0.9])\n",
+    "plt.ylim(z[0], z[-1]);"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "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.5"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}

docs/transient/02examples/index.rst

@@ -59,4 +59,16 @@ Pathlines
     :caption: Pathlines
     :hidden:
 
-    pathline_trace
+    pathline_trace
+
+Head in leaky layers
+--------------------
+
+- :doc:`head_in_vertical`
+
+.. toctree::
+    :maxdepth: 2
+    :caption: Head in leaky layers
+    :hidden:
+
+    head_in_vertical

timflow/transient/aquifer.py

@@ -24,6 +24,7 @@ def __init__(
         Saq,
         Sll,
         leffaq,
+        leffll,
         poraq,
         porll,
         ltype,
@@ -49,6 +50,7 @@ def __init__(
         self.Sll = np.atleast_1d(Sll).astype(float)
         self.Sll[self.Sll < 1e-20] = 1e-20  # Cannot be zero
         self.leffaq = np.atleast_1d(leffaq).astype(float)
+        self.leffll = np.atleast_1d(leffll).astype(float)
         self.poraq = np.atleast_1d(poraq).astype(float)
         self.porll = np.atleast_1d(porll).astype(float)
         self.ltype = np.atleast_1d(ltype)
@@ -113,6 +115,8 @@ def initialize(self):
         # Compute Saq and Sll
         self.Scoefaq = self.Saq * self.Haq
         self.Scoefll = self.Sll * self.Hll
+        # Compute kappa (hydraulic conductivity of leaky layer)
+        self.kappa = self.Hll / self.c
         if (self.topboundary == "con") and self.phreatictop:
             self.Scoefaq[0] = self.Scoefaq[0] / self.Haq[0]
         elif (self.topboundary == "lea") and self.phreatictop:
@@ -131,10 +135,15 @@ def initialize(self):
         bmat = np.diag(np.ones(self.naq))
         self.a = np.zeros((self.model.npval, len(self.c)), dtype=complex)
         self.b = np.zeros((self.model.npval, len(self.c)), dtype=complex)
+        self.d = np.zeros(self.model.npval, dtype=complex)
+        self.alpha = np.zeros((len(self.c), self.model.npval), dtype=complex)
+
         for i in range(self.model.npval):
-            w, v, a, b = self.compute_lab_eigvec(self.model.p[i])
+            w, v, a, b, dzero = self.compute_lab_eigvec(self.model.p[i])
             self.a[i] = a
             self.b[i] = b
+            self.d[i] = dzero
+            self.alpha[:, i] = np.sqrt(self.model.p[i] * self.Sll / self.kappa)
             # Eigenvectors are columns of v
             self.eigval[:, i] = w
             self.eigvec[:, :, i] = v
@@ -151,6 +160,7 @@ def initialize(self):
     def compute_lab_eigvec(self, p, returnA=False, B=None):
         sqrtpSc = np.sqrt(p * self.Scoefll * self.c)
         a, b = np.zeros_like(sqrtpSc), np.zeros_like(sqrtpSc)
+        dzero = 0j
         small = np.abs(sqrtpSc) < 200
         a[small] = sqrtpSc[small] / np.tanh(sqrtpSc[small])
         b[small] = sqrtpSc[small] / np.sinh(sqrtpSc[small])
@@ -185,7 +195,7 @@ def compute_lab_eigvec(self, p, returnA=False, B=None):
         index = np.argsort(abs(w))[::-1]
         w = w[index]
         v = v[:, index]
-        return w, v, a, b
+        return w, v, a, b, dzero
 
     def head_to_potential(self, h, layers):
         return h * self.Tcol[layers]
@@ -268,6 +278,7 @@ def __init__(
         Saq,
         Sll,
         leffaq,
+        leffll,
         poraq,
         porll,
         ltype,
@@ -286,6 +297,7 @@ def __init__(
             Saq,
             Sll,
             leffaq,
+            leffll,
             poraq,
             porll,
             ltype,

timflow/transient/aquifer_parameters.py

@@ -14,6 +14,7 @@ def param_maq(
     Saq=[0.001],
     Sll=[0],
     leffaq=0,
+    leffll=0,
     poraq=[0.3],
     porll=[0.3],
     topboundary="conf",
@@ -28,6 +29,7 @@ def param_maq(
     Sll = np.atleast_1d(Sll).astype(float)
     porll = np.atleast_1d(porll).astype(float)
     leffaq = np.atleast_1d(leffaq).astype(float)
+    leffll = np.atleast_1d(leffll).astype(float)
     H = z[:-1] - z[1:]
     assert np.all(H >= 0), (
         "Error: Not all layer thicknesses are" + " non-negative" + str(H)
@@ -93,20 +95,23 @@ def param_maq(
             porll = porll * np.ones(naq)
         if len(leffaq) == 1:
             leffaq = leffaq * np.ones(naq)
+        if len(leffll) == 1:
+            leffll = leffll * np.ones(naq)
         assert len(kaq) == naq, "Error: Length of kaq needs to be " + str(naq)
         assert len(Saq) == naq, "Error: Length of Saq needs to be " + str(naq)
         assert len(poraq) == naq, "Error: Length of poraq needs to be " + str(naq)
         assert len(c) == naq, "Error: Length of c needs to be " + str(naq)
         assert len(Sll) == naq, "Error: Length of Sll needs to be " + str(naq)
         assert len(porll) == naq, "Error: Length of porll needs to be " + str(naq)
         assert len(leffaq) == naq, "Error: Length of leffaq needs to be " + str(naq)
+        assert len(leffll) == naq, "Error: Length of leffll needs to be " + str(naq)
         Haq = H[1::2]
         Hll = H[::2]
         # layertype
         nlayers = len(z) - 1
         ltype = np.array(nlayers * ["a"])
         ltype[0::2] = "l"
-    return kaq, Haq, Hll, c, Saq, Sll, leffaq, poraq, porll, ltype
+    return kaq, Haq, Hll, c, Saq, Sll, leffaq, leffll, poraq, porll, ltype
 
 
 def param_3d(
@@ -115,6 +120,7 @@ def param_3d(
     Saq=[0.001],
     kzoverkh=1,
     leffaq=0,
+    leffll=0,
     poraq=0.3,
     phreatictop=False,
     topboundary="conf",
@@ -138,6 +144,9 @@ def param_3d(
     leffaq = np.atleast_1d(leffaq).astype(float)
     if len(leffaq) == 1:
         leffaq = leffaq * np.ones(naq)
+    leffll = np.atleast_1d(leffll).astype(float)
+    if len(leffll) == 1:
+        leffll = leffll * np.ones(naq)
     poraq = np.atleast_1d(poraq).astype(float)
     if len(poraq) == 1:
         poraq = poraq * np.ones(naq)
@@ -161,4 +170,4 @@ def param_3d(
         porll[0] = toppor
         ltype = np.hstack(("l", ltype))
         z = np.hstack((z[0] + topthick, z))
-    return kaq, Haq, Hll, c, Saq, Sll, leffaq, poraq, porll, ltype, z
+    return kaq, Haq, Hll, c, Saq, Sll, leffaq, leffll, poraq, porll, ltype, z