diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst
index 525f330ef..b315052d4 100644
--- a/docs/api/changelog.rst
+++ b/docs/api/changelog.rst
@@ -51,6 +51,14 @@ Added
   :class:`imod.mf6.HorizontalFlowBarrierHydraulicCharacteristic`.
 - :class:`imod.mf6.LayeredWell` to specify wells directly to layers instead
   assigning them with filter depths.
+- :func:`imod.prepare.cleanup_drn`, :func:`imod.prepare.cleanup_ghb`,
+  :func:`imod.prepare.cleanup_riv`. These are utility functions to clean up
+  drainage, general head boundaries, and rivers, respectively.
+- :meth:`imod.mf6.Drainage.cleanup`,
+  :meth:`imod.mf6.GeneralHeadboundary.cleanup`, :meth:`imod.mf6.River.cleanup`
+  convenience methods to call the corresponding cleanup utility functions with
+  the appropriate arguments.
+
 
 Removed
 ~~~~~~~
diff --git a/docs/api/mf6.rst b/docs/api/mf6.rst
index aaf9a1a44..f7046786d 100644
--- a/docs/api/mf6.rst
+++ b/docs/api/mf6.rst
@@ -69,8 +69,14 @@ Flow Packages
     Buoyancy
     ConstantHead
     Drainage
+    Drainage.mask
+    Drainage.regrid_like
+    Drainage.cleanup
     Evapotranspiration
     GeneralHeadBoundary
+    GeneralHeadBoundary.mask
+    GeneralHeadBoundary.regrid_like
+    GeneralHeadBoundary.cleanup
     HorizontalFlowBarrierHydraulicCharacteristic
     HorizontalFlowBarrierMultiplier
     HorizontalFlowBarrierResistance
@@ -83,6 +89,9 @@ Flow Packages
     NodePropertyFlow
     Recharge
     River
+    River.mask
+    River.regrid_like
+    River.cleanup
     SpecificStorage
     StorageCoefficient
     UnsaturatedZoneFlow
diff --git a/docs/api/prepare.rst b/docs/api/prepare.rst
index 55ede360d..b8aed4a19 100644
--- a/docs/api/prepare.rst
+++ b/docs/api/prepare.rst
@@ -47,3 +47,7 @@ Prepare model input
     distribute_drn_conductance
     distribute_ghb_conductance
     distribute_riv_conductance
+
+    cleanup_drn
+    cleanup_ghb
+    cleanup_riv
diff --git a/imod/mf6/drn.py b/imod/mf6/drn.py
index aef93f278..7cc54b371 100644
--- a/imod/mf6/drn.py
+++ b/imod/mf6/drn.py
@@ -4,9 +4,10 @@
 
 import numpy as np
 
-from imod.logging import init_log_decorator
+from imod.logging import init_log_decorator, standard_log_decorator
 from imod.mf6.boundary_condition import BoundaryCondition
 from imod.mf6.dis import StructuredDiscretization
+from imod.mf6.disv import VerticesDiscretization
 from imod.mf6.interfaces.iregridpackage import IRegridPackage
 from imod.mf6.npf import NodePropertyFlow
 from imod.mf6.regrid.regrid_schemes import DrainageRegridMethod
@@ -15,6 +16,7 @@
     _regrid_package_data,
 )
 from imod.mf6.validation import BOUNDARY_DIMS_SCHEMA, CONC_DIMS_SCHEMA
+from imod.prepare.cleanup import cleanup_drn
 from imod.prepare.topsystem.allocation import ALLOCATION_OPTION, allocate_drn_cells
 from imod.prepare.topsystem.conductance import (
     DISTRIBUTING_OPTION,
@@ -163,6 +165,20 @@ def _validate(self, schemata, **kwargs):
 
         return errors
 
+    @standard_log_decorator()
+    def cleanup(self, dis: StructuredDiscretization | VerticesDiscretization) -> None:
+        """
+        Clean up package inplace. This method calls
+        :func:`imod.prepare.cleanup.cleanup_drn`, see documentation of that
+        function for details on cleanup.
+
+        dis: imod.mf6.StructuredDiscretization | imod.mf6.VerticesDiscretization
+            Model discretization package.
+        """
+        dis_dict = {"idomain": dis.dataset["idomain"]}
+        cleaned_dict = self._call_func_on_grids(cleanup_drn, dis_dict)
+        super().__init__(cleaned_dict)
+
     @classmethod
     def from_imod5_data(
         cls,
diff --git a/imod/mf6/ghb.py b/imod/mf6/ghb.py
index 287627fe5..e40d32b86 100644
--- a/imod/mf6/ghb.py
+++ b/imod/mf6/ghb.py
@@ -4,8 +4,10 @@
 
 import numpy as np
 
-from imod.logging import init_log_decorator
+from imod.logging import init_log_decorator, standard_log_decorator
 from imod.mf6.boundary_condition import BoundaryCondition
+from imod.mf6.dis import StructuredDiscretization
+from imod.mf6.disv import VerticesDiscretization
 from imod.mf6.interfaces.iregridpackage import IRegridPackage
 from imod.mf6.npf import NodePropertyFlow
 from imod.mf6.regrid.regrid_schemes import (
@@ -14,6 +16,7 @@
 )
 from imod.mf6.utilities.regrid import RegridderWeightsCache, _regrid_package_data
 from imod.mf6.validation import BOUNDARY_DIMS_SCHEMA, CONC_DIMS_SCHEMA
+from imod.prepare.cleanup import cleanup_ghb
 from imod.prepare.topsystem.allocation import ALLOCATION_OPTION, allocate_ghb_cells
 from imod.prepare.topsystem.conductance import (
     DISTRIBUTING_OPTION,
@@ -165,6 +168,20 @@ def _validate(self, schemata, **kwargs):
 
         return errors
 
+    @standard_log_decorator()
+    def cleanup(self, dis: StructuredDiscretization | VerticesDiscretization) -> None:
+        """
+        Clean up package inplace. This method calls
+        :func:`imod.prepare.cleanup.cleanup_ghb`, see documentation of that
+        function for details on cleanup.
+
+        dis: imod.mf6.StructuredDiscretization | imod.mf6.VerticesDiscretization
+            Model discretization package.
+        """
+        dis_dict = {"idomain": dis.dataset["idomain"]}
+        cleaned_dict = self._call_func_on_grids(cleanup_ghb, dis_dict)
+        super().__init__(cleaned_dict)
+
     @classmethod
     def from_imod5_data(
         cls,
diff --git a/imod/mf6/package.py b/imod/mf6/package.py
index 0b81a18d1..ab2fbbc7a 100644
--- a/imod/mf6/package.py
+++ b/imod/mf6/package.py
@@ -4,7 +4,7 @@
 import pathlib
 from collections import defaultdict
 from copy import deepcopy
-from typing import Any, Dict, List, Mapping, Optional, Tuple, Union
+from typing import Any, Callable, Dict, List, Mapping, Optional, Tuple, Union
 
 import cftime
 import jinja2
@@ -81,6 +81,9 @@ def sel(self):
             f"{type(self).__name__}(**{self._pkg_id}.dataset.sel(**selection))"
         )
 
+    def cleanup(self, dis: Any):
+        raise NotImplementedError("Method not implemented for this package.")
+
     @staticmethod
     def _valid(value: Any) -> bool:
         return _is_valid(value)
@@ -630,6 +633,28 @@ def get_non_grid_data(self, grid_names: list[str]) -> dict[str, Any]:
                 result[name] = self.dataset[name].values[()]
         return result
 
+    def _call_func_on_grids(
+        self, func: Callable, dis: dict
+    ) -> dict[str, GridDataArray]:
+        """
+        Call function on dictionary of grids and merge settings back into
+        dictionary.
+
+        Parameters
+        ----------
+        func: Callable
+            Function to call on all grids
+        """
+        grid_varnames = list(self._write_schemata.keys())
+        grids = {
+            varname: self.dataset[varname]
+            for varname in grid_varnames
+            if varname in self.dataset.keys()
+        }
+        cleaned_grids = func(**dis, **grids)
+        settings = self.get_non_grid_data(grid_varnames)
+        return cleaned_grids | settings
+
     def is_splitting_supported(self) -> bool:
         return True
 
diff --git a/imod/mf6/riv.py b/imod/mf6/riv.py
index d91610651..0787f387c 100644
--- a/imod/mf6/riv.py
+++ b/imod/mf6/riv.py
@@ -6,10 +6,11 @@
 import xarray as xr
 
 from imod import logging
-from imod.logging import init_log_decorator
+from imod.logging import init_log_decorator, standard_log_decorator
 from imod.logging.loglevel import LogLevel
 from imod.mf6.boundary_condition import BoundaryCondition
 from imod.mf6.dis import StructuredDiscretization
+from imod.mf6.disv import VerticesDiscretization
 from imod.mf6.drn import Drainage
 from imod.mf6.interfaces.iregridpackage import IRegridPackage
 from imod.mf6.regrid.regrid_schemes import RiverRegridMethod
@@ -18,6 +19,7 @@
     _regrid_package_data,
 )
 from imod.mf6.validation import BOUNDARY_DIMS_SCHEMA, CONC_DIMS_SCHEMA
+from imod.prepare.cleanup import cleanup_riv
 from imod.prepare.topsystem.allocation import ALLOCATION_OPTION, allocate_riv_cells
 from imod.prepare.topsystem.conductance import (
     DISTRIBUTING_OPTION,
@@ -184,6 +186,20 @@ def _validate(self, schemata, **kwargs):
 
         return errors
 
+    @standard_log_decorator()
+    def cleanup(self, dis: StructuredDiscretization | VerticesDiscretization) -> None:
+        """
+        Clean up package inplace. This method calls
+        :func:`imod.prepare.cleanup.cleanup_riv`, see documentation of that
+        function for details on cleanup.
+
+        dis: imod.mf6.StructuredDiscretization | imod.mf6.VerticesDiscretization
+            Model discretization package.
+        """
+        dis_dict = {"idomain": dis.dataset["idomain"], "bottom": dis.dataset["bottom"]}
+        cleaned_dict = self._call_func_on_grids(cleanup_riv, dis_dict)
+        super().__init__(cleaned_dict)
+
     @classmethod
     def from_imod5_data(
         cls,
diff --git a/imod/mf6/wel.py b/imod/mf6/wel.py
index 5e705307d..b4072fd4d 100644
--- a/imod/mf6/wel.py
+++ b/imod/mf6/wel.py
@@ -30,6 +30,7 @@
 from imod.mf6.validation import validation_pkg_error_message
 from imod.mf6.write_context import WriteContext
 from imod.prepare import assign_wells
+from imod.prepare.cleanup import cleanup_wel
 from imod.prepare.layer import create_layered_top
 from imod.schemata import (
     AllValueSchema,
@@ -971,6 +972,9 @@ def _validate_imod5_depth_information(
             logger.log(loglevel=LogLevel.ERROR, message=log_msg, additional_depth=2)
             raise ValueError(log_msg)
 
+    def cleanup(self):
+        self.dataset = cleanup_wel(self.dataset)
+
 
 class LayeredWell(GridAgnosticWell):
     """
diff --git a/imod/prepare/__init__.py b/imod/prepare/__init__.py
index caf5fc202..2ea86c3cb 100644
--- a/imod/prepare/__init__.py
+++ b/imod/prepare/__init__.py
@@ -14,6 +14,7 @@
 """
 
 from imod.prepare import hfb, spatial, subsoil, surface_water
+from imod.prepare.cleanup import cleanup_drn, cleanup_ghb, cleanup_riv
 from imod.prepare.hfb import (
     linestring_to_square_zpolygons,
     linestring_to_trapezoid_zpolygons,
diff --git a/imod/prepare/cleanup.py b/imod/prepare/cleanup.py
new file mode 100644
index 000000000..0d6edf10d
--- /dev/null
+++ b/imod/prepare/cleanup.py
@@ -0,0 +1,219 @@
+"""Cleanup utilities"""
+
+from enum import Enum
+from typing import Optional
+
+import xarray as xr
+
+from imod.mf6.utilities.mask import mask_arrays
+from imod.schemata import scalar_None
+from imod.typing import GridDataArray, GridDataset
+
+
+class AlignLevelsMode(Enum):
+    TOPDOWN = 0
+    BOTTOMUP = 1
+
+
+def align_nodata(grids: dict[str, xr.DataArray]) -> dict[str, xr.DataArray]:
+    return mask_arrays(grids)
+
+
+def align_interface_levels(
+    top: GridDataArray,
+    bottom: GridDataArray,
+    method: AlignLevelsMode = AlignLevelsMode.TOPDOWN,
+) -> tuple[GridDataArray, GridDataArray]:
+    to_align = top < bottom
+
+    match method:
+        case AlignLevelsMode.BOTTOMUP:
+            return top.where(~to_align, bottom), bottom
+        case AlignLevelsMode.TOPDOWN:
+            return top, bottom.where(~to_align, top)
+        case _:
+            raise TypeError(f"Unmatched case for method, got {method}")
+
+
+def _cleanup_robin_boundary(
+    idomain: GridDataArray, grids: dict[str, GridDataArray]
+) -> dict[str, GridDataArray]:
+    """Cleanup robin boundary condition (i.e. bc with conductance)"""
+    active = idomain == 1
+    # Deactivate conductance cells outside active domain; this nodata
+    # inconsistency will be aligned in the final call to align_nodata
+    conductance = grids["conductance"].where(active)
+    concentration = grids["concentration"]
+    # Make conductance cells with erronous values inactive
+    grids["conductance"] = conductance.where(conductance > 0.0)
+    # Clip negative concentration cells to 0.0
+    if (concentration is not None) and not scalar_None(concentration):
+        grids["concentration"] = concentration.clip(min=0.0)
+    else:
+        grids.pop("concentration")
+
+    # Align nodata
+    return align_nodata(grids)
+
+
+def cleanup_riv(
+    idomain: GridDataArray,
+    bottom: GridDataArray,
+    stage: GridDataArray,
+    conductance: GridDataArray,
+    bottom_elevation: GridDataArray,
+    concentration: Optional[GridDataArray] = None,
+) -> dict[str, GridDataArray]:
+    """
+    Clean up river data, fixes some common mistakes causing ValidationErrors by
+    doing the following:
+
+    - Cells where conductance <= 0 are deactivated.
+    - Cells where concentration < 0 are set to 0.0.
+    - Cells outside active domain (idomain==1) are removed.
+    - Align NoData: If one variable has an inactive cell in one cell, ensure
+      this cell is deactivated for all variables.
+    - River bottom elevations below model bottom of a layer are set to model
+      bottom of that layer.
+    - River bottom elevations which exceed river stage are lowered to river
+      stage.
+
+    Parameters
+    ----------
+    idomain: xarray.DataArray | xugrid.UgridDataArray
+        MODFLOW 6 model domain. idomain==1 is considered active domain.
+    bottom: xarray.DataArray | xugrid.UgridDataArray
+        Grid with`model bottoms
+    stage: xarray.DataArray | xugrid.UgridDataArray
+        Grid with river stages
+    conductance: xarray.DataArray | xugrid.UgridDataArray
+        Grid with conductances
+    bottom_elevation: xarray.DataArray | xugrid.UgridDataArray
+        Grid with river bottom elevations
+    concentration: xarray.DataArray | xugrid.UgridDataArray, optional
+        Optional grid with concentrations
+
+    Returns
+    -------
+    dict[str, xarray.DataArray | xugrid.UgridDataArray]
+        Dict of cleaned up grids. Has keys: "stage", "conductance",
+        "bottom_elevation", "concentration".
+    """
+    # Output dict
+    output_dict = {
+        "stage": stage,
+        "conductance": conductance,
+        "bottom_elevation": bottom_elevation,
+        "concentration": concentration,
+    }
+    output_dict = _cleanup_robin_boundary(idomain, output_dict)
+    if (output_dict["stage"] < bottom).any():
+        raise ValueError(
+            "River stage below bottom of model layer, cannot fix this. "
+            "Probably rivers are assigned to the wrong layer, you can reallocate "
+            "river data to model layers with: "
+            "``imod.prepare.topsystem.allocate_riv_cells``."
+        )
+    # Ensure bottom elevation above model bottom
+    output_dict["bottom_elevation"], _ = align_interface_levels(
+        output_dict["bottom_elevation"], bottom, AlignLevelsMode.BOTTOMUP
+    )
+    # Ensure stage above bottom_elevation
+    output_dict["stage"], output_dict["bottom_elevation"] = align_interface_levels(
+        output_dict["stage"], output_dict["bottom_elevation"], AlignLevelsMode.TOPDOWN
+    )
+    return output_dict
+
+
+def cleanup_drn(
+    idomain: GridDataArray,
+    elevation: GridDataArray,
+    conductance: GridDataArray,
+    concentration: Optional[GridDataArray] = None,
+) -> dict[str, GridDataArray]:
+    """
+    Clean up drain data, fixes some common mistakes causing ValidationErrors by
+    doing the following:
+
+    - Cells where conductance <= 0 are deactivated.
+    - Cells where concentration < 0 are set to 0.0.
+    - Cells outside active domain (idomain==1) are removed.
+    - Align NoData: If one variable has an inactive cell in one cell, ensure
+      this cell is deactivated for all variables.
+
+    Parameters
+    ----------
+    idomain: xarray.DataArray | xugrid.UgridDataArray
+        MODFLOW 6 model domain. idomain==1 is considered active domain.
+    elevation: xarray.DataArray | xugrid.UgridDataArray
+        Grid with drain elevations
+    conductance: xarray.DataArray | xugrid.UgridDataArray
+        Grid with conductances
+    concentration: xarray.DataArray | xugrid.UgridDataArray, optional
+        Optional grid with concentrations
+
+    Returns
+    -------
+    dict[str, xarray.DataArray | xugrid.UgridDataArray]
+        Dict of cleaned up grids. Has keys: "elevation", "conductance",
+        "concentration".
+    """
+    # Output dict
+    output_dict = {
+        "elevation": elevation,
+        "conductance": conductance,
+        "concentration": concentration,
+    }
+    return _cleanup_robin_boundary(idomain, output_dict)
+
+
+def cleanup_ghb(
+    idomain: GridDataArray,
+    head: GridDataArray,
+    conductance: GridDataArray,
+    concentration: Optional[GridDataArray] = None,
+) -> dict[str, GridDataArray]:
+    """
+    Clean up general head boundary data, fixes some common mistakes causing
+    ValidationErrors by doing the following:
+
+    - Cells where conductance <= 0 are deactivated.
+    - Cells where concentration < 0 are set to 0.0.
+    - Cells outside active domain (idomain==1) are removed.
+    - Align NoData: If one variable has an inactive cell in one cell, ensure
+      this cell is deactivated for all variables.
+
+    Parameters
+    ----------
+    idomain: xarray.DataArray | xugrid.UgridDataArray
+        MODFLOW 6 model domain. idomain==1 is considered active domain.
+    head: xarray.DataArray | xugrid.UgridDataArray
+        Grid with heads
+    conductance: xarray.DataArray | xugrid.UgridDataArray
+        Grid with conductances
+    concentration: xarray.DataArray | xugrid.UgridDataArray, optional
+        Optional grid with concentrations
+
+    Returns
+    -------
+    dict[str, xarray.DataArray | xugrid.UgridDataArray]
+        Dict of cleaned up grids. Has keys: "head", "conductance",
+        "concentration".
+    """
+    # Output dict
+    output_dict = {
+        "head": head,
+        "conductance": conductance,
+        "concentration": concentration,
+    }
+    return _cleanup_robin_boundary(idomain, output_dict)
+
+
+def cleanup_wel(wel_ds: GridDataset):
+    """
+    Clean up wells
+
+    - Removes wells where the screen bottom elevation exceeds screen top.
+    """
+    deactivate = wel_ds["screen_top"] < wel_ds["screen_bottom"]
+    return wel_ds.where(~deactivate, drop=True)
diff --git a/imod/tests/test_mf6/test_mf6_riv.py b/imod/tests/test_mf6/test_mf6_riv.py
index d48cab220..d57ef157b 100644
--- a/imod/tests/test_mf6/test_mf6_riv.py
+++ b/imod/tests/test_mf6/test_mf6_riv.py
@@ -12,14 +12,19 @@
 
 import imod
 from imod.mf6.dis import StructuredDiscretization
+from imod.mf6.disv import VerticesDiscretization
 from imod.mf6.write_context import WriteContext
 from imod.prepare.topsystem.allocation import ALLOCATION_OPTION
 from imod.prepare.topsystem.conductance import DISTRIBUTING_OPTION
 from imod.schemata import ValidationError
 from imod.typing.grid import ones_like, zeros_like
 
+TYPE_DIS_PKG = {
+    xu.UgridDataArray: VerticesDiscretization,
+    xr.DataArray: StructuredDiscretization,
+}
+
 
-@pytest.fixture(scope="function")
 def make_da():
     x = [5.0, 15.0, 25.0]
     y = [25.0, 15.0, 5.0]
@@ -34,9 +39,8 @@ def make_da():
     )
 
 
-@pytest.fixture(scope="function")
-def dis_dict(make_da):
-    da = make_da
+def dis_dict():
+    da = make_da()
     bottom = da - xr.DataArray(
         data=[1.5, 2.5], dims=("layer",), coords={"layer": [2, 3]}
     )
@@ -44,16 +48,15 @@ def dis_dict(make_da):
     return {"idomain": da.astype(int), "top": da.sel(layer=2), "bottom": bottom}
 
 
-@pytest.fixture(scope="function")
-def riv_dict(make_da):
-    da = make_da
+def riv_dict():
+    da = make_da()
     da[:, 1, 1] = np.nan
 
     bottom = da - xr.DataArray(
         data=[1.0, 2.0], dims=("layer",), coords={"layer": [2, 3]}
     )
 
-    return {"stage": da, "conductance": da, "bottom_elevation": bottom}
+    return {"stage": da, "conductance": da.copy(), "bottom_elevation": bottom}
 
 
 def make_dict_unstructured(d):
@@ -61,19 +64,19 @@ def make_dict_unstructured(d):
 
 
 class RivCases:
-    def case_structured(self, riv_dict):
-        return riv_dict
+    def case_structured(self):
+        return riv_dict()
 
-    def case_unstructured(self, riv_dict):
-        return make_dict_unstructured(riv_dict)
+    def case_unstructured(self):
+        return make_dict_unstructured(riv_dict())
 
 
 class RivDisCases:
-    def case_structured(self, riv_dict, dis_dict):
-        return riv_dict, dis_dict
+    def case_structured(self):
+        return riv_dict(), dis_dict()
 
-    def case_unstructured(self, riv_dict, dis_dict):
-        return make_dict_unstructured(riv_dict), make_dict_unstructured(dis_dict)
+    def case_unstructured(self):
+        return make_dict_unstructured(riv_dict()), make_dict_unstructured(dis_dict())
 
 
 @parametrize_with_cases("riv_data", cases=RivCases)
@@ -118,19 +121,32 @@ def test_all_nan(riv_data, dis_data):
     errors = river._validate(river._write_schemata, **dis_data)
 
     assert len(errors) == 1
+    assert "stage" in errors.keys()
 
-    for var, var_errors in errors.items():
-        assert var == "stage"
+
+@parametrize_with_cases("riv_data,dis_data", cases=RivDisCases)
+def test_validate_inconsistent_nan(riv_data, dis_data):
+    riv_data["stage"][..., 2] = np.nan
+    river = imod.mf6.River(**riv_data)
+
+    errors = river._validate(river._write_schemata, **dis_data)
+
+    assert len(errors) == 2
+    assert "bottom_elevation" in errors.keys()
+    assert "conductance" in errors.keys()
 
 
 @parametrize_with_cases("riv_data,dis_data", cases=RivDisCases)
-def test_inconsistent_nan(riv_data, dis_data):
+def test_cleanup_inconsistent_nan(riv_data, dis_data):
     riv_data["stage"][..., 2] = np.nan
     river = imod.mf6.River(**riv_data)
+    type_grid = type(riv_data["stage"])
+    dis_pkg = TYPE_DIS_PKG[type_grid](**dis_data)
 
+    river.cleanup(dis_pkg)
     errors = river._validate(river._write_schemata, **dis_data)
 
-    assert len(errors) == 1
+    assert len(errors) == 0
 
 
 @parametrize_with_cases("riv_data,dis_data", cases=RivDisCases)
@@ -210,11 +226,11 @@ def test_check_dimsize_zero():
 
 
 @parametrize_with_cases("riv_data,dis_data", cases=RivDisCases)
-def test_check_zero_conductance(riv_data, dis_data):
+def test_validate_zero_conductance(riv_data, dis_data):
     """
-    Test for zero conductance
+    Test for validation zero conductance
     """
-    riv_data["conductance"] = riv_data["conductance"] * 0.0
+    riv_data["conductance"][..., 2] = 0.0
 
     river = imod.mf6.River(**riv_data)
 
@@ -226,9 +242,25 @@ def test_check_zero_conductance(riv_data, dis_data):
 
 
 @parametrize_with_cases("riv_data,dis_data", cases=RivDisCases)
-def test_check_bottom_above_stage(riv_data, dis_data):
+def test_cleanup_zero_conductance(riv_data, dis_data):
     """
-    Check that river bottom is not above stage.
+    Cleanup zero conductance
+    """
+    riv_data["conductance"][..., 2] = 0.0
+    type_grid = type(riv_data["stage"])
+    dis_pkg = TYPE_DIS_PKG[type_grid](**dis_data)
+
+    river = imod.mf6.River(**riv_data)
+    river.cleanup(dis_pkg)
+
+    errors = river._validate(river._write_schemata, **dis_data)
+    assert len(errors) == 0
+
+
+@parametrize_with_cases("riv_data,dis_data", cases=RivDisCases)
+def test_validate_bottom_above_stage(riv_data, dis_data):
+    """
+    Validate that river bottom is not above stage.
     """
 
     riv_data["bottom_elevation"] = riv_data["bottom_elevation"] + 10.0
@@ -238,8 +270,26 @@ def test_check_bottom_above_stage(riv_data, dis_data):
     errors = river._validate(river._write_schemata, **dis_data)
 
     assert len(errors) == 1
-    for var, var_errors in errors.items():
-        assert var == "stage"
+    assert "stage" in errors.keys()
+
+
+@parametrize_with_cases("riv_data,dis_data", cases=RivDisCases)
+def test_cleanup_bottom_above_stage(riv_data, dis_data):
+    """
+    Cleanup river bottom above stage.
+    """
+
+    riv_data["bottom_elevation"] = riv_data["bottom_elevation"] + 10.0
+    type_grid = type(riv_data["stage"])
+    dis_pkg = TYPE_DIS_PKG[type_grid](**dis_data)
+
+    river = imod.mf6.River(**riv_data)
+    river.cleanup(dis_pkg)
+
+    errors = river._validate(river._write_schemata, **dis_data)
+
+    assert len(errors) == 0
+    assert river.dataset["bottom_elevation"].equals(river.dataset["stage"])
 
 
 @parametrize_with_cases("riv_data,dis_data", cases=RivDisCases)
@@ -280,12 +330,12 @@ def test_check_boundary_outside_active_domain(riv_data, dis_data):
     assert len(errors) == 1
 
 
-def test_check_dim_monotonicity(riv_dict):
+def test_check_dim_monotonicity():
     """
     Test if dimensions are monotonically increasing or, in case of the y coord,
     decreasing
     """
-    riv_ds = xr.merge([riv_dict])
+    riv_ds = xr.merge([riv_dict()])
 
     message = textwrap.dedent(
         """
@@ -327,19 +377,19 @@ def test_check_dim_monotonicity(riv_dict):
         imod.mf6.River(**riv_ds.sel(layer=slice(None, None, -1)))
 
 
-def test_validate_false(riv_dict):
+def test_validate_false():
     """
     Test turning off validation
     """
 
-    riv_ds = xr.merge([riv_dict])
+    riv_ds = xr.merge([riv_dict()])
 
     imod.mf6.River(validate=False, **riv_ds.sel(layer=slice(None, None, -1)))
 
 
 @pytest.mark.usefixtures("concentration_fc")
-def test_render_concentration(riv_dict, concentration_fc):
-    riv_ds = xr.merge([riv_dict])
+def test_render_concentration(concentration_fc):
+    riv_ds = xr.merge([riv_dict()])
 
     concentration = concentration_fc.sel(
         layer=[2, 3], time=np.datetime64("2000-01-01"), drop=True
diff --git a/imod/tests/test_prepare/test_cleanup.py b/imod/tests/test_prepare/test_cleanup.py
new file mode 100644
index 000000000..1802ca22d
--- /dev/null
+++ b/imod/tests/test_prepare/test_cleanup.py
@@ -0,0 +1,178 @@
+from typing import Callable
+
+import numpy as np
+import pytest
+import xugrid as xu
+from pytest_cases import parametrize, parametrize_with_cases
+
+from imod.prepare.cleanup import cleanup_drn, cleanup_ghb, cleanup_riv
+from imod.tests.test_mf6.test_mf6_riv import RivDisCases
+from imod.typing import GridDataArray
+
+
+def _first(grid: GridDataArray):
+    """
+    helper function to get first value, regardless of unstructured or
+    structured grid."""
+    return grid.values.ravel()[0]
+
+
+def _first_index(grid: GridDataArray) -> tuple:
+    if isinstance(grid, xu.UgridDataArray):
+        return (0, 0)
+    else:
+        return (0, 0, 0)
+
+
+def _rename_data_dict(data: dict, func: Callable):
+    renamed = data.copy()
+    to_rename = _RENAME_DICT[func]
+    for src, dst in to_rename.items():
+        mv_data = renamed.pop(src)
+        if dst is not None:
+            renamed[dst] = mv_data
+    return renamed
+
+
+def _prepare_dis_dict(dis_dict: dict, func: Callable):
+    """Keep required dis args for specific cleanup functions"""
+    keep_vars = _KEEP_FROM_DIS_DICT[func]
+    return {var: dis_dict[var] for var in keep_vars}
+
+
+_RENAME_DICT = {
+    cleanup_riv: {},
+    cleanup_drn: {"stage": "elevation", "bottom_elevation": None},
+    cleanup_ghb: {"stage": "head", "bottom_elevation": None},
+}
+
+_KEEP_FROM_DIS_DICT = {
+    cleanup_riv: ["idomain", "bottom"],
+    cleanup_drn: ["idomain"],
+    cleanup_ghb: ["idomain"],
+}
+
+
+@parametrize_with_cases("riv_data, dis_data", cases=RivDisCases)
+@parametrize("cleanup_func", [cleanup_drn, cleanup_ghb, cleanup_riv])
+def test_cleanup__align_nodata(riv_data: dict, dis_data: dict, cleanup_func: Callable):
+    dis_dict = _prepare_dis_dict(dis_data, cleanup_func)
+    data_dict = _rename_data_dict(riv_data, cleanup_func)
+    # Assure conductance not modified by previous tests.
+    np.testing.assert_equal(_first(data_dict["conductance"]), 1.0)
+    idx = _first_index(data_dict["conductance"])
+    # Arrange: Deactivate one cell
+    first_key = next(iter(data_dict.keys()))
+    data_dict[first_key][idx] = np.nan
+    # Act
+    data_cleaned = cleanup_func(**dis_dict, **data_dict)
+    # Assert
+    for key in data_cleaned.keys():
+        np.testing.assert_equal(_first(data_cleaned[key][idx]), np.nan)
+
+
+@parametrize_with_cases("riv_data, dis_data", cases=RivDisCases)
+@parametrize("cleanup_func", [cleanup_drn, cleanup_ghb, cleanup_riv])
+def test_cleanup__zero_conductance(
+    riv_data: dict, dis_data: dict, cleanup_func: Callable
+):
+    dis_dict = _prepare_dis_dict(dis_data, cleanup_func)
+    data_dict = _rename_data_dict(riv_data, cleanup_func)
+    # Assure conductance not modified by previous tests.
+    np.testing.assert_equal(_first(data_dict["conductance"]), 1.0)
+    idx = _first_index(data_dict["conductance"])
+    # Arrange: Deactivate one cell
+    data_dict["conductance"][idx] = 0.0
+    # Act
+    data_cleaned = cleanup_func(**dis_dict, **data_dict)
+    # Assert
+    for key in data_cleaned.keys():
+        np.testing.assert_equal(_first(data_cleaned[key][idx]), np.nan)
+
+
+@parametrize_with_cases("riv_data, dis_data", cases=RivDisCases)
+@parametrize("cleanup_func", [cleanup_drn, cleanup_ghb, cleanup_riv])
+def test_cleanup__negative_concentration(
+    riv_data: dict, dis_data: dict, cleanup_func: Callable
+):
+    dis_dict = _prepare_dis_dict(dis_data, cleanup_func)
+    data_dict = _rename_data_dict(riv_data, cleanup_func)
+    first_key = next(iter(data_dict.keys()))
+    # Create concentration data
+    data_dict["concentration"] = data_dict[first_key].copy()
+    # Assure conductance not modified by previous tests.
+    np.testing.assert_equal(_first(data_dict["conductance"]), 1.0)
+    idx = _first_index(data_dict["conductance"])
+    # Arrange: Deactivate one cell
+    data_dict["concentration"][idx] = -10.0
+    # Act
+    data_cleaned = cleanup_func(**dis_dict, **data_dict)
+    # Assert
+    np.testing.assert_equal(_first(data_cleaned["concentration"]), 0.0)
+
+
+@parametrize_with_cases("riv_data, dis_data", cases=RivDisCases)
+@parametrize("cleanup_func", [cleanup_drn, cleanup_ghb, cleanup_riv])
+def test_cleanup__outside_active_domain(
+    riv_data: dict, dis_data: dict, cleanup_func: Callable
+):
+    dis_dict = _prepare_dis_dict(dis_data, cleanup_func)
+    data_dict = _rename_data_dict(riv_data, cleanup_func)
+    # Assure conductance not modified by previous tests.
+    np.testing.assert_equal(_first(data_dict["conductance"]), 1.0)
+    idx = _first_index(data_dict["conductance"])
+    # Arrange: Deactivate one cell
+    dis_dict["idomain"][idx] = 0.0
+    # Act
+    data_cleaned = cleanup_func(**dis_dict, **data_dict)
+    # Assert
+    for key in data_cleaned.keys():
+        np.testing.assert_equal(_first(data_cleaned[key][idx]), np.nan)
+
+
+@parametrize_with_cases("riv_data, dis_data", cases=RivDisCases)
+def test_cleanup_riv__fix_bottom_elevation_to_bottom(riv_data: dict, dis_data: dict):
+    dis_dict = _prepare_dis_dict(dis_data, cleanup_riv)
+    # Arrange: Set bottom elevation model layer bottom
+    riv_data["bottom_elevation"] -= 3.0
+    # Assure conductance not modified by previous tests.
+    np.testing.assert_equal(_first(riv_data["conductance"]), 1.0)
+    # Act
+    riv_data_cleaned = cleanup_riv(**dis_dict, **riv_data)
+    # Assert
+    # Account for cells inactive river cells.
+    riv_active = riv_data_cleaned["stage"].notnull()
+    expected = dis_dict["bottom"].where(riv_active)
+
+    np.testing.assert_equal(
+        riv_data_cleaned["bottom_elevation"].values, expected.values
+    )
+
+
+@parametrize_with_cases("riv_data, dis_data", cases=RivDisCases)
+def test_cleanup_riv__fix_bottom_elevation_to_stage(riv_data: dict, dis_data: dict):
+    dis_dict = _prepare_dis_dict(dis_data, cleanup_riv)
+    # Arrange: Set bottom elevation above stage
+    riv_data["bottom_elevation"] += 3.0
+    # Assure conductance not modified by previous tests.
+    np.testing.assert_equal(_first(riv_data["conductance"]), 1.0)
+    # Act
+    riv_data_cleaned = cleanup_riv(**dis_dict, **riv_data)
+    # Assert
+    np.testing.assert_equal(
+        riv_data_cleaned["bottom_elevation"].values, riv_data_cleaned["stage"].values
+    )
+
+
+@parametrize_with_cases("riv_data, dis_data", cases=RivDisCases)
+def test_cleanup_riv__raise_error(riv_data: dict, dis_data: dict):
+    """
+    Test if error raised when stage below model layer bottom and see if user is
+    guided to the right prepare function.
+    """
+    dis_dict = _prepare_dis_dict(dis_data, cleanup_riv)
+    # Arrange: Set bottom elevation above stage
+    riv_data["stage"] -= 10.0
+    # Act
+    with pytest.raises(ValueError, match="imod.prepare.topsystem.allocate_riv_cells"):
+        cleanup_riv(**dis_dict, **riv_data)