diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst
index 629ff6f17..703a719b7 100644
--- a/docs/api/changelog.rst
+++ b/docs/api/changelog.rst
@@ -47,8 +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 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
 ~~~~~~~
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
diff --git a/imod/tests/test_prepare/test_topsystem.py b/imod/tests/test_prepare/test_topsystem.py
index a16a7b016..65295c7ff 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
     )
@@ -309,6 +351,29 @@ 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_",
@@ -499,6 +564,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_",