From c3558b76a4dd2dd75fc1a47d9bf81ff712ff754b Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Tue, 22 Sep 2026 14:30:48 +0200
Subject: [PATCH 01/29] Add mask_packages method

---
 docs/api/changelog.rst           |  4 +++
 docs/api/mf6.rst                 |  2 ++
 imod/common/interfaces/imodel.py | 14 +++++++++-
 imod/common/utilities/mask.py    | 13 +++++----
 imod/mf6/model.py                | 46 ++++++++++++++++++++++++++++----
 5 files changed, 68 insertions(+), 11 deletions(-)

diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst
index a5b40a49a..6f02e203c 100644
--- a/docs/api/changelog.rst
+++ b/docs/api/changelog.rst
@@ -18,6 +18,10 @@ Added
   :meth:`imod.msw.SprinklingPoints.from_imod5_data`.
 - :class:`imod.mf6.LayeredWell.from_imod5_cap_data` now also supports loading
   wells from IPF files in an iMOD5 CAP dataset.
+- Added :meth:`imod.mf6.GroundwaterFlowModel.mask_packages` and
+  :meth:`imod.mf6.GroundwaterTransportModel.mask_packages` to mask specific
+  packages of a groundwater flow model and a groundwater transport model
+  respectively.
 
 Fixed
 ~~~~~
diff --git a/docs/api/mf6.rst b/docs/api/mf6.rst
index 32fb67246..dda794721 100644
--- a/docs/api/mf6.rst
+++ b/docs/api/mf6.rst
@@ -46,6 +46,7 @@ Model objects & methods
     Modflow6Simulation.set_validation_settings
     GroundwaterFlowModel
     GroundwaterFlowModel.mask_all_packages
+    GroundwaterFlowModel.mask_packages
     GroundwaterFlowModel.prepare_wel_for_mf6
     GroundwaterFlowModel.regrid_like
     GroundwaterFlowModel.dump
@@ -60,6 +61,7 @@ Model objects & methods
     GroundwaterFlowModel.get_diskey
     GroundwaterTransportModel
     GroundwaterTransportModel.mask_all_packages
+    GroundwaterTransportModel.mask_packages
     GroundwaterTransportModel.dump
     GroundwaterTransportModel.clip_box
     GroundwaterTransportModel.regrid_like
diff --git a/imod/common/interfaces/imodel.py b/imod/common/interfaces/imodel.py
index 2fe91e0a4..0b30e0a6e 100644
--- a/imod/common/interfaces/imodel.py
+++ b/imod/common/interfaces/imodel.py
@@ -16,9 +16,21 @@ class IModel(IDict):
     def mask_all_packages(self, mask: GridDataArray):
         raise NotImplementedError
 
+    @abstractmethod
+    def mask_packages(
+        self,
+        package_names: list[str],
+        mask: GridDataArray,
+        ignore_time_purge_empty: bool = False,
+    ):
+        raise NotImplementedError
+
     @abstractmethod
     def purge_empty_packages(
-        self, model_name: Optional[str] = "", ignore_time: bool = False
+        self,
+        model_name: Optional[str] = "",
+        ignore_time: bool = False,
+        package_names: list[str] | None = None,
     ) -> None:
         raise NotImplementedError
 
diff --git a/imod/common/utilities/mask.py b/imod/common/utilities/mask.py
index 21f689f27..eba69d86f 100644
--- a/imod/common/utilities/mask.py
+++ b/imod/common/utilities/mask.py
@@ -62,15 +62,18 @@ def mask_all_models(
             )
 
 
-def mask_all_packages(
+def mask_packages(
     model: IModel,
+    package_names: list[str],
     mask: GridDataArray,
     ignore_time_purge_empty: bool = False,
-):
+) -> None:
     _validate_coords_mask(mask)
-    for pkgname, pkg in model.items():
-        model[pkgname] = pkg.mask(mask)
-    model.purge_empty_packages(ignore_time=ignore_time_purge_empty)
+    for pkgname in package_names:
+        model[pkgname] = model[pkgname].mask(mask)
+    model.purge_empty_packages(
+        ignore_time=ignore_time_purge_empty, package_names=package_names
+    )
 
 
 def mask_package(package: IPackage, mask: GridDataArray) -> IPackage:
diff --git a/imod/mf6/model.py b/imod/mf6/model.py
index 2a931732f..54208cefe 100644
--- a/imod/mf6/model.py
+++ b/imod/mf6/model.py
@@ -22,7 +22,7 @@
 from imod.common.statusinfo import NestedStatusInfo, StatusInfo, StatusInfoBase
 from imod.common.utilities.clip import clip_box_dataset
 from imod.common.utilities.dump_model import dump_model
-from imod.common.utilities.mask import mask_all_packages
+from imod.common.utilities.mask import mask_packages
 from imod.common.utilities.regrid import _regrid_like
 from imod.common.utilities.schemata import (
     concatenate_schemata_dicts,
@@ -930,11 +930,41 @@ def mask_all_packages(
             Whether to ignore time dimension when purging empty packages. Can
             improve performance when masking models with many time steps.
         """
+        package_names = list(self.keys())
+        mask_packages(self, package_names, mask, ignore_time_purge_empty)
 
-        mask_all_packages(self, mask, ignore_time_purge_empty)
+    def mask_packages(
+        self,
+        package_names: list[str],
+        mask: GridDataArray,
+        ignore_time_purge_empty: bool = False,
+    ) -> None:
+        """
+        This function applies a mask to packages in a model. The mask must
+        be presented as an idomain-like integer array that has 0 (inactive) or
+        <0 (vertical passthrough) values in filtered cells and >0 in active
+        cells.
+        Masking will overwrite idomain with the mask where the mask is <=0.
+        Where the mask is >0, the original value of idomain will be kept. Masking
+        will update the packages accordingly, blanking their input where needed,
+        and is therefore not a reversible operation.
+
+        Parameters
+        ----------
+        mask: xr.DataArray, xu.UgridDataArray of ints
+            idomain-like integer array. >0 sets cells to active, 0 sets cells to inactive,
+            <0 sets cells to vertical passthrough
+        ignore_time_purge_empty: bool, default False
+            Whether to ignore time dimension when purging empty packages. Can
+            improve performance when masking models with many time steps.
+        """
+        mask_packages(self, package_names, mask, ignore_time_purge_empty)
 
     def purge_empty_packages(
-        self, model_name: Optional[str] = "", ignore_time: bool = False
+        self,
+        model_name: Optional[str] = "",
+        ignore_time: bool = False,
+        package_names: Optional[list[str]] = None,
     ) -> None:
         """
         This method removes empty packages from the model in place.
@@ -948,11 +978,17 @@ def purge_empty_packages(
             timesteps. If True, packages are considered empty if they have no
             data at the first time step. The latter can increase performance
             considerably.
+        package_names: list[str], optional
+            List of package names to check for emptiness. If None, all packages
+            are checked.
         """
+        if package_names is None:
+            package_names = list(self.keys())
+
         empty_packages = [
             package_name
-            for package_name, package in self.items()
-            if package.is_empty(ignore_time=ignore_time)
+            for package_name in package_names
+            if self[package_name].is_empty(ignore_time=ignore_time)
         ]
         logger.info(
             f"packages: {empty_packages} removed in {model_name}, because all empty"

From 1c8ffb0e6d3ac9f0f6bc141fa4cb7aeddf0fe379 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Tue, 22 Sep 2026 15:30:32 +0200
Subject: [PATCH 02/29] Add mask_topsystem_packages utility function

---
 imod/mf6/model_gwf.py                 |  8 ++++++
 imod/mf6/utilities/imod5_converter.py | 37 ++++++++++++++++++++++++++-
 2 files changed, 44 insertions(+), 1 deletion(-)

diff --git a/imod/mf6/model_gwf.py b/imod/mf6/model_gwf.py
index 545f948da..7e5d764af 100644
--- a/imod/mf6/model_gwf.py
+++ b/imod/mf6/model_gwf.py
@@ -37,6 +37,7 @@
 from imod.mf6.riv import River
 from imod.mf6.sto import StorageCoefficient
 from imod.mf6.utilities.chd_concat import concat_layered_chd_packages
+from imod.mf6.utilities.imod5_converter import mask_topsystem_packages
 from imod.mf6.validation_settings import ValidationSettings
 from imod.mf6.wel import LayeredWell, Well
 from imod.prepare.topsystem.default_allocation_methods import (
@@ -453,4 +454,11 @@ def from_imod5_data(
             for key, chd_package in chd_packages.items():
                 result[key] = chd_package
 
+        # Mask all topsystem packages where IBOUND == -1
+        mask_topsystem_packages(
+            imod5_data,
+            result,
+            cast(ConstantHeadRegridMethod, regridder_types.get("topsystem_mask")),
+            regrid_cache,
+        )
         return result
diff --git a/imod/mf6/utilities/imod5_converter.py b/imod/mf6/utilities/imod5_converter.py
index e27c21581..eb3a19bb6 100644
--- a/imod/mf6/utilities/imod5_converter.py
+++ b/imod/mf6/utilities/imod5_converter.py
@@ -4,10 +4,12 @@
 import pandas as pd
 import xarray as xr
 
+from imod.common.interfaces.imodel import IModel
 from imod.common.interfaces.iregridpackage import IRegridPackage
 from imod.common.utilities.dataclass_type import DataclassType
 from imod.common.utilities.regrid import _regrid_package_data, regrid_imod5_cap_data
 from imod.mf6.package import Package
+from imod.mf6.regrid.regrid_schemes import ConstantHeadRegridMethod
 from imod.typing import GridDataArray, GridDataDict, Imod5DataDict
 from imod.typing.grid import full_like
 from imod.util.regrid import RegridderWeightsCache
@@ -140,7 +142,7 @@ def well_from_imod5_cap_data(
 
 
 def regrid_imod5_pkg_data(
-    cls: type[Package],
+    cls: Optional[type[Package]],
     imod5_pkg_data: GridDataDict,
     target_dis: Package,
     regridder_types: Optional[DataclassType],
@@ -150,6 +152,11 @@ def regrid_imod5_pkg_data(
     Regrid iMOD5 package data to target idomain. Optionally get regrid methods
     from class if not provided.
     """
+    if cls is None and regridder_types is None:
+        raise ValueError(
+            "Either cls or regridder_types must be provided for regridding."
+        )
+
     target_idomain = target_dis.dataset["idomain"]
 
     # set up regridder methods
@@ -175,3 +182,31 @@ def chd_cells_from_imod5_data(
     head = head.where(target_idomain > 0)
 
     return {"head": head}
+
+
+def mask_topsystem_packages(
+    imod5_data: Imod5DataDict,
+    model: IModel,
+    regridder_types: ConstantHeadRegridMethod,
+    regrid_cache: RegridderWeightsCache,
+) -> None:
+    """
+    Mask all top system packages where IBOUND == -1.
+    """
+    from imod.mf6.topsystem import TopSystemBoundaryCondition
+
+    ibound = imod5_data["bnd"]["ibound"]
+    regridded_ibound = regrid_imod5_pkg_data(
+        cls=None,
+        imod5_pkg_data={"ibound": ibound},
+        target_dis=model["dis"],
+        regridder_types=regridder_types,
+        regrid_cache=regrid_cache,
+    )["ibound"]
+    mask = regridded_ibound == -1
+
+    topsystem_packages = [
+        key for key, pkg in model.items() if isinstance(pkg, TopSystemBoundaryCondition)
+    ]
+    for key in topsystem_packages:
+        model[key].mask(mask)

From be94e0f32a30f88d1daa18aded8fcef6a5b137bb Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Wed, 23 Sep 2026 12:36:15 +0200
Subject: [PATCH 03/29] Provide proper mask and fix mypy issues

---
 imod/mf6/utilities/imod5_converter.py | 27 +++++++++++++++++----------
 1 file changed, 17 insertions(+), 10 deletions(-)

diff --git a/imod/mf6/utilities/imod5_converter.py b/imod/mf6/utilities/imod5_converter.py
index eb3a19bb6..31995448d 100644
--- a/imod/mf6/utilities/imod5_converter.py
+++ b/imod/mf6/utilities/imod5_converter.py
@@ -1,4 +1,4 @@
-from typing import Optional, Union
+from typing import Optional, Union, cast
 
 import numpy as np
 import pandas as pd
@@ -152,16 +152,18 @@ def regrid_imod5_pkg_data(
     Regrid iMOD5 package data to target idomain. Optionally get regrid methods
     from class if not provided.
     """
-    if cls is None and regridder_types is None:
+    if (cls is None) and (regridder_types is None):
         raise ValueError(
             "Either cls or regridder_types must be provided for regridding."
         )
+    # set up regridder methods
+    elif (cls is not None) and (regridder_types is None):  # check cls not None for mypy
+        regridder_types = cls.get_regrid_methods()
+    # For mypy to succeed
+    regridder_types = cast(DataclassType, regridder_types)
 
     target_idomain = target_dis.dataset["idomain"]
 
-    # set up regridder methods
-    if regridder_types is None:
-        regridder_types = cls.get_regrid_methods()
     # regrid the input data
     regridded_pkg_data = _regrid_package_data(
         imod5_pkg_data, target_idomain, regridder_types, regrid_cache, {}
@@ -185,16 +187,21 @@ def chd_cells_from_imod5_data(
 
 
 def mask_topsystem_packages(
-    imod5_data: Imod5DataDict,
+    imod5_data: dict[str, dict[str, GridDataArray]],
     model: IModel,
-    regridder_types: ConstantHeadRegridMethod,
+    regridder_types: Optional[ConstantHeadRegridMethod],
     regrid_cache: RegridderWeightsCache,
 ) -> None:
     """
-    Mask all top system packages where IBOUND == -1.
+    Mask all top system packages where IBOUND < 0. These locations are assigned
+    a constant head.
     """
+    # Import here to avoid circular import issues
     from imod.mf6.topsystem import TopSystemBoundaryCondition
 
+    if regridder_types is None:
+        regridder_types = ConstantHeadRegridMethod()
+
     ibound = imod5_data["bnd"]["ibound"]
     regridded_ibound = regrid_imod5_pkg_data(
         cls=None,
@@ -203,10 +210,10 @@ def mask_topsystem_packages(
         regridder_types=regridder_types,
         regrid_cache=regrid_cache,
     )["ibound"]
-    mask = regridded_ibound == -1
+    is_active = regridded_ibound >= 0
 
     topsystem_packages = [
         key for key, pkg in model.items() if isinstance(pkg, TopSystemBoundaryCondition)
     ]
     for key in topsystem_packages:
-        model[key].mask(mask)
+        model[key] = model[key].mask(is_active)

From 69bd29c852791bad0e35ab188afefdec32759991 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Wed, 23 Sep 2026 12:37:35 +0200
Subject: [PATCH 04/29] Add test

---
 imod/tests/test_mf6/test_mf6_simulation.py | 37 ++++++++++++++++++++++
 1 file changed, 37 insertions(+)

diff --git a/imod/tests/test_mf6/test_mf6_simulation.py b/imod/tests/test_mf6/test_mf6_simulation.py
index e980b0eb9..3c1c1a718 100644
--- a/imod/tests/test_mf6/test_mf6_simulation.py
+++ b/imod/tests/test_mf6/test_mf6_simulation.py
@@ -610,6 +610,43 @@ def test_import_from_imod5(imod5_dataset, tmp_path):
     assert simulation._validation_context.strict_well_validation is False
 
 
+@pytest.mark.unittest_jit
+def test_import_from_imod5_mask_topsystem(imod5_dataset):
+    """Test importing from imod5 masks the top system packages"""
+    # Arrange
+    imod5_data = imod5_dataset[0]
+    period_data = imod5_dataset[1]
+
+    datelist = pd.date_range(start="1/1/1989", end="1/1/2013", freq="W")
+    # Act
+    simulation = Modflow6Simulation.from_imod5_data(
+        imod5_data,
+        period_data,
+        datelist,
+        SimulationAllocationOptions,
+        SimulationDistributingOptions,
+    )
+    # Assert
+    ibound = imod5_data["bnd"]["ibound"].isel(layer=0, drop=True)
+    topsystem_mask = ibound < 0
+    topsystem_keys = [
+        "rch",
+        "drn-1",
+        "drn-2",
+        "riv-1riv",
+        "riv-1drn",
+        "riv-2riv",
+        "riv-2drn",
+    ]
+    for key in topsystem_keys:
+        topsystem_pkg = simulation["imported_model"][key]
+        # Take first grid var
+        gridded_var = topsystem_pkg.dataset[topsystem_pkg._period_data[0]]
+        # True wherever topsystem is masked but gridded_var still has a value
+        bad = gridded_var.notnull() & topsystem_mask
+        assert not bad.any().item()
+
+
 @pytest.mark.unittest_jit
 def test_import_from_imod5__custom_name(imod5_dataset):
     imod5_data = imod5_dataset[0]

From d0970c3d0facba5bf1873840237472c359351ef3 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Wed, 23 Sep 2026 12:37:42 +0200
Subject: [PATCH 05/29] Fix docstring

---
 imod/mf6/model_gwf.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/imod/mf6/model_gwf.py b/imod/mf6/model_gwf.py
index 7e5d764af..0faf210eb 100644
--- a/imod/mf6/model_gwf.py
+++ b/imod/mf6/model_gwf.py
@@ -454,7 +454,7 @@ def from_imod5_data(
             for key, chd_package in chd_packages.items():
                 result[key] = chd_package
 
-        # Mask all topsystem packages where IBOUND == -1
+        # Mask all topsystem packages where IBOUND < 0
         mask_topsystem_packages(
             imod5_data,
             result,

From c8cbc5f2b029531da0f874518e2fe802f29c9355 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Wed, 23 Sep 2026 13:04:11 +0200
Subject: [PATCH 06/29] Move mask topsystem function to separate utility and
 add ITopSystemBoundaryCondition interface.

---
 imod/common/interfaces/itopsystembc.py | 19 +++++++++++++++++++
 imod/mf6/model_gwf.py                  |  4 ++--
 imod/mf6/topsystem.py                  |  5 ++++-
 imod/mf6/utilities/imod5_converter.py  | 11 +++--------
 imod/mf6/utilities/mask.py             | 25 +++++++++++++++++++++++++
 5 files changed, 53 insertions(+), 11 deletions(-)
 create mode 100644 imod/common/interfaces/itopsystembc.py
 create mode 100644 imod/mf6/utilities/mask.py

diff --git a/imod/common/interfaces/itopsystembc.py b/imod/common/interfaces/itopsystembc.py
new file mode 100644
index 000000000..e202064c4
--- /dev/null
+++ b/imod/common/interfaces/itopsystembc.py
@@ -0,0 +1,19 @@
+from abc import abstractmethod
+
+from imod.common.interfaces.ipackage import IPackage
+from imod.typing import GridDataDict, GridDataset
+
+
+class ITopSystemBoundaryCondition(IPackage):
+    """
+    Interface for top system boundary condition packages in MODFLOW 6.
+    """
+
+    @classmethod
+    @abstractmethod
+    def aggregate_layers(cls, dataset: GridDataset) -> GridDataDict:
+        raise NotImplementedError
+
+    @abstractmethod
+    def reallocate(self):
+        raise NotImplementedError
diff --git a/imod/mf6/model_gwf.py b/imod/mf6/model_gwf.py
index 0faf210eb..5296bdc62 100644
--- a/imod/mf6/model_gwf.py
+++ b/imod/mf6/model_gwf.py
@@ -37,7 +37,7 @@
 from imod.mf6.riv import River
 from imod.mf6.sto import StorageCoefficient
 from imod.mf6.utilities.chd_concat import concat_layered_chd_packages
-from imod.mf6.utilities.imod5_converter import mask_topsystem_packages
+from imod.mf6.utilities.imod5_converter import mask_topsystem_packages_with_ibound
 from imod.mf6.validation_settings import ValidationSettings
 from imod.mf6.wel import LayeredWell, Well
 from imod.prepare.topsystem.default_allocation_methods import (
@@ -455,7 +455,7 @@ def from_imod5_data(
                 result[key] = chd_package
 
         # Mask all topsystem packages where IBOUND < 0
-        mask_topsystem_packages(
+        mask_topsystem_packages_with_ibound(
             imod5_data,
             result,
             cast(ConstantHeadRegridMethod, regridder_types.get("topsystem_mask")),
diff --git a/imod/mf6/topsystem.py b/imod/mf6/topsystem.py
index 14d67e0f3..692fa1d7a 100644
--- a/imod/mf6/topsystem.py
+++ b/imod/mf6/topsystem.py
@@ -3,6 +3,7 @@
 from dataclasses import asdict
 from typing import Optional, Self, cast
 
+from imod.common.interfaces.itopsystembc import ITopSystemBoundaryCondition
 from imod.common.utilities.dataclass_type import DataclassType
 from imod.mf6.aggregate.aggregate_schemes import EmptyAggregationMethod
 from imod.mf6.boundary_condition import BoundaryCondition
@@ -41,7 +42,9 @@ def _handle_reallocate_arguments(
     return allocation_option, distributing_option
 
 
-class TopSystemBoundaryCondition(BoundaryCondition, abc.ABC):
+class TopSystemBoundaryCondition(
+    BoundaryCondition, ITopSystemBoundaryCondition, abc.ABC
+):
     """
     Base class to add some extra functionality for topsystem packages, such as
     RCH, DRN, RIV, and GHB.
diff --git a/imod/mf6/utilities/imod5_converter.py b/imod/mf6/utilities/imod5_converter.py
index 31995448d..0da531a2a 100644
--- a/imod/mf6/utilities/imod5_converter.py
+++ b/imod/mf6/utilities/imod5_converter.py
@@ -10,6 +10,7 @@
 from imod.common.utilities.regrid import _regrid_package_data, regrid_imod5_cap_data
 from imod.mf6.package import Package
 from imod.mf6.regrid.regrid_schemes import ConstantHeadRegridMethod
+from imod.mf6.utilities.mask import mask_topsystem
 from imod.typing import GridDataArray, GridDataDict, Imod5DataDict
 from imod.typing.grid import full_like
 from imod.util.regrid import RegridderWeightsCache
@@ -186,7 +187,7 @@ def chd_cells_from_imod5_data(
     return {"head": head}
 
 
-def mask_topsystem_packages(
+def mask_topsystem_packages_with_ibound(
     imod5_data: dict[str, dict[str, GridDataArray]],
     model: IModel,
     regridder_types: Optional[ConstantHeadRegridMethod],
@@ -196,8 +197,6 @@ def mask_topsystem_packages(
     Mask all top system packages where IBOUND < 0. These locations are assigned
     a constant head.
     """
-    # Import here to avoid circular import issues
-    from imod.mf6.topsystem import TopSystemBoundaryCondition
 
     if regridder_types is None:
         regridder_types = ConstantHeadRegridMethod()
@@ -212,8 +211,4 @@ def mask_topsystem_packages(
     )["ibound"]
     is_active = regridded_ibound >= 0
 
-    topsystem_packages = [
-        key for key, pkg in model.items() if isinstance(pkg, TopSystemBoundaryCondition)
-    ]
-    for key in topsystem_packages:
-        model[key] = model[key].mask(is_active)
+    mask_topsystem(model, is_active)
diff --git a/imod/mf6/utilities/mask.py b/imod/mf6/utilities/mask.py
new file mode 100644
index 000000000..574422d48
--- /dev/null
+++ b/imod/mf6/utilities/mask.py
@@ -0,0 +1,25 @@
+from imod.common.interfaces.imodel import IModel
+from imod.common.interfaces.itopsystembc import ITopSystemBoundaryCondition
+from imod.typing import GridDataArray
+
+
+def mask_topsystem(model: IModel, is_active: GridDataArray):
+    """
+    Mask all top system packages in the model inplace with a boolean mask
+    indicating active cells.
+
+    Parameters
+    ----------
+    model : IModel
+        The MODFLOW 6 model containing top system packages.
+    is_active : GridDataArray
+        A boolean array indicating active cells. Top system packages will be masked
+        where this array is False.
+    """
+    topsystem_packages = [
+        key
+        for key, pkg in model.items()
+        if isinstance(pkg, ITopSystemBoundaryCondition)
+    ]
+    for key in topsystem_packages:
+        model[key] = model[key].mask(is_active)

From 29204bd381952fd4a8ecf89282453f1bc4961980 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Wed, 23 Sep 2026 13:09:28 +0200
Subject: [PATCH 07/29] Remove method

---
 imod/common/interfaces/itopsystembc.py | 4 ----
 1 file changed, 4 deletions(-)

diff --git a/imod/common/interfaces/itopsystembc.py b/imod/common/interfaces/itopsystembc.py
index e202064c4..de541ec09 100644
--- a/imod/common/interfaces/itopsystembc.py
+++ b/imod/common/interfaces/itopsystembc.py
@@ -13,7 +13,3 @@ class ITopSystemBoundaryCondition(IPackage):
     @abstractmethod
     def aggregate_layers(cls, dataset: GridDataset) -> GridDataDict:
         raise NotImplementedError
-
-    @abstractmethod
-    def reallocate(self):
-        raise NotImplementedError

From a4ba60d0ba03cfcc3f3deb338d8d2cc02f5d6970 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Wed, 23 Sep 2026 14:02:25 +0200
Subject: [PATCH 08/29] Also clip topsystems for clip_box when
 states_for_boundary are provided

---
 imod/mf6/model.py                     |  8 ++++++++
 imod/tests/test_mf6/test_ex01_twri.py | 11 +++++++++++
 2 files changed, 19 insertions(+)

diff --git a/imod/mf6/model.py b/imod/mf6/model.py
index 54208cefe..70f21325f 100644
--- a/imod/mf6/model.py
+++ b/imod/mf6/model.py
@@ -43,11 +43,13 @@
     StateType,
     create_clipped_boundary,
 )
+from imod.mf6.utilities.mask import mask_topsystem
 from imod.mf6.utilities.mf6hfb import merge_hfb_packages
 from imod.mf6.validation_settings import ValidationSettings
 from imod.mf6.wel import GridAgnosticWell
 from imod.mf6.write_context import WriteContext
 from imod.schemata import SchemataDict, ValidationError
+from imod.select.grid import active_grid_boundary_xy
 from imod.typing import GridDataArray
 from imod.util.regrid import RegridderWeightsCache
 
@@ -813,6 +815,12 @@ def clip_box(
         if clipped_boundary_condition is not None:
             clipped[pkg_name] = clipped_boundary_condition
 
+            # Clip topsystem packages where the active grid boundary has
+            # changed.
+            _, _, idomain_clipped = clipped._get_domain_geometry()
+            active_bounds_clipped = active_grid_boundary_xy(idomain_clipped > 0)
+            mask_topsystem(clipped, ~active_bounds_clipped)
+
         clipped.purge_empty_packages(ignore_time=ignore_time_purge_empty)
 
         return clipped
diff --git a/imod/tests/test_mf6/test_ex01_twri.py b/imod/tests/test_mf6/test_ex01_twri.py
index 9cd2663fa..54482791b 100644
--- a/imod/tests/test_mf6/test_ex01_twri.py
+++ b/imod/tests/test_mf6/test_ex01_twri.py
@@ -577,6 +577,17 @@ def test_slice_and_run_with_state(transient_twri_model_extended, tmp_path):
     np_array = clipped_boundary["head"].values
     assert (np_array == 1.23).sum() == 33
 
+    # Test that topsystem packages are masked.
+    topsystem_keys = ["rch", "drn"]
+    topsystem_mask = clipped_boundary["head"].notnull().compute()
+    for key in topsystem_keys:
+        topsystem_pkg = clipped_simulation["GWF_1"][key]
+        # Take first grid var
+        gridded_var = topsystem_pkg.dataset[topsystem_pkg._period_data[0]].compute()
+        # True wherever topsystem is masked but gridded_var still has a value
+        bad = gridded_var.notnull() & topsystem_mask
+        assert not bad.any().item()
+
 
 @pytest.mark.skipif(sys.version_info < (3, 7), reason="capture_output added in 3.7")
 def test_slice_and_run_purge_empty_package(transient_twri_model, tmp_path):

From 433cd1000bc93a64dc72fcd14d983378678996ef Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Wed, 23 Sep 2026 14:10:42 +0200
Subject: [PATCH 09/29] Add test for masking the topsystem

---
 .../test_mf6/test_utilities/test_mask.py      | 19 +++++++++++++++++++
 1 file changed, 19 insertions(+)
 create mode 100644 imod/tests/test_mf6/test_utilities/test_mask.py

diff --git a/imod/tests/test_mf6/test_utilities/test_mask.py b/imod/tests/test_mf6/test_utilities/test_mask.py
new file mode 100644
index 000000000..19275469a
--- /dev/null
+++ b/imod/tests/test_mf6/test_utilities/test_mask.py
@@ -0,0 +1,19 @@
+from imod.mf6.utilities.mask import mask_topsystem
+from imod.typing.grid import zeros_like
+
+
+def test_mask_topsystem(twri_model):
+    """
+    Test the mask_topsystem utility function by deactivating all cells in the
+    grid.
+    """
+    # Arrange
+    gwf_model = twri_model["GWF_1"]
+    mask = zeros_like(gwf_model.domain)
+    # Act
+    mask_topsystem(gwf_model, mask)
+    # Assert
+    for key in ["rch", "drn"]:
+        pkg = gwf_model[key]
+        gridded_var = pkg.dataset[pkg._period_data[0]].compute()
+        assert not gridded_var.notnull().any().item()

From 32e6772ae0023214e5733f6faf04080d61fca09d Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Wed, 23 Sep 2026 14:13:13 +0200
Subject: [PATCH 10/29] Return None in type annotation

---
 imod/mf6/utilities/mask.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/imod/mf6/utilities/mask.py b/imod/mf6/utilities/mask.py
index 574422d48..3de4f422d 100644
--- a/imod/mf6/utilities/mask.py
+++ b/imod/mf6/utilities/mask.py
@@ -3,7 +3,7 @@
 from imod.typing import GridDataArray
 
 
-def mask_topsystem(model: IModel, is_active: GridDataArray):
+def mask_topsystem(model: IModel, is_active: GridDataArray) -> None:
     """
     Mask all top system packages in the model inplace with a boolean mask
     indicating active cells.

From b3a59fc32140863a541ab27e4f0361a7f8056bf6 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Wed, 23 Sep 2026 14:25:14 +0200
Subject: [PATCH 11/29] Use mask_packages method

---
 imod/mf6/utilities/mask.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/imod/mf6/utilities/mask.py b/imod/mf6/utilities/mask.py
index 3de4f422d..99a969d9d 100644
--- a/imod/mf6/utilities/mask.py
+++ b/imod/mf6/utilities/mask.py
@@ -21,5 +21,4 @@ def mask_topsystem(model: IModel, is_active: GridDataArray) -> None:
         for key, pkg in model.items()
         if isinstance(pkg, ITopSystemBoundaryCondition)
     ]
-    for key in topsystem_packages:
-        model[key] = model[key].mask(is_active)
+    model.mask_packages(topsystem_packages, is_active)

From 698f36e43ca45fc466dfb7a23420b26d766eab01 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Wed, 23 Sep 2026 15:07:17 +0200
Subject: [PATCH 12/29] Regrid iMOD5 IBOUND data when regridding cap data and
 also mask where ibound is not active

---
 imod/common/utilities/regrid.py       | 4 ++++
 imod/mf6/rch.py                       | 8 +++++---
 imod/msw/grid_data.py                 | 3 ++-
 imod/msw/regrid/regrid_schemes.py     | 1 +
 imod/msw/utilities/imod5_converter.py | 4 +++-
 imod/typing/__init__.py               | 1 +
 6 files changed, 16 insertions(+), 5 deletions(-)

diff --git a/imod/common/utilities/regrid.py b/imod/common/utilities/regrid.py
index aaf28abb9..f7c389ebf 100644
--- a/imod/common/utilities/regrid.py
+++ b/imod/common/utilities/regrid.py
@@ -492,9 +492,13 @@ def regrid_imod5_cap_data(
     cap_data_regridded = _regrid_package_data(
         imod5_cap_no_layer["cap"], target_grid, regridder_types, regrid_cache
     )
+    bnd_data_regridded = _regrid_package_data(
+        imod5_data["bnd"], target_grid, regridder_types, regrid_cache
+    )
     extra_paths = imod5_data["extra"]["paths"]
     imod5_regridded: Imod5DataDict = {
         "cap": cap_data_regridded,
+        "bnd": bnd_data_regridded,
         "extra": {"paths": extra_paths},
     }
     return imod5_regridded
diff --git a/imod/mf6/rch.py b/imod/mf6/rch.py
index 1933f480e..b221b451a 100644
--- a/imod/mf6/rch.py
+++ b/imod/mf6/rch.py
@@ -298,12 +298,14 @@ def from_imod5_cap_data(
         used to couple MODFLOW6 to MetaSWAP models. Active cells will have a
         recharge rate of 0.0.
         """
-        cap_data = regrid_imod5_cap_data(
+        imod5_data_regridded = regrid_imod5_cap_data(
             imod5_data, target_dis, regridder_types, regrid_cache
-        )["cap"]
+        )
+        cap_data = imod5_data_regridded["cap"]
+        bnd_data = imod5_data_regridded["bnd"]
 
         msw_area = get_cell_area_from_imod5_data(cap_data)
-        msw_active = is_msw_active_cell(target_dis, cap_data, msw_area)
+        msw_active = is_msw_active_cell(target_dis, cap_data, msw_area, bnd_data)
         active = msw_active.all
 
         data = {}
diff --git a/imod/msw/grid_data.py b/imod/msw/grid_data.py
index 24fbece43..a15fe8a47 100644
--- a/imod/msw/grid_data.py
+++ b/imod/msw/grid_data.py
@@ -202,6 +202,7 @@ def from_imod5_data(
             as aggregated over subunits.
         """
         imod5_cap = imod5_data["cap"]
+        imod5_bnd = imod5_data["bnd"]
 
         data = {}
         data["area"] = get_cell_area_from_imod5_data(imod5_cap)
@@ -210,7 +211,7 @@ def from_imod5_data(
         data["surface_elevation"] = imod5_cap["surface_elevation"]
         data["soil_physical_unit"] = imod5_cap["soil_physical_unit"].astype(int)
 
-        msw_active = is_msw_active_cell(target_dis, imod5_cap, data["area"])
+        msw_active = is_msw_active_cell(target_dis, imod5_cap, data["area"], imod5_bnd)
         data_active = mask_and_broadcast_pkg_data(cls, data, msw_active)
         data_active["active"] = msw_active.all
         return cls(**data_active), msw_active
diff --git a/imod/msw/regrid/regrid_schemes.py b/imod/msw/regrid/regrid_schemes.py
index 30547abfc..a7e73207a 100644
--- a/imod/msw/regrid/regrid_schemes.py
+++ b/imod/msw/regrid/regrid_schemes.py
@@ -73,6 +73,7 @@ class CapDataRegridMethod(DataclassType):
     steering_location: RegridVarType = (RegridderType.OVERLAP, "mode")
     plot_drainage_level: RegridVarType = (RegridderType.OVERLAP, "mean")
     plot_drainage_resistance: RegridVarType = (RegridderType.OVERLAP, "mean")
+    ibound: RegridVarType = (RegridderType.OVERLAP, "mode")
 
 
 @dataclass(config=_CONFIG)
diff --git a/imod/msw/utilities/imod5_converter.py b/imod/msw/utilities/imod5_converter.py
index 17e2923e3..772f11277 100644
--- a/imod/msw/utilities/imod5_converter.py
+++ b/imod/msw/utilities/imod5_converter.py
@@ -96,6 +96,7 @@ def is_msw_active_cell(
     target_dis: StructuredDiscretization,
     imod5_cap: GridDataDict,
     msw_area: GridDataArray,
+    imod5_bnd: GridDataArray,
 ) -> MetaSwapActive:
     """
     Return grid of cells that are active in the coupled computation, based on
@@ -113,7 +114,8 @@ def is_msw_active_cell(
         Cells active per subunit
     """
     mf6_top_active = target_dis["idomain"].isel(layer=0, drop=True)
-    subunit_active = (imod5_cap["boundary"] > 0) & (msw_area > 0) & (mf6_top_active > 0)
+    imod5_active = (imod5_bnd["ibound"] > 0) & (imod5_cap["boundary"] > 0)
+    subunit_active = imod5_active & (msw_area > 0) & (mf6_top_active > 0)
     active = subunit_active.any(dim="subunit")
     return MetaSwapActive(active, subunit_active)
 
diff --git a/imod/typing/__init__.py b/imod/typing/__init__.py
index 1c2a4885c..250d9b69a 100644
--- a/imod/typing/__init__.py
+++ b/imod/typing/__init__.py
@@ -33,6 +33,7 @@ class DropVarsType(TypedDict, total=False):
 
 
 class Imod5DataDict(TypedDict, total=False):
+    bnd: GridDataDict
     cap: GridDataDict
     extra: dict[str, list[str]]
 

From 03aabc0578cd197e1a0f731049632488104e0c37 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Wed, 23 Sep 2026 15:09:01 +0200
Subject: [PATCH 13/29] Also drop bnd layer

---
 imod/util/dims.py | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/imod/util/dims.py b/imod/util/dims.py
index c9228ee39..de3c1fd01 100644
--- a/imod/util/dims.py
+++ b/imod/util/dims.py
@@ -57,4 +57,8 @@ def _drop_layer_if_dataarray(
 
 def drop_layer_dim_cap_data(imod5_data: Imod5DataDict) -> Imod5DataDict:
     cap_data = imod5_data["cap"]
-    return {"cap": {key: _drop_layer_if_dataarray(da) for key, da in cap_data.items()}}
+    bnd_data = imod5_data["bnd"]
+    return {
+        "cap": {key: _drop_layer_if_dataarray(da) for key, da in cap_data.items()},
+        "bnd": {key: _drop_layer_if_dataarray(da) for key, da in bnd_data.items()},
+    }

From 3171d54129c34848b20ad54b848ce5bd679cd409 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Wed, 23 Sep 2026 16:02:37 +0200
Subject: [PATCH 14/29] Call correct var an improve varname

---
 imod/common/utilities/regrid.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/imod/common/utilities/regrid.py b/imod/common/utilities/regrid.py
index f7c389ebf..fd264378a 100644
--- a/imod/common/utilities/regrid.py
+++ b/imod/common/utilities/regrid.py
@@ -486,14 +486,14 @@ def regrid_imod5_cap_data(
     and ``imod.mf6.Recharge.from_imod5_cap_data``.
     """
     # Drop layer coords
-    imod5_cap_no_layer = drop_layer_dim_cap_data(imod5_data)
+    imod5_no_layer = drop_layer_dim_cap_data(imod5_data)
     target_grid = target_dis.dataset["idomain"].isel(layer=0, drop=True)
     # Regrid the input data
     cap_data_regridded = _regrid_package_data(
-        imod5_cap_no_layer["cap"], target_grid, regridder_types, regrid_cache
+        imod5_no_layer["cap"], target_grid, regridder_types, regrid_cache
     )
     bnd_data_regridded = _regrid_package_data(
-        imod5_data["bnd"], target_grid, regridder_types, regrid_cache
+        imod5_no_layer["bnd"], target_grid, regridder_types, regrid_cache
     )
     extra_paths = imod5_data["extra"]["paths"]
     imod5_regridded: Imod5DataDict = {

From ac10b19cce1f7ec478f39a0febed8b0d44d8cd75 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Wed, 23 Sep 2026 16:02:46 +0200
Subject: [PATCH 15/29] Add docstring

---
 imod/msw/utilities/imod5_converter.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/imod/msw/utilities/imod5_converter.py b/imod/msw/utilities/imod5_converter.py
index 772f11277..7f18a4acf 100644
--- a/imod/msw/utilities/imod5_converter.py
+++ b/imod/msw/utilities/imod5_converter.py
@@ -114,6 +114,8 @@ def is_msw_active_cell(
         Cells active per subunit
     """
     mf6_top_active = target_dis["idomain"].isel(layer=0, drop=True)
+    # Where IBOUND = -1, there also shouldn't be any active cells in the CAP
+    # boundary array.
     imod5_active = (imod5_bnd["ibound"] > 0) & (imod5_cap["boundary"] > 0)
     subunit_active = imod5_active & (msw_area > 0) & (mf6_top_active > 0)
     active = subunit_active.any(dim="subunit")

From 65cd181622298d2a859a5f9844f7db1dee42b061 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Wed, 23 Sep 2026 16:03:41 +0200
Subject: [PATCH 16/29] Add ibound to test fixture and expand tests to test for
 cell inactivity

---
 imod/tests/fixtures/msw_imod5_cap_fixture.py | 13 +++++++++++++
 imod/tests/test_msw/test_grid_data.py        | 17 +++++++++++++++--
 2 files changed, 28 insertions(+), 2 deletions(-)

diff --git a/imod/tests/fixtures/msw_imod5_cap_fixture.py b/imod/tests/fixtures/msw_imod5_cap_fixture.py
index 28c53335b..76bb895f6 100644
--- a/imod/tests/fixtures/msw_imod5_cap_fixture.py
+++ b/imod/tests/fixtures/msw_imod5_cap_fixture.py
@@ -232,6 +232,19 @@ def imod5_cap_data() -> GridDataDict:
         ),
         **da_kwargs
     )
+    d2 = {}
+    d2["ibound"] = xr.DataArray(
+        np.array([
+            [
+                [1, 1, 1],
+                [0, 0, 0],
+                [1, 1, 0],
+            ]],
+        dtype=int),
+        **da_kwargs
+    )
     # fmt: on
     imod5_data["cap"] = d
+    imod5_data["bnd"] = d2
+
     return imod5_data
diff --git a/imod/tests/test_msw/test_grid_data.py b/imod/tests/test_msw/test_grid_data.py
index f38559d14..c3b40d657 100644
--- a/imod/tests/test_msw/test_grid_data.py
+++ b/imod/tests/test_msw/test_grid_data.py
@@ -474,7 +474,12 @@ def test_from_imod5_data(grid_data_dict: dict[str, xr.DataArray]):
     cap_data["soil_physical_unit"] = xr.ones_like(like, dtype=int)
     cap_data["active"] = xr.ones_like(like, dtype=bool)
 
-    imod5_data = {"cap": cap_data}
+    ibound = xr.ones_like(like, dtype=int)
+    ibound[0, 0] = -1
+
+    bnd_data = {}
+    bnd_data["ibound"] = ibound
+    imod5_data = {"cap": cap_data, "bnd": bnd_data}
 
     layer = xr.DataArray([1, 1], coords={"layer": [1, 2]}, dims=("layer",))
     idomain = layer * xr.ones_like(like, dtype=int)
@@ -484,7 +489,15 @@ def test_from_imod5_data(grid_data_dict: dict[str, xr.DataArray]):
 
     griddata, _ = GridData.from_imod5_data(imod5_data, target_dis=dis)
     expected_rootzone_depth = cap_data["rootzone_thickness"] * 0.01
+    expected_rootzone_depth[0, 0] = np.nan
     xr.testing.assert_allclose(
         expected_rootzone_depth, griddata["rootzone_depth"].sel(subunit=0, drop=True)
     )
-    assert (griddata["landuse"].sel(subunit=1, drop=True) == 18).all()
+    # Test if all cells in subunit = 1 set to "urban" landuse code
+    np.testing.assert_array_equal(
+        np.unique(griddata["landuse"].sel(subunit=1, drop=True)), [0, 18]
+    )
+    # Test if cell where IBOUND == -1 is set to inactive (landuse = 0)
+    np.testing.assert_array_equal(griddata["landuse"][:, 0, 0], [0, 0])
+    np.testing.assert_array_equal(griddata["rootzone_depth"][:, 0, 0], [np.nan, np.nan])
+    np.testing.assert_array_equal(griddata["active"][0, 0], [False, False])

From 9441fd40c915eaf25520d1ab15b648887a45f8c3 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Wed, 23 Sep 2026 16:14:24 +0200
Subject: [PATCH 17/29] Update changelog

---
 docs/api/changelog.rst | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst
index 6f02e203c..1ecc7c0d5 100644
--- a/docs/api/changelog.rst
+++ b/docs/api/changelog.rst
@@ -29,6 +29,12 @@ Fixed
 - Fixed resampling in :meth:`imod.mf6.Well.from_imod5_data` and
   :meth:`imod.mf6.LayeredWell.from_imod5_data` when simulation timesteps precede
   the first well timestep.
+- :meth:`imod.mf6.GroundwaterFlowModel.from_imod5_data` now masks cells in
+  topsystem packages (:class:`imod.mf6.River`,
+  :class:`imod.mf6.GeneralHeadBoundary`, :class:`imod.mf6.Drainage`,
+  :class:`imod.mf6.Recharge`) where IBOUND is less than 0.
+- :meth:`imod.msw.MetaSwapModel.from_imod5_data` now masks cells where IBOUND is
+  less than 0.
 
 Changed
 ~~~~~~~
@@ -36,6 +42,11 @@ Changed
 - Deprecated :class:`imod.msw.Sprinkling` in favor of
   :class:`imod.msw.SprinklingGrid`. Call :class:`imod.msw.SprinklingGrid` to get
   the same behavior as you were used to.
+- If ``states_for_boundary`` is provided to
+  :meth:`imod.mf6.GroundwaterFlowModel.clip_box`, topsystem packages
+  (:class:`imod.mf6.River`, :class:`imod.mf6.GeneralHeadBoundary`,
+  :class:`imod.mf6.Drainage`, :class:`imod.mf6.Recharge`) will also be masked
+  where constant head cells are placed. 
 
 [1.1.0] - 2026-08-03
 --------------------

From e6e7a0365a21d1b7c09ff538ba0ea22244e0f08d Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Wed, 23 Sep 2026 17:00:00 +0200
Subject: [PATCH 18/29] Rename to avoid duplicate test module names

---
 .../test_utilities/{test_mask.py => test_mf6_mask_util.py}        | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 rename imod/tests/test_mf6/test_utilities/{test_mask.py => test_mf6_mask_util.py} (100%)

diff --git a/imod/tests/test_mf6/test_utilities/test_mask.py b/imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py
similarity index 100%
rename from imod/tests/test_mf6/test_utilities/test_mask.py
rename to imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py

From fb1b855660fb95d5f2e98dc13d74c298d955721b Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Wed, 23 Sep 2026 17:03:16 +0200
Subject: [PATCH 19/29] Also rename msw mask util test module

---
 .../test_utilities/{test_mask.py => test_msw_mask_util.py}        | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 rename imod/tests/test_msw/test_utilities/{test_mask.py => test_msw_mask_util.py} (100%)

diff --git a/imod/tests/test_msw/test_utilities/test_mask.py b/imod/tests/test_msw/test_utilities/test_msw_mask_util.py
similarity index 100%
rename from imod/tests/test_msw/test_utilities/test_mask.py
rename to imod/tests/test_msw/test_utilities/test_msw_mask_util.py

From 52b6ea42c5af451ce3c8b9000dd340559d616c62 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Wed, 23 Sep 2026 17:18:33 +0200
Subject: [PATCH 20/29] Fix and expand mf6 mask tests

---
 .../test_utilities/test_mf6_mask_util.py      | 31 ++++++++++++++++---
 1 file changed, 26 insertions(+), 5 deletions(-)

diff --git a/imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py b/imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py
index 19275469a..8993fa7b3 100644
--- a/imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py
+++ b/imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py
@@ -1,6 +1,6 @@
 from imod.mf6.utilities.mask import mask_topsystem
-from imod.typing.grid import zeros_like
-
+from imod.typing.grid import zeros_like, ones_like
+import numpy as np
 
 def test_mask_topsystem(twri_model):
     """
@@ -9,11 +9,32 @@ def test_mask_topsystem(twri_model):
     """
     # Arrange
     gwf_model = twri_model["GWF_1"]
-    mask = zeros_like(gwf_model.domain)
+    is_active = ones_like(gwf_model.domain)
+    # Mask first cell
+    is_active[0, 0, 0] = 0
+
     # Act
-    mask_topsystem(gwf_model, mask)
+    mask_topsystem(gwf_model, is_active)
     # Assert
     for key in ["rch", "drn"]:
         pkg = gwf_model[key]
         gridded_var = pkg.dataset[pkg._period_data[0]].compute()
-        assert not gridded_var.notnull().any().item()
+        first_cell = gridded_var.data.ravel()[0]
+        assert np.isnan(first_cell).item()
+
+
+
+def test_mask_topsystem__all_removed(twri_model):
+    """
+    Test the mask_topsystem utility function by deactivating all cells in the
+    grid.
+    """
+    # Arrange
+    gwf_model = twri_model["GWF_1"]
+    is_active = zeros_like(gwf_model.domain)
+    # Act
+    mask_topsystem(gwf_model, is_active)
+    # Assert
+    for key in ["rch", "drn"]:
+        assert key not in gwf_model.keys()
+

From fc7d39dcc4659a86f6b7a1bd210034355046e35b Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Wed, 23 Sep 2026 17:27:35 +0200
Subject: [PATCH 21/29] Update mock setup

---
 imod/tests/test_mf6/test_mf6_model.py | 22 ++++++++++++++++++++++
 1 file changed, 22 insertions(+)

diff --git a/imod/tests/test_mf6/test_mf6_model.py b/imod/tests/test_mf6/test_mf6_model.py
index b7a3ecd4c..cea2b47ee 100644
--- a/imod/tests/test_mf6/test_mf6_model.py
+++ b/imod/tests/test_mf6/test_mf6_model.py
@@ -263,9 +263,20 @@ def test_clip_box_with_state_for_boundary(
         # Arrange.
         state_for_boundary = MagicMock(spec_set=UgridDataArray)
 
+        idomain = xr.DataArray(
+            np.ones((1, 2, 2), dtype=np.int32), dims=("layer", "y", "x")
+        )
+        top = xr.DataArray(np.ones((2, 2), dtype=np.float64), dims=("y", "x"))
+        bottom = xr.DataArray(np.array([-1.0], dtype=np.float64), dims=("layer",))
+
         discretization_mock = MagicMock(spec_set=Package)
         discretization_mock._pkg_id = "dis"
         discretization_mock.clip_box.return_value = discretization_mock
+        discretization_mock.__getitem__.side_effect = {
+            "idomain": idomain,
+            "top": top,
+            "bottom": bottom,
+        }.__getitem__
 
         clipped_boundary_mock = MagicMock(spec_set=pkg_type)
         clipped_boundary_mock.is_empty.return_value = False
@@ -305,10 +316,21 @@ def test_clip_box_with_unassigned_boundaries_in_original_model(
         # Arrange.
         state_for_boundary = MagicMock(spec_set=UgridDataArray)
 
+        idomain = xr.DataArray(
+            np.ones((1, 2, 2), dtype=np.int32), dims=("layer", "y", "x")
+        )
+        top = xr.DataArray(np.ones((2, 2), dtype=np.float64), dims=("y", "x"))
+        bottom = xr.DataArray(np.array([-1.0], dtype=np.float64), dims=("layer",))
+
         discretization_mock = MagicMock(spec_set=Package)
         discretization_mock._pkg_id = "dis"
         discretization_mock.is_empty.side_effect = [False, False]
         discretization_mock.clip_box.return_value = discretization_mock
+        discretization_mock.__getitem__.side_effect = {
+            "idomain": idomain,
+            "top": top,
+            "bottom": bottom,
+        }.__getitem__
 
         constant_boundary_mock = MagicMock(spec_set=pkg_type)
         constant_boundary_mock.is_empty.side_effect = [False, False]

From 4e060446d94f8390b43b27f0502ec64fac319712 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Wed, 23 Sep 2026 17:28:13 +0200
Subject: [PATCH 22/29] format

---
 imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py b/imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py
index 8993fa7b3..eded019f6 100644
--- a/imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py
+++ b/imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py
@@ -1,7 +1,9 @@
-from imod.mf6.utilities.mask import mask_topsystem
-from imod.typing.grid import zeros_like, ones_like
 import numpy as np
 
+from imod.mf6.utilities.mask import mask_topsystem
+from imod.typing.grid import ones_like, zeros_like
+
+
 def test_mask_topsystem(twri_model):
     """
     Test the mask_topsystem utility function by deactivating all cells in the
@@ -23,7 +25,6 @@ def test_mask_topsystem(twri_model):
         assert np.isnan(first_cell).item()
 
 
-
 def test_mask_topsystem__all_removed(twri_model):
     """
     Test the mask_topsystem utility function by deactivating all cells in the
@@ -37,4 +38,3 @@ def test_mask_topsystem__all_removed(twri_model):
     # Assert
     for key in ["rch", "drn"]:
         assert key not in gwf_model.keys()
-

From cb51c400c64bc6dec428d39a671b2d1e02dc8229 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Thu, 24 Sep 2026 11:51:58 +0200
Subject: [PATCH 23/29] Add ibound to regrid schemes where it was missing and
 slightly improve docstring

---
 imod/mf6/regrid/regrid_schemes.py | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/imod/mf6/regrid/regrid_schemes.py b/imod/mf6/regrid/regrid_schemes.py
index f2a957316..1348488ff 100644
--- a/imod/mf6/regrid/regrid_schemes.py
+++ b/imod/mf6/regrid/regrid_schemes.py
@@ -456,13 +456,15 @@ class StorageCoefficientRegridMethod(DataclassType):
 class CapDataRechargeRegridMethod(DataclassType):
     """
     Object containing regridder methods for CAP data for the
-    :class:`imod.mf6.Recharge.from_imod5_cap_data` method. This contains regridder
-    methods for only the relevant CAP variables for the recharge package.
+    :class:`imod.mf6.Recharge.from_imod5_cap_data` method. This contains
+    regridder methods for only the relevant iMOD5 CAP and BND variables for the
+    recharge package.
     """
 
     boundary: RegridVarType = (RegridderType.OVERLAP, "mode")
     wetted_area: RegridVarType = (RegridderType.RELATIVEOVERLAP, "conductance")
     urban_area: RegridVarType = (RegridderType.RELATIVEOVERLAP, "conductance")
+    ibound: RegridVarType = (RegridderType.OVERLAP, "mode")
 
 
 @dataclass(config=_CONFIG)
@@ -470,8 +472,10 @@ class CapDataWellRegridMethod(DataclassType):
     """
     Object containing regridder methods for CAP data for the
     :class:`imod.mf6.LayeredWell.from_imod5_cap_data` method. This contains
-    regridder methods for only the relevant CAP variables for the well package.
+    regridder methods for only the relevant iMOD5 CAP and BND variables for the
+    well package.
     """
 
     artificial_recharge: RegridVarType = (RegridderType.OVERLAP, "mean")
     artificial_recharge_layer: RegridVarType = (RegridderType.OVERLAP, "mode")
+    ibound: RegridVarType = (RegridderType.OVERLAP, "mode")

From 4776d41dfde075654c7327f09576f7f09b0c8cb4 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Thu, 24 Sep 2026 11:52:18 +0200
Subject: [PATCH 24/29] Include bnd ibound data in test fixtures where missing.

---
 imod/tests/fixtures/imod5_cap_data.py | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/imod/tests/fixtures/imod5_cap_data.py b/imod/tests/fixtures/imod5_cap_data.py
index 8e9b29a9b..3330a4099 100644
--- a/imod/tests/fixtures/imod5_cap_data.py
+++ b/imod/tests/fixtures/imod5_cap_data.py
@@ -82,8 +82,8 @@ def cap_data_sprinkling_grid() -> Imod5DataDict:
         "artificial_recharge_layer": layer,
         "artificial_recharge_capacity": xr.DataArray(25.0),
     }
-
-    return {"cap": cap_data, "extra": {"paths": ["path1", "path2"]}}
+    bnd_data = {"ibound": zeros_grid(n) + 1}
+    return {"cap": cap_data, "bnd": bnd_data, "extra": {"paths": ["path1", "path2"]}}
 
 
 @pytest.fixture(scope="function")
@@ -105,8 +105,9 @@ def cap_data_sprinkling_grid__big() -> Imod5DataDict:
         "artificial_recharge_layer": layer,
         "artificial_recharge_capacity": xr.DataArray(25.0),
     }
+    bnd_data = {"ibound": zeros_dask_grid(n) + 1}
 
-    return {"cap": cap_data, "extra": {"paths": ["path1", "path2"]}}
+    return {"cap": cap_data, "bnd": bnd_data, "extra": {"paths": ["path1", "path2"]}}
 
 
 @pytest.fixture(scope="function")

From 65ec47d31a2acfc76fc3a70c8f4f06aa5b3b9c0a Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Thu, 24 Sep 2026 14:52:52 +0200
Subject: [PATCH 25/29] Include ignore_time_purge_empty in mask_topsystem
 calls. Fix creation of mask for topsystem in clip_box.

---
 imod/mf6/model.py                     | 25 +++++++++++++++----------
 imod/mf6/model_gwf.py                 |  1 +
 imod/mf6/utilities/imod5_converter.py |  3 ++-
 imod/mf6/utilities/mask.py            |  8 ++++++--
 4 files changed, 24 insertions(+), 13 deletions(-)

diff --git a/imod/mf6/model.py b/imod/mf6/model.py
index 70f21325f..4f3400665 100644
--- a/imod/mf6/model.py
+++ b/imod/mf6/model.py
@@ -49,7 +49,6 @@
 from imod.mf6.wel import GridAgnosticWell
 from imod.mf6.write_context import WriteContext
 from imod.schemata import SchemataDict, ValidationError
-from imod.select.grid import active_grid_boundary_xy
 from imod.typing import GridDataArray
 from imod.util.regrid import RegridderWeightsCache
 
@@ -810,18 +809,24 @@ def clip_box(
         clipped_boundary_condition = _create_boundary_condition_clipped_boundary(
             self, clipped, state_for_boundary, clip_box_args
         )
-        state_pkg_id = self._boundary_state_pkg_type._pkg_id
-        pkg_name = f"{state_pkg_id}_clipped"
         if clipped_boundary_condition is not None:
-            clipped[pkg_name] = clipped_boundary_condition
+            # Assign clipped boundary condition package
+            state_pkg_id = self._boundary_state_pkg_type._pkg_id
+            pkg_name = f"{state_pkg_id}_clipped"
 
-            # Clip topsystem packages where the active grid boundary has
-            # changed.
-            _, _, idomain_clipped = clipped._get_domain_geometry()
-            active_bounds_clipped = active_grid_boundary_xy(idomain_clipped > 0)
-            mask_topsystem(clipped, ~active_bounds_clipped)
+            clipped[pkg_name] = clipped_boundary_condition
 
-        clipped.purge_empty_packages(ignore_time=ignore_time_purge_empty)
+            # Mask topsystem packages where the state boundary cells have been
+            # added.
+            state_varname = clipped_boundary_condition._period_data[0]
+            state_var = clipped_boundary_condition.dataset[state_varname].isel(
+                time=0, missing_dims="ignore"
+            )
+            not_added_bc = np.isnan(state_var)
+            # Purge empty packages called by the mask_topsystem function
+            mask_topsystem(clipped, not_added_bc, ignore_time_purge_empty)
+        else:
+            clipped.purge_empty_packages(ignore_time=ignore_time_purge_empty)
 
         return clipped
 
diff --git a/imod/mf6/model_gwf.py b/imod/mf6/model_gwf.py
index 5296bdc62..906d1e4bb 100644
--- a/imod/mf6/model_gwf.py
+++ b/imod/mf6/model_gwf.py
@@ -460,5 +460,6 @@ def from_imod5_data(
             result,
             cast(ConstantHeadRegridMethod, regridder_types.get("topsystem_mask")),
             regrid_cache,
+            ignore_time_purge_empty=True,
         )
         return result
diff --git a/imod/mf6/utilities/imod5_converter.py b/imod/mf6/utilities/imod5_converter.py
index 0da531a2a..81c5e84d7 100644
--- a/imod/mf6/utilities/imod5_converter.py
+++ b/imod/mf6/utilities/imod5_converter.py
@@ -192,6 +192,7 @@ def mask_topsystem_packages_with_ibound(
     model: IModel,
     regridder_types: Optional[ConstantHeadRegridMethod],
     regrid_cache: RegridderWeightsCache,
+    ignore_time_purge_empty: bool,
 ) -> None:
     """
     Mask all top system packages where IBOUND < 0. These locations are assigned
@@ -211,4 +212,4 @@ def mask_topsystem_packages_with_ibound(
     )["ibound"]
     is_active = regridded_ibound >= 0
 
-    mask_topsystem(model, is_active)
+    mask_topsystem(model, is_active, ignore_time_purge_empty)
diff --git a/imod/mf6/utilities/mask.py b/imod/mf6/utilities/mask.py
index 99a969d9d..d734d5518 100644
--- a/imod/mf6/utilities/mask.py
+++ b/imod/mf6/utilities/mask.py
@@ -3,7 +3,9 @@
 from imod.typing import GridDataArray
 
 
-def mask_topsystem(model: IModel, is_active: GridDataArray) -> None:
+def mask_topsystem(
+    model: IModel, is_active: GridDataArray, ignore_time_purge_empty: bool
+) -> None:
     """
     Mask all top system packages in the model inplace with a boolean mask
     indicating active cells.
@@ -15,10 +17,12 @@ def mask_topsystem(model: IModel, is_active: GridDataArray) -> None:
     is_active : GridDataArray
         A boolean array indicating active cells. Top system packages will be masked
         where this array is False.
+    ignore_time_purge_empty : bool
+        If True, ignore the time dimension when masking the packages.
     """
     topsystem_packages = [
         key
         for key, pkg in model.items()
         if isinstance(pkg, ITopSystemBoundaryCondition)
     ]
-    model.mask_packages(topsystem_packages, is_active)
+    model.mask_packages(topsystem_packages, is_active, ignore_time_purge_empty)

From 50f71e5f5965e7a8950b99ccc9f245ea80a435f8 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Thu, 24 Sep 2026 15:16:39 +0200
Subject: [PATCH 26/29] Refactor: Move boundary condition creator utiltity
 functions from model module to their respective

---
 imod/common/interfaces/imodel.py              |   6 +
 imod/mf6/model.py                             |  90 +-------------
 imod/mf6/utilities/clipped_bc_creator.py      | 113 +++++++++++++++++-
 ..._mf6_clipped_boundary_condition_creator.py |   6 +-
 4 files changed, 122 insertions(+), 93 deletions(-)

diff --git a/imod/common/interfaces/imodel.py b/imod/common/interfaces/imodel.py
index 0b30e0a6e..fd143feb7 100644
--- a/imod/common/interfaces/imodel.py
+++ b/imod/common/interfaces/imodel.py
@@ -2,6 +2,7 @@
 from typing import Any, Optional, Tuple
 
 from imod.common.interfaces.idict import IDict
+from imod.common.interfaces.ipackage import IPackage
 from imod.common.statusinfo import StatusInfoBase
 from imod.mf6.validation_settings import ValidationSettings
 from imod.typing import GridDataArray
@@ -68,3 +69,8 @@ def _is_splitting_supported(self) -> Tuple[bool, str]:
     @abstractmethod
     def _is_clipping_supported(self) -> Tuple[bool, str]:
         raise NotImplementedError
+
+    @property
+    @abstractmethod
+    def _boundary_state_pkg_type(self) -> type[IPackage]:
+        raise NotImplementedError
diff --git a/imod/mf6/model.py b/imod/mf6/model.py
index 4f3400665..83bfa725d 100644
--- a/imod/mf6/model.py
+++ b/imod/mf6/model.py
@@ -20,7 +20,6 @@
 from imod.common.interfaces.imodel import IModel
 from imod.common.serializer import EngineType
 from imod.common.statusinfo import NestedStatusInfo, StatusInfo, StatusInfoBase
-from imod.common.utilities.clip import clip_box_dataset
 from imod.common.utilities.dump_model import dump_model
 from imod.common.utilities.mask import mask_packages
 from imod.common.utilities.regrid import _regrid_like
@@ -40,8 +39,7 @@
 from imod.mf6.riv import River
 from imod.mf6.utilities.clipped_bc_creator import (
     StateClassType,
-    StateType,
-    create_clipped_boundary,
+    create_boundary_condition_clipped_boundary,
 )
 from imod.mf6.utilities.mask import mask_topsystem
 from imod.mf6.utilities.mf6hfb import merge_hfb_packages
@@ -63,90 +61,6 @@ def pkg_has_cleanup(pkg: Package):
     return any(isinstance(pkg, pkgtype) for pkgtype in PKGTYPES_WITH_CLEANUP)
 
 
-def _create_boundary_condition_for_unassigned_boundary(
-    model: Modflow6Model,
-    state_for_boundary: Optional[GridDataArray],
-    additional_boundaries: list[Optional[StateType]] = [None],
-) -> Optional[StateType]:
-    if state_for_boundary is None:
-        return None
-
-    pkg_type = model._boundary_state_pkg_type
-    constant_state_packages = [
-        pkg for _, pkg in model.items() if isinstance(pkg, pkg_type)
-    ]
-
-    filtered_boundaries: list[StateType] = [
-        item for item in additional_boundaries or [] if item is not None
-    ]
-
-    constant_state_packages.extend(filtered_boundaries)
-
-    return create_clipped_boundary(
-        model.domain, state_for_boundary, constant_state_packages, pkg_type
-    )
-
-
-def _create_boundary_condition_clipped_boundary(
-    original_model: Modflow6Model,
-    clipped_model: Modflow6Model,
-    state_for_boundary: Optional[GridDataArray],
-    clip_box_args: tuple[Any, ...],
-) -> Optional[StateType]:
-    # Create temporary boundary condition for the original model boundary. This
-    # is used later to see which boundaries can be ignored as they were already
-    # present in the original model. We want to just end up with the boundary
-    # created by the clip.
-    unassigned_boundary_original_domain = (
-        _create_boundary_condition_for_unassigned_boundary(
-            original_model, state_for_boundary
-        )
-    )
-    # Clip the unassigned boundary to the clipped model's domain, required to
-    # avoid topological errors later.
-    if unassigned_boundary_original_domain is not None:
-        unassigned_boundary_clipped = unassigned_boundary_original_domain.clip_box(
-            *clip_box_args
-        )
-    else:
-        unassigned_boundary_clipped = None
-
-    if state_for_boundary is not None:
-        # Clip box as dataset, temporarily add variable name to convert to
-        # dataset, then turn back into DataArray.
-        varname = original_model._boundary_state_pkg_type._period_data[0]
-        state_for_boundary = state_for_boundary.to_dataset(name=varname)
-        state_for_boundary_clipped = clip_box_dataset(
-            state_for_boundary, *clip_box_args
-        )[varname]
-    else:
-        state_for_boundary_clipped = None
-
-    bc_constant_pkg = _create_boundary_condition_for_unassigned_boundary(
-        clipped_model, state_for_boundary_clipped, [unassigned_boundary_clipped]
-    )
-
-    # Remove all indices before first timestep of state_for_clipped_boundary.
-    # This to prevent empty dataarrays unnecessarily being made for these
-    # indices, which can lead to them to be removed when purging empty packages
-    # with ignore_time=True. Unfortunately, this is needs to be handled here and
-    # not in _create_boundary_condition_for_unassigned_boundary, as otherwise
-    # this function is called twice which could result in broadcasting errors in
-    # the second call if the time domain of state_for_boundary and assigned
-    # packages have no overlap.
-    if (
-        (state_for_boundary is not None)
-        and (state_for_boundary.indexes.get("time") is not None)
-        and (bc_constant_pkg is not None)
-    ):
-        start_time = state_for_boundary.indexes["time"][0]
-        bc_constant_pkg.dataset = bc_constant_pkg.dataset.sel(
-            time=slice(start_time, None)
-        )
-
-    return bc_constant_pkg
-
-
 class Modflow6Model(collections.UserDict[str, Package], IModel, abc.ABC):
     _mandatory_packages: tuple[str, ...] = ()
     _init_schemata: SchemataDict = {}
@@ -806,7 +720,7 @@ def clip_box(
             *clip_box_args,
         )
 
-        clipped_boundary_condition = _create_boundary_condition_clipped_boundary(
+        clipped_boundary_condition = create_boundary_condition_clipped_boundary(
             self, clipped, state_for_boundary, clip_box_args
         )
         if clipped_boundary_condition is not None:
diff --git a/imod/mf6/utilities/clipped_bc_creator.py b/imod/mf6/utilities/clipped_bc_creator.py
index ea6f0f415..c666ae261 100644
--- a/imod/mf6/utilities/clipped_bc_creator.py
+++ b/imod/mf6/utilities/clipped_bc_creator.py
@@ -1,7 +1,9 @@
-from typing import Optional, Tuple, TypeAlias
+from typing import Any, Optional, Tuple, TypeAlias, cast
 
 import xarray as xr
 
+from imod.common.interfaces.imodel import IModel
+from imod.common.utilities.clip import clip_box_dataset
 from imod.mf6 import ConstantConcentration, ConstantHead
 from imod.select.grid import active_grid_boundary_xy
 from imod.typing import GridDataArray
@@ -117,7 +119,7 @@ def _create_clipped_boundary_state(
     return state_for_clipped_boundary.where(unassigned_grid_boundaries)
 
 
-def create_clipped_boundary(
+def _create_clipped_boundary_pkg(
     idomain: GridDataArray,
     state_for_clipped_boundary: GridDataArray,
     original_constant_head_boundaries: list[StateType],
@@ -152,3 +154,110 @@ def create_clipped_boundary(
     )
 
     return pkg_type(constant_state, print_input=True, print_flows=True, save_flows=True)
+
+
+def _create_boundary_condition_for_unassigned_boundary(
+    model: IModel,
+    state_for_boundary: Optional[GridDataArray],
+    additional_boundaries: list[Optional[StateType]] = [None],
+) -> Optional[StateType]:
+    if state_for_boundary is None:
+        return None
+
+    pkg_type = cast(StateClassType, model._boundary_state_pkg_type)
+    constant_state_packages = [
+        pkg for _, pkg in model.items() if isinstance(pkg, pkg_type)
+    ]
+
+    filtered_boundaries: list[StateType] = [
+        item for item in additional_boundaries or [] if item is not None
+    ]
+
+    constant_state_packages.extend(filtered_boundaries)
+
+    return _create_clipped_boundary_pkg(
+        model.domain, state_for_boundary, constant_state_packages, pkg_type
+    )
+
+
+def create_boundary_condition_clipped_boundary(
+    original_model: IModel,
+    clipped_model: IModel,
+    state_for_boundary: Optional[GridDataArray],
+    clip_box_args: tuple[Any, ...],
+) -> Optional[StateType]:
+    """
+    Create a clipped boundary condition for a given state in the clipped model.
+    The function takes the original model as a reference to determine where
+    boundary conditions should NOT be placed, then applies this information to
+    create the boundary condition in the clipped model.
+
+    Parameters
+    ----------
+    original_model : IModel
+        The original model containing the unassigned boundary condition.
+    clipped_model : IModel
+        The clipped model where the boundary condition will be applied.
+    state_for_boundary : Optional[GridDataArray]
+        The state array for the boundary condition.
+    clip_box_args : tuple[Any, ...]
+        Arguments defining the clipping box.
+
+    Returns
+    -------
+    Optional[StateType]
+        The clipped boundary condition package, or None if no boundary condition is created.
+    """
+    # Create temporary boundary condition for the original model boundary. This
+    # is used later to see which boundaries can be ignored as they were already
+    # present in the original model. We want to just end up with the boundary
+    # created by the clip.
+    unassigned_boundary_original_domain = (
+        _create_boundary_condition_for_unassigned_boundary(
+            original_model, state_for_boundary
+        )
+    )
+    # Clip the unassigned boundary to the clipped model's domain, required to
+    # avoid topological errors later.
+    if unassigned_boundary_original_domain is not None:
+        unassigned_boundary_clipped = unassigned_boundary_original_domain.clip_box(
+            *clip_box_args
+        )
+    else:
+        unassigned_boundary_clipped = None
+
+    if state_for_boundary is not None:
+        # Clip box as dataset, temporarily add variable name to convert to
+        # dataset, then turn back into DataArray.
+        state_cls = cast(StateClassType, original_model._boundary_state_pkg_type)
+        varname = state_cls._period_data[0]
+        state_for_boundary = state_for_boundary.to_dataset(name=varname)
+        state_for_boundary_clipped = clip_box_dataset(
+            state_for_boundary, *clip_box_args
+        )[varname]
+    else:
+        state_for_boundary_clipped = None
+
+    bc_constant_pkg = _create_boundary_condition_for_unassigned_boundary(
+        clipped_model, state_for_boundary_clipped, [unassigned_boundary_clipped]
+    )
+
+    # Remove all indices before first timestep of state_for_clipped_boundary.
+    # This to prevent empty dataarrays unnecessarily being made for these
+    # indices, which can lead to them to be removed when purging empty packages
+    # with ignore_time=True. Unfortunately, this is needs to be handled here and
+    # not in _create_boundary_condition_for_unassigned_boundary, as otherwise
+    # this function is called twice which could result in broadcasting errors in
+    # the second call if the time domain of state_for_boundary and assigned
+    # packages have no overlap.
+    if (
+        (state_for_boundary is not None)
+        and (state_for_boundary.indexes.get("time") is not None)
+        and (bc_constant_pkg is not None)
+    ):
+        start_time = state_for_boundary.indexes["time"][0]
+        bc_constant_pkg.dataset = bc_constant_pkg.dataset.sel(
+            time=slice(start_time, None)
+        )
+
+    return bc_constant_pkg
diff --git a/imod/tests/test_mf6/test_utilities/test_mf6_clipped_boundary_condition_creator.py b/imod/tests/test_mf6/test_utilities/test_mf6_clipped_boundary_condition_creator.py
index aa599ebca..2be514194 100644
--- a/imod/tests/test_mf6/test_utilities/test_mf6_clipped_boundary_condition_creator.py
+++ b/imod/tests/test_mf6/test_utilities/test_mf6_clipped_boundary_condition_creator.py
@@ -5,7 +5,7 @@
 
 from imod.mf6 import ConstantHead
 from imod.mf6.utilities.clipped_bc_creator import (
-    create_clipped_boundary,
+    _create_clipped_boundary_pkg,
 )
 from imod.select.grid import grid_boundary_xy
 
@@ -55,7 +55,7 @@ def test_create_different_n_clipped_cells(self, circle_dis, n_clipped_cells):
         )
 
         # Act.
-        constant_head_pkg_clipped_domain = create_clipped_boundary(
+        constant_head_pkg_clipped_domain = _create_clipped_boundary_pkg(
             idomain,
             clipped_boundary_values,
             [reduced_boundary_constant_head_pkg],
@@ -104,7 +104,7 @@ def test_create_different_dis(self, dis, grid_data_array, request):
         )
 
         # Act.
-        constant_head_pkg_clipped_domain = create_clipped_boundary(
+        constant_head_pkg_clipped_domain = _create_clipped_boundary_pkg(
             idomain,
             clipped_boundary_values,
             [reduced_boundary_constant_head_pkg],

From d8b191ae8912ed1871afbdc28d6c51627bb0b8fb Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Thu, 24 Sep 2026 16:06:09 +0200
Subject: [PATCH 27/29] Update mocking framework

---
 imod/tests/test_mf6/test_mf6_model.py | 34 ++++++++++++++++++++-------
 1 file changed, 26 insertions(+), 8 deletions(-)

diff --git a/imod/tests/test_mf6/test_mf6_model.py b/imod/tests/test_mf6/test_mf6_model.py
index cea2b47ee..42f19e85b 100644
--- a/imod/tests/test_mf6/test_mf6_model.py
+++ b/imod/tests/test_mf6/test_mf6_model.py
@@ -102,12 +102,18 @@ def test_circle_roundtrip(circle_model, tmp_path):
     roundtrip(circle_model["GWF_1"], tmp_path)
 
 
+class ConcreteModflow6Model(Modflow6Model):
+    """Concrete implementation of the abstract Modflow6Model for testing purposes."""
+
+    _boundary_state_pkg_type = ConstantHead
+
+
 class TestModel:
     def test_write_valid_model_without_error(self, tmpdir_factory):
         # Arrange.
         tmp_path = tmpdir_factory.mktemp("TestSimulation")
         model_name = "Test model"
-        model = Modflow6Model()
+        model = ConcreteModflow6Model()
         # create write context
         validation_context = ValidationSettings()
         write_context = WriteContext(tmp_path)
@@ -135,7 +141,7 @@ def test_write_without_dis_pkg_return_error(self, tmpdir_factory):
         # Arrange.
         tmp_path = tmpdir_factory.mktemp("TestSimulation")
         model_name = "Test model"
-        model = Modflow6Model()
+        model = ConcreteModflow6Model()
         # create write context
         validation_context = ValidationSettings()
         write_context = WriteContext(tmp_path)
@@ -158,7 +164,7 @@ def test_write_with_invalid_pkg_returns_error(self, tmpdir_factory):
         # Arrange.
         tmp_path = tmpdir_factory.mktemp("TestSimulation")
         model_name = "Test model"
-        model = Modflow6Model()
+        model = ConcreteModflow6Model()
         # create write context
         validation_context = ValidationSettings()
         write_context = WriteContext(tmp_path)
@@ -192,7 +198,7 @@ def test_write_with_two_invalid_pkg_returns_two_errors(self, tmpdir_factory):
         validation_context = ValidationSettings()
         write_context = WriteContext(simulation_directory=tmp_path)
 
-        model = Modflow6Model()
+        model = ConcreteModflow6Model()
 
         discretization_mock = MagicMock(spec_set=Package)
         discretization_mock._pkg_id = "dis"
@@ -249,7 +255,8 @@ def test_clip_box_without_state_for_boundary(self, model_type, pkg_type):
         pkg_id = pkg_type._pkg_id
         assert f"{pkg_id}_clipped" not in clipped
 
-    @mock.patch("imod.mf6.model.create_clipped_boundary")
+    @mock.patch("imod.mf6.model.mask_topsystem")
+    @mock.patch("imod.mf6.utilities.clipped_bc_creator._create_clipped_boundary_pkg")
     @pytest.mark.parametrize(
         "model_type, pkg_type",
         [
@@ -258,7 +265,11 @@ def test_clip_box_without_state_for_boundary(self, model_type, pkg_type):
         ],
     )
     def test_clip_box_with_state_for_boundary(
-        self, create_clipped_boundary_mock, model_type, pkg_type
+        self,
+        create_clipped_boundary_mock,
+        mask_topsystem_mock,
+        model_type,
+        pkg_type,
     ):
         # Arrange.
         state_for_boundary = MagicMock(spec_set=UgridDataArray)
@@ -301,8 +312,10 @@ def test_clip_box_with_state_for_boundary(
             [],
             pkg_type,
         )
+        mask_topsystem_mock.assert_called_once()
 
-    @mock.patch("imod.mf6.model.create_clipped_boundary")
+    @mock.patch("imod.mf6.model.mask_topsystem")
+    @mock.patch("imod.mf6.utilities.clipped_bc_creator._create_clipped_boundary_pkg")
     @pytest.mark.parametrize(
         "model_type, pkg_type",
         [
@@ -311,7 +324,11 @@ def test_clip_box_with_state_for_boundary(
         ],
     )
     def test_clip_box_with_unassigned_boundaries_in_original_model(
-        self, create_clipped_boundary_mock, model_type, pkg_type
+        self,
+        create_clipped_boundary_mock,
+        mask_topsystem_mock,
+        model_type,
+        pkg_type,
     ):
         # Arrange.
         state_for_boundary = MagicMock(spec_set=UgridDataArray)
@@ -361,6 +378,7 @@ def test_clip_box_with_unassigned_boundaries_in_original_model(
             [constant_boundary_mock, unassigned_original_constant_boundary.clip_box()],
             pkg_type,
         )
+        mask_topsystem_mock.assert_called_once()
 
 
 class TestGroundwaterFlowModel:

From 0f41608e526be423db4eb7e5570132d4d2d8ccd5 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Thu, 24 Sep 2026 16:06:21 +0200
Subject: [PATCH 28/29] Update missing args

---
 imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py b/imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py
index eded019f6..167d6d1d4 100644
--- a/imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py
+++ b/imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py
@@ -16,7 +16,7 @@ def test_mask_topsystem(twri_model):
     is_active[0, 0, 0] = 0
 
     # Act
-    mask_topsystem(gwf_model, is_active)
+    mask_topsystem(gwf_model, is_active, True)
     # Assert
     for key in ["rch", "drn"]:
         pkg = gwf_model[key]
@@ -34,7 +34,7 @@ def test_mask_topsystem__all_removed(twri_model):
     gwf_model = twri_model["GWF_1"]
     is_active = zeros_like(gwf_model.domain)
     # Act
-    mask_topsystem(gwf_model, is_active)
+    mask_topsystem(gwf_model, is_active, True)
     # Assert
     for key in ["rch", "drn"]:
         assert key not in gwf_model.keys()

From 497f94a55b4563df15667732e74ad0e29b4372bd Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Thu, 24 Sep 2026 16:12:17 +0200
Subject: [PATCH 29/29] Drop time coord and add docstring

---
 imod/mf6/model.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/imod/mf6/model.py b/imod/mf6/model.py
index 83bfa725d..6d783fd40 100644
--- a/imod/mf6/model.py
+++ b/imod/mf6/model.py
@@ -733,8 +733,10 @@ def clip_box(
             # Mask topsystem packages where the state boundary cells have been
             # added.
             state_varname = clipped_boundary_condition._period_data[0]
+            # Select the state variable for the first time step as mask.
+            # Its location will be constant through time.
             state_var = clipped_boundary_condition.dataset[state_varname].isel(
-                time=0, missing_dims="ignore"
+                time=0, missing_dims="ignore", drop=True
             )
             not_added_bc = np.isnan(state_var)
             # Purge empty packages called by the mask_topsystem function