Bugfix for issue #1901

landlab/components/space/space_large_scale_eroder.py

@@ -496,12 +496,6 @@ def _calc_erosion_rates(self):
             1.0 - np.exp(-omega_br_over_sp_crit)
         )
 
-        # Do not allow for the formation of potholes (addition v2)
-        r = self._grid.at_node["flow__receiver_node"]
-        br_e_max = br - br[r]
-        br_e_max[br_e_max < 0] = 0
-        self._br_erosion_term = np.minimum(self._br_erosion_term, br_e_max)
-
         self._Es = self._sed_erosion_term * (1.0 - np.exp(-H / self._H_star))
         self._Er = self._br_erosion_term * np.exp(-H / self._H_star)
 

tests/components/space/test_space_large_scale_eroder.py

@@ -1018,6 +1018,116 @@ def test_matches_bedrock_alluvial_solution_PF():
     # %%
 
 
+# %%
+@pytest.mark.slow
+@pytest.mark.skipif(not with_richdem, reason="richdem is not installed")
+def test_matches_bedrock_alluvial_solution_PF_extended_range():
+    """
+    Test that model matches the bedrock-alluvial analytical solution
+    for slope/area relationship at steady state:
+    S=((U * v_s * (1 - F_f)) / (K_sed * A^m) + U / (K_br * A^m))^(1/n).
+
+    Also test that the soil depth everywhere matches the bedrock-alluvial
+    analytical solution at steady state:
+    H = -H_star * ln(1 - (v_s / (K_sed / (K_br * (1 - F_f)) + v_s))).
+    """
+    # %%
+    # set up a 5x5 grid with one open outlet node and low initial elevations.
+    nr = 5
+    nc = 5
+    mg = RasterModelGrid((nr, nc), xy_spacing=10.0)
+
+    z = mg.add_zeros("topographic__elevation", at="node")
+    br = mg.add_zeros("bedrock__elevation", at="node")
+    soil = mg.add_zeros("soil__depth", at="node")
+
+    np.random.seed(seed=5000)
+    mg["node"]["topographic__elevation"] += (
+        mg.node_y / 100000 + mg.node_x / 100000 + np.random.rand(len(mg.node_y)) / 10000
+    )
+
+    mg.set_watershed_boundary_condition_outlet_id(
+        0, mg["node"]["topographic__elevation"], -9999.0
+    )
+    soil[:] += 0.0  # initial condition of no soil depth.
+    br[:] = z[:]
+    z[:] += soil[:]
+
+    # Create a D8 flow handler
+    fa = PriorityFloodFlowRouter(
+        mg, surface="topographic__elevation", flow_metric="D8", suppress_out=True
+    )
+    fa.run_one_step()
+
+    # Parameter values for detachment-limited test
+    K_br = 0.005
+    K_sed = 0.01
+    U = 0.0001
+    dt = 10.0
+    F_f = 0.  
+    m_sp = 0.5
+    n_sp = 1.0
+    v_s = 1.0
+    H_star = 1.0
+
+    # Instantiate the SpaceLargeScaleEroder component...
+    sp = SpaceLargeScaleEroder(
+        mg,
+        K_sed=K_sed,
+        K_br=K_br,
+        F_f=F_f,
+        phi=0.0,
+        H_star=H_star,
+        v_s=v_s,
+        m_sp=m_sp,
+        n_sp=n_sp,
+        sp_crit_sed=0,
+        sp_crit_br=0,
+    )
+
+    # ... and run it to steady state (10000x1-year timesteps).
+    for _ in range(10000):
+        fa.run_one_step()
+        sp.run_one_step(dt=dt)
+        br[mg.core_nodes] += U * dt  # m
+        soil[0] = 0.0  # enforce 0 soil depth at boundary to keep lowering steady
+        z[:] = br[:] + soil[:]
+
+    # compare numerical and analytical slope solutions
+    num_slope = mg.at_node["topographic__steepest_slope"][mg.core_nodes]
+    analytical_slope = np.power(
+        (
+            (U * v_s * (1 - F_f))
+            / (K_sed * np.power(mg.at_node["drainage_area"][mg.core_nodes], m_sp))
+        )
+        + (U / (K_br * np.power(mg.at_node["drainage_area"][mg.core_nodes], m_sp))),
+        1.0 / n_sp,
+    )
+
+    # test for match with analytical slope-area relationship
+    testing.assert_array_almost_equal(
+        num_slope,
+        analytical_slope,
+        decimal=8,
+        err_msg="SpaceLargeScaleEroder bedrock-alluvial slope-area test failed",
+        verbose=True,
+    )
+
+    # compare numerical and analytical soil depth solutions
+    num_h = mg.at_node["soil__depth"][mg.core_nodes]
+    analytical_h = -H_star * np.log(1 - (v_s / (K_sed / (K_br * (1 - F_f)) + v_s)))
+
+    # test for match with analytical sediment depth
+    testing.assert_array_almost_equal(
+        num_h,
+        analytical_h,
+        decimal=5,
+        err_msg="SpaceLargeScaleEroder bedrock-alluvial soil thickness test failed",
+        verbose=True,
+    )
+
+
+#%%
 @pytest.mark.slow
 def test_MassBalance():
     # %%