From d958102bcc636e095f32a08f79aafcb981616bce Mon Sep 17 00:00:00 2001
From: Mark Bakker <markbak@gmail.com>
Date: Wed, 18 Mar 2026 16:22:46 +0100
Subject: [PATCH 01/12] compute kappa and alpha

---
 timflow/transient/aquifer.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/timflow/transient/aquifer.py b/timflow/transient/aquifer.py
index 87acc8d..451cb15 100644
--- a/timflow/transient/aquifer.py
+++ b/timflow/transient/aquifer.py
@@ -113,6 +113,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 +133,13 @@ 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.alpha = np.zeros((self.model.npval, len(self.c)), dtype=complex)
+        
         for i in range(self.model.npval):
             w, v, a, b = self.compute_lab_eigvec(self.model.p[i])
             self.a[i] = a
             self.b[i] = b
+            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

From cb2ef3ceda6e7408da600482617ed68f5461dc40 Mon Sep 17 00:00:00 2001
From: Mark Bakker <markbak@gmail.com>
Date: Wed, 18 Mar 2026 21:34:42 +0100
Subject: [PATCH 02/12] start of computing head in ll

---
 timflow/transient/aquifer.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/timflow/transient/aquifer.py b/timflow/transient/aquifer.py
index 451cb15..27b98c2 100644
--- a/timflow/transient/aquifer.py
+++ b/timflow/transient/aquifer.py
@@ -133,13 +133,13 @@ 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.alpha = np.zeros((self.model.npval, len(self.c)), 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])
             self.a[i] = a
             self.b[i] = b
-            self.alpha[i] = np.sqrt(self.model.p[i] * self.Sll / self.kappa)
+            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

From fe77b5a37087261521eae3725ed8b1504fd16035 Mon Sep 17 00:00:00 2001
From: Mark Bakker <markbak@gmail.com>
Date: Wed, 18 Mar 2026 21:36:11 +0100
Subject: [PATCH 03/12] one more step

---
 timflow/transient/model.py | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/timflow/transient/model.py b/timflow/transient/model.py
index f0664f4..eb83c9a 100644
--- a/timflow/transient/model.py
+++ b/timflow/transient/model.py
@@ -352,6 +352,17 @@ def head(self, x, y, t, layers=None, aq=None, derivative=0, neglect_steady=False
                 h += htimml[:, np.newaxis]
         return h
 
+    def headll(self, x, y, t, ll_layer, aq=None):
+        if aq is None:
+            aq = self.aq.find_aquifer_data(x, y)
+        if ll_layer == 0:
+            print('not implemented yet for ll_layer = 0')
+            return
+        head = self.potential(x, y, t, returnphi=True) / aq.T[np.newaxis, :, np.newaxis]
+        eta = np.zeros_like(pot) # first row is leaky layer on top
+        #for i in range(1, aq.naq):
+        #    eta[:, 1:] = head[:, :-1] * np.sinh(alpha
+
     def velocompold(self, x, y, z, t, aq=None, layer_ltype=[0, 0]):
         # implemented for one layer
         if aq is None:

From 37dcb5f62ea68645e3f39bf14fb8fab49fed3b98 Mon Sep 17 00:00:00 2001
From: Mark Bakker <markbak@gmail.com>
Date: Thu, 19 Mar 2026 21:54:48 +0100
Subject: [PATCH 04/12] working on head in leaky layer

---
 timflow/transient/model.py | 33 +++++++++++++++++++++++++++------
 1 file changed, 27 insertions(+), 6 deletions(-)

diff --git a/timflow/transient/model.py b/timflow/transient/model.py
index eb83c9a..97ef803 100644
--- a/timflow/transient/model.py
+++ b/timflow/transient/model.py
@@ -352,16 +352,37 @@ def head(self, x, y, t, layers=None, aq=None, derivative=0, neglect_steady=False
                 h += htimml[:, np.newaxis]
         return h
 
-    def headll(self, x, y, t, ll_layer, aq=None):
+    def headll(self, x, y, z, t, aq=None, returneta=False):
         if aq is None:
             aq = self.aq.find_aquifer_data(x, y)
-        if ll_layer == 0:
+        lay, ltype, _ = aq.findlayer(z)
+        assert ltype == 'l'
+        headbar = self.potential(x, y, t, returnphi=True) / aq.T[np.newaxis, :, np.newaxis]
+        if lay == 0:
             print('not implemented yet for ll_layer = 0')
             return
-        head = self.potential(x, y, t, returnphi=True) / aq.T[np.newaxis, :, np.newaxis]
-        eta = np.zeros_like(pot) # first row is leaky layer on top
-        #for i in range(1, aq.naq):
-        #    eta[:, 1:] = head[:, :-1] * np.sinh(alpha
+        else:
+            eta = (headbar[:, lay - 1] * np.sinh(aq.alpha[lay] * (z - aq.zaqtop[lay])) +
+                   headbar[:, lay] * np.sinh(aq.alpha[lay] * (aq.zaqbot[lay - 1] - z))) / np.sinh(aq.alpha[lay] * aq.Hll[lay])
+        eta = eta[:, np.newaxis, :]
+        if returneta:
+            return eta
+        time = np.atleast_1d(t)
+        nlayers = 1
+        rv = invlapcomp(
+            time,
+            eta,
+            self.npint,
+            self.M,
+            self.tintervals,
+            self.enumber,
+            self.etstart,
+            self.ebc,
+            nlayers,
+        )
+        return rv
+
+        
 
     def velocompold(self, x, y, z, t, aq=None, layer_ltype=[0, 0]):
         # implemented for one layer

From b9cd0ed3d9baccfb574d2d55d3fcb9ffa9f9bbbd Mon Sep 17 00:00:00 2001
From: Mark Bakker <markbak@gmail.com>
Date: Fri, 20 Mar 2026 19:10:10 +0100
Subject: [PATCH 05/12] updated hlll

Can now compute heads for array of time and z. Only need to incorporate varying head on top of aquifer. Should probably be renamed headalongvertical
---
 timflow/transient/model.py | 69 ++++++++++++++++++++++++++------------
 1 file changed, 48 insertions(+), 21 deletions(-)

diff --git a/timflow/transient/model.py b/timflow/transient/model.py
index 97ef803..affcdc8 100644
--- a/timflow/transient/model.py
+++ b/timflow/transient/model.py
@@ -353,37 +353,64 @@ def head(self, x, y, t, layers=None, aq=None, derivative=0, neglect_steady=False
         return h
 
     def headll(self, x, y, z, t, aq=None, returneta=False):
+        """Head at x, y, t where t can be multiple times.
+
+        Parameters
+        ----------
+        x : float
+        y : float
+        z : float or array
+        t : list or array
+            times for which grid is returned
+
+        Returns
+        -------
+        h : array size `ntimes, nz`
+        """
+        z = np.atleast_1d(z)
+        t = np.atleast_1d(t)
+        rv = np.zeros((len(z), len(t)))
         if aq is None:
             aq = self.aq.find_aquifer_data(x, y)
-        lay, ltype, _ = aq.findlayer(z)
-        assert ltype == 'l'
-        headbar = self.potential(x, y, t, returnphi=True) / aq.T[np.newaxis, :, np.newaxis]
-        if lay == 0:
-            print('not implemented yet for ll_layer = 0')
-            return
-        else:
-            eta = (headbar[:, lay - 1] * np.sinh(aq.alpha[lay] * (z - aq.zaqtop[lay])) +
-                   headbar[:, lay] * np.sinh(aq.alpha[lay] * (aq.zaqbot[lay - 1] - z))) / np.sinh(aq.alpha[lay] * aq.Hll[lay])
-        eta = eta[:, np.newaxis, :]
-        if returneta:
-            return eta
-        time = np.atleast_1d(t)
-        nlayers = 1
-        rv = invlapcomp(
-            time,
-            eta,
+        headbar = self.potential(x, y, t, aq=aq, returnphi=True) / aq.T[np.newaxis, :, np.newaxis]
+        haq = invlapcomp(
+            t,
+            headbar,
             self.npint,
             self.M,
             self.tintervals,
             self.enumber,
             self.etstart,
             self.ebc,
-            nlayers,
-        )
+            aq.nlayers,
+            )
+        for iz in range(len(z)):       
+            lay, ltype, _ = aq.findlayer(z[iz])
+            if ltype == 'a':
+                rv[iz] = haq[lay]
+            elif ltype == 'l':
+                if lay == 0:
+                    rv[iz] = np.nan
+                else:
+                    eta = (headbar[:, lay - 1] * np.sinh(aq.alpha[lay] * (z[iz] - aq.zaqtop[lay])) +
+                           headbar[:, lay] * np.sinh(aq.alpha[lay] * (aq.zaqbot[lay - 1] - z[iz]))) / np.sinh(aq.alpha[lay] * aq.Hll[lay])
+                    eta = eta[:, np.newaxis, :]
+                    if returneta:
+                        return eta
+                    nlayers = 1
+                    rv[iz] = invlapcomp(
+                        t,
+                        eta,
+                        self.npint,
+                        self.M,
+                        self.tintervals,
+                        self.enumber,
+                        self.etstart,
+                        self.ebc,
+                        nlayers,
+                    )
         return rv
 
-        
-
     def velocompold(self, x, y, z, t, aq=None, layer_ltype=[0, 0]):
         # implemented for one layer
         if aq is None:

From 3b30d0a33096e1028c0ba8039428fbf7b15bbd1c Mon Sep 17 00:00:00 2001
From: Mark Bakker <markbak@gmail.com>
Date: Mon, 23 Mar 2026 13:23:46 +0100
Subject: [PATCH 06/12] add head in  leaky layer on top

also renamed function to headinvertical and added example
---
 .../02examples/head_in_vertical.ipynb         | 170 ++++++++++++++++++
 timflow/transient/aquifer.py                  |   2 +-
 timflow/transient/model.py                    |  67 ++++---
 3 files changed, 209 insertions(+), 30 deletions(-)
 create mode 100755 docs/transient/02examples/head_in_vertical.ipynb

diff --git a/docs/transient/02examples/head_in_vertical.ipynb b/docs/transient/02examples/head_in_vertical.ipynb
new file mode 100755
index 0000000..9d8e46f
--- /dev/null
+++ b/docs/transient/02examples/head_in_vertical.ipynb
@@ -0,0 +1,170 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Multi-layer example"
+   ]
+  },
+  {
+   "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": [
+    "plt.plot(h[:, 0], z)\n",
+    "plt.plot(h[:, 1], 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])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Semi-onfined 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",
+    "h = ml.headinvertical(10, 10, z, [0.1, 1])\n",
+    "plt.plot(h[:, 0], z)\n",
+    "plt.plot(h[:, 1], 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])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "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
+}
diff --git a/timflow/transient/aquifer.py b/timflow/transient/aquifer.py
index 27b98c2..a649828 100644
--- a/timflow/transient/aquifer.py
+++ b/timflow/transient/aquifer.py
@@ -134,7 +134,7 @@ def initialize(self):
         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.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])
             self.a[i] = a
diff --git a/timflow/transient/model.py b/timflow/transient/model.py
index affcdc8..fbe0fa4 100644
--- a/timflow/transient/model.py
+++ b/timflow/transient/model.py
@@ -352,16 +352,15 @@ def head(self, x, y, t, layers=None, aq=None, derivative=0, neglect_steady=False
                 h += htimml[:, np.newaxis]
         return h
 
-    def headll(self, x, y, z, t, aq=None, returneta=False):
-        """Head at x, y, t where t can be multiple times.
+    def headinvertical(self, x, y, z, t, aq=None, returneta=False):
+        """Head along vertical line at x, y, z, t.
 
         Parameters
         ----------
         x : float
         y : float
-        z : float or array
-        t : list or array
-            times for which grid is returned
+        z : float, list, or array
+        t : float, list, or array
 
         Returns
         -------
@@ -372,7 +371,10 @@ def headll(self, x, y, z, t, aq=None, returneta=False):
         rv = np.zeros((len(z), len(t)))
         if aq is None:
             aq = self.aq.find_aquifer_data(x, y)
-        headbar = self.potential(x, y, t, aq=aq, returnphi=True) / aq.T[np.newaxis, :, np.newaxis]
+        headbar = (
+            self.potential(x, y, t, aq=aq, returnphi=True)
+            / aq.T[np.newaxis, :, np.newaxis]
+        )
         haq = invlapcomp(
             t,
             headbar,
@@ -382,33 +384,40 @@ def headll(self, x, y, z, t, aq=None, returneta=False):
             self.enumber,
             self.etstart,
             self.ebc,
-            aq.nlayers,
-            )
-        for iz in range(len(z)):       
+            aq.naq,
+        )
+        for iz in range(len(z)):
             lay, ltype, _ = aq.findlayer(z[iz])
-            if ltype == 'a':
+            if ltype == "a":
                 rv[iz] = haq[lay]
-            elif ltype == 'l':
+            elif ltype == "l":
                 if lay == 0:
-                    rv[iz] = np.nan
-                else:
-                    eta = (headbar[:, lay - 1] * np.sinh(aq.alpha[lay] * (z[iz] - aq.zaqtop[lay])) +
-                           headbar[:, lay] * np.sinh(aq.alpha[lay] * (aq.zaqbot[lay - 1] - z[iz]))) / np.sinh(aq.alpha[lay] * aq.Hll[lay])
-                    eta = eta[:, np.newaxis, :]
-                    if returneta:
-                        return eta
-                    nlayers = 1
-                    rv[iz] = invlapcomp(
-                        t,
-                        eta,
-                        self.npint,
-                        self.M,
-                        self.tintervals,
-                        self.enumber,
-                        self.etstart,
-                        self.ebc,
-                        nlayers,
+                    eta = (
+                        headbar[:, lay]
+                        * np.sinh(aq.alpha[lay] * (aq.z[lay] - z[iz]))
+                        / np.sinh(aq.alpha[lay] * aq.Hll[lay])
                     )
+                else:
+                    eta = (
+                        headbar[:, lay - 1]
+                        * np.sinh(aq.alpha[lay] * (z[iz] - aq.zaqtop[lay]))
+                        + headbar[:, lay]
+                        * np.sinh(aq.alpha[lay] * (aq.zaqbot[lay - 1] - z[iz]))
+                    ) / np.sinh(aq.alpha[lay] * aq.Hll[lay])
+                eta = eta[:, np.newaxis, :]
+                if returneta:
+                    return eta
+                rv[iz] = invlapcomp(
+                    t,
+                    eta,
+                    self.npint,
+                    self.M,
+                    self.tintervals,
+                    self.enumber,
+                    self.etstart,
+                    self.ebc,
+                    1,
+                )
         return rv
 
     def velocompold(self, x, y, z, t, aq=None, layer_ltype=[0, 0]):

From 3fe40f8e50b75b659d31c9cb8d45fae9063e0c21 Mon Sep 17 00:00:00 2001
From: Mark Bakker <markbak@gmail.com>
Date: Mon, 23 Mar 2026 16:14:27 +0100
Subject: [PATCH 07/12] added leaky top boundary for headinvertical

---
 .../02examples/head_in_vertical.ipynb         | 62 ++++++++++++++++---
 timflow/transient/aquifer.py                  |  7 ++-
 timflow/transient/model.py                    | 17 +++--
 3 files changed, 70 insertions(+), 16 deletions(-)

diff --git a/docs/transient/02examples/head_in_vertical.ipynb b/docs/transient/02examples/head_in_vertical.ipynb
index 9d8e46f..3465c53 100755
--- a/docs/transient/02examples/head_in_vertical.ipynb
+++ b/docs/transient/02examples/head_in_vertical.ipynb
@@ -90,17 +90,22 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plt.plot(h[:, 0], z)\n",
-    "plt.plot(h[:, 1], z)\n",
+    "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])"
+    "    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-onfined top"
+    "## Semi-confined top"
    ]
   },
   {
@@ -131,11 +136,20 @@
    "outputs": [],
    "source": [
     "z = np.linspace(0, 30, 100)\n",
-    "h = ml.headinvertical(10, 10, z, [0.1, 1])\n",
-    "plt.plot(h[:, 0], z)\n",
-    "plt.plot(h[:, 1], z)\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])"
+    "    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"
    ]
   },
   {
@@ -143,7 +157,37 @@
    "execution_count": null,
    "metadata": {},
    "outputs": [],
-   "source": []
+   "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": {
diff --git a/timflow/transient/aquifer.py b/timflow/transient/aquifer.py
index a649828..d5906ac 100644
--- a/timflow/transient/aquifer.py
+++ b/timflow/transient/aquifer.py
@@ -133,12 +133,14 @@ 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
@@ -156,6 +158,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])
@@ -190,7 +193,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]
diff --git a/timflow/transient/model.py b/timflow/transient/model.py
index fbe0fa4..70e1034 100644
--- a/timflow/transient/model.py
+++ b/timflow/transient/model.py
@@ -392,11 +392,18 @@ def headinvertical(self, x, y, z, t, aq=None, returneta=False):
                 rv[iz] = haq[lay]
             elif ltype == "l":
                 if lay == 0:
-                    eta = (
-                        headbar[:, lay]
-                        * np.sinh(aq.alpha[lay] * (aq.z[lay] - z[iz]))
-                        / np.sinh(aq.alpha[lay] * aq.Hll[lay])
-                    )
+                    if aq.topboundary == "sem":
+                        eta = (
+                            headbar[:, lay]
+                            * np.sinh(aq.alpha[lay] * (aq.z[lay] - z[iz]))
+                            / np.sinh(aq.alpha[lay] * aq.Hll[lay])
+                        )
+                    elif aq.topboundary == "lea":
+                        eta = (
+                            headbar[:, lay]
+                            * np.cosh(aq.alpha[lay] * (aq.z[lay] - z[iz]))
+                            / np.cosh(aq.alpha[lay] * aq.Hll[lay])
+                        )
                 else:
                     eta = (
                         headbar[:, lay - 1]

From 434087af4229546fe008d6917025e485f9ac7745 Mon Sep 17 00:00:00 2001
From: Mark Bakker <markbak@gmail.com>
Date: Mon, 23 Mar 2026 20:49:06 +0100
Subject: [PATCH 08/12] Update docs

---
 docs/transient/02examples/head_in_vertical.ipynb |  2 +-
 docs/transient/02examples/index.rst              | 14 +++++++++++++-
 2 files changed, 14 insertions(+), 2 deletions(-)

diff --git a/docs/transient/02examples/head_in_vertical.ipynb b/docs/transient/02examples/head_in_vertical.ipynb
index 3465c53..dcc6bbb 100755
--- a/docs/transient/02examples/head_in_vertical.ipynb
+++ b/docs/transient/02examples/head_in_vertical.ipynb
@@ -4,7 +4,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# Multi-layer example"
+    "# Head propagation in leaky layers"
    ]
   },
   {
diff --git a/docs/transient/02examples/index.rst b/docs/transient/02examples/index.rst
index 372a360..291b846 100644
--- a/docs/transient/02examples/index.rst
+++ b/docs/transient/02examples/index.rst
@@ -59,4 +59,16 @@ Pathlines
     :caption: Pathlines
     :hidden:
 
-    pathline_trace
\ No newline at end of file
+    pathline_trace
+
+Head in leaky layers
+--------------------
+
+- :doc:`head_in_vertical`
+
+.. toctree::
+    :maxdepth: 2
+    :caption: Head in leaky layers
+    :hidden:
+
+    head_in_vertical
\ No newline at end of file

From e3fab292b03d1dd5ae128978766a645791bf73d8 Mon Sep 17 00:00:00 2001
From: Mark Bakker <markbak@gmail.com>
Date: Mon, 23 Mar 2026 20:49:27 +0100
Subject: [PATCH 09/12] Start with head in leaky layers because of loading eff

---
 timflow/transient/model.py | 21 +++++++++++++++++++++
 1 file changed, 21 insertions(+)

diff --git a/timflow/transient/model.py b/timflow/transient/model.py
index 70e1034..59ea9a4 100644
--- a/timflow/transient/model.py
+++ b/timflow/transient/model.py
@@ -22,6 +22,7 @@
     compute_laplace_parameters_numba,
     invlap,
     invlapcomp,
+    invlapgen,
 )
 from timflow.transient.plots import PlotTransient
 
@@ -427,6 +428,26 @@ def headinvertical(self, x, y, z, t, aq=None, returneta=False):
                 )
         return rv
 
+    def headll(self, z, t, aq, returneta=False):
+        """Head in leaky layer caused by loading efficiency
+        """
+        z = np.atleast_1d(z)
+        t = np.atleast_1d(t)
+        rv = np.zeros((len(z), len(t)))
+        headbar = -1 / self.p
+        for iz in range(len(z)):
+            lay, ltype, _ = aq.findlayer(z[iz])
+            eta = (
+                headbar
+                * np.sinh(aq.alpha[lay] * (z[iz] - aq.zaqtop[lay]))
+                + headbar
+                * np.sinh(aq.alpha[lay] * (aq.zaqbot[lay - 1] - z[iz]))
+                ) / np.sinh(aq.alpha[lay] * aq.Hll[lay])
+            if returneta:
+                return eta
+            rv[iz] = invlapgen(t, eta, self.M, self.tintervals, np.array([0.0]), np.array([1.0]))
+        return rv
+
     def velocompold(self, x, y, z, t, aq=None, layer_ltype=[0, 0]):
         # implemented for one layer
         if aq is None:

From 1ec64b416ac36097384ee1a16b5ea0277762f31b Mon Sep 17 00:00:00 2001
From: Mark Bakker <markbak@gmail.com>
Date: Mon, 23 Mar 2026 20:52:49 +0100
Subject: [PATCH 10/12] ruff

---
 timflow/transient/model.py | 15 +++++++--------
 1 file changed, 7 insertions(+), 8 deletions(-)

diff --git a/timflow/transient/model.py b/timflow/transient/model.py
index 59ea9a4..1e6ab50 100644
--- a/timflow/transient/model.py
+++ b/timflow/transient/model.py
@@ -429,8 +429,7 @@ def headinvertical(self, x, y, z, t, aq=None, returneta=False):
         return rv
 
     def headll(self, z, t, aq, returneta=False):
-        """Head in leaky layer caused by loading efficiency
-        """
+        """Head in leaky layer caused by loading efficiency."""
         z = np.atleast_1d(z)
         t = np.atleast_1d(t)
         rv = np.zeros((len(z), len(t)))
@@ -438,14 +437,14 @@ def headll(self, z, t, aq, returneta=False):
         for iz in range(len(z)):
             lay, ltype, _ = aq.findlayer(z[iz])
             eta = (
-                headbar
-                * np.sinh(aq.alpha[lay] * (z[iz] - aq.zaqtop[lay]))
-                + headbar
-                * np.sinh(aq.alpha[lay] * (aq.zaqbot[lay - 1] - z[iz]))
-                ) / np.sinh(aq.alpha[lay] * aq.Hll[lay])
+                headbar * np.sinh(aq.alpha[lay] * (z[iz] - aq.zaqtop[lay]))
+                + headbar * np.sinh(aq.alpha[lay] * (aq.zaqbot[lay - 1] - z[iz]))
+            ) / np.sinh(aq.alpha[lay] * aq.Hll[lay])
             if returneta:
                 return eta
-            rv[iz] = invlapgen(t, eta, self.M, self.tintervals, np.array([0.0]), np.array([1.0]))
+            rv[iz] = invlapgen(
+                t, eta, self.M, self.tintervals, np.array([0.0]), np.array([1.0])
+            )
         return rv
 
     def velocompold(self, x, y, z, t, aq=None, layer_ltype=[0, 0]):

From a9473f18aeb83ce5904af7cdd04d43435e1d2b5c Mon Sep 17 00:00:00 2001
From: Mark Bakker <markbak@gmail.com>
Date: Tue, 24 Mar 2026 09:53:24 +0100
Subject: [PATCH 11/12] added leffll

---
 timflow/transient/aquifer.py            |  4 ++++
 timflow/transient/aquifer_parameters.py | 13 +++++++++--
 timflow/transient/circinhom.py          | 24 ++++++++++++++++----
 timflow/transient/inhom1d.py            | 29 ++++++++++++++++++++++---
 timflow/transient/model.py              | 29 ++++++++++++++++++++++---
 5 files changed, 87 insertions(+), 12 deletions(-)

diff --git a/timflow/transient/aquifer.py b/timflow/transient/aquifer.py
index d5906ac..c9a815c 100644
--- a/timflow/transient/aquifer.py
+++ b/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)
@@ -276,6 +278,7 @@ def __init__(
         Saq,
         Sll,
         leffaq,
+        leffll,
         poraq,
         porll,
         ltype,
@@ -294,6 +297,7 @@ def __init__(
             Saq,
             Sll,
             leffaq,
+            leffll,
             poraq,
             porll,
             ltype,
diff --git a/timflow/transient/aquifer_parameters.py b/timflow/transient/aquifer_parameters.py
index 4b7ea38..bcfcba1 100644
--- a/timflow/transient/aquifer_parameters.py
+++ b/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,6 +95,8 @@ 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)
@@ -100,13 +104,14 @@ def param_maq(
         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
diff --git a/timflow/transient/circinhom.py b/timflow/transient/circinhom.py
index 9a2bf4b..ff7ba48 100644
--- a/timflow/transient/circinhom.py
+++ b/timflow/transient/circinhom.py
@@ -64,8 +64,14 @@ def __init__(
         topboundary="imp",
         phreatictop=False,
     ):
-        kaq, Haq, Hll, c, Saq, Sll = param_maq(
-            kaq, z, c, Saq, Sll, topboundary, phreatictop
+        kaq, Haq, Hll, c, Saq, Sll, _, _, _, _, _ = param_maq(
+            kaq=kaq,
+            z=z,
+            c=c,
+            Saq=Saq,
+            Sll=Sll,
+            topboundary=topboundary,
+            phreatictop=phreatictop,
         )
         CircInhomData.__init__(
             self, model, x0, y0, R, kaq, Haq, c, Saq, Sll, topboundary, phreatictop
@@ -89,8 +95,18 @@ def __init__(
         topthick=0,
         topSll=0,
     ):
-        kaq, Haq, Hll, c, Saq, Sll = param_3d(
-            kaq, z, Saq, kzoverkh, phreatictop, topboundary, topres, topthick, topSll
+        kaq, Haq, Hll, c, Saq, Sll, _, _, _, _, _, z = param_3d(
+            kaq=kaq,
+            z=z,
+            Saq=Saq,
+            kzoverkh=kzoverkh,
+            poraq=0.3,
+            phreatictop=phreatictop,
+            topboundary=topboundary,
+            topres=topres,
+            topthick=topthick,
+            topSll=topSll,
+            toppor=0.3,
         )
         CircInhomData.__init__(self, model, x0, y0, R, kaq, Haq, c, Saq, Sll, "imp")
 
diff --git a/timflow/transient/inhom1d.py b/timflow/transient/inhom1d.py
index e36f01b..c241619 100644
--- a/timflow/transient/inhom1d.py
+++ b/timflow/transient/inhom1d.py
@@ -46,6 +46,8 @@ class Xsection(AquiferData):
         Specific storage of the leaky layers.
     leffaq : array
         loading efficiency of the aquifer
+    leffll : array
+        loading efficiency of the leaky layer
     poraq : array
         Porosities of the aquifers.
     porll : array
@@ -86,6 +88,7 @@ def __init__(
         Saq,
         Sll,
         leffaq,
+        leffll,
         poraq,
         porll,
         ltype,
@@ -107,6 +110,7 @@ def __init__(
             Saq,
             Sll,
             leffaq,
+            leffll,
             poraq,
             porll,
             ltype,
@@ -390,6 +394,8 @@ class XsectionMaq(Xsection):
         Specific storage of the leaky layers.
     leffaq : array
         loading efficiency of the aquifer
+    leffll : array
+        loading efficiency of the leaky layer
     poraq : array
         Porosities of the aquifers.
     porll : array
@@ -420,6 +426,7 @@ def __init__(
         Saq=0.001,
         Sll=0,
         leffaq=0,
+        leffll=0,
         poraq=0.3,
         porll=0.3,
         topboundary="conf",
@@ -428,8 +435,18 @@ def __init__(
         tsandN=None,
         name=None,
     ):
-        kaq, Haq, Hll, c, Saq, Sll, leffaq, poraq, porll, ltype = param_maq(
-            kaq, z, c, Saq, Sll, leffaq, poraq, porll, topboundary, phreatictop
+        kaq, Haq, Hll, c, Saq, Sll, leffaq, leffll, poraq, porll, ltype = param_maq(
+            kaq,
+            z,
+            c,
+            Saq,
+            Sll,
+            leffaq,
+            leffll,
+            poraq,
+            porll,
+            topboundary,
+            phreatictop,
         )
         super().__init__(
             model,
@@ -443,6 +460,7 @@ def __init__(
             Saq,
             Sll,
             leffaq,
+            leffll,
             poraq,
             porll,
             ltype,
@@ -479,6 +497,8 @@ class Xsection3D(Xsection):
         conductivity.
     leffaq : array
         Loading efficiency
+    leffll : array
+        loading efficiency of the leaky layer
     poraq : array
         Porosities of the aquifers.
     topboundary : str
@@ -514,6 +534,7 @@ def __init__(
         Saq=0.001,
         kzoverkh=0.1,
         leffaq=0,
+        leffll=0,
         poraq=0.3,
         topboundary="conf",
         phreatictop=False,
@@ -525,12 +546,13 @@ def __init__(
         tsandN=None,
         name=None,
     ):
-        kaq, Haq, Hll, c, Saq, Sll, leffaq, poraq, porll, ltype, z = param_3d(
+        kaq, Haq, Hll, c, Saq, Sll, leffaq, leffll, poraq, porll, ltype, z = param_3d(
             kaq,
             z,
             Saq,
             kzoverkh,
             leffaq,
+            leffll,
             poraq,
             phreatictop,
             topboundary,
@@ -551,6 +573,7 @@ def __init__(
             Saq,
             Sll,
             leffaq,
+            leffll,
             poraq,
             porll,
             ltype,
diff --git a/timflow/transient/model.py b/timflow/transient/model.py
index 1e6ab50..f665ccd 100644
--- a/timflow/transient/model.py
+++ b/timflow/transient/model.py
@@ -38,6 +38,7 @@ def __init__(
         Saq=[1e-4, 1e-4],
         Sll=[0],
         leffaq=0,
+        leffll=0,
         poraq=0.3,
         porll=0.3,
         ltype=["a", "a"],
@@ -71,6 +72,7 @@ def __init__(
             Saq,
             Sll,
             leffaq,
+            leffll,
             poraq,
             porll,
             ltype,
@@ -833,6 +835,9 @@ class ModelMaq(TimModel):
     leffaq : float, array or list
         loading efficiency of the aquifer
         only used when topboundary='semi' and hstar varies with time
+    leffll : float, array or list
+        loading efficiency of the leaky layer
+        only used when topboundary='semi' and hstar varies with time
     topboundary : string, 'conf' or 'semi' (default is 'conf')
         indicating whether the top is confined ('conf') or
         semi-confined ('semi')
@@ -862,6 +867,7 @@ def __init__(
         Saq=[0.001],
         Sll=[0],
         leffaq=0,
+        leffll=0,
         poraq=[0.3],
         porll=[0.3],
         topboundary="conf",
@@ -873,8 +879,18 @@ def __init__(
         steady=None,
     ):
         self.storeinput(inspect.currentframe())
-        kaq, Haq, Hll, c, Saq, Sll, leffaq, poraq, porll, ltype = param_maq(
-            kaq, z, c, Saq, Sll, leffaq, poraq, porll, topboundary, phreatictop
+        kaq, Haq, Hll, c, Saq, Sll, leffaq, leffll, poraq, porll, ltype = param_maq(
+            kaq,
+            z,
+            c,
+            Saq,
+            Sll,
+            leffaq,
+            leffll,
+            poraq,
+            porll,
+            topboundary,
+            phreatictop,
         )
         super().__init__(
             kaq,
@@ -885,6 +901,7 @@ def __init__(
             Saq,
             Sll,
             leffaq,
+            leffll,
             poraq,
             porll,
             ltype,
@@ -931,6 +948,9 @@ class Model3D(TimModel):
     leffaq : float, array or list
         loading efficiency of the aquifer
         only used when topboundary='semi' and hstar varies with time
+    leffll : float, array or list
+        loading efficiency of the leaky layer
+        only used when topboundary='semi' and hstar varies with time
     topboundary : string, 'conf' or 'semi' (default is 'conf')
         indicating whether the top is confined ('conf') or
         semi-confined ('semi').
@@ -964,6 +984,7 @@ def __init__(
         Saq=0.001,
         kzoverkh=0.1,
         leffaq=0,
+        leffll=0,
         poraq=0.3,
         topboundary="conf",
         phreatictop=False,
@@ -979,12 +1000,13 @@ def __init__(
     ):
         """Z must have the length of the number of layers + 1."""
         self.storeinput(inspect.currentframe())
-        kaq, Haq, Hll, c, Saq, Sll, leffaq, poraq, porll, ltype, z = param_3d(
+        kaq, Haq, Hll, c, Saq, Sll, leffaq, leffll, poraq, porll, ltype, z = param_3d(
             kaq,
             z,
             Saq,
             kzoverkh,
             leffaq,
+            leffll,
             poraq,
             phreatictop,
             topboundary,
@@ -1002,6 +1024,7 @@ def __init__(
             Saq,
             Sll,
             leffaq,
+            leffll,
             poraq,
             porll,
             ltype,

From 0ada6581308f5840eab922a5bbb5f1988a460c33 Mon Sep 17 00:00:00 2001
From: Mark Bakker <markbak@gmail.com>
Date: Wed, 25 Mar 2026 10:10:27 +0100
Subject: [PATCH 12/12] added effect in top clay layer

add effect in top clay layer of head variation on top of semi-confined aquifer
---
 timflow/transient/aquifer.py | 36 ++++++++++++++++++++++++++++++++++++
 timflow/transient/model.py   |  2 ++
 2 files changed, 38 insertions(+)

diff --git a/timflow/transient/aquifer.py b/timflow/transient/aquifer.py
index c9a815c..2b01fd6 100644
--- a/timflow/transient/aquifer.py
+++ b/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__(
@@ -265,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__(
diff --git a/timflow/transient/model.py b/timflow/transient/model.py
index f665ccd..c3184be 100644
--- a/timflow/transient/model.py
+++ b/timflow/transient/model.py
@@ -428,6 +428,8 @@ def headinvertical(self, x, y, z, t, aq=None, returneta=False):
                     self.ebc,
                     1,
                 )
+        if aq.topboundary[:3] == "sem":
+            rv += aq.headsemitoplayer(x, y, z, t)
         return rv
 
     def headll(self, z, t, aq, returneta=False):