diff --git a/docs/transient/02examples/head_in_vertical.ipynb b/docs/transient/02examples/head_in_vertical.ipynb
new file mode 100755
index 0000000..dcc6bbb
--- /dev/null
+++ b/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
+}
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
diff --git a/timflow/transient/aquifer.py b/timflow/transient/aquifer.py
index 87acc8d..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__(
@@ -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,
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 f0664f4..c3184be 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
 
@@ -37,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"],
@@ -70,6 +72,7 @@ def __init__(
             Saq,
             Sll,
             leffaq,
+            leffll,
             poraq,
             porll,
             ltype,
@@ -352,6 +355,102 @@ def head(self, x, y, t, layers=None, aq=None, derivative=0, neglect_steady=False
                 h += htimml[:, np.newaxis]
         return h
 
+    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, list, or array
+        t : float, list, or array
+
+        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)
+        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,
+            aq.naq,
+        )
+        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:
+                    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]
+                        * 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,
+                )
+        if aq.topboundary[:3] == "sem":
+            rv += aq.headsemitoplayer(x, y, z, t)
+        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:
@@ -738,6 +837,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')
@@ -767,6 +869,7 @@ def __init__(
         Saq=[0.001],
         Sll=[0],
         leffaq=0,
+        leffll=0,
         poraq=[0.3],
         porll=[0.3],
         topboundary="conf",
@@ -778,8 +881,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,
@@ -790,6 +903,7 @@ def __init__(
             Saq,
             Sll,
             leffaq,
+            leffll,
             poraq,
             porll,
             ltype,
@@ -836,6 +950,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').
@@ -869,6 +986,7 @@ def __init__(
         Saq=0.001,
         kzoverkh=0.1,
         leffaq=0,
+        leffll=0,
         poraq=0.3,
         topboundary="conf",
         phreatictop=False,
@@ -884,12 +1002,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,
@@ -907,6 +1026,7 @@ def __init__(
             Saq,
             Sll,
             leffaq,
+            leffll,
             poraq,
             porll,
             ltype,