From 582e4e312e9bebb7db9db8e074acb67865e6aa95 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Tue, 16 Dec 2025 16:51:08 +0100
Subject: [PATCH 01/14] Add test to clip and add state at boundary from file

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

diff --git a/imod/tests/test_mf6/test_circle.py b/imod/tests/test_mf6/test_circle.py
index 4b0ac79cf..a3b6e10f5 100644
--- a/imod/tests/test_mf6/test_circle.py
+++ b/imod/tests/test_mf6/test_circle.py
@@ -10,6 +10,7 @@
 from imod.logging import LoggerType, LogLevel
 from imod.mf6.validation_settings import ValidationSettings
 from imod.mf6.write_context import WriteContext
+from imod.prepare.partition import create_partition_labels
 
 
 def test_simulation_write_and_run(circle_model, tmp_path):
@@ -243,3 +244,83 @@ def test_simulation_clip_and_state_at_boundary(circle_model_transport, tmp_path)
     )
     assert head_half.shape == (52, 2, 108)
     assert concentration_half.shape == (52, 2, 108)
+
+
+def test_simulation_clip_and_state_at_boundary__from_file(
+    circle_model_transport, tmp_path
+):
+    # Arrange
+    simulation = circle_model_transport
+    idomain = simulation["GWF_1"]["disv"]["idomain"]
+
+    partition_labels = create_partition_labels(idomain, npartitions=2)
+
+    simulation_split = simulation.split(partition_labels)
+    simulation_split.write(tmp_path / "full")
+    simulation_split.run()
+    head = simulation_split.open_head().reindex_like(idomain)
+    concentration = simulation_split.open_concentration().reindex_like(idomain)
+
+    simulation.dump(tmp_path / "dumped")
+    simulated_from_dump = imod.mf6.Modflow6Simulation.from_file(
+        tmp_path / "dumped" / "circle.toml"
+    )
+
+    paths_head = (tmp_path / "full").glob("**/*.hds")
+    paths_conc = (tmp_path / "full").glob("**/*.ucn")
+    paths_grb = (tmp_path / "full").glob("**/disv.disv.grb")
+
+    heads = []
+    concs = []
+    tstart = simulation["time_discretization"]["time"].isel(time=0).values
+    for path_head, path_conc, path_grb in zip(paths_head, paths_conc, paths_grb):
+        head_part = imod.mf6.open_hds(
+            path_head,
+            path_grb,
+            simulation_start_time=tstart,
+        )
+        heads.append(head_part)
+        conc_part = imod.mf6.open_conc(
+            path_conc,
+            path_grb,
+            simulation_start_time=tstart,
+        )
+        concs.append(conc_part)
+
+    head = xu.merge_partitions(heads)
+    conc = xu.merge_partitions(concs)
+    head = head.reindex_like(idomain)
+    conc = conc.reindex_like(idomain)
+
+    states_for_boundary = {
+        "GWF_1": head["head"].isel(time=-1, drop=True),
+        "transport": concentration["concentration"].isel(time=-1, drop=True),
+    }
+    # Act
+    half_simulation = simulated_from_dump.clip_box(
+        x_max=0.1, states_for_boundary=states_for_boundary
+    )
+    # Assert
+    # Test if model dims halved
+    idomain_half = half_simulation["GWF_1"]["disv"]["idomain"]
+    dim = idomain_half.grid.face_dimension
+    np.testing.assert_array_equal(idomain_half.sizes[dim] / idomain.sizes[dim], 0.5)
+    n_conc = (
+        half_simulation["transport"]["cnc_clipped"]["concentration"]
+        .notnull()
+        .sum()
+        .compute()
+    )
+    n_head = half_simulation["GWF_1"]["chd_clipped"]["head"].notnull().sum().compute()
+    assert n_conc == 24
+    assert n_head == 20
+    # Test if model runs
+    half_simulation.write(tmp_path / "half")
+    half_simulation.run()
+    # Test if the clipped model output has the correct dimension
+    head_half = half_simulation.open_head().compute().reindex_like(idomain_half)
+    concentration_half = (
+        half_simulation.open_concentration().compute().reindex_like(idomain_half)
+    )
+    assert head_half.shape == (52, 2, 108)
+    assert concentration_half.shape == (52, 2, 108)

From 47298c8e456502d9225e78e1bf142bcb8970732b Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Thu, 18 Dec 2025 16:57:37 +0100
Subject: [PATCH 02/14] Add test with transient_chd

---
 imod/tests/test_mf6/test_circle.py | 47 ++++++++++++++++++++++++++++--
 1 file changed, 44 insertions(+), 3 deletions(-)

diff --git a/imod/tests/test_mf6/test_circle.py b/imod/tests/test_mf6/test_circle.py
index a3b6e10f5..dd731cfbe 100644
--- a/imod/tests/test_mf6/test_circle.py
+++ b/imod/tests/test_mf6/test_circle.py
@@ -246,6 +246,49 @@ def test_simulation_clip_and_state_at_boundary(circle_model_transport, tmp_path)
     assert concentration_half.shape == (52, 2, 108)
 
 
+def test_simulation_clip_and_state_at_boundary__transient_chd(circle_model_transport, tmp_path):
+    # Arrange
+    simulation = circle_model_transport
+    idomain = simulation["GWF_1"]["disv"]["idomain"].compute()
+    time = simulation["time_discretization"].dataset.coords["time"]
+
+    simulation["GWF_1"]["chd"].dataset["head"] = simulation["GWF_1"]["chd"].dataset["head"].expand_dims(time=time.values)
+
+    simulation.write(tmp_path / "full")
+    simulation.run()
+    head = simulation.open_head().compute().reindex_like(idomain)
+    concentration = simulation.open_concentration().compute().reindex_like(idomain)
+
+    states_for_boundary = {
+        "GWF_1": head.isel(time=-1, drop=True),
+        "transport": concentration.isel(time=-1, drop=True),
+    }
+    # Act
+    half_simulation = simulation.clip_box(
+        x_max=0.1, states_for_boundary=states_for_boundary
+    )
+    # Assert
+    # Test if model dims halved
+    idomain_half = half_simulation["GWF_1"]["disv"]["idomain"]
+    dim = idomain_half.grid.face_dimension
+    np.testing.assert_array_equal(idomain_half.sizes[dim] / idomain.sizes[dim], 0.5)
+    assert (
+        half_simulation["transport"]["cnc_clipped"]["concentration"].notnull().sum()
+        == 20
+    )
+    assert half_simulation["GWF_1"]["chd_clipped"]["head"].notnull().sum() == 20
+    # Test if model runs
+    half_simulation.write(tmp_path / "half")
+    half_simulation.run()
+    # Test if the clipped model output has the correct dimension
+    head_half = half_simulation.open_head().compute().reindex_like(idomain_half)
+    concentration_half = (
+        half_simulation.open_concentration().compute().reindex_like(idomain_half)
+    )
+    assert head_half.shape == (52, 2, 108)
+    assert concentration_half.shape == (52, 2, 108)
+
+
 def test_simulation_clip_and_state_at_boundary__from_file(
     circle_model_transport, tmp_path
 ):
@@ -258,8 +301,6 @@ def test_simulation_clip_and_state_at_boundary__from_file(
     simulation_split = simulation.split(partition_labels)
     simulation_split.write(tmp_path / "full")
     simulation_split.run()
-    head = simulation_split.open_head().reindex_like(idomain)
-    concentration = simulation_split.open_concentration().reindex_like(idomain)
 
     simulation.dump(tmp_path / "dumped")
     simulated_from_dump = imod.mf6.Modflow6Simulation.from_file(
@@ -294,7 +335,7 @@ def test_simulation_clip_and_state_at_boundary__from_file(
 
     states_for_boundary = {
         "GWF_1": head["head"].isel(time=-1, drop=True),
-        "transport": concentration["concentration"].isel(time=-1, drop=True),
+        "transport": conc["concentration"].isel(time=-1, drop=True),
     }
     # Act
     half_simulation = simulated_from_dump.clip_box(

From 217d484d685ebacb58f204f92cd10f78b905d393 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Thu, 18 Dec 2025 17:44:04 +0100
Subject: [PATCH 03/14] Fix cnc template for stress periods

---
 imod/templates/mf6/gwf-cnc.j2 | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/imod/templates/mf6/gwf-cnc.j2 b/imod/templates/mf6/gwf-cnc.j2
index 4ee8634b5..9de057ae4 100644
--- a/imod/templates/mf6/gwf-cnc.j2
+++ b/imod/templates/mf6/gwf-cnc.j2
@@ -25,4 +25,5 @@ end dimensions
 
 {% for i, path in periods.items() %}begin period {{i}}
   open/close {{path}}{% if binary %} (binary){% endif %}
-end period{% endfor %}
\ No newline at end of file
+end period
+{% endfor %}
\ No newline at end of file

From 517a3f4465b0983b100838c1589104639b4242d1 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Thu, 18 Dec 2025 18:04:20 +0100
Subject: [PATCH 04/14] Enforce dimension order before creating clipped
 boundary package

---
 imod/mf6/utilities/clipped_bc_creator.py | 22 ++++++++++++++++++----
 1 file changed, 18 insertions(+), 4 deletions(-)

diff --git a/imod/mf6/utilities/clipped_bc_creator.py b/imod/mf6/utilities/clipped_bc_creator.py
index 30662ae82..144e274eb 100644
--- a/imod/mf6/utilities/clipped_bc_creator.py
+++ b/imod/mf6/utilities/clipped_bc_creator.py
@@ -3,6 +3,7 @@
 from imod.mf6 import ConstantConcentration, ConstantHead
 from imod.select.grid import active_grid_boundary_xy
 from imod.typing import GridDataArray
+from imod.util.dims import enforced_dim_order
 
 StateType: TypeAlias = ConstantHead | ConstantConcentration
 StateClassType: TypeAlias = type[ConstantHead] | type[ConstantConcentration]
@@ -24,6 +25,19 @@ def _find_unassigned_grid_boundaries(
     return unassigned_grid_boundaries
 
 
+@enforced_dim_order
+def _create_clipped_boundary_state(
+    idomain: GridDataArray,
+    state_for_clipped_boundary: GridDataArray,
+    original_constant_head_boundaries: list[StateType],
+):
+    """"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(
+        active_grid_boundary, original_constant_head_boundaries
+    )
+    return state_for_clipped_boundary.where(unassigned_grid_boundaries)
+
 def create_clipped_boundary(
     idomain: GridDataArray,
     state_for_clipped_boundary: GridDataArray,
@@ -52,10 +66,10 @@ def create_clipped_boundary(
         packages
 
     """
-    active_grid_boundary = active_grid_boundary_xy(idomain > 0)
-    unassigned_grid_boundaries = _find_unassigned_grid_boundaries(
-        active_grid_boundary, original_constant_head_boundaries
+    constant_state = _create_clipped_boundary_state(
+        idomain,
+        state_for_clipped_boundary,
+        original_constant_head_boundaries,
     )
-    constant_state = state_for_clipped_boundary.where(unassigned_grid_boundaries)
 
     return pkg_type(constant_state, print_input=True, print_flows=True, save_flows=True)

From 859cbfe47fd5db70fddf0fe261a8cd13e39a7df1 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Thu, 18 Dec 2025 18:05:07 +0100
Subject: [PATCH 05/14] Add tests

---
 imod/tests/test_mf6/test_circle.py | 84 ++++++++++--------------------
 1 file changed, 27 insertions(+), 57 deletions(-)

diff --git a/imod/tests/test_mf6/test_circle.py b/imod/tests/test_mf6/test_circle.py
index dd731cfbe..fc5424084 100644
--- a/imod/tests/test_mf6/test_circle.py
+++ b/imod/tests/test_mf6/test_circle.py
@@ -245,19 +245,19 @@ def test_simulation_clip_and_state_at_boundary(circle_model_transport, tmp_path)
     assert head_half.shape == (52, 2, 108)
     assert concentration_half.shape == (52, 2, 108)
 
-
-def test_simulation_clip_and_state_at_boundary__transient_chd(circle_model_transport, tmp_path):
+def test_simulation_clip_and_constant_state_at_boundary__transient_chd(circle_model_transport, tmp_path):
     # Arrange
     simulation = circle_model_transport
     idomain = simulation["GWF_1"]["disv"]["idomain"].compute()
-    time = simulation["time_discretization"].dataset.coords["time"]
+    time = simulation["time_discretization"].dataset.coords["time"].values
 
-    simulation["GWF_1"]["chd"].dataset["head"] = simulation["GWF_1"]["chd"].dataset["head"].expand_dims(time=time.values)
+    simulation["GWF_1"]["chd"].dataset["head"] = simulation["GWF_1"]["chd"].dataset["head"].expand_dims(time=time)
 
     simulation.write(tmp_path / "full")
     simulation.run()
-    head = simulation.open_head().compute().reindex_like(idomain)
-    concentration = simulation.open_concentration().compute().reindex_like(idomain)
+    # 
+    head = simulation.open_head(simulation_start_time=time[0]).compute().reindex_like(idomain)
+    concentration = simulation.open_concentration(simulation_start_time=time[0]).compute().reindex_like(idomain)
 
     states_for_boundary = {
         "GWF_1": head.isel(time=-1, drop=True),
@@ -276,7 +276,7 @@ def test_simulation_clip_and_state_at_boundary__transient_chd(circle_model_trans
         half_simulation["transport"]["cnc_clipped"]["concentration"].notnull().sum()
         == 20
     )
-    assert half_simulation["GWF_1"]["chd_clipped"]["head"].notnull().sum() == 20
+    assert half_simulation["GWF_1"]["chd_clipped"]["head"].notnull().sum() == 20 * 52
     # Test if model runs
     half_simulation.write(tmp_path / "half")
     half_simulation.run()
@@ -289,56 +289,26 @@ def test_simulation_clip_and_state_at_boundary__transient_chd(circle_model_trans
     assert concentration_half.shape == (52, 2, 108)
 
 
-def test_simulation_clip_and_state_at_boundary__from_file(
-    circle_model_transport, tmp_path
-):
+def test_simulation_clip_and_transient_state_at_boundary__transient_chd(circle_model_transport, tmp_path):
     # Arrange
     simulation = circle_model_transport
-    idomain = simulation["GWF_1"]["disv"]["idomain"]
-
-    partition_labels = create_partition_labels(idomain, npartitions=2)
-
-    simulation_split = simulation.split(partition_labels)
-    simulation_split.write(tmp_path / "full")
-    simulation_split.run()
+    idomain = simulation["GWF_1"]["disv"]["idomain"].compute()
+    time = simulation["time_discretization"].dataset.coords["time"].values
 
-    simulation.dump(tmp_path / "dumped")
-    simulated_from_dump = imod.mf6.Modflow6Simulation.from_file(
-        tmp_path / "dumped" / "circle.toml"
-    )
+    simulation["GWF_1"]["chd"].dataset["head"] = simulation["GWF_1"]["chd"].dataset["head"].expand_dims(time=time)
 
-    paths_head = (tmp_path / "full").glob("**/*.hds")
-    paths_conc = (tmp_path / "full").glob("**/*.ucn")
-    paths_grb = (tmp_path / "full").glob("**/disv.disv.grb")
-
-    heads = []
-    concs = []
-    tstart = simulation["time_discretization"]["time"].isel(time=0).values
-    for path_head, path_conc, path_grb in zip(paths_head, paths_conc, paths_grb):
-        head_part = imod.mf6.open_hds(
-            path_head,
-            path_grb,
-            simulation_start_time=tstart,
-        )
-        heads.append(head_part)
-        conc_part = imod.mf6.open_conc(
-            path_conc,
-            path_grb,
-            simulation_start_time=tstart,
-        )
-        concs.append(conc_part)
-
-    head = xu.merge_partitions(heads)
-    conc = xu.merge_partitions(concs)
-    head = head.reindex_like(idomain)
-    conc = conc.reindex_like(idomain)
+    simulation.write(tmp_path / "full")
+    simulation.run()
+    # 
+    head = simulation.open_head(simulation_start_time=time[0]).compute().reindex_like(idomain)
+    concentration = simulation.open_concentration(simulation_start_time=time[0]).compute().reindex_like(idomain)
 
     states_for_boundary = {
-        "GWF_1": head["head"].isel(time=-1, drop=True),
-        "transport": conc["concentration"].isel(time=-1, drop=True),
+        "GWF_1": head.isel(time=slice(12, None)),
+        "transport": concentration.isel(time=slice(12, None)),
     }
     # Act
-    half_simulation = simulated_from_dump.clip_box(
+    half_simulation = simulation.clip_box(
         x_max=0.1, states_for_boundary=states_for_boundary
     )
     # Assert
@@ -346,15 +316,14 @@ def test_simulation_clip_and_state_at_boundary__from_file(
     idomain_half = half_simulation["GWF_1"]["disv"]["idomain"]
     dim = idomain_half.grid.face_dimension
     np.testing.assert_array_equal(idomain_half.sizes[dim] / idomain.sizes[dim], 0.5)
-    n_conc = (
-        half_simulation["transport"]["cnc_clipped"]["concentration"]
-        .notnull()
-        .sum()
-        .compute()
+    # Conc data for 40 time steps (from 12 to 52)
+    assert (
+        half_simulation["transport"]["cnc_clipped"]["concentration"].notnull().sum()
+        == 20 * 40
     )
-    n_head = half_simulation["GWF_1"]["chd_clipped"]["head"].notnull().sum().compute()
-    assert n_conc == 24
-    assert n_head == 20
+    # Head data for 39 time steps (from 12 to 51) as the last time step of
+    # results is not present in the chd package already present in the model.
+    assert half_simulation["GWF_1"]["chd_clipped"]["head"].notnull().sum() == 20 * 39
     # Test if model runs
     half_simulation.write(tmp_path / "half")
     half_simulation.run()
@@ -365,3 +334,4 @@ def test_simulation_clip_and_state_at_boundary__from_file(
     )
     assert head_half.shape == (52, 2, 108)
     assert concentration_half.shape == (52, 2, 108)
+

From 312a022dd0e75e0b6e2f035facb1555aa2e09df6 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Fri, 19 Dec 2025 10:45:13 +0100
Subject: [PATCH 06/14] format

---
 imod/mf6/utilities/clipped_bc_creator.py |  3 +-
 imod/tests/test_mf6/test_circle.py       | 47 ++++++++++++++++++------
 2 files changed, 37 insertions(+), 13 deletions(-)

diff --git a/imod/mf6/utilities/clipped_bc_creator.py b/imod/mf6/utilities/clipped_bc_creator.py
index 144e274eb..9713b31bb 100644
--- a/imod/mf6/utilities/clipped_bc_creator.py
+++ b/imod/mf6/utilities/clipped_bc_creator.py
@@ -31,13 +31,14 @@ def _create_clipped_boundary_state(
     state_for_clipped_boundary: GridDataArray,
     original_constant_head_boundaries: list[StateType],
 ):
-    """"Helper function to make sure dimension order is enforced"""
+    """ "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(
         active_grid_boundary, original_constant_head_boundaries
     )
     return state_for_clipped_boundary.where(unassigned_grid_boundaries)
 
+
 def create_clipped_boundary(
     idomain: GridDataArray,
     state_for_clipped_boundary: GridDataArray,
diff --git a/imod/tests/test_mf6/test_circle.py b/imod/tests/test_mf6/test_circle.py
index fc5424084..d853acff2 100644
--- a/imod/tests/test_mf6/test_circle.py
+++ b/imod/tests/test_mf6/test_circle.py
@@ -10,7 +10,6 @@
 from imod.logging import LoggerType, LogLevel
 from imod.mf6.validation_settings import ValidationSettings
 from imod.mf6.write_context import WriteContext
-from imod.prepare.partition import create_partition_labels
 
 
 def test_simulation_write_and_run(circle_model, tmp_path):
@@ -245,19 +244,32 @@ def test_simulation_clip_and_state_at_boundary(circle_model_transport, tmp_path)
     assert head_half.shape == (52, 2, 108)
     assert concentration_half.shape == (52, 2, 108)
 
-def test_simulation_clip_and_constant_state_at_boundary__transient_chd(circle_model_transport, tmp_path):
+
+def test_simulation_clip_and_constant_state_at_boundary__transient_chd(
+    circle_model_transport, tmp_path
+):
     # Arrange
     simulation = circle_model_transport
     idomain = simulation["GWF_1"]["disv"]["idomain"].compute()
     time = simulation["time_discretization"].dataset.coords["time"].values
 
-    simulation["GWF_1"]["chd"].dataset["head"] = simulation["GWF_1"]["chd"].dataset["head"].expand_dims(time=time)
+    simulation["GWF_1"]["chd"].dataset["head"] = (
+        simulation["GWF_1"]["chd"].dataset["head"].expand_dims(time=time)
+    )
 
     simulation.write(tmp_path / "full")
     simulation.run()
-    # 
-    head = simulation.open_head(simulation_start_time=time[0]).compute().reindex_like(idomain)
-    concentration = simulation.open_concentration(simulation_start_time=time[0]).compute().reindex_like(idomain)
+    #
+    head = (
+        simulation.open_head(simulation_start_time=time[0])
+        .compute()
+        .reindex_like(idomain)
+    )
+    concentration = (
+        simulation.open_concentration(simulation_start_time=time[0])
+        .compute()
+        .reindex_like(idomain)
+    )
 
     states_for_boundary = {
         "GWF_1": head.isel(time=-1, drop=True),
@@ -289,19 +301,31 @@ def test_simulation_clip_and_constant_state_at_boundary__transient_chd(circle_mo
     assert concentration_half.shape == (52, 2, 108)
 
 
-def test_simulation_clip_and_transient_state_at_boundary__transient_chd(circle_model_transport, tmp_path):
+def test_simulation_clip_and_transient_state_at_boundary__transient_chd(
+    circle_model_transport, tmp_path
+):
     # Arrange
     simulation = circle_model_transport
     idomain = simulation["GWF_1"]["disv"]["idomain"].compute()
     time = simulation["time_discretization"].dataset.coords["time"].values
 
-    simulation["GWF_1"]["chd"].dataset["head"] = simulation["GWF_1"]["chd"].dataset["head"].expand_dims(time=time)
+    simulation["GWF_1"]["chd"].dataset["head"] = (
+        simulation["GWF_1"]["chd"].dataset["head"].expand_dims(time=time)
+    )
 
     simulation.write(tmp_path / "full")
     simulation.run()
-    # 
-    head = simulation.open_head(simulation_start_time=time[0]).compute().reindex_like(idomain)
-    concentration = simulation.open_concentration(simulation_start_time=time[0]).compute().reindex_like(idomain)
+    #
+    head = (
+        simulation.open_head(simulation_start_time=time[0])
+        .compute()
+        .reindex_like(idomain)
+    )
+    concentration = (
+        simulation.open_concentration(simulation_start_time=time[0])
+        .compute()
+        .reindex_like(idomain)
+    )
 
     states_for_boundary = {
         "GWF_1": head.isel(time=slice(12, None)),
@@ -334,4 +358,3 @@ def test_simulation_clip_and_transient_state_at_boundary__transient_chd(circle_m
     )
     assert head_half.shape == (52, 2, 108)
     assert concentration_half.shape == (52, 2, 108)
-

From 6a5db959f43935ebcbe46346f04e48f3b1766544 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Fri, 19 Dec 2025 17:26:47 +0100
Subject: [PATCH 07/14] Fix broadcasting trouble in clipped_bc_creator.py

---
 imod/mf6/model.py                        | 24 ++++++-
 imod/mf6/utilities/clipped_bc_creator.py | 84 +++++++++++++++++++++++-
 imod/tests/test_mf6/test_circle.py       | 22 +++----
 3 files changed, 113 insertions(+), 17 deletions(-)

diff --git a/imod/mf6/model.py b/imod/mf6/model.py
index ee0de6da6..3fdd91dfc 100644
--- a/imod/mf6/model.py
+++ b/imod/mf6/model.py
@@ -65,7 +65,7 @@ def _create_boundary_condition_for_unassigned_boundary(
     model: Modflow6Model,
     state_for_boundary: Optional[GridDataArray],
     additional_boundaries: Optional[list[StateClassType]] = None,
-):
+) -> Optional[StateClassType]:
     if state_for_boundary is None:
         return None
 
@@ -90,7 +90,7 @@ def _create_boundary_condition_clipped_boundary(
     clipped_model: Modflow6Model,
     state_for_boundary: Optional[GridDataArray],
     clip_box_args: tuple,
-):
+) -> Optional[StateClassType]:
     # 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
@@ -120,10 +120,28 @@ def _create_boundary_condition_clipped_boundary(
     else:
         state_for_boundary_clipped = None
 
-    return _create_boundary_condition_for_unassigned_boundary(
+    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 and twice and it could result in broadcasting
+    # errors in the second call otherwise 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
+    ):
+        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, IModel, abc.ABC):
     _mandatory_packages: tuple[str, ...] = ()
diff --git a/imod/mf6/utilities/clipped_bc_creator.py b/imod/mf6/utilities/clipped_bc_creator.py
index 9713b31bb..8e3182334 100644
--- a/imod/mf6/utilities/clipped_bc_creator.py
+++ b/imod/mf6/utilities/clipped_bc_creator.py
@@ -1,4 +1,6 @@
-from typing import TypeAlias
+from typing import Optional, Tuple, TypeAlias
+
+import xarray as xr
 
 from imod.mf6 import ConstantConcentration, ConstantHead
 from imod.select.grid import active_grid_boundary_xy
@@ -13,7 +15,7 @@ def _find_unassigned_grid_boundaries(
     active_grid_boundary: GridDataArray,
     boundary_conditions: list[StateType],
 ) -> GridDataArray:
-    unassigned_grid_boundaries = active_grid_boundary
+    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
         # "concentration".
@@ -25,17 +27,93 @@ def _find_unassigned_grid_boundaries(
     return unassigned_grid_boundaries
 
 
+def _align_time_indexes_boundaries(
+    state_for_clipped_boundary: GridDataArray,
+    unassigned_grid_boundaries: GridDataArray,
+) -> Tuple[GridDataArray, GridDataArray]:
+    """
+    Create an outer time index for aligning boundaries. Furthermore deal with
+    cases where one or both boundaries don't have a time dimension. In a graphic
+    way, we want this:
+
+    State
+       a-----b-----c
+    Unassigned
+    d-----e
+
+    Needs to align to:
+    d---a--e--b-----c
+    """
+    index_unassigned = unassigned_grid_boundaries.indexes
+    index_state = state_for_clipped_boundary.indexes
+    if "time" in index_unassigned and "time" in index_state:
+        return index_unassigned["time"].join(index_state["time"], how="outer")
+    elif "time" in index_unassigned:
+        return index_unassigned["time"]
+    elif "time" in index_state:
+        return index_state["time"]
+    else:
+        return None
+
+
+def _align_boundaries(
+    state_for_clipped_boundary: GridDataArray,
+    unassigned_grid_boundaries: Optional[GridDataArray],
+) -> Tuple[GridDataArray, Optional[GridDataArray]]:
+    """
+    Customly align the state_for_clipped_boundary and unassigned grid boundaries.
+    - "layer" dimension requires outer alignment
+    - "time" requires reindexing with ffill
+    - planar coordinates are expected to be aligned already
+    """
+    # Align dimensions
+    if unassigned_grid_boundaries is not None:
+        # Align along layer dimension with outer join, xarray API only supports
+        # excluding dims, not specifying dims to align along, so we have to do
+        # it this way.
+        dims_to_exclude = set(state_for_clipped_boundary.dims) | set(
+            unassigned_grid_boundaries.dims
+        )
+        dims_to_exclude.remove("layer")
+        state_for_clipped_boundary, unassigned_grid_boundaries = xr.align(
+            state_for_clipped_boundary,
+            unassigned_grid_boundaries,
+            join="outer",
+            exclude=dims_to_exclude,
+        )
+        # Align along time dimension by finding the outer time indexes and then
+        # reindexing with a ffill.
+        outer_time_index = _align_time_indexes_boundaries(
+            state_for_clipped_boundary, unassigned_grid_boundaries
+        )
+        if "time" in state_for_clipped_boundary.indexes:
+            state_for_clipped_boundary = state_for_clipped_boundary.reindex(
+                {"time": outer_time_index}, method="ffill"
+            )
+        if "time" in unassigned_grid_boundaries.indexes:
+            unassigned_grid_boundaries = unassigned_grid_boundaries.reindex(
+                {"time": outer_time_index}, method="ffill"
+            )
+
+    return state_for_clipped_boundary, unassigned_grid_boundaries
+
+
 @enforced_dim_order
 def _create_clipped_boundary_state(
     idomain: GridDataArray,
     state_for_clipped_boundary: GridDataArray,
     original_constant_head_boundaries: list[StateType],
 ):
-    """ "Helper function to make sure dimension order is enforced"""
+    """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(
         active_grid_boundary, original_constant_head_boundaries
     )
+
+    state_for_clipped_boundary, unassigned_grid_boundaries = _align_boundaries(
+        state_for_clipped_boundary, unassigned_grid_boundaries
+    )
+
     return state_for_clipped_boundary.where(unassigned_grid_boundaries)
 
 
diff --git a/imod/tests/test_mf6/test_circle.py b/imod/tests/test_mf6/test_circle.py
index d853acff2..e3e93887e 100644
--- a/imod/tests/test_mf6/test_circle.py
+++ b/imod/tests/test_mf6/test_circle.py
@@ -310,7 +310,7 @@ def test_simulation_clip_and_transient_state_at_boundary__transient_chd(
     time = simulation["time_discretization"].dataset.coords["time"].values
 
     simulation["GWF_1"]["chd"].dataset["head"] = (
-        simulation["GWF_1"]["chd"].dataset["head"].expand_dims(time=time)
+        simulation["GWF_1"]["chd"].dataset["head"].expand_dims(time=time[:5])
     )
 
     simulation.write(tmp_path / "full")
@@ -328,8 +328,8 @@ def test_simulation_clip_and_transient_state_at_boundary__transient_chd(
     )
 
     states_for_boundary = {
-        "GWF_1": head.isel(time=slice(12, None)),
-        "transport": concentration.isel(time=slice(12, None)),
+        "GWF_1": head.isel(time=slice(11, -11)),
+        "transport": concentration.isel(time=slice(11, -11)),
     }
     # Act
     half_simulation = simulation.clip_box(
@@ -340,14 +340,14 @@ def test_simulation_clip_and_transient_state_at_boundary__transient_chd(
     idomain_half = half_simulation["GWF_1"]["disv"]["idomain"]
     dim = idomain_half.grid.face_dimension
     np.testing.assert_array_equal(idomain_half.sizes[dim] / idomain.sizes[dim], 0.5)
-    # Conc data for 40 time steps (from 12 to 52)
-    assert (
-        half_simulation["transport"]["cnc_clipped"]["concentration"].notnull().sum()
-        == 20 * 40
-    )
-    # Head data for 39 time steps (from 12 to 51) as the last time step of
-    # results is not present in the chd package already present in the model.
-    assert half_simulation["GWF_1"]["chd_clipped"]["head"].notnull().sum() == 20 * 39
+    # Transport data for just 30 time steps (from 12 up to 41)
+    conc_clipped = half_simulation["transport"]["cnc_clipped"]["concentration"]
+    assert conc_clipped.sizes["time"] == 30
+    assert conc_clipped.notnull().sum() == 20 * 30
+    # Head data for chd for just 30 time steps (from 12 up to 41)
+    head_clipped = half_simulation["GWF_1"]["chd_clipped"]["head"]
+    assert head_clipped.sizes["time"] == 30
+    assert head_clipped.notnull().sum() == 20 * 30
     # Test if model runs
     half_simulation.write(tmp_path / "half")
     half_simulation.run()

From 05f8b01adb0d05c698ea62028c238e670acf58fe Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Fri, 19 Dec 2025 17:35:24 +0100
Subject: [PATCH 08/14] Fix mypy errors

---
 imod/mf6/model.py                        | 13 ++++++++-----
 imod/mf6/utilities/clipped_bc_creator.py |  2 +-
 2 files changed, 9 insertions(+), 6 deletions(-)

diff --git a/imod/mf6/model.py b/imod/mf6/model.py
index 3fdd91dfc..67cd81760 100644
--- a/imod/mf6/model.py
+++ b/imod/mf6/model.py
@@ -40,6 +40,7 @@
 from imod.mf6.riv import River
 from imod.mf6.utilities.clipped_bc_creator import (
     StateClassType,
+    StateType,
     create_clipped_boundary,
 )
 from imod.mf6.utilities.mf6hfb import merge_hfb_packages
@@ -64,8 +65,8 @@ def pkg_has_cleanup(pkg: Package):
 def _create_boundary_condition_for_unassigned_boundary(
     model: Modflow6Model,
     state_for_boundary: Optional[GridDataArray],
-    additional_boundaries: Optional[list[StateClassType]] = None,
-) -> Optional[StateClassType]:
+    additional_boundaries: list[Optional[StateType]] = [None],
+) -> Optional[StateType]:
     if state_for_boundary is None:
         return None
 
@@ -90,7 +91,7 @@ def _create_boundary_condition_clipped_boundary(
     clipped_model: Modflow6Model,
     state_for_boundary: Optional[GridDataArray],
     clip_box_args: tuple,
-) -> Optional[StateClassType]:
+) -> 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
@@ -132,8 +133,10 @@ def _create_boundary_condition_clipped_boundary(
     # this function is called and twice and it could result in broadcasting
     # errors in the second call otherwise 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
+    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(
diff --git a/imod/mf6/utilities/clipped_bc_creator.py b/imod/mf6/utilities/clipped_bc_creator.py
index 8e3182334..ea6f0f415 100644
--- a/imod/mf6/utilities/clipped_bc_creator.py
+++ b/imod/mf6/utilities/clipped_bc_creator.py
@@ -30,7 +30,7 @@ def _find_unassigned_grid_boundaries(
 def _align_time_indexes_boundaries(
     state_for_clipped_boundary: GridDataArray,
     unassigned_grid_boundaries: GridDataArray,
-) -> Tuple[GridDataArray, GridDataArray]:
+) -> Optional[GridDataArray]:
     """
     Create an outer time index for aligning boundaries. Furthermore deal with
     cases where one or both boundaries don't have a time dimension. In a graphic

From e64ef49b8abde42954fd9764f89851dc9a9610a9 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Fri, 19 Dec 2025 17:38:02 +0100
Subject: [PATCH 09/14] Fix typos

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

diff --git a/imod/mf6/model.py b/imod/mf6/model.py
index 67cd81760..3af1bc219 100644
--- a/imod/mf6/model.py
+++ b/imod/mf6/model.py
@@ -130,9 +130,9 @@ def _create_boundary_condition_clipped_boundary(
     # 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 and twice and it could result in broadcasting
-    # errors in the second call otherwise if the time domain of
-    # state_for_boundary and assigned packages have no overlap.
+    # 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)

From cc8849c9883e396ac4d78132a374d5e92b903c99 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Fri, 19 Dec 2025 17:53:23 +0100
Subject: [PATCH 10/14] Clean up tests somewhat

---
 imod/tests/test_mf6/test_circle.py | 98 +++++++++++++-----------------
 1 file changed, 42 insertions(+), 56 deletions(-)

diff --git a/imod/tests/test_mf6/test_circle.py b/imod/tests/test_mf6/test_circle.py
index e3e93887e..4c761214c 100644
--- a/imod/tests/test_mf6/test_circle.py
+++ b/imod/tests/test_mf6/test_circle.py
@@ -205,15 +205,30 @@ def test_simulation_write_and_run_transport_vsc(circle_model_transport_vsc, tmp_
     assert concentration.shape == (52, 2, 216)
 
 
+def run_simulation(simulation, modeldir):
+    idomain = simulation["GWF_1"]["disv"]["idomain"].compute()
+    time = simulation["time_discretization"].dataset.coords["time"].values
+    simulation.write(modeldir)
+    simulation.run()
+    head = (
+        simulation.open_head(simulation_start_time=time[0])
+        .compute()
+        .reindex_like(idomain)
+    )
+    concentration = (
+        simulation.open_concentration()
+        .compute(simulation_start_time=time[0])
+        .reindex_like(idomain)
+    )
+    return head, concentration
+
+
 def test_simulation_clip_and_state_at_boundary(circle_model_transport, tmp_path):
     # Arrange
     simulation = circle_model_transport
     idomain = simulation["GWF_1"]["disv"]["idomain"].compute()
 
-    simulation.write(tmp_path / "full")
-    simulation.run()
-    head = simulation.open_head().compute().reindex_like(idomain)
-    concentration = simulation.open_concentration().compute().reindex_like(idomain)
+    head, concentration = run_simulation(simulation, tmp_path / "full")
 
     states_for_boundary = {
         "GWF_1": head.isel(time=-1, drop=True),
@@ -233,14 +248,9 @@ def test_simulation_clip_and_state_at_boundary(circle_model_transport, tmp_path)
         == 20
     )
     assert half_simulation["GWF_1"]["chd_clipped"]["head"].notnull().sum() == 20
+
     # Test if model runs
-    half_simulation.write(tmp_path / "half")
-    half_simulation.run()
-    # Test if the clipped model output has the correct dimension
-    head_half = half_simulation.open_head().compute().reindex_like(idomain_half)
-    concentration_half = (
-        half_simulation.open_concentration().compute().reindex_like(idomain_half)
-    )
+    head_half, concentration_half = run_simulation(half_simulation, tmp_path / "half")
     assert head_half.shape == (52, 2, 108)
     assert concentration_half.shape == (52, 2, 108)
 
@@ -257,19 +267,7 @@ def test_simulation_clip_and_constant_state_at_boundary__transient_chd(
         simulation["GWF_1"]["chd"].dataset["head"].expand_dims(time=time)
     )
 
-    simulation.write(tmp_path / "full")
-    simulation.run()
-    #
-    head = (
-        simulation.open_head(simulation_start_time=time[0])
-        .compute()
-        .reindex_like(idomain)
-    )
-    concentration = (
-        simulation.open_concentration(simulation_start_time=time[0])
-        .compute()
-        .reindex_like(idomain)
-    )
+    head, concentration = run_simulation(simulation, tmp_path / "full")
 
     states_for_boundary = {
         "GWF_1": head.isel(time=-1, drop=True),
@@ -284,19 +282,17 @@ def test_simulation_clip_and_constant_state_at_boundary__transient_chd(
     idomain_half = half_simulation["GWF_1"]["disv"]["idomain"]
     dim = idomain_half.grid.face_dimension
     np.testing.assert_array_equal(idomain_half.sizes[dim] / idomain.sizes[dim], 0.5)
-    assert (
-        half_simulation["transport"]["cnc_clipped"]["concentration"].notnull().sum()
-        == 20
-    )
-    assert half_simulation["GWF_1"]["chd_clipped"]["head"].notnull().sum() == 20 * 52
+    # Test if the clipped boundary conditions have the correct dimension and values
+    conc_clipped = half_simulation["transport"]["cnc_clipped"]["concentration"]
+    assert conc_clipped.sizes["layer"] == 2
+    assert "time" not in conc_clipped.dims
+    assert conc_clipped.notnull().sum() == 20
+    head_clipped = half_simulation["GWF_1"]["chd_clipped"]["head"]
+    assert head_clipped.sizes["layer"] == 2
+    assert head_clipped.sizes["time"] == 52
+    assert head_clipped.notnull().sum() == 20 * 52
     # Test if model runs
-    half_simulation.write(tmp_path / "half")
-    half_simulation.run()
-    # Test if the clipped model output has the correct dimension
-    head_half = half_simulation.open_head().compute().reindex_like(idomain_half)
-    concentration_half = (
-        half_simulation.open_concentration().compute().reindex_like(idomain_half)
-    )
+    head_half, concentration_half = run_simulation(half_simulation, tmp_path / "half")
     assert head_half.shape == (52, 2, 108)
     assert concentration_half.shape == (52, 2, 108)
 
@@ -304,6 +300,11 @@ def test_simulation_clip_and_constant_state_at_boundary__transient_chd(
 def test_simulation_clip_and_transient_state_at_boundary__transient_chd(
     circle_model_transport, tmp_path
 ):
+    """
+    Test with a transient chd boundary condition and transient states for
+    boundary conditions. The time steps of the chd package are outside of the
+    time domain of the states_for_boundary.
+    """
     # Arrange
     simulation = circle_model_transport
     idomain = simulation["GWF_1"]["disv"]["idomain"].compute()
@@ -313,19 +314,7 @@ def test_simulation_clip_and_transient_state_at_boundary__transient_chd(
         simulation["GWF_1"]["chd"].dataset["head"].expand_dims(time=time[:5])
     )
 
-    simulation.write(tmp_path / "full")
-    simulation.run()
-    #
-    head = (
-        simulation.open_head(simulation_start_time=time[0])
-        .compute()
-        .reindex_like(idomain)
-    )
-    concentration = (
-        simulation.open_concentration(simulation_start_time=time[0])
-        .compute()
-        .reindex_like(idomain)
-    )
+    head, concentration = run_simulation(simulation, tmp_path / "full")
 
     states_for_boundary = {
         "GWF_1": head.isel(time=slice(11, -11)),
@@ -340,21 +329,18 @@ def test_simulation_clip_and_transient_state_at_boundary__transient_chd(
     idomain_half = half_simulation["GWF_1"]["disv"]["idomain"]
     dim = idomain_half.grid.face_dimension
     np.testing.assert_array_equal(idomain_half.sizes[dim] / idomain.sizes[dim], 0.5)
+    # Test if the clipped boundary conditions have the correct dimension and values
     # Transport data for just 30 time steps (from 12 up to 41)
     conc_clipped = half_simulation["transport"]["cnc_clipped"]["concentration"]
+    assert conc_clipped.sizes["layer"] == 2
     assert conc_clipped.sizes["time"] == 30
     assert conc_clipped.notnull().sum() == 20 * 30
     # Head data for chd for just 30 time steps (from 12 up to 41)
     head_clipped = half_simulation["GWF_1"]["chd_clipped"]["head"]
+    assert head_clipped.sizes["layer"] == 2
     assert head_clipped.sizes["time"] == 30
     assert head_clipped.notnull().sum() == 20 * 30
     # Test if model runs
-    half_simulation.write(tmp_path / "half")
-    half_simulation.run()
-    # Test if the clipped model output has the correct dimension
-    head_half = half_simulation.open_head().compute().reindex_like(idomain_half)
-    concentration_half = (
-        half_simulation.open_concentration().compute().reindex_like(idomain_half)
-    )
+    head_half, concentration_half = run_simulation(half_simulation, tmp_path / "half")
     assert head_half.shape == (52, 2, 108)
     assert concentration_half.shape == (52, 2, 108)

From 9321f431d48e606b0f7fe673b4a06612be621cf2 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Fri, 19 Dec 2025 18:00:29 +0100
Subject: [PATCH 11/14] Update changelog

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

diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst
index 82a3f5d96..12618b401 100644
--- a/docs/api/changelog.rst
+++ b/docs/api/changelog.rst
@@ -15,6 +15,17 @@ Fixed
 - Fixed bug where :class:`imod.mf6.Evapotranspiration` package would write files
   to binary, which could not be parsed by MODFLOW 6 when ``proportion_depth``
   and ``proportion_rate`` were provided without segments.
+- Bug where :class:`imod.mf6.ConstantConcentration` package could not be written
+  for multiple timesteps.
+- Bug where :meth:`imod.mf6.Modflow6Simulation.clip_box` would drop layers if
+  ``states_for_boundary`` were provided and the model already contained a
+  :class:`imod.mf6.ConstantHead` or :class:`imod.mf6.ConstantConcentration` with
+  less layers.
+- Bug where :meth:`imod.mf6.Modflow6Simulation.clip_box` would not properly
+  align timesteps and forward fill data if ``states_for_boundary`` were provided
+  and the model already contained a :class:`imod.mf6.ConstantHead` or
+  :class:`imod.mf6.ConstantConcentration`, both with timesteps, which were
+  unaligned.
 
 Changed
 ~~~~~~~

From 77017582426da3334664b06ac58f13646bbb3af0 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Fri, 19 Dec 2025 18:04:07 +0100
Subject: [PATCH 12/14] Update changelog

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

diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst
index 12618b401..87f2eee6f 100644
--- a/docs/api/changelog.rst
+++ b/docs/api/changelog.rst
@@ -17,6 +17,10 @@ Fixed
   and ``proportion_rate`` were provided without segments.
 - Bug where :class:`imod.mf6.ConstantConcentration` package could not be written
   for multiple timesteps.
+- Bug where :meth:`imod.mf6.Modflow6Simulation.clip_box` where a ValidationError
+  was thrown when clipping a model with a :class:`imod.mf6.ConstantHead` or
+  :class:`imod.mf6.ConstantConcentration` package with a ``time`` dimension and
+  providing ``states_for_boundary``.
 - Bug where :meth:`imod.mf6.Modflow6Simulation.clip_box` would drop layers if
   ``states_for_boundary`` were provided and the model already contained a
   :class:`imod.mf6.ConstantHead` or :class:`imod.mf6.ConstantConcentration` with

From 6f4014d21a4f0817bbfc567583d1140c3f7f2070 Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Tue, 6 Jan 2026 14:21:04 +0100
Subject: [PATCH 13/14] Add test for rendering transient cnc and actually test
 that test_write_period_data writes the right amount of data

---
 imod/tests/test_mf6/test_mf6_cnc.py | 69 +++++++++++++++++++++++++----
 1 file changed, 61 insertions(+), 8 deletions(-)

diff --git a/imod/tests/test_mf6/test_mf6_cnc.py b/imod/tests/test_mf6/test_mf6_cnc.py
index 3caf6d048..31fbce973 100644
--- a/imod/tests/test_mf6/test_mf6_cnc.py
+++ b/imod/tests/test_mf6/test_mf6_cnc.py
@@ -74,16 +74,29 @@ def concentration_transient():
     idomain = xr.DataArray(np.ones(shape), coords=coords, dims=dims)
 
     # Constant cocnentration
-    concentration = xr.full_like(idomain, np.nan)
+    concentration = xr.full_like(idomain, np.nan).sel(layer=[1, 2])
     concentration[...] = np.nan
     concentration[..., 0] = 0.0
 
-    return concentration
+    time_da = xr.DataArray(
+        [1.0, 2.0],
+        coords={"time": globaltimes[:-1]},
+        dims=("time",),
+    )
 
+    return time_da * concentration
 
-def test_render(concentration_steadystate):
+
+def test_render_steady_state(concentration_steadystate):
     directory = pathlib.Path("mymodel")
-    globaltimes = np.array(["2000-01-01"], dtype="datetime64[ns]")
+    globaltimes = np.array(
+        [
+            "2000-01-01",
+            "2000-01-02",
+            "2000-01-03",
+        ],
+        dtype="datetime64[ns]",
+    )
 
     cnc = imod.mf6.ConstantConcentration(
         concentration_steadystate, print_input=True, print_flows=True, save_flows=True
@@ -104,7 +117,47 @@ def test_render(concentration_steadystate):
 
         begin period 1
           open/close mymodel/cnc/cnc.bin (binary)
-        end period"""
+        end period
+        """
+    )
+    assert actual == expected
+
+
+def test_render_transient(concentration_transient):
+    directory = pathlib.Path("mymodel")
+    globaltimes = np.array(
+        [
+            "2000-01-01",
+            "2000-01-02",
+            "2000-01-03",
+        ],
+        dtype="datetime64[ns]",
+    )
+
+    cnc = imod.mf6.ConstantConcentration(
+        concentration_transient, print_input=True, print_flows=True, save_flows=True
+    )
+    actual = cnc._render(directory, "cnc", globaltimes, True)
+
+    expected = textwrap.dedent(
+        """\
+        begin options
+          print_input
+          print_flows
+          save_flows
+        end options
+
+        begin dimensions
+          maxbound 30
+        end dimensions
+
+        begin period 1
+          open/close mymodel/cnc/cnc-0.bin (binary)
+        end period
+        begin period 2
+          open/close mymodel/cnc/cnc-1.bin (binary)
+        end period
+        """
     )
     assert actual == expected
 
@@ -118,7 +171,7 @@ def test_write_period_data(concentration_transient):
         ],
         dtype="datetime64[ns]",
     )
-    concentration_transient[:] = 2
+    concentration_transient += 999  # to avoid having zeros in the data
     cnc = imod.mf6.ConstantConcentration(
         concentration_transient,
         print_input=True,
@@ -131,5 +184,5 @@ def test_write_period_data(concentration_transient):
         with open(output_dir + "/cnc/cnc-0.dat", "r") as f:
             data = f.read()
             assert (
-                data.count("2") == 1080
-            )  # the number 2 is in the concentration data, and in the cell indices.
+                data.count("999") == 30
+            )  # the number 999 is in the concentration data, and in the cell indices.

From 4974aa90ec837661aaa94018b22f70fbe85b478c Mon Sep 17 00:00:00 2001
From: JoerivanEngelen <joerivanengelen@hotmail.com>
Date: Tue, 6 Jan 2026 15:46:34 +0100
Subject: [PATCH 14/14] deal with comment

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

diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst
index c851dfa00..c9df74458 100644
--- a/docs/api/changelog.rst
+++ b/docs/api/changelog.rst
@@ -24,17 +24,17 @@ Fixed
 - Fixed bug where :class:`imod.mf6.Evapotranspiration` package would write files
   to binary, which could not be parsed by MODFLOW 6 when ``proportion_depth``
   and ``proportion_rate`` were provided without segments.
-- Bug where :class:`imod.mf6.ConstantConcentration` package could not be written
+- Fixed bug where :class:`imod.mf6.ConstantConcentration` package could not be written
   for multiple timesteps.
-- Bug where :meth:`imod.mf6.Modflow6Simulation.clip_box` where a ValidationError
+- Fixed bug where :meth:`imod.mf6.Modflow6Simulation.clip_box` where a ValidationError
   was thrown when clipping a model with a :class:`imod.mf6.ConstantHead` or
   :class:`imod.mf6.ConstantConcentration` package with a ``time`` dimension and
   providing ``states_for_boundary``.
-- Bug where :meth:`imod.mf6.Modflow6Simulation.clip_box` would drop layers if
+- Fixed bug where :meth:`imod.mf6.Modflow6Simulation.clip_box` would drop layers if
   ``states_for_boundary`` were provided and the model already contained a
   :class:`imod.mf6.ConstantHead` or :class:`imod.mf6.ConstantConcentration` with
   less layers.
-- Bug where :meth:`imod.mf6.Modflow6Simulation.clip_box` would not properly
+- Fixed bug where :meth:`imod.mf6.Modflow6Simulation.clip_box` would not properly
   align timesteps and forward fill data if ``states_for_boundary`` were provided
   and the model already contained a :class:`imod.mf6.ConstantHead` or
   :class:`imod.mf6.ConstantConcentration`, both with timesteps, which were