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

@@ -11,6 +11,9 @@
 import numpy as np
 import pandas as pd
 
+from timflow.transient.invlapnumba import invlapcomp
+from timflow.transient.stripareasink import HstarXsection
+
 
 class AquiferData:
     def __init__(
@@ -24,6 +27,7 @@ def __init__(
         Saq,
         Sll,
         leffaq,
+        leffll,
         poraq,
         porll,
         ltype,
@@ -49,6 +53,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 +118,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 +138,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 +163,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 +198,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]
@@ -255,6 +268,39 @@ def summary(self):
         summary.loc[:, "layer"] = self.layernumber
         return summary  # .set_index("layer")
 
+    def headsemitoplayer(self, x, y, z, t):
+        z = np.atleast_1d(z)
+        t = np.atleast_1d(t)
+        rv = np.zeros((len(z), len(t)))
+        # find if there is HstarXsection
+        for e in self.model.elementlist:
+            if isinstance(e, HstarXsection):
+                if e.aq == self:
+                    if x >= e.x1 and x <= e.x2:
+                        break
+        else:
+            return rv
+        htop = np.zeros((1, self.model.npval), dtype=complex)
+        htop[0] = 1 / self.model.p
+        for iz in range(len(z)):
+            if z[iz] <= self.z[0] and z[iz] >= self.zaqtop[0]:
+                eta = (
+                    htop * np.sinh(self.alpha[0] * (z[iz] - self.zaqtop[0]))
+                ) / np.sinh(self.alpha[0] * self.Hll[0])
+                eta = eta[:, np.newaxis, :]
+                rv[iz] = invlapcomp(
+                    t,
+                    eta,
+                    self.model.npint,
+                    self.model.M,
+                    self.model.tintervals,
+                    np.array(len(e.tstart) * [0]),
+                    e.tstart,
+                    e.bc,
+                    1,
+                )
+        return rv
+
 
 class Aquifer(AquiferData):
     def __init__(
@@ -268,6 +314,7 @@ def __init__(
         Saq,
         Sll,
         leffaq,
+        leffll,
         poraq,
         porll,
         ltype,
@@ -286,6 +333,7 @@ def __init__(
             Saq,
             Sll,
             leffaq,
+            leffll,
             poraq,
             porll,
             ltype,