diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst
index 6a0f7b718..7592a886d 100644
--- a/docs/api/changelog.rst
+++ b/docs/api/changelog.rst
@@ -74,6 +74,12 @@ Fixed
   settings were specified as dask array, which could happen when loading a model
   lazily with :meth:`imod.mf6.Modflow6Simulation.from_file` and not
   computing the data before writing.
+- Fixed bug where ``concentration`` variables were needlessly loaded into
+  memory. Affected :class:`imod.mf6.River`, :class:`imod.mf6.Drain`,
+  :class:`imod.mf6.ConstantHeadBoundary`, :class:`imod.mf6.Recharge`,
+  :class:`imod.mf6.Well` and :class:`imod.mf6.GeneralHeadBoundary`.
+- Fix bug where ``maxbound`` of the ``.wel`` file computed by
+  :class:`imod.mf6.Mf6Wel` was twice or thrice too large.
 
 Changed
 ~~~~~~~
diff --git a/imod/common/utilities/clip.py b/imod/common/utilities/clip.py
index da1862a5c..2a303a3ec 100644
--- a/imod/common/utilities/clip.py
+++ b/imod/common/utilities/clip.py
@@ -20,7 +20,7 @@
     vertical_polygons_to_zbound_linestrings,
 )
 from imod.common.utilities.mask import mask_da
-from imod.common.utilities.value_filters import is_valid
+from imod.common.utilities.value_filters import is_scalar_nan
 from imod.typing import GeoDataFrameType, GridDataArray, GridDataset
 from imod.typing.grid import bounding_polygon, is_spatial_grid
 from imod.util.imports import MissingOptionalModule
@@ -430,7 +430,7 @@ def clip_time_slice(
             time_end=time_end,
         )
 
-        if "repeat_stress" in selection.data_vars and is_valid(
+        if "repeat_stress" in selection.data_vars and is_scalar_nan(
             selection["repeat_stress"].values[()]
         ):
             repeat_indexer, repeat_stress = clip_repeat_stress(
diff --git a/imod/common/utilities/regrid.py b/imod/common/utilities/regrid.py
index a3d5995b3..f9ab40d0a 100644
--- a/imod/common/utilities/regrid.py
+++ b/imod/common/utilities/regrid.py
@@ -19,7 +19,7 @@
 from imod.common.utilities.clip import clip_by_grid
 from imod.common.utilities.dataclass_type import DataclassType, EmptyRegridMethod
 from imod.common.utilities.dtype import is_integer
-from imod.common.utilities.value_filters import enforce_scalar, is_valid
+from imod.common.utilities.value_filters import enforce_scalar, is_scalar_nan
 from imod.typing import Imod5DataDict
 from imod.typing.grid import (
     GridDataArray,
@@ -73,7 +73,7 @@ def _regrid_array(
 
     # skip regridding for scalar arrays with no valid values (such as "None")
     scalar_da: bool = is_scalar(da)
-    if scalar_da and not is_valid(enforce_scalar(da)):
+    if scalar_da and not is_scalar_nan(enforce_scalar(da)):
         return None
 
     # the dataarray might be a scalar. If it is, then it does not need regridding.
diff --git a/imod/common/utilities/value_filters.py b/imod/common/utilities/value_filters.py
index ea7163d22..ad890c40b 100644
--- a/imod/common/utilities/value_filters.py
+++ b/imod/common/utilities/value_filters.py
@@ -1,4 +1,3 @@
-import numbers
 from typing import Any
 
 import numpy as np
@@ -8,18 +7,7 @@
 from imod.typing import GridDataArray, GridDataset
 
 
-def is_scalar_nan(da: GridDataArray):
-    """
-    Test if is_scalar_nan, carefully avoid loading grids in memory
-    """
-    scalar_data: bool = is_scalar(da)
-    if scalar_data:
-        stripped_value = da.to_numpy()[()]
-        return isinstance(stripped_value, numbers.Real) and np.isnan(stripped_value)  # type: ignore[call-overload]
-    return False
-
-
-def is_valid(value: Any) -> bool:
+def is_scalar_nan(value: Any) -> bool:
     """
     Filters values that are None, False, or a numpy.bool_ False.
     Needs to be this specific, since 0.0 and 0 are valid values, but are
diff --git a/imod/mf6/boundary_condition.py b/imod/mf6/boundary_condition.py
index d8861c8a6..882ee10df 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.utilities.value_filters import enforce_scalar
 from imod.mf6.auxiliary_variables import (
     expand_transient_auxiliary_variables,
     get_variable_names,
@@ -92,7 +93,12 @@ def _max_active_n(self):
         Determine the maximum active number of cells that are active
         during a stress period.
         """
-        da = self.dataset[self._get_period_varnames()[0]]
+        if not hasattr(self, "_period_data") or len(self._period_data) == 0:
+            raise ValueError("No period data variables defined for this package.")
+        # Get the last period data element to work around the first variable
+        # (cellid) of the Mf6HorizontalFlowBarrier having an additional
+        # dimension.
+        da = self.dataset[self._period_data[-1]]
         if "time" in da.coords:
             nmax = int(da.groupby("time").count(xr.ALL_DIMS).max())
         else:
@@ -206,12 +212,13 @@ def _get_unfiltered_pkg_options(
         options = copy(predefined_options)
 
         if not_options is None:
-            not_options = self._get_period_varnames()
+            other_exceptions = ["concentration", "repeat_stress"]
+            not_options = self._get_period_varnames() + other_exceptions
 
         for varname in self.dataset.data_vars.keys():  # pylint:disable=no-member
             if varname in not_options:
                 continue
-            v = self.dataset[varname].values[()]
+            v = enforce_scalar(self.dataset[varname])
             options[str(varname)] = v
         return options
 
@@ -243,12 +250,12 @@ def _render(self, directory, pkgname, globaltimes, binary):
         """Render fills in the template only, doesn't write binary data"""
         d = {"binary": binary}
         bin_ds = self._get_bin_ds()
-        d["periods"] = self._period_paths(
-            directory, pkgname, globaltimes, bin_ds, binary
-        )
         # construct the rest (dict for render)
         d = self._get_pkg_options(d)
         d["maxbound"] = self._max_active_n()
+        d["periods"] = self._period_paths(
+            directory, pkgname, globaltimes, bin_ds, binary
+        )
 
         if (hasattr(self, "_auxiliary_data")) and (names := get_variable_names(self)):
             d["auxiliary"] = names
@@ -298,7 +305,9 @@ def _write(
 
     def _get_period_varnames(self) -> list[str]:
         """
-        Get variable names for transient data of this package.
+        Get variable names for transient data of this package. These will be the
+        variables that are written to binary files and referenced in the block
+        file.
 
         Returns
         -------
@@ -313,7 +322,7 @@ def _get_period_varnames(self) -> list[str]:
 
         >>> river = imod.mf6.River.from_file("river_with_concentration.nc")
         >>> river._get_period_varnames()
-        >>> # prints: ['stage', 'conductance', 'bottom_elevation', 'concentration']
+        >>> # prints: ['stage', 'conductance', 'bottom_elevation', 'species1', 'species2']
         """
         result = []
         if hasattr(self, "_period_data"):
diff --git a/imod/mf6/mf6_hfb_adapter.py b/imod/mf6/mf6_hfb_adapter.py
index d922264fc..399907437 100644
--- a/imod/mf6/mf6_hfb_adapter.py
+++ b/imod/mf6/mf6_hfb_adapter.py
@@ -102,7 +102,7 @@ class Mf6HorizontalFlowBarrier(BoundaryCondition):
     _write_schemata = {
         "idomain": (AnyValueSchema(">", 0),),
     }
-    _period_data = ("hydraulic_characteristic",)
+    _period_data = ("cell_id1", "cell_id2", "layer", "hydraulic_characteristic")
     _keyword_map = {}
     _template = BoundaryCondition._initialize_template(_pkg_id)
 
@@ -124,13 +124,6 @@ def __init__(
         }
         super().__init__(dict_dataset)
 
-    def _get_bin_ds(self):
-        bin_ds = self.dataset[
-            ["cell_id1", "cell_id2", "layer", "hydraulic_characteristic"]
-        ]
-
-        return bin_ds
-
     def _ds_to_arrdict(self, ds):
         arrdict = super()._ds_to_arrdict(ds)
         arrdict["cell_dims1"] = ds.coords["cell_dims1"].values
diff --git a/imod/mf6/package.py b/imod/mf6/package.py
index 54898acf7..c7d3e1764 100644
--- a/imod/mf6/package.py
+++ b/imod/mf6/package.py
@@ -37,7 +37,7 @@
     validate_schemata_dict,
     validate_with_error_message,
 )
-from imod.common.utilities.value_filters import is_valid
+from imod.common.utilities.value_filters import is_scalar_nan
 from imod.common.utilities.version import prepend_content_with_version_info
 from imod.logging import standard_log_decorator
 from imod.mf6.auxiliary_variables import (
@@ -87,7 +87,7 @@ def __init__(self, allargs: Mapping[str, GridDataArray | float | int | bool | st
 
     @staticmethod
     def _valid(value: Any) -> bool:
-        return is_valid(value)
+        return is_scalar_nan(value)
 
     @staticmethod
     def _number_format(dtype: type):
diff --git a/imod/mf6/uzf.py b/imod/mf6/uzf.py
index a9336a4d7..a4d936c07 100644
--- a/imod/mf6/uzf.py
+++ b/imod/mf6/uzf.py
@@ -422,14 +422,17 @@ def _render(self, directory, pkgname, globaltimes, binary):
         d["periods"] = self._period_paths(
             directory, pkgname, globaltimes, bin_ds, binary=False
         )
+        # List variables that are not in the period data or package data, but
+        # are needed for rendering the template and are not options.
+        vars_from_init = ["iuzno", "ivertcon", "landflag", "stress_period_active"]
         not_options = (
-            list(self._period_data) + list(self._package_data) + ["iuzno" + "ivertcon"]
+            list(self._period_data) + list(self._package_data) + vars_from_init
         )
         d = self._get_pkg_options(d, not_options=not_options)
         path = directory / pkgname / f"{self._pkg_id}-pkgdata.dat"
         d["packagedata"] = path.as_posix()
         # max uzf-cells for which time period data will be supplied
-        d["nuzfcells"] = np.count_nonzero(np.isfinite(d["landflag"]))
+        d["nuzfcells"] = np.count_nonzero(np.isfinite(self.dataset["landflag"]))
         return self._template.render(d)
 
     def _to_struct_array(self, arrdict, layer):
diff --git a/imod/tests/fixtures/mf6_flow_with_transport_fixture.py b/imod/tests/fixtures/mf6_flow_with_transport_fixture.py
index 2f9d312d7..974710216 100644
--- a/imod/tests/fixtures/mf6_flow_with_transport_fixture.py
+++ b/imod/tests/fixtures/mf6_flow_with_transport_fixture.py
@@ -83,6 +83,7 @@ def conductance_fc():
 
     # Constant head
     conductance = xr.full_like(idomain, np.nan)
+    conductance[:, 0, 7, 7:9] = 1.0
     return conductance
 
 
diff --git a/imod/tests/test_mf6/test_ex01_twri.py b/imod/tests/test_mf6/test_ex01_twri.py
index 71b24c480..9cd2663fa 100644
--- a/imod/tests/test_mf6/test_ex01_twri.py
+++ b/imod/tests/test_mf6/test_ex01_twri.py
@@ -319,7 +319,7 @@ def test_wel_render(twri_model_extended, tmp_path):
             end options
 
             begin dimensions
-              maxbound 45
+              maxbound 15
             end dimensions
 
             begin period 1
diff --git a/imod/tests/test_mf6/test_mf6_wel_lowlvl.py b/imod/tests/test_mf6/test_mf6_wel_lowlvl.py
index 4f961305c..7e4ea7957 100644
--- a/imod/tests/test_mf6/test_mf6_wel_lowlvl.py
+++ b/imod/tests/test_mf6/test_mf6_wel_lowlvl.py
@@ -149,7 +149,7 @@ def test_mf6wel_render__transient(
     end options
 
     begin dimensions
-      maxbound 36
+      maxbound 12
     end dimensions
 
     begin period 1