From 99634da96117d7f5ab7b0ee53ac35b24db09cb34 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Tue, 17 Mar 2026 11:20:56 +0100 Subject: [PATCH 1/7] Add test for edge case where stage equals bottom_elevation equals model bottom --- imod/tests/test_prepare/test_topsystem.py | 46 ++++++++++++++++++++++- 1 file changed, 44 insertions(+), 2 deletions(-) diff --git a/imod/tests/test_prepare/test_topsystem.py b/imod/tests/test_prepare/test_topsystem.py index a16a7b016..a19cf1125 100644 --- a/imod/tests/test_prepare/test_topsystem.py +++ b/imod/tests/test_prepare/test_topsystem.py @@ -246,8 +246,50 @@ def test_riv_allocation__elevation_above_surface_level( def test_riv_allocation__stage_equals_bottom_elevation( active, top, bottom, stage, bottom_elevation, option, expected_riv, expected_drn ): - # Put elevations a lot above surface level. Need to be allocated to first - # layer. + # Bottom elevation equals stage here. + actual_riv_da, actual_drn_da = allocate_riv_cells( + option, active, top, bottom, stage, stage + ) + + # Override expected values + if option is ALLOCATION_OPTION.first_active_to_elevation: + expected_riv = [True, True, False, False] + elif option is not ALLOCATION_OPTION.at_first_active: + expected_riv = [False, True, False, False] + if expected_drn: + expected_drn = [True, False, False, False] + + actual_riv = take_nth_layer_column(actual_riv_da, 0) + empty_riv = take_nth_layer_column(actual_riv_da, 1) + + if actual_drn_da is None: + actual_drn = None + empty_drn = None + else: + actual_drn = take_nth_layer_column(actual_drn_da, 0) + empty_drn = take_nth_layer_column(actual_drn_da, 1) + + np.testing.assert_equal(actual_riv, expected_riv) + np.testing.assert_equal(actual_drn, expected_drn) + assert np.all(~empty_riv) + if empty_drn is not None: + assert np.all(~empty_drn) + + +@parametrize_with_cases( + argnames="active,top,bottom,stage,bottom_elevation", + prefix="riv_", +) +@parametrize_with_cases( + argnames="option,expected_riv,expected_drn", prefix="allocation_", has_tag="riv" +) +def test_riv_allocation__stage_equals_bottom_elevation_equals_bottom( + active, top, bottom, stage, bottom_elevation, option, expected_riv, expected_drn +): + # Set bottom in layer 2 to stage, take first value (stage is equal everywhere.) + bottom.loc[bottom.coords["layer"] == 2] = stage.values.ravel()[0] + + # Bottom elevation equals stage here. actual_riv_da, actual_drn_da = allocate_riv_cells( option, active, top, bottom, stage, stage ) From 5b6999ea61d8118f24e0efef4e329231acccff6a Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Tue, 17 Mar 2026 11:29:09 +0100 Subject: [PATCH 2/7] Deal with edge case where river stage equals model bottom. --- imod/prepare/topsystem/allocation.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/imod/prepare/topsystem/allocation.py b/imod/prepare/topsystem/allocation.py index 116a580fb..b74bec2f9 100644 --- a/imod/prepare/topsystem/allocation.py +++ b/imod/prepare/topsystem/allocation.py @@ -321,7 +321,7 @@ def get_above_lower_bound(bottom_elevation: GridDataArray, top_layered: GridData bottom_elevation grid or in the first layer. """ top_layer_label = {"layer": min(top_layered.coords["layer"])} - is_above_lower_bound = bottom_elevation <= top_layered + is_above_lower_bound = bottom_elevation < top_layered # Bottom elevation above top surface is allowed, so these are set to True # regardless. is_above_lower_bound.loc[top_layer_label] = ~bottom_elevation.isnull() @@ -366,7 +366,7 @@ def _allocate_cells__stage_to_riv_bot( top_layered = _enforce_layered_top(top, bottom) is_above_lower_bound = get_above_lower_bound(bottom_elevation, top_layered) - is_below_upper_bound = stage > bottom + is_below_upper_bound = stage >= bottom riv_cells = is_below_upper_bound & is_above_lower_bound return riv_cells, None @@ -465,7 +465,7 @@ def _allocate_cells__stage_to_riv_bot_drn_above( layer = active.coords["layer"] is_below_upper_bound = upper_active_layer <= layer is_below_upper_bound_and_active = is_below_upper_bound & active - drn_cells = is_below_upper_bound_and_active & (bottom >= stage) + drn_cells = is_below_upper_bound_and_active & (bottom > stage) riv_cells = (is_below_upper_bound_and_active & is_above_lower_bound) != drn_cells return riv_cells, drn_cells From 53e3e10be5fd9b73445208bbde16029bf643d822 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Tue, 17 Mar 2026 12:44:27 +0100 Subject: [PATCH 3/7] Add test for distributing conductances to model bottom --- imod/tests/test_prepare/test_topsystem.py | 37 +++++++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/imod/tests/test_prepare/test_topsystem.py b/imod/tests/test_prepare/test_topsystem.py index a19cf1125..c7bf56be6 100644 --- a/imod/tests/test_prepare/test_topsystem.py +++ b/imod/tests/test_prepare/test_topsystem.py @@ -541,6 +541,43 @@ def test_distribute_riv_conductance__stage_equal_to_bottom_elevation( np.testing.assert_equal(actual, expected) +@parametrize_with_cases( + argnames="active,top,bottom,stage,bottom_elevation", + prefix="riv_", +) +@parametrize_with_cases( + argnames="option,allocated_layer,_", prefix="distribution_", has_tag="riv" +) +def test_distribute_riv_conductance__stage_equal_to_bottom_elevation_equal_to_bottom( + active, top, bottom, stage, bottom_elevation, option, allocated_layer, _ +): + # Set bottom in layer 2 to stage, take first value (stage is equal everywhere.) + bottom.loc[bottom.coords["layer"] == 2] = stage.values.ravel()[0] + + allocated_layer.data = np.array([False, True, False, False]) + expected = [np.nan, 1.0, np.nan, np.nan] + allocated = enforce_dim_order(active & allocated_layer) + k = xr.DataArray( + [2.0, 2.0, 1.0, 1.0], coords={"layer": [1, 2, 3, 4]}, dims=("layer",) + ) + + conductance = zeros_like(bottom_elevation) + 1.0 + + actual_da = distribute_riv_conductance( + option, + allocated, + conductance, + top, + bottom, + k, + stage, + stage, + ) + actual = take_nth_layer_column(actual_da, 0) + + np.testing.assert_equal(actual, expected) + + @parametrize_with_cases( argnames="active,top,bottom,elevation", prefix="ghb_", From 04be21f5fca9ed4009421d58d2c83b8f38bb5d10 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Tue, 17 Mar 2026 15:19:52 +0100 Subject: [PATCH 4/7] Add test where drain equals model bottom --- imod/tests/test_prepare/test_topsystem.py | 24 +++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/imod/tests/test_prepare/test_topsystem.py b/imod/tests/test_prepare/test_topsystem.py index c7bf56be6..6101158cf 100644 --- a/imod/tests/test_prepare/test_topsystem.py +++ b/imod/tests/test_prepare/test_topsystem.py @@ -351,6 +351,30 @@ def test_drn_allocation__elevation_above_surface_level( assert np.all(~empty) +@parametrize_with_cases( + argnames="active,top,bottom,drn_elevation", + prefix="drn_", +) +@parametrize_with_cases( + argnames="option,expected,_", prefix="allocation_", has_tag="drn" +) +def test_drn_allocation__elevation_equal_to_bottom( + active, top, bottom, drn_elevation, option, expected, _ + ): + # Set bottom in layer 3 to drain elevation, take first value (drain + # elevation is equal everywhere.) + bottom.loc[bottom.coords["layer"] == 3] = drn_elevation.values.ravel()[0] + + actual_da = allocate_drn_cells(option, active, top, bottom, drn_elevation) + + actual = take_nth_layer_column(actual_da, 0) + empty = take_nth_layer_column(actual_da, 1) + + np.testing.assert_equal(actual, expected) + assert np.all(~empty) + + + @parametrize_with_cases( argnames="active,top,bottom,head", prefix="ghb_", From fc618fc7927b9483c92970c5f1c7961d38735c87 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Tue, 17 Mar 2026 15:36:20 +0100 Subject: [PATCH 5/7] Update changelog --- docs/api/changelog.rst | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst index 2fc2228bf..717f15a92 100644 --- a/docs/api/changelog.rst +++ b/docs/api/changelog.rst @@ -45,6 +45,12 @@ Fixed produced incorrect results when point water heads were below elevation levels for unstructured grids. - Support pandas 3.0. +- Fixed edge case where allocation of :class:`imod.mf6.River` package with the + ``stage_to_riv_bot`` or ``stage_to_riv_bot_drn_above`` option of + :func:`imod.prepare.ALLOCATION_OPTION` would assign river cells to the wrong + layer when was exactly equal to a bottom of a layer in the model + discretization, which would these cells to be dropped when distributing + conductances later. Changed From b502b93b94d7a4cf22f159e1e55d1b3c43d6cb20 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Tue, 17 Mar 2026 15:38:09 +0100 Subject: [PATCH 6/7] format --- imod/tests/test_prepare/test_topsystem.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/imod/tests/test_prepare/test_topsystem.py b/imod/tests/test_prepare/test_topsystem.py index 6101158cf..65295c7ff 100644 --- a/imod/tests/test_prepare/test_topsystem.py +++ b/imod/tests/test_prepare/test_topsystem.py @@ -359,8 +359,8 @@ def test_drn_allocation__elevation_above_surface_level( argnames="option,expected,_", prefix="allocation_", has_tag="drn" ) def test_drn_allocation__elevation_equal_to_bottom( - active, top, bottom, drn_elevation, option, expected, _ - ): + active, top, bottom, drn_elevation, option, expected, _ +): # Set bottom in layer 3 to drain elevation, take first value (drain # elevation is equal everywhere.) bottom.loc[bottom.coords["layer"] == 3] = drn_elevation.values.ravel()[0] @@ -374,7 +374,6 @@ def test_drn_allocation__elevation_equal_to_bottom( assert np.all(~empty) - @parametrize_with_cases( argnames="active,top,bottom,head", prefix="ghb_", From 135cbac2ae507e377f79ab70672e1395ad82b56d Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Thu, 19 Mar 2026 08:23:08 +0100 Subject: [PATCH 7/7] Fix typos in changelog --- docs/api/changelog.rst | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst index 5538fb90b..703a719b7 100644 --- a/docs/api/changelog.rst +++ b/docs/api/changelog.rst @@ -47,13 +47,14 @@ Fixed for unstructured grids. - Support pandas 3.0. - :class:`imod.msw.IdfMapping` when model clipping is applied, the global - row/column indices are converted to local indices, as writen in idf_svat.inp. + row/column indices are converted to local indices, as written in + ``idf_svat.inp``. - Fixed edge case where allocation of :class:`imod.mf6.River` package with the ``stage_to_riv_bot`` or ``stage_to_riv_bot_drn_above`` option of :func:`imod.prepare.ALLOCATION_OPTION` would assign river cells to the wrong - layer when was exactly equal to a bottom of a layer in the model - discretization, which would cause these cells to be dropped when distributing - conductances later. + layer, when the stage and bottom_elevation were exactly equal to the bottom of + a layer in the model discretization, which would cause these cells to be + dropped when distributing conductances later. Changed ~~~~~~~