diff --git a/docs/transient/05benchmarks/index.rst b/docs/transient/05benchmarks/index.rst
index 4dd43d9f..19363bc5 100644
--- a/docs/transient/05benchmarks/index.rst
+++ b/docs/transient/05benchmarks/index.rst
@@ -18,6 +18,7 @@ timflow solutions or understand the numerical accuracy of specific features.
 - :doc:`synthetic_test_calibrate`
 - :doc:`synthetic_calibrate_2aquifers`
 - :doc:`river1d`
+- :doc:`river1d_2layer_mf6.ipynb`
 
 .. toctree::
     :maxdepth: 3
@@ -34,4 +35,5 @@ timflow solutions or understand the numerical accuracy of specific features.
     validation_tidal_wave_with_Bruggeman
     synthetic_test_calibrate
     synthetic_calibrate_2aquifers
-    river1d
\ No newline at end of file
+    river1d
+    river1d_2layer_mf6.ipynb
\ No newline at end of file
diff --git a/timflow/plots/plots.py b/timflow/plots/plots.py
index 569fda8f..56235c88 100644
--- a/timflow/plots/plots.py
+++ b/timflow/plots/plots.py
@@ -652,7 +652,11 @@ def _xsection_leaky_layer_params(
             # Transient: resistance c and storage Sll
             ssfmt = ".2e"
             cstr = f"$c$ = {self._ml.aq.c[lli]:{fmt}}"
-            sstr = f"$S_s$ = {self._ml.aq.Sll[lli]:{ssfmt}}"
+            Slli = self._ml.aq.Sll[lli]
+            if Slli > 1e-20:
+                sstr = f"$S_s$ = {Slli:{ssfmt}}"
+            else:
+                sstr = "$S_s$ = 0.0"
             if units is not None:
                 c_unitstr = f" {units['c']}" if "c" in units else ""
                 # Prefer Sll unit; fall back to Saq for compatibility.
@@ -661,6 +665,8 @@ def _xsection_leaky_layer_params(
                 c_unitstr = ""
                 ss_unitstr = ""
             paramtxt = cstr + c_unitstr + sep + sstr + ss_unitstr
+            if hasattr(self._ml.aq, "leffll") and self._ml.aq.leffll[lli] != 0:
+                paramtxt += f"{sep}$\\beta$ = {self._ml.aq.leffll[lli]:{fmt}}"
 
         ax.text(
             r0 + 0.75 * r if labels else r0 + 0.5 * r,
@@ -731,6 +737,8 @@ def _xsection_aquifer_params(
                 paramtxt += f"{sep}$S$ = {self._ml.aq.Saq[aqi]:{fmt}}"
             else:
                 paramtxt += f"{sep}$S_s$ = {self._ml.aq.Saq[aqi]:{ssfmt}}" + ss_unitstr
+            if hasattr(self._ml.aq, "leffaq") and self._ml.aq.leffaq[aqi] != 0:
+                paramtxt += f"{sep}$\\beta$ = {self._ml.aq.leffaq[aqi]:{fmt}}"
 
         ax.text(
             r0 + 0.75 * r if labels else r0 + 0.5 * r,
diff --git a/timflow/transient/aquifer_parameters.py b/timflow/transient/aquifer_parameters.py
index 9ccd0e09..7b0c67aa 100644
--- a/timflow/transient/aquifer_parameters.py
+++ b/timflow/transient/aquifer_parameters.py
@@ -50,6 +50,8 @@ def param_maq(
             Sll = Sll * np.ones(naq - 1)
         if len(porll) == 1:
             porll = porll * np.ones(naq - 1)
+        leffaq = np.zeros(naq)  # always 0.0, regardless of user input
+        leffll = np.zeros(naq)  # always 0.0, regardless of user input
         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)
diff --git a/timflow/transient/inhom1d.py b/timflow/transient/inhom1d.py
index b58fe5e3..3942a477 100644
--- a/timflow/transient/inhom1d.py
+++ b/timflow/transient/inhom1d.py
@@ -47,8 +47,10 @@ class Xsection(AquiferData):
         Specific storage of the leaky layers.
     leffaq : array
         loading efficiency of the aquifer
+        only used when topboundary='semi' and hstar varies with time
     leffll : array
         loading efficiency of the leaky layer
+        only used when topboundary='semi' and hstar varies with time
     poraq : array
         Porosities of the aquifers.
     porll : array
@@ -338,13 +340,19 @@ def plot(
                     )
                 if params:
                     cstr = f"$c$ = {self.c[lli]:{fmt}}"
-                    sstr = f"$S_s$ = {self.Sll[lli]:{ssfmt}}"
+                    Slli = self.Sll[lli]
+                    if Slli > 1e-20:
+                        sstr = f"$S_s$ = {Slli:{ssfmt}}"
+                    else:
+                        sstr = "$S_s$ = 0.0"
                     cstr_with_unit = cstr + c_unitstr
                     sstr_with_unit = sstr + ss_unitstr
                     if sep == "\n":
                         paramtxt = cstr_with_unit + sep + sstr_with_unit
                     else:
                         paramtxt = cstr_with_unit + sep + sstr_with_unit
+                    if self.leffll[lli] != 0.0:
+                        paramtxt += f"{sep}$\\beta$ = {self.leffll[lli]:{fmt}}"
                     ax.text(
                         r0 + 0.75 * r if labels else r0 + 0.5 * r,
                         np.mean(self.z[i : i + 2]),
@@ -380,6 +388,8 @@ def plot(
                     paramtxt = khstr + kh_unitstr + "\n" + sstr + ss_unitstr
                 else:
                     paramtxt = khstr + kh_unitstr + sep + sstr + ss_unitstr
+                if self.leffaq[aqi] != 0.0:
+                    paramtxt += f"{sep}$\\beta$ = {self.leffaq[aqi]:{fmt}}"
                 ax.text(
                     r0 + 0.75 * r if labels else r0 + 0.5 * r,
                     np.mean(self.z[i : i + 2]),
@@ -429,8 +439,10 @@ class XsectionMaq(Xsection):
         Specific storage of the leaky layers.
     leffaq : array
         loading efficiency of the aquifer
+        only used when topboundary='semi' and hstar varies with time
     leffll : array
         loading efficiency of the leaky layer
+        only used when topboundary='semi' and hstar varies with time
     poraq : array
         Porosities of the aquifers.
     porll : array
@@ -531,9 +543,11 @@ class Xsection3D(Xsection):
         Ratio of vertical hydraulic conductivity to horizontal hydraulic
         conductivity.
     leffaq : array
-        Loading efficiency
+        loading efficiency of the aquifer
+        only used when topboundary='semi' and hstar varies with time
     leffll : array
         loading efficiency of the leaky layer
+        only used when topboundary='semi' and hstar varies with time
     poraq : array
         Porosities of the aquifers.
     topboundary : string, 'confined', 'phreatic', or 'semi' (default is 'conf')