From efbaed1f4e5b71b59f20e9cf0f09f9ad2e44d2b8 Mon Sep 17 00:00:00 2001
From: Claire Donnelly <claire.donnelly@deltares.nl>
Date: Thu, 23 Jul 2026 15:09:40 +0200
Subject: [PATCH 1/4] Create new user guide for assigning wells

---
 ...ured-grids.py => 04-unstructured-grids.py} |   0
 ...zy-evaluation.py => 05-lazy-evaluation.py} |   0
 ...etization.py => 06-time-discretization.py} |   0
 .../{08-regridding.py => 07-regridding.py}    |   0
 examples/user-guide/08-assign-wells.py        | 201 ++++++++++++++++++
 5 files changed, 201 insertions(+)
 rename examples/user-guide/{05-unstructured-grids.py => 04-unstructured-grids.py} (100%)
 rename examples/user-guide/{06-lazy-evaluation.py => 05-lazy-evaluation.py} (100%)
 rename examples/user-guide/{07-time-discretization.py => 06-time-discretization.py} (100%)
 rename examples/user-guide/{08-regridding.py => 07-regridding.py} (100%)
 create mode 100644 examples/user-guide/08-assign-wells.py

diff --git a/examples/user-guide/05-unstructured-grids.py b/examples/user-guide/04-unstructured-grids.py
similarity index 100%
rename from examples/user-guide/05-unstructured-grids.py
rename to examples/user-guide/04-unstructured-grids.py
diff --git a/examples/user-guide/06-lazy-evaluation.py b/examples/user-guide/05-lazy-evaluation.py
similarity index 100%
rename from examples/user-guide/06-lazy-evaluation.py
rename to examples/user-guide/05-lazy-evaluation.py
diff --git a/examples/user-guide/07-time-discretization.py b/examples/user-guide/06-time-discretization.py
similarity index 100%
rename from examples/user-guide/07-time-discretization.py
rename to examples/user-guide/06-time-discretization.py
diff --git a/examples/user-guide/08-regridding.py b/examples/user-guide/07-regridding.py
similarity index 100%
rename from examples/user-guide/08-regridding.py
rename to examples/user-guide/07-regridding.py
diff --git a/examples/user-guide/08-assign-wells.py b/examples/user-guide/08-assign-wells.py
new file mode 100644
index 000000000..3e9caf832
--- /dev/null
+++ b/examples/user-guide/08-assign-wells.py
@@ -0,0 +1,201 @@
+"""
+Assigning Wells to Model Layers
+==========
+
+iMOD Python provides two grid-agnostic well classes:
+:class:`imod.mf6.Well` and :class:`imod.mf6.LayeredWell` 
+to build MODFLOW 6 well package input.
+
+Use :class:`imod.mf6.Well` when the physical top and bottom of a well screen
+are known. During conversion, iMOD Python intersects each screen with model
+layers and distributes the specified rate over eligible cells in proportion to
+transmissivity.
+
+Use :class:`imod.mf6.LayeredWell` when the target model layer is already known
+for every well record. In that case, the supplied layer and rate are kept as
+provided.
+
+In both cases, wells in inactive cells are removed during conversion to a 
+MODFLOW 6 well package.
+
+"""
+
+# sphinx_gallery_thumbnail_number = -1
+
+# %%
+# Example data
+# ------------
+#
+# Let's load the data first. We have a layer model containing a basic
+# hydrogeological schemitization of our model, so the tops and bottoms of model
+# layers, the hydraulic conductivity (k), and which cells are active (idomain=1)
+# or vertical passthrough (idomain=-1).
+
+import imod
+
+layer_model = imod.data.hondsrug_layermodel_topsystem()
+
+layer_model
+
+# %%
+
+# %%
+#
+# Let's extract the layer model data into separate variables for convenience.
+idomain = layer_model["idomain"]
+top = layer_model["top"]
+bottom = layer_model["bottom"]
+k = layer_model["k"]
+# %%
+
+# %%
+# 
+# Let's define a cross-section line through the model, and some well locations along that line.
+from shapely.geometry import LineString
+geometry = LineString([[238725, 560000], [242000, 563500]])
+
+x=[239380.0, 240362.5, 241345.0]
+y=[560700.0, 561750.0, 562800.0]
+rate=[-10.0, -25.0, -15.0]
+
+# %%
+
+# %%
+# Create a Well object
+# ------------
+#
+# Now that we have the model data and well data, we can create a :class:`imod.mf6.Well` object
+# and convert it to a MODFLOW 6 well package input.
+
+screen_based = imod.mf6.Well(
+       x=x,
+       y=y,
+       screen_top=[6.0, 7.0, 6.0],
+       screen_bottom=[5.0, 6.5, 4.5],
+       rate=rate,
+)
+screen_based_mf6 = screen_based.to_mf6_pkg(idomain, top, bottom, k)
+
+screen_based_mf6["cellid"]
+
+# %%
+
+# %%
+# 
+# Let's plot the top elevation of the model on a map, with the well locations and cross-section overlaid. 
+# You can see we have a ridge roughly the centre of the model, sided by two low-lying areas.
+
+import numpy as np
+import geopandas as gpd
+
+overlays = [
+    {"gdf": gpd.GeoDataFrame(geometry=[geometry]), "edgecolor": "black", "linewidth": 3}
+]
+
+fig, ax = imod.visualize.plot_map(
+    layer_model["top"].sel(layer=1), "viridis", np.linspace(1, 20, 11), overlays
+)
+
+ax.scatter(screen_based.x, screen_based.y, c="red", s=60, marker="o")
+
+# %%
+
+# %%
+# 
+# We can also visualise the well location per model layer, with respect to the hydraulic conductivity.
+
+from matplotlib import pyplot as plt
+
+cellid = screen_based_mf6["cellid"]
+screen_layers = cellid.sel(dim_cellid="layer").values.astype(int)
+well_x = cellid["x"].values
+well_y = cellid["y"].values
+
+unique_layers = np.unique(screen_layers)
+color_levels = np.linspace(float(k.min()), float(k.max()), 11)
+
+fig, axes = plt.subplots(
+    3,
+    2,
+    figsize=(12, 7),
+    constrained_layout=True,
+)
+fig.suptitle("Well locations vs hydraulic conductivity", fontsize=16)
+axes = axes.ravel()
+
+for ax, layer in zip(axes, unique_layers):
+    layer_mask = screen_layers == layer
+
+    imod.visualize.plot_map(
+        k.sel(layer=int(layer)),
+        "viridis",
+        color_levels,
+        fig=fig,
+        ax=ax,
+    )
+    ax.scatter(well_x[layer_mask], well_y[layer_mask], c="red", s=60, marker="o")
+    ax.set_title(f"Model layer {int(layer)}")
+
+for ax in axes[len(unique_layers):]:
+    ax.set_visible(False)
+
+# %%
+
+# %%
+# Create a Layered Well object
+# ------------
+#
+# Now we can follow a similar process to create a :class:`imod.mf6.LayeredWell` object
+# and convert it to a MODFLOW 6 well package input. The process is the similar, except we 
+# specify the target model layer for each well.
+layer_based = imod.mf6.LayeredWell(
+    x=x,
+    y=y,
+    layer=[6, 7, 6],
+    rate=rate,
+)
+layer_based_mf6 = layer_based.to_mf6_pkg(idomain, top, bottom, k)
+
+layer_based_mf6["cellid"]
+
+# %%
+
+# %%
+# To visualise the difference between the two well types, 
+# we can plot the wells on a cross-section of the model. 
+import xarray as xr
+from shapely.geometry import Point
+from matplotlib import pyplot as plt
+
+layer_grid = layer_model.layer * xr.ones_like(layer_model["top"])
+layer_grid.coords["top"] = layer_model["top"]
+layer_grid.coords["bottom"] = layer_model["bottom"]
+xsection_layer_nr = imod.select.cross_section_linestring(layer_grid, geometry)
+
+fig, axes = plt.subplots(1, 2, figsize=(14, 6), constrained_layout=True)
+
+# Well
+well_df = screen_based.dataset[["x", "y", "screen_top", "screen_bottom"]].to_dataframe().reset_index(drop=True)
+well_df["s"] = [geometry.project(Point(x, y)) for x, y in zip(well_df["x"], well_df["y"])]
+
+imod.visualize.cross_section(xsection_layer_nr, "tab20", np.arange(21), fig=fig, ax=axes[0])
+for _, row in well_df.iterrows():
+    axes[0].vlines(
+        row["s"],
+        row["screen_bottom"],
+        row["screen_top"],
+        color="black",
+        linewidth=3,
+    )
+    axes[0].scatter(row["s"], row["screen_top"], color="black", marker="1", s=100, linewidths=1.6)
+    axes[0].scatter(row["s"], row["screen_bottom"], color="black", marker="2", s=100, linewidths=1.6)
+axes[0].set_title("Well on layer cross section", fontsize=16)
+
+# LayeredWell
+layered_well_df = layer_based.dataset[["x", "y", "layer"]].to_dataframe().reset_index(drop=True)
+layered_well_df["s"] = [geometry.project(Point(x, y)) for x, y in zip(layered_well_df["x"], layered_well_df["y"])]
+
+imod.visualize.cross_section(xsection_layer_nr, "tab20", np.arange(21), fig=fig, ax=axes[1])
+for _, row in layered_well_df.iterrows():
+    axes[1].scatter(row["s"], row["layer"], color="black", marker="x", s=120, linewidths=1.8)
+axes[1].set_title("LayeredWell on layer cross section", fontsize=16)

From 20f034ba015e94a1ebade2dce19281099d332902 Mon Sep 17 00:00:00 2001
From: Claire Donnelly <claire.donnelly@deltares.nl>
Date: Thu, 23 Jul 2026 15:10:15 +0200
Subject: [PATCH 2/4] Better explain the difference between Well and
 LayeredWell

---
 imod/mf6/wel.py | 31 ++++++++++++++++++++++++++++---
 1 file changed, 28 insertions(+), 3 deletions(-)

diff --git a/imod/mf6/wel.py b/imod/mf6/wel.py
index e39af906a..04a96d900 100644
--- a/imod/mf6/wel.py
+++ b/imod/mf6/wel.py
@@ -446,9 +446,17 @@ def to_mf6_pkg(
         Write package to Modflow 6 package.
 
         Based on the model grid and top and bottoms, cellids are determined.
-        When well screens hit multiple layers, groundwater extractions are
-        distributed based on layer transmissivities. Wells located in inactive
-        cells are removed.
+
+        For :class:`Well`, the screen interval is split over
+        intersected cells, and the rate is distributed by layer
+        transmissivity when multiple layers are intersected. Wells located in
+        inactive cells are removed.
+
+        For :class:`LayeredWell`, the provided layer and rate are
+        used directly, without screen splitting, screen-depth allocation, or
+        transmissivity-based rate distribution. If the provided layer is
+        inactive, the well is removed.
+        
 
         Note
         ----
@@ -768,6 +776,15 @@ class Well(GridAgnosticWell):
     """
     Agnostic WEL package, which accepts x, y and a top and bottom of the well screens.
 
+    Use :class:`Well` when input is defined by screen elevations rather than
+    by model layer. During conversion to a MODFLOW 6 package, the screen
+    interval is intersected with model layers and the specified rate is
+    distributed over eligible cells in proportion to transmissivity (based on
+    horizontal hydraulic conductivity and screen-overlap thickness).
+
+    Use :class:`LayeredWell` when the target layer is already known and each
+    rate must remain assigned to that layer.
+
     This package can be written to any provided model grid.
     Any number of WEL Packages can be specified for a single groundwater flow model.
     https://water.usgs.gov/water-resources/software/MODFLOW-6/mf6io_6.0.4.pdf#page=63
@@ -1169,6 +1186,14 @@ class LayeredWell(GridAgnosticWell):
     """
     Agnostic WEL package, which accepts x, y and layers.
 
+    Use :class:`LayeredWell` when input is defined by model layer. During
+    conversion to a MODFLOW 6 package, the supplied layer and rate are
+    retained; no layer inference from screen elevations or
+    transmissivity-based rate redistribution is applied.
+
+    Use :class:`Well` instead when input is defined by a physical well screen
+    that can intersect one or more model layers.
+
     This package can be written to any provided model grid, given that it has
     enough layers. Any number of WEL Packages can be specified for a single
     groundwater flow model.

From ab83906d261b6534434c1acb41276d0c9d2bd263 Mon Sep 17 00:00:00 2001
From: Claire Donnelly <claire.donnelly@deltares.nl>
Date: Thu, 23 Jul 2026 15:11:34 +0200
Subject: [PATCH 3/4] Fix formatting and ignore mydask.png

---
 .gitignore                             |  2 +-
 examples/user-guide/08-assign-wells.py | 80 ++++++++++++++++----------
 imod/mf6/wel.py                        |  2 +-
 3 files changed, 53 insertions(+), 31 deletions(-)

diff --git a/.gitignore b/.gitignore
index c8ef7b959..4a2412a13 100644
--- a/.gitignore
+++ b/.gitignore
@@ -142,6 +142,6 @@ examples/data
 # pixi environments
 .pixi
 
-/imod/tests/mydask.png
+**/mydask.png
 /imod/tests/*_report.xml
 docs/sg_execution_times.rst
diff --git a/examples/user-guide/08-assign-wells.py b/examples/user-guide/08-assign-wells.py
index 3e9caf832..b1f725529 100644
--- a/examples/user-guide/08-assign-wells.py
+++ b/examples/user-guide/08-assign-wells.py
@@ -3,7 +3,7 @@
 ==========
 
 iMOD Python provides two grid-agnostic well classes:
-:class:`imod.mf6.Well` and :class:`imod.mf6.LayeredWell` 
+:class:`imod.mf6.Well` and :class:`imod.mf6.LayeredWell`
 to build MODFLOW 6 well package input.
 
 Use :class:`imod.mf6.Well` when the physical top and bottom of a well screen
@@ -15,7 +15,7 @@
 for every well record. In that case, the supplied layer and rate are kept as
 provided.
 
-In both cases, wells in inactive cells are removed during conversion to a 
+In both cases, wells in inactive cells are removed during conversion to a
 MODFLOW 6 well package.
 
 """
@@ -49,14 +49,15 @@
 # %%
 
 # %%
-# 
+#
 # Let's define a cross-section line through the model, and some well locations along that line.
 from shapely.geometry import LineString
+
 geometry = LineString([[238725, 560000], [242000, 563500]])
 
-x=[239380.0, 240362.5, 241345.0]
-y=[560700.0, 561750.0, 562800.0]
-rate=[-10.0, -25.0, -15.0]
+x = [239380.0, 240362.5, 241345.0]
+y = [560700.0, 561750.0, 562800.0]
+rate = [-10.0, -25.0, -15.0]
 
 # %%
 
@@ -68,11 +69,11 @@
 # and convert it to a MODFLOW 6 well package input.
 
 screen_based = imod.mf6.Well(
-       x=x,
-       y=y,
-       screen_top=[6.0, 7.0, 6.0],
-       screen_bottom=[5.0, 6.5, 4.5],
-       rate=rate,
+    x=x,
+    y=y,
+    screen_top=[6.0, 7.0, 6.0],
+    screen_bottom=[5.0, 6.5, 4.5],
+    rate=rate,
 )
 screen_based_mf6 = screen_based.to_mf6_pkg(idomain, top, bottom, k)
 
@@ -81,12 +82,12 @@
 # %%
 
 # %%
-# 
-# Let's plot the top elevation of the model on a map, with the well locations and cross-section overlaid. 
+#
+# Let's plot the top elevation of the model on a map, with the well locations and cross-section overlaid.
 # You can see we have a ridge roughly the centre of the model, sided by two low-lying areas.
 
-import numpy as np
 import geopandas as gpd
+import numpy as np
 
 overlays = [
     {"gdf": gpd.GeoDataFrame(geometry=[geometry]), "edgecolor": "black", "linewidth": 3}
@@ -101,7 +102,7 @@
 # %%
 
 # %%
-# 
+#
 # We can also visualise the well location per model layer, with respect to the hydraulic conductivity.
 
 from matplotlib import pyplot as plt
@@ -136,7 +137,7 @@
     ax.scatter(well_x[layer_mask], well_y[layer_mask], c="red", s=60, marker="o")
     ax.set_title(f"Model layer {int(layer)}")
 
-for ax in axes[len(unique_layers):]:
+for ax in axes[len(unique_layers) :]:
     ax.set_visible(False)
 
 # %%
@@ -146,7 +147,7 @@
 # ------------
 #
 # Now we can follow a similar process to create a :class:`imod.mf6.LayeredWell` object
-# and convert it to a MODFLOW 6 well package input. The process is the similar, except we 
+# and convert it to a MODFLOW 6 well package input. The process is the similar, except we
 # specify the target model layer for each well.
 layer_based = imod.mf6.LayeredWell(
     x=x,
@@ -161,11 +162,11 @@
 # %%
 
 # %%
-# To visualise the difference between the two well types, 
-# we can plot the wells on a cross-section of the model. 
+# To visualise the difference between the two well types,
+# we can plot the wells on a cross-section of the model.
 import xarray as xr
-from shapely.geometry import Point
 from matplotlib import pyplot as plt
+from shapely.geometry import Point
 
 layer_grid = layer_model.layer * xr.ones_like(layer_model["top"])
 layer_grid.coords["top"] = layer_model["top"]
@@ -175,10 +176,18 @@
 fig, axes = plt.subplots(1, 2, figsize=(14, 6), constrained_layout=True)
 
 # Well
-well_df = screen_based.dataset[["x", "y", "screen_top", "screen_bottom"]].to_dataframe().reset_index(drop=True)
-well_df["s"] = [geometry.project(Point(x, y)) for x, y in zip(well_df["x"], well_df["y"])]
+well_df = (
+    screen_based.dataset[["x", "y", "screen_top", "screen_bottom"]]
+    .to_dataframe()
+    .reset_index(drop=True)
+)
+well_df["s"] = [
+    geometry.project(Point(x, y)) for x, y in zip(well_df["x"], well_df["y"])
+]
 
-imod.visualize.cross_section(xsection_layer_nr, "tab20", np.arange(21), fig=fig, ax=axes[0])
+imod.visualize.cross_section(
+    xsection_layer_nr, "tab20", np.arange(21), fig=fig, ax=axes[0]
+)
 for _, row in well_df.iterrows():
     axes[0].vlines(
         row["s"],
@@ -187,15 +196,28 @@
         color="black",
         linewidth=3,
     )
-    axes[0].scatter(row["s"], row["screen_top"], color="black", marker="1", s=100, linewidths=1.6)
-    axes[0].scatter(row["s"], row["screen_bottom"], color="black", marker="2", s=100, linewidths=1.6)
+    axes[0].scatter(
+        row["s"], row["screen_top"], color="black", marker="1", s=100, linewidths=1.6
+    )
+    axes[0].scatter(
+        row["s"], row["screen_bottom"], color="black", marker="2", s=100, linewidths=1.6
+    )
 axes[0].set_title("Well on layer cross section", fontsize=16)
 
 # LayeredWell
-layered_well_df = layer_based.dataset[["x", "y", "layer"]].to_dataframe().reset_index(drop=True)
-layered_well_df["s"] = [geometry.project(Point(x, y)) for x, y in zip(layered_well_df["x"], layered_well_df["y"])]
+layered_well_df = (
+    layer_based.dataset[["x", "y", "layer"]].to_dataframe().reset_index(drop=True)
+)
+layered_well_df["s"] = [
+    geometry.project(Point(x, y))
+    for x, y in zip(layered_well_df["x"], layered_well_df["y"])
+]
 
-imod.visualize.cross_section(xsection_layer_nr, "tab20", np.arange(21), fig=fig, ax=axes[1])
+imod.visualize.cross_section(
+    xsection_layer_nr, "tab20", np.arange(21), fig=fig, ax=axes[1]
+)
 for _, row in layered_well_df.iterrows():
-    axes[1].scatter(row["s"], row["layer"], color="black", marker="x", s=120, linewidths=1.8)
+    axes[1].scatter(
+        row["s"], row["layer"], color="black", marker="x", s=120, linewidths=1.8
+    )
 axes[1].set_title("LayeredWell on layer cross section", fontsize=16)
diff --git a/imod/mf6/wel.py b/imod/mf6/wel.py
index 04a96d900..2bacf9056 100644
--- a/imod/mf6/wel.py
+++ b/imod/mf6/wel.py
@@ -456,7 +456,7 @@ def to_mf6_pkg(
         used directly, without screen splitting, screen-depth allocation, or
         transmissivity-based rate distribution. If the provided layer is
         inactive, the well is removed.
-        
+
 
         Note
         ----

From 2203b66bfb62917c1f8348697285b7e08d9bd6c0 Mon Sep 17 00:00:00 2001
From: Claire Donnelly <claire.donnelly@deltares.nl>
Date: Tue, 28 Jul 2026 14:04:16 +0200
Subject: [PATCH 4/4] Make user guide easier to follow and apply other review
 suggestions

---
 examples/user-guide/08-assign-wells.py | 67 ++++++++++++++++++--------
 imod/mf6/wel.py                        |  5 +-
 2 files changed, 52 insertions(+), 20 deletions(-)

diff --git a/examples/user-guide/08-assign-wells.py b/examples/user-guide/08-assign-wells.py
index b1f725529..72034b7b5 100644
--- a/examples/user-guide/08-assign-wells.py
+++ b/examples/user-guide/08-assign-wells.py
@@ -11,9 +11,8 @@
 layers and distributes the specified rate over eligible cells in proportion to
 transmissivity.
 
-Use :class:`imod.mf6.LayeredWell` when the target model layer is already known
-for every well record. In that case, the supplied layer and rate are kept as
-provided.
+Use :class:`imod.mf6.LayeredWell` for direct control over well assignment to model layers.
+In that case, the supplied layer and rate are kept as provided.
 
 In both cases, wells in inactive cells are removed during conversion to a
 MODFLOW 6 well package.
@@ -50,15 +49,15 @@
 
 # %%
 #
-# Let's define a cross-section line through the model, and some well locations along that line.
+# Let's define some well locations, and then draw a cross-section line through them.
 from shapely.geometry import LineString
 
-geometry = LineString([[238725, 560000], [242000, 563500]])
-
 x = [239380.0, 240362.5, 241345.0]
 y = [560700.0, 561750.0, 562800.0]
 rate = [-10.0, -25.0, -15.0]
 
+geometry = LineString([[238725, 560000], [242000, 563500]])
+
 # %%
 
 # %%
@@ -68,11 +67,15 @@
 # Now that we have the model data and well data, we can create a :class:`imod.mf6.Well` object
 # and convert it to a MODFLOW 6 well package input.
 
+# Define the top and bottom elevations of the well screen for each well.
+screen_top = [6.0, 7.0, 6.0]
+screen_bottom = [5.0, 6.5, 4.5]
+
 screen_based = imod.mf6.Well(
     x=x,
     y=y,
-    screen_top=[6.0, 7.0, 6.0],
-    screen_bottom=[5.0, 6.5, 4.5],
+    screen_top=screen_top,
+    screen_bottom=screen_bottom,
     rate=rate,
 )
 screen_based_mf6 = screen_based.to_mf6_pkg(idomain, top, bottom, k)
@@ -149,10 +152,14 @@
 # Now we can follow a similar process to create a :class:`imod.mf6.LayeredWell` object
 # and convert it to a MODFLOW 6 well package input. The process is the similar, except we
 # specify the target model layer for each well.
+
+# Assign the wells to model layers directly, instead of using screen top and bottom elevations.
+layer = [6, 7, 6]
+
 layer_based = imod.mf6.LayeredWell(
     x=x,
     y=y,
-    layer=[6, 7, 6],
+    layer=layer,
     rate=rate,
 )
 layer_based_mf6 = layer_based.to_mf6_pkg(idomain, top, bottom, k)
@@ -168,56 +175,78 @@
 from matplotlib import pyplot as plt
 from shapely.geometry import Point
 
+# Create a grid containing layer numbers and add top/bottom elevations as coordinates
 layer_grid = layer_model.layer * xr.ones_like(layer_model["top"])
 layer_grid.coords["top"] = layer_model["top"]
 layer_grid.coords["bottom"] = layer_model["bottom"]
-xsection_layer_nr = imod.select.cross_section_linestring(layer_grid, geometry)
 
-fig, axes = plt.subplots(1, 2, figsize=(14, 6), constrained_layout=True)
+# Extract a cross-section along the specified geometry line
+xsection_layer_nr = imod.select.cross_section_linestring(layer_grid, geometry)
 
-# Well
+# Prepare the screen-based well data for visualization
 well_df = (
     screen_based.dataset[["x", "y", "screen_top", "screen_bottom"]]
     .to_dataframe()
     .reset_index(drop=True)
 )
-well_df["s"] = [
+# Project well locations onto the cross-section line to get their position along the line
+well_df["position_along_line"] = [
     geometry.project(Point(x, y)) for x, y in zip(well_df["x"], well_df["y"])
 ]
 
+# Create subplots
+fig, axes = plt.subplots(1, 2, figsize=(14, 6), constrained_layout=True)
+
+# Plot the screen-based well data on the first subplot
 imod.visualize.cross_section(
     xsection_layer_nr, "tab20", np.arange(21), fig=fig, ax=axes[0]
 )
 for _, row in well_df.iterrows():
     axes[0].vlines(
-        row["s"],
+        row["position_along_line"],
         row["screen_bottom"],
         row["screen_top"],
         color="black",
         linewidth=3,
     )
     axes[0].scatter(
-        row["s"], row["screen_top"], color="black", marker="1", s=100, linewidths=1.6
+        row["position_along_line"],
+        row["screen_top"],
+        color="black",
+        marker="1",
+        s=100,
+        linewidths=1.6,
     )
     axes[0].scatter(
-        row["s"], row["screen_bottom"], color="black", marker="2", s=100, linewidths=1.6
+        row["position_along_line"],
+        row["screen_bottom"],
+        color="black",
+        marker="2",
+        s=100,
+        linewidths=1.6,
     )
 axes[0].set_title("Well on layer cross section", fontsize=16)
 
-# LayeredWell
+# Prepare LayeredWell data for visualization
 layered_well_df = (
     layer_based.dataset[["x", "y", "layer"]].to_dataframe().reset_index(drop=True)
 )
-layered_well_df["s"] = [
+layered_well_df["position_along_line"] = [
     geometry.project(Point(x, y))
     for x, y in zip(layered_well_df["x"], layered_well_df["y"])
 ]
 
+# Plot the LayeredWell data on the second subplot
 imod.visualize.cross_section(
     xsection_layer_nr, "tab20", np.arange(21), fig=fig, ax=axes[1]
 )
 for _, row in layered_well_df.iterrows():
     axes[1].scatter(
-        row["s"], row["layer"], color="black", marker="x", s=120, linewidths=1.8
+        row["position_along_line"],
+        row["layer"],
+        color="black",
+        marker="x",
+        s=120,
+        linewidths=1.8,
     )
 axes[1].set_title("LayeredWell on layer cross section", fontsize=16)
diff --git a/imod/mf6/wel.py b/imod/mf6/wel.py
index 2bacf9056..3319c5744 100644
--- a/imod/mf6/wel.py
+++ b/imod/mf6/wel.py
@@ -450,7 +450,8 @@ def to_mf6_pkg(
         For :class:`Well`, the screen interval is split over
         intersected cells, and the rate is distributed by layer
         transmissivity when multiple layers are intersected. Wells located in
-        inactive cells are removed.
+        inactive cells, or in a layer less than the ``minimum_k`` or
+        ``minimum_thickness`` are removed.
 
         For :class:`LayeredWell`, the provided layer and rate are
         used directly, without screen splitting, screen-depth allocation, or
@@ -781,6 +782,8 @@ class Well(GridAgnosticWell):
     interval is intersected with model layers and the specified rate is
     distributed over eligible cells in proportion to transmissivity (based on
     horizontal hydraulic conductivity and screen-overlap thickness).
+    Wells assigned to layers that fall below the minimum_thickness or minimum_k
+    threshold are dropped.
 
     Use :class:`LayeredWell` when the target layer is already known and each
     rate must remain assigned to that layer.