diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst
index a5b40a49a..1ecc7c0d5 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
 ~~~~~
@@ -25,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
 ~~~~~~~
@@ -32,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
 --------------------
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/iboundarycondition.py b/imod/common/interfaces/iboundarycondition.py
new file mode 100644
index 000000000..7db2f9db2
--- /dev/null
+++ b/imod/common/interfaces/iboundarycondition.py
@@ -0,0 +1,12 @@
+from abc import ABC
+from typing import ClassVar
+
+from imod.common.interfaces.ipackage import IPackage
+
+
+class IBoundaryCondition(IPackage, ABC):
+    """
+    Interface for boundary condition packages with stress-period variables.
+    """
+
+    _period_data: ClassVar[tuple[str, ...]]
diff --git a/imod/common/interfaces/imodel.py b/imod/common/interfaces/imodel.py
index 2fe91e0a4..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
@@ -16,9 +17,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
 
@@ -56,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/common/interfaces/itopsystembc.py b/imod/common/interfaces/itopsystembc.py
new file mode 100644
index 000000000..de541ec09
--- /dev/null
+++ b/imod/common/interfaces/itopsystembc.py
@@ -0,0 +1,15 @@
+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
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/common/utilities/regrid.py b/imod/common/utilities/regrid.py
index aaf28abb9..94ab76f75 100644
--- a/imod/common/utilities/regrid.py
+++ b/imod/common/utilities/regrid.py
@@ -27,7 +27,7 @@
     is_unstructured,
     ones_like,
 )
-from imod.util.dims import drop_layer_dim_cap_data, enforced_dim_order
+from imod.util.dims import drop_layer_dim_cap_and_bnd_data, enforced_dim_order
 from imod.util.regrid import (
     RegridderType,
     RegridderWeightsCache,
@@ -471,7 +471,7 @@ def _get_regridding_domain(
     return new_idomain
 
 
-def regrid_imod5_cap_data(
+def regrid_imod5_cap_and_bnd_data(
     imod5_data: Imod5DataDict,
     target_dis: IRegridPackage,
     regridder_types: DataclassType,
@@ -486,15 +486,19 @@ 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_and_bnd_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_no_layer["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/boundary_condition.py b/imod/mf6/boundary_condition.py
index 882ee10df..468c18dc6 100644
--- a/imod/mf6/boundary_condition.py
+++ b/imod/mf6/boundary_condition.py
@@ -7,6 +7,7 @@
 import xarray as xr
 import xugrid as xu
 
+from imod.common.interfaces.iboundarycondition import IBoundaryCondition
 from imod.common.utilities.value_filters import enforce_scalar
 from imod.mf6.auxiliary_variables import (
     expand_transient_auxiliary_variables,
@@ -57,7 +58,7 @@ def _disv_recarr(arrdict, layer, notnull):
     return recarr
 
 
-class BoundaryCondition(Package, abc.ABC):
+class BoundaryCondition(Package, IBoundaryCondition, abc.ABC):
     """
     BoundaryCondition is used to share methods for specific stress packages
     with a time component.
@@ -324,7 +325,7 @@ def _get_period_varnames(self) -> list[str]:
         >>> river._get_period_varnames()
         >>> # prints: ['stage', 'conductance', 'bottom_elevation', 'species1', 'species2']
         """
-        result = []
+        result: list[str] = []
         if hasattr(self, "_period_data"):
             result.extend(self._period_data)
         if hasattr(self, "_optional_data"):
diff --git a/imod/mf6/model.py b/imod/mf6/model.py
index 2a931732f..74cf37198 100644
--- a/imod/mf6/model.py
+++ b/imod/mf6/model.py
@@ -20,9 +20,8 @@
 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_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,
@@ -40,9 +39,9 @@
 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_where_bc
 from imod.mf6.utilities.mf6hfb import merge_hfb_packages
 from imod.mf6.validation_settings import ValidationSettings
 from imod.mf6.wel import GridAgnosticWell
@@ -62,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 = {}
@@ -805,15 +720,21 @@ 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
         )
-        state_pkg_id = self._boundary_state_pkg_type._pkg_id
-        pkg_name = f"{state_pkg_id}_clipped"
         if clipped_boundary_condition is not None:
+            # Assign clipped boundary condition package
+            state_pkg_id = self._boundary_state_pkg_type._pkg_id
+            pkg_name = f"{state_pkg_id}_clipped"
+
             clipped[pkg_name] = clipped_boundary_condition
 
-        clipped.purge_empty_packages(ignore_time=ignore_time_purge_empty)
+            mask_topsystem_where_bc(
+                clipped, clipped_boundary_condition, ignore_time_purge_empty
+            )
+        else:
+            clipped.purge_empty_packages(ignore_time=ignore_time_purge_empty)
 
         return clipped
 
@@ -930,11 +851,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)
+
+    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.
 
-        mask_all_packages(self, mask, ignore_time_purge_empty)
+        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 +899,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"
diff --git a/imod/mf6/model_gwf.py b/imod/mf6/model_gwf.py
index 545f948da..906d1e4bb 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_with_ibound
 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,12 @@ def from_imod5_data(
             for key, chd_package in chd_packages.items():
                 result[key] = chd_package
 
+        # Mask all topsystem packages where IBOUND < 0
+        mask_topsystem_packages_with_ibound(
+            imod5_data,
+            result,
+            cast(ConstantHeadRegridMethod, regridder_types.get("topsystem_mask")),
+            regrid_cache,
+            ignore_time_purge_empty=True,
+        )
         return result
diff --git a/imod/mf6/rch.py b/imod/mf6/rch.py
index 1933f480e..731676776 100644
--- a/imod/mf6/rch.py
+++ b/imod/mf6/rch.py
@@ -6,7 +6,7 @@
 
 from imod.common.interfaces.iregridpackage import IRegridPackage
 from imod.common.utilities.dataclass_type import DataclassType
-from imod.common.utilities.regrid import regrid_imod5_cap_data
+from imod.common.utilities.regrid import regrid_imod5_cap_and_bnd_data
 from imod.logging import init_log_decorator
 from imod.mf6.aggregate.aggregate_schemes import RechargeAggregationMethod
 from imod.mf6.dis import StructuredDiscretization, VerticesDiscretization
@@ -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_and_bnd_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/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")
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/clipped_bc_creator.py b/imod/mf6/utilities/clipped_bc_creator.py
index ea6f0f415..b8456cf64 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
@@ -15,6 +17,9 @@ def _find_unassigned_grid_boundaries(
     active_grid_boundary: GridDataArray,
     boundary_conditions: list[StateType],
 ) -> GridDataArray:
+    """
+    Find the grid boundaries that have not been assigned any boundary conditions.
+    """
     unassigned_grid_boundaries = active_grid_boundary.copy()
     for boundary_condition in boundary_conditions:
         # Fetch variable name from the first boundary condition, can be "head" or
@@ -103,7 +108,7 @@ def _create_clipped_boundary_state(
     idomain: GridDataArray,
     state_for_clipped_boundary: GridDataArray,
     original_constant_head_boundaries: list[StateType],
-):
+) -> GridDataArray:
     """Helper function to make sure dimension order is enforced"""
     active_grid_boundary = active_grid_boundary_xy(idomain > 0)
     unassigned_grid_boundaries = _find_unassigned_grid_boundaries(
@@ -117,7 +122,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 +157,118 @@ 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]:
+    """
+    Create a boundary condition for the unnassigned boundary cells of a model.
+
+    Constant state packages are collected from the model and any additional
+    boundaries provided will be added to this list. Cells which are not covered
+    by any of these packages will have a new boundary condition created for
+    them.
+    """
+    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 boundary condition for the clipped model.
+
+    The function takes the original model as a reference to determine where
+    boundary conditions should NOT be placed.
+
+    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 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/mf6/utilities/imod5_converter.py b/imod/mf6/utilities/imod5_converter.py
index e27c21581..0eac5e611 100644
--- a/imod/mf6/utilities/imod5_converter.py
+++ b/imod/mf6/utilities/imod5_converter.py
@@ -1,13 +1,19 @@
-from typing import Optional, Union
+from typing import Optional, Union, cast
 
 import numpy as np
 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.common.utilities.regrid import (
+    _regrid_package_data,
+    regrid_imod5_cap_and_bnd_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
@@ -16,6 +22,22 @@
 def convert_ibound_to_idomain(
     ibound: xr.DataArray, thickness: xr.DataArray
 ) -> xr.DataArray:
+    """
+    Convert IBOUND array to IDOMAIN array. IBOUND -1 will be set to 1 in
+    IDOMAIN. When the thickness is <= 0, IDOMAIN will be set to -1.
+
+    Parameters
+    ----------
+    ibound : xr.DataArray
+        The IBOUND array from iMOD5.
+    thickness : xr.DataArray
+        The thickness array of the model layers.
+
+    Returns
+    -------
+    xr.DataArray
+        The corresponding IDOMAIN array.
+    """
     # Convert IBOUND to IDOMAIN
     # -1 to 1, these will have to be filled with
     # CHD cells.
@@ -121,6 +143,26 @@ def well_from_imod5_cap_data(
         ``"artificial_recharge_capacity"`` is ignored as the abstraction
         capacity is already defined in the point data. This is an ``n:1``
         mapping: multiple grid cells can map to one well.
+
+    Parameters
+    ----------
+    imod5_data : Imod5DataDict
+        The iMOD5 data containing the "cap" package with abstraction
+        information.
+    target_dis : Optional[IRegridPackage]
+        The target discretization package for regridding the data. Required if
+        the data is in grid format (IDF).
+    regridder_types : DataclassType
+        The regrid methods to use for regridding the data.
+    regrid_cache : RegridderWeightsCache
+        Cache for storing regridder weights to speed up repeated regridding
+        operations.
+
+    Returns
+    -------
+    dict[str, np.ndarray]
+        A dictionary containing well information extracted from the iMOD5 cap
+        data.
     """
     cap_data = imod5_data["cap"]
     has_ipf_well = isinstance(cap_data["artificial_recharge_layer"], pd.DataFrame)
@@ -133,14 +175,14 @@ def well_from_imod5_cap_data(
                 "target_dis must be provided when converting iMOD5 cap data "
                 "from grids (IDF)"
             )
-        cap_data_regridded = regrid_imod5_cap_data(
+        cap_data_regridded = regrid_imod5_cap_and_bnd_data(
             imod5_data, target_dis, regridder_types, regrid_cache
         )["cap"]
         return _well_from_imod5_cap_grid_data(cap_data_regridded)
 
 
 def regrid_imod5_pkg_data(
-    cls: type[Package],
+    pkg_type: Optional[type[Package]],
     imod5_pkg_data: GridDataDict,
     target_dis: Package,
     regridder_types: Optional[DataclassType],
@@ -149,12 +191,42 @@ def regrid_imod5_pkg_data(
     """
     Regrid iMOD5 package data to target idomain. Optionally get regrid methods
     from class if not provided.
+
+    Parameters
+    ----------
+    pkg_type:
+        The type of the package being regridded. This is used to determine the
+        appropriate regrid methods if regridder_types is not provided.
+    imod5_pkg_data:
+        The iMOD5 package data to be regridded.
+    target_dis:
+        The target discretization package containing the idomain to regrid to.
+    regridder_types:
+        Optional regrid methods to use for regridding. If not provided, they
+        will be obtained from the pkg_type.
+    regrid_cache:
+        Cache for storing regridder weights to speed up repeated regridding
+        operations.
+
+    Returns
+    -------
+    GridDataDict
+        The regridded iMOD5 package data.
     """
+    if (pkg_type is None) and (regridder_types is None):
+        raise ValueError(
+            "Either pkg_type or regridder_types must be provided for regridding."
+        )
+    # set up regridder methods
+    elif (pkg_type is not None) and (
+        regridder_types is None
+    ):  # check pkg_type not None for mypy
+        regridder_types = pkg_type.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, {}
@@ -165,6 +237,22 @@ def regrid_imod5_pkg_data(
 def chd_cells_from_imod5_data(
     imod5_pkg_data: GridDataDict, target_idomain: GridDataArray
 ) -> GridDataDict:
+    """
+    Get CHD cells from iMOD5 package data based on IBOUND and target idomain.
+
+    Parameters
+    ----------
+    imod5_pkg_data:
+        The iMOD5 package data containing "head" and "ibound".
+    target_idomain:
+        The target idomain to filter active cells.
+
+    Returns
+    -------
+    GridDataDict
+        The filtered CHD cells with "head" values where IBOUND < 0 and target
+        idomain > 0.
+    """
     head = imod5_pkg_data["head"]
     ibound = imod5_pkg_data["ibound"]
 
@@ -175,3 +263,51 @@ def chd_cells_from_imod5_data(
     head = head.where(target_idomain > 0)
 
     return {"head": head}
+
+
+def mask_topsystem_packages_with_ibound(
+    imod5_data: dict[str, dict[str, GridDataArray]],
+    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
+    a constant head.
+
+    Parameters
+    ----------
+    imod5_data:
+        The iMOD5 data containing the "bnd" package with "ibound".
+    model:
+        The target MODFLOW 6 model.
+    regridder_types:
+        Optional regrid methods to use for regridding. If not provided, default
+        methods will be used.
+    regrid_cache:
+        Cache for storing regridder weights to speed up repeated regridding
+        operations.
+    ignore_time_purge_empty:
+        Flag indicating whether to ignore time when purging empty cells.
+
+    Returns
+    -------
+    None
+        The function modifies the top system packages in the model in place.
+    """
+
+    if regridder_types is None:
+        regridder_types = ConstantHeadRegridMethod()
+
+    ibound = imod5_data["bnd"]["ibound"]
+    regridded_ibound = regrid_imod5_pkg_data(
+        pkg_type=None,
+        imod5_pkg_data={"ibound": ibound},
+        target_dis=model["dis"],
+        regridder_types=regridder_types,
+        regrid_cache=regrid_cache,
+    )["ibound"]
+    is_active = regridded_ibound >= 0
+
+    mask_topsystem(model, is_active, ignore_time_purge_empty)
diff --git a/imod/mf6/utilities/mask.py b/imod/mf6/utilities/mask.py
new file mode 100644
index 000000000..d2d6f13c8
--- /dev/null
+++ b/imod/mf6/utilities/mask.py
@@ -0,0 +1,61 @@
+import numpy as np
+
+from imod.common.interfaces.iboundarycondition import IBoundaryCondition
+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, ignore_time_purge_empty: bool
+) -> None:
+    """
+    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.
+    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, ignore_time_purge_empty)
+
+
+def mask_topsystem_where_bc(
+    model: IModel, boundary_condition: IBoundaryCondition, ignore_time: bool
+) -> None:
+    """
+    Mask all top system packages in the model inplace where the boundary
+    condition is active.
+
+    Parameters
+    ----------
+    model : IModel
+        The MODFLOW 6 model containing top system packages.
+    boundary_condition : IBoundaryCondition
+        The boundary condition package used to determine which cells are not
+        added.
+    ignore_time : bool
+        If True, ignore the time dimension when masking the packages and set the
+        mask to where the first time step of the boundary condition contains
+        active cells. Else, aggregate over all time steps to determine the mask.
+    """
+    state_varname = boundary_condition._period_data[0]
+    state_var = boundary_condition.dataset[state_varname]
+    if "time" in state_var.dims:
+        if ignore_time:
+            state_var = state_var.isel(time=0, drop=True)
+        else:
+            state_var = state_var.min(dim="time")
+    not_added_bc = np.isnan(state_var)
+    mask_topsystem(model, not_added_bc, ignore_time)
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/model.py b/imod/msw/model.py
index 2a4081296..19832db6e 100644
--- a/imod/msw/model.py
+++ b/imod/msw/model.py
@@ -19,7 +19,7 @@
 from imod.common.utilities.clip import clip_by_grid
 from imod.common.utilities.dump_model import dump_model
 from imod.common.utilities.partitioninfo import create_partition_info
-from imod.common.utilities.regrid import regrid_imod5_cap_data
+from imod.common.utilities.regrid import regrid_imod5_cap_and_bnd_data
 from imod.common.utilities.version import prepend_content_with_version_info
 from imod.mf6.dis import StructuredDiscretization
 from imod.mf6.mf6_wel_adapter import Mf6Wel
@@ -823,7 +823,7 @@ def from_imod5_data(
         parasim_settings = read_para_sim(path_to_parasim)
         unsa_svat_path = cast(str, parasim_settings["unsa_svat_path"])
         # Regrid iMOD5 CAP data to target discretization.
-        imod5_regridded = regrid_imod5_cap_data(
+        imod5_regridded = regrid_imod5_cap_and_bnd_data(
             imod5_data, target_dis, regridder_types, regrid_cache
         )
         # Test with regridded data instead of masked, as masking broadcasts
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..7f18a4acf 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,10 @@ 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)
+    # 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")
     return MetaSwapActive(active, subunit_active)
 
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")
diff --git a/imod/tests/fixtures/msw_imod5_cap_fixture.py b/imod/tests/fixtures/msw_imod5_cap_fixture.py
index 28c53335b..3d46a3bda 100644
--- a/imod/tests/fixtures/msw_imod5_cap_fixture.py
+++ b/imod/tests/fixtures/msw_imod5_cap_fixture.py
@@ -19,10 +19,10 @@ def imod5_cap_data() -> GridDataDict:
     da_kwargs["coords"] = {"layer": layer, "y": y, "x": x, "dx": dx, "dy": dy}
 
     imod5_data = {}
-    d = {}
+    d_cap = {}
 
     # fmt: off
-    d["boundary"] = xr.DataArray(
+    d_cap["boundary"] = xr.DataArray(
         np.array([
             [
                 [1, 1, 1],
@@ -32,7 +32,7 @@ def imod5_cap_data() -> GridDataDict:
         dtype=int),
         **da_kwargs
     )
-    d["landuse"] = xr.DataArray(
+    d_cap["landuse"] = xr.DataArray(
         np.array([
             [
                 [1, 2, 3],
@@ -42,7 +42,7 @@ def imod5_cap_data() -> GridDataDict:
         dtype=int),
         **da_kwargs
     )
-    d["rootzone_thickness"] = xr.DataArray(
+    d_cap["rootzone_thickness"] = xr.DataArray(
         np.array([
             [
                 [0.1, 0.1, 0.1],
@@ -52,7 +52,7 @@ def imod5_cap_data() -> GridDataDict:
         ),
         **da_kwargs
     )
-    d["soil_physical_unit"] = xr.DataArray(
+    d_cap["soil_physical_unit"] = xr.DataArray(
         np.array([
             [
                 [1, 1, 1],
@@ -62,7 +62,7 @@ def imod5_cap_data() -> GridDataDict:
         dtype=int),
         **da_kwargs
     )
-    d["surface_elevation"] = xr.DataArray(
+    d_cap["surface_elevation"] = xr.DataArray(
         np.array([
             [
                 [1.1, 1.2, 1.3],
@@ -72,7 +72,7 @@ def imod5_cap_data() -> GridDataDict:
         ),
         **da_kwargs
     )
-    d["artificial_recharge"] = xr.DataArray(
+    d_cap["artificial_recharge"] = xr.DataArray(
         np.array([
             [
                 [0.2, 0.2, 0.2],
@@ -82,7 +82,7 @@ def imod5_cap_data() -> GridDataDict:
         ),
         **da_kwargs
     )
-    d["artificial_recharge_layer"] = xr.DataArray(
+    d_cap["artificial_recharge_layer"] = xr.DataArray(
         np.array([
             [
                 [1, 2, 3],
@@ -92,7 +92,7 @@ def imod5_cap_data() -> GridDataDict:
         dtype=int),
         **da_kwargs
     )
-    d["artificial_recharge_capacity"] = xr.DataArray(
+    d_cap["artificial_recharge_capacity"] = xr.DataArray(
         np.array([
             [
                 [0.4, 0.4, 0.4],
@@ -102,7 +102,7 @@ def imod5_cap_data() -> GridDataDict:
         ),
         **da_kwargs
     )
-    d["wetted_area"] = xr.DataArray(
+    d_cap["wetted_area"] = xr.DataArray(
         np.array([
             [
                 [0.1, 0.1, 0.1],
@@ -112,7 +112,7 @@ def imod5_cap_data() -> GridDataDict:
         ),
         **da_kwargs
     )
-    d["urban_area"] = xr.DataArray(
+    d_cap["urban_area"] = xr.DataArray(
         np.array([
             [
                 [0.2, 0.2, 0.2],
@@ -122,7 +122,7 @@ def imod5_cap_data() -> GridDataDict:
         ),
         **da_kwargs
     )
-    d["urban_ponding_depth"] = xr.DataArray(
+    d_cap["urban_ponding_depth"] = xr.DataArray(
         np.array([
             [
                 [2.2, 2.2, 2.2],
@@ -132,7 +132,7 @@ def imod5_cap_data() -> GridDataDict:
         ),
         **da_kwargs
     )
-    d["rural_ponding_depth"] = xr.DataArray(
+    d_cap["rural_ponding_depth"] = xr.DataArray(
         np.array([
             [
                 [1.2, 1.2, 1.2],
@@ -142,7 +142,7 @@ def imod5_cap_data() -> GridDataDict:
         ),
         **da_kwargs
     )
-    d["urban_runoff_resistance"] = xr.DataArray(
+    d_cap["urban_runoff_resistance"] = xr.DataArray(
         np.array([
             [
                 [3.2, 3.2, 3.2],
@@ -152,7 +152,7 @@ def imod5_cap_data() -> GridDataDict:
         ),
         **da_kwargs
     )
-    d["rural_runoff_resistance"] = xr.DataArray(
+    d_cap["rural_runoff_resistance"] = xr.DataArray(
         np.array([
             [
                 [3.6, 3.6, 3.6],
@@ -162,7 +162,7 @@ def imod5_cap_data() -> GridDataDict:
         ),
         **da_kwargs
     )
-    d["urban_runon_resistance"] = xr.DataArray(
+    d_cap["urban_runon_resistance"] = xr.DataArray(
         np.array([
             [
                 [5.2, 5.2, 5.2],
@@ -172,7 +172,7 @@ def imod5_cap_data() -> GridDataDict:
         ),
         **da_kwargs
     )
-    d["rural_runon_resistance"] = xr.DataArray(
+    d_cap["rural_runon_resistance"] = xr.DataArray(
         np.array([
             [
                 [5.6, 5.6, 5.6],
@@ -182,7 +182,7 @@ def imod5_cap_data() -> GridDataDict:
         ),
         **da_kwargs
     )
-    d["urban_infiltration_capacity"] = xr.DataArray(
+    d_cap["urban_infiltration_capacity"] = xr.DataArray(
         np.array([
             [
                 [10.2, 10.2, 10.2],
@@ -192,7 +192,7 @@ def imod5_cap_data() -> GridDataDict:
         ),
         **da_kwargs
     )
-    d["rural_infiltration_capacity"] = xr.DataArray(
+    d_cap["rural_infiltration_capacity"] = xr.DataArray(
         np.array([
             [
                 [20.2, 20.2, 20.2],
@@ -202,7 +202,7 @@ def imod5_cap_data() -> GridDataDict:
         ),
         **da_kwargs
     )
-    d["perched_water_table_level"]= xr.DataArray(
+    d_cap["perched_water_table_level"]= xr.DataArray(
         np.array([
             [
                 [2.0, 2.0, 2.0],
@@ -212,7 +212,7 @@ def imod5_cap_data() -> GridDataDict:
         ),
         **da_kwargs
     )
-    d["soil_moisture_fraction"]= xr.DataArray(
+    d_cap["soil_moisture_fraction"]= xr.DataArray(
         np.array([
             [
                 [1.5, 1.5, 1.5],
@@ -222,7 +222,7 @@ def imod5_cap_data() -> GridDataDict:
         ),
         **da_kwargs
     )
-    d["conductivity_factor"]= xr.DataArray(
+    d_cap["conductivity_factor"]= xr.DataArray(
         np.array([
             [
                 [2.5, 2.5, 2.5],
@@ -232,6 +232,19 @@ def imod5_cap_data() -> GridDataDict:
         ),
         **da_kwargs
     )
+    d_bnd = {}
+    d_bnd["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["cap"] = d_cap
+    imod5_data["bnd"] = d_bnd
+
     return imod5_data
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):
diff --git a/imod/tests/test_mf6/test_mf6_model.py b/imod/tests/test_mf6/test_mf6_model.py
index b7a3ecd4c..fc371b88c 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_where_bc")
+    @mock.patch("imod.mf6.utilities.clipped_bc_creator._create_clipped_boundary_pkg")
     @pytest.mark.parametrize(
         "model_type, pkg_type",
         [
@@ -258,14 +265,29 @@ 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)
 
+        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
@@ -290,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_where_bc")
+    @mock.patch("imod.mf6.utilities.clipped_bc_creator._create_clipped_boundary_pkg")
     @pytest.mark.parametrize(
         "model_type, pkg_type",
         [
@@ -300,15 +324,30 @@ 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)
 
+        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]
@@ -339,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:
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]
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],
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
new file mode 100644
index 000000000..9edf9c549
--- /dev/null
+++ b/imod/tests/test_mf6/test_utilities/test_mf6_mask_util.py
@@ -0,0 +1,104 @@
+import numpy as np
+import pandas as pd
+import pytest_cases
+import xarray as xr
+
+from imod.mf6.utilities.mask import mask_topsystem, mask_topsystem_where_bc
+from imod.typing.grid import ones_like, zeros_like
+
+
+def test_mask_topsystem(twri_model):
+    """
+    Test the mask_topsystem utility function by deactivating the first cell in
+    the grid.
+    """
+    # Arrange
+    gwf_model = twri_model["GWF_1"]
+    is_active = ones_like(gwf_model.domain)
+    # Mask first cell
+    is_active[0, 0, 0] = 0
+
+    # Act
+    mask_topsystem(gwf_model, is_active, True)
+    # Assert
+    for key in ["rch", "drn"]:
+        pkg = gwf_model[key]
+        gridded_var = pkg.dataset[pkg._period_data[0]].compute()
+        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, True)
+    # Assert
+    for key in ["rch", "drn"]:
+        assert key not in gwf_model.keys()
+
+
+def test_mask_topsystem_where_bc__no_time(twri_model):
+    """
+    Test the mask_topsystem utility function with boundary conditions by deactivating the first cell in the grid.
+    """
+    # Arrange
+    gwf_model = twri_model["GWF_1"]
+    bc_condition = gwf_model["chd"]
+
+    # Test if fixture still as expected, first cell should be active in the
+    # boundary condition
+    chd_active = np.isnan(bc_condition.dataset["head"])
+    assert chd_active[0, 0, 0].item() is False
+    assert chd_active[0, 0, 1].item() is True
+    # Should be all active
+    rch_active_prior = np.isnan(gwf_model["rch"].dataset["rate"])
+    assert rch_active_prior[0, 0].item() is False
+    assert rch_active_prior[0, 1].item() is False
+
+    # Act
+    mask_topsystem_where_bc(gwf_model, bc_condition, True)
+
+    # Assert
+    rch_active_post = np.isnan(gwf_model["rch"].dataset["rate"])
+    assert rch_active_post[0, 0].item() is True
+    assert rch_active_post[0, 1].item() is False
+
+
+@pytest_cases.parametrize("ignore_time", [True, False])
+def test_mask_topsystem_where_bc__time(twri_model, ignore_time):
+    """
+    Test the mask_topsystem utility function with boundary conditions by deactivating the first cell in the grid.
+    """
+    # Arrange
+    gwf_model = twri_model["GWF_1"]
+    bc_condition = gwf_model["chd"]
+
+    # Make package transient
+    time_coord = pd.date_range("2000-01-01", periods=3)
+    time_da = xr.DataArray([1, 1, 1], dims=["time"], coords={"time": time_coord})
+    ds = bc_condition.dataset.copy()
+    bc_condition.dataset["head"] = time_da * ds["head"]
+
+    # Test if fixture still as expected, first cell should be active in the
+    # boundary condition
+    chd_active = np.isnan(bc_condition.dataset["head"])
+    assert chd_active[1, 0, 0, 0].item() is False
+    assert chd_active[1, 0, 0, 1].item() is True
+    # Should be all active
+    rch_active_prior = np.isnan(gwf_model["rch"].dataset["rate"])
+    assert rch_active_prior[0, 0].item() is False
+    assert rch_active_prior[0, 1].item() is False
+
+    # Act
+    mask_topsystem_where_bc(gwf_model, bc_condition, ignore_time)
+
+    # Assert
+    rch_active_post = np.isnan(gwf_model["rch"].dataset["rate"])
+    assert rch_active_post[0, 0].item() is True
+    assert rch_active_post[0, 1].item() is False
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])
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
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]]
 
diff --git a/imod/util/dims.py b/imod/util/dims.py
index c9228ee39..535ef6cb1 100644
--- a/imod/util/dims.py
+++ b/imod/util/dims.py
@@ -55,6 +55,10 @@ def _drop_layer_if_dataarray(
     return _drop_layer_dim_and_coord(da)
 
 
-def drop_layer_dim_cap_data(imod5_data: Imod5DataDict) -> Imod5DataDict:
+def drop_layer_dim_cap_and_bnd_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()},
+    }