diff --git a/news/1931.misc b/news/1931.misc
new file mode 100644
index 0000000000..624243aedb
--- /dev/null
+++ b/news/1931.misc
@@ -0,0 +1,2 @@
+
+Provide effective sediment erosion and deposition rates for 'SpaceLargeScaleEroder'
diff --git a/notebooks/tutorials/landscape_evolution/space/SPACE_large_scale_eroder_user_guide_and_examples.ipynb b/notebooks/tutorials/landscape_evolution/space/SPACE_large_scale_eroder_user_guide_and_examples.ipynb
index 6ac527e710..8ce59da30a 100644
--- a/notebooks/tutorials/landscape_evolution/space/SPACE_large_scale_eroder_user_guide_and_examples.ipynb
+++ b/notebooks/tutorials/landscape_evolution/space/SPACE_large_scale_eroder_user_guide_and_examples.ipynb
@@ -113,11 +113,10 @@
"source": [
"import matplotlib.pyplot as plt # For plotting results; optional\n",
"import numpy as np\n",
+ "from tqdm.notebook import trange\n",
"\n",
- "from landlab import RasterModelGrid # Grid utility\n",
- "from landlab import imshow_grid # For plotting results; optional\n",
- "from landlab.components import SpaceLargeScaleEroder # SpaceLargeScaleEroder model\n",
- "from landlab.components import PriorityFloodFlowRouter"
+ "from landlab import RasterModelGrid\n",
+ "from landlab.components import PriorityFloodFlowRouter, SpaceLargeScaleEroder"
]
},
{
@@ -126,7 +125,7 @@
"source": [
"Two Landlab components are essential to running the SPACE model: the model itself, and the `PriorityFloodFlowRouter`, which calculates drainage pathways, topographic slopes, and surface water discharge across the grid. \n",
"\n",
- "In addition to the relevant process components, some Landlab utilities are required to generate the model grid (in this example `RasterModelGrid`) and to visualize output (`imshow_grid`). Note that while it is possible to visualize output through functionality in other libraries (e.g., matplotlib), `imshow_grid` provides a simple way to generate 2-D maps of model variables.\n",
+ "In addition to the relevant process components, some Landlab utilities are required to generate the model grid (in this example `RasterModelGrid`).\n",
"\n",
"Most Landlab functionality requires the Numpy package for scientific computing in python. The matplotlib plotting library has also been imported to aid visualization of results."
]
@@ -154,44 +153,26 @@
},
"outputs": [],
"source": [
- "# Set grid parameters\n",
- "num_rows = 20\n",
- "num_columns = 20\n",
- "node_spacing = 100.0\n",
+ "rng = np.random.default_rng(seed=5000)\n",
"\n",
- "# track sediment flux at the node adjacent to the outlet at lower-left\n",
- "node_next_to_outlet = num_columns + 1\n",
- "\n",
- "# Instantiate model grid\n",
- "mg = RasterModelGrid((num_rows, num_columns), node_spacing)\n",
- "# add field ’topographic elevation’ to the grid\n",
- "mg.add_zeros(\"node\", \"topographic__elevation\")\n",
- "# set constant random seed for consistent topographic roughness\n",
- "np.random.seed(seed=5000)\n",
+ "mg = RasterModelGrid((20, 20), xy_spacing=100.0)\n",
"\n",
- "# Create initial model topography:\n",
"# plane tilted towards the lower−left corner\n",
- "topo = mg.node_y / 100000.0 + mg.node_x / 100000.0\n",
- "\n",
- "# add topographic roughness\n",
- "random_noise = (\n",
- " np.random.rand(len(mg.node_y)) / 1000.0\n",
- ") # impose topography values on model grid\n",
- "mg[\"node\"][\"topographic__elevation\"] += topo + random_noise\n",
- "\n",
- "# add field 'soil__depth' to the grid\n",
- "mg.add_zeros(\"node\", \"soil__depth\")\n",
+ "mg.at_node[\"bedrock__elevation\"] = (\n",
+ " mg.y_of_node / 100000.0\n",
+ " + mg.x_of_node / 100000.0\n",
+ " + rng.uniform(0.0, 1e-3, size=mg.number_of_nodes)\n",
+ ")\n",
"\n",
"# Set 2 m of initial soil depth at core nodes\n",
+ "mg.add_zeros(\"soil__depth\", at=\"node\")\n",
"mg.at_node[\"soil__depth\"][mg.core_nodes] = 2.0 # meters\n",
"\n",
- "# Add field 'bedrock__elevation' to the grid\n",
- "mg.add_zeros(\"bedrock__elevation\", at=\"node\")\n",
- "\n",
"# Sum 'soil__depth' and 'bedrock__elevation'\n",
"# to yield 'topographic elevation'\n",
- "mg.at_node[\"bedrock__elevation\"][:] = mg.at_node[\"topographic__elevation\"]\n",
- "mg.at_node[\"topographic__elevation\"][:] += mg.at_node[\"soil__depth\"]"
+ "mg.at_node[\"topographic__elevation\"] = (\n",
+ " mg.at_node[\"bedrock__elevation\"] + mg.at_node[\"soil__depth\"]\n",
+ ")"
]
},
{
@@ -218,10 +199,19 @@
"\n",
"# Set lower-left (southwest) corner as an open boundary\n",
"mg.set_watershed_boundary_condition_outlet_id(\n",
- " 0, mg[\"node\"][\"topographic__elevation\"], -9999.0\n",
+ " 0, mg.at_node[\"topographic__elevation\"], -9999.0\n",
")"
]
},
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "mg.imshow(\"topographic__elevation\", cmap=\"terrain\")"
+ ]
+ },
{
"cell_type": "markdown",
"metadata": {},
@@ -288,36 +278,22 @@
},
"outputs": [],
"source": [
- "# Set model timestep\n",
- "timestep = 1.0 # years\n",
+ "n_steps = 500\n",
+ "dt = 1.0 # years\n",
"\n",
- "# Set elapsed time to zero\n",
- "elapsed_time = 0.0 # years\n",
- "\n",
- "# Set timestep count to zero\n",
- "count = 0\n",
- "\n",
- "# Set model run time\n",
- "run_time = 500.0 # years\n",
- "\n",
- "# Array to save sediment flux values\n",
- "sed_flux = np.zeros(int(run_time // timestep))\n",
+ "# track sediment flux at the node adjacent to the outlet at lower-left\n",
+ "node_next_to_outlet = mg.shape[0] + 1\n",
+ "sed_flux = []\n",
"\n",
- "while elapsed_time < run_time: # time units of years\n",
+ "for time in trange(n_steps):\n",
" # Run the flow router\n",
" fr.run_one_step()\n",
"\n",
" # Run SPACE for one time step\n",
- " sp.run_one_step(dt=timestep)\n",
+ " sp.run_one_step(dt=dt)\n",
"\n",
" # Save sediment flux value to array\n",
- " sed_flux[count] = mg.at_node[\"sediment__flux\"][node_next_to_outlet]\n",
- "\n",
- " # Add to value of elapsed time\n",
- " elapsed_time += timestep\n",
- "\n",
- " # Increase timestep count\n",
- " count += 1"
+ " sed_flux.append(mg.at_node[\"sediment__flux\"][node_next_to_outlet])"
]
},
{
@@ -344,15 +320,7 @@
},
"outputs": [],
"source": [
- "# Instantiate figure\n",
- "fig = plt.figure()\n",
- "\n",
- "# Instantiate subplot\n",
- "plot = plt.subplot()\n",
- "\n",
- "# Show sediment flux map\n",
- "imshow_grid(\n",
- " mg,\n",
+ "mg.imshow(\n",
" \"sediment__flux\",\n",
" plot_name=\"Sediment flux\",\n",
" var_name=\"Sediment flux\",\n",
@@ -362,18 +330,20 @@
")"
]
},
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Elevation map"
+ ]
+ },
{
"cell_type": "code",
"execution_count": null,
- "metadata": {
- "jupyter": {
- "outputs_hidden": false
- }
- },
+ "metadata": {},
"outputs": [],
"source": [
- "# Export figure to image\n",
- "fig.savefig(\"sediment_flux_map.eps\")"
+ "mg.imshow(\"topographic__elevation\", cmap=\"terrain\")"
]
},
{
@@ -395,18 +365,10 @@
},
"outputs": [],
"source": [
- "# Instantiate figure\n",
- "fig = plt.figure()\n",
- "\n",
- "# Instantiate subplot\n",
- "sedfluxplot = plt.subplot()\n",
- "\n",
- "# Plot data\n",
- "sedfluxplot.plot(np.arange(500), sed_flux, color=\"k\", linewidth=3.0)\n",
+ "plt.plot(np.arange(len(sed_flux)) * dt, sed_flux, color=\"k\", linewidth=3.0)\n",
"\n",
- "# Add axis labels\n",
- "sedfluxplot.set_xlabel(\"Time [yr]\")\n",
- "sedfluxplot.set_ylabel(r\"Sediment flux [m$^3$/yr]\")"
+ "plt.gca().set_xlabel(\"Time [yr]\")\n",
+ "plt.gca().set_ylabel(r\"Sediment flux [m$^3$/yr]\")"
]
},
{
@@ -426,7 +388,7 @@
],
"metadata": {
"kernelspec": {
- "display_name": "Python 3",
+ "display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
@@ -440,7 +402,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
- "version": "3.9.1"
+ "version": "3.12.5"
}
},
"nbformat": 4,
diff --git a/src/landlab/components/space/ext/calc_sequential_ero_depo.pyx b/src/landlab/components/space/ext/calc_sequential_ero_depo.pyx
index f201ddeabd..caf0adc4ce 100644
--- a/src/landlab/components/space/ext/calc_sequential_ero_depo.pyx
+++ b/src/landlab/components/space/ext/calc_sequential_ero_depo.pyx
@@ -25,14 +25,15 @@ def _sequential_ero_depo(
const cython.floating [:] sed_erosion_term,
const cython.floating [:] bed_erosion_term,
const cython.floating [:] K_sed,
- double v,
- double phi,
- double F_f,
- double H_star,
- double dt,
- double thickness_lim,
+ cython.floating [:] ero_sed_effective,
+ cython.floating [:] depo_effective,
+ const double v,
+ const double phi,
+ const double F_f,
+ const double H_star,
+ const double dt,
+ const double thickness_lim,
):
-
"""Calculate and qs and qs_in."""
# define internal variables
cdef unsigned int node_id
@@ -107,4 +108,10 @@ def _sequential_ero_depo(
br[node_id] += -dt * ero_bed
vol_SSY_riv += F_f*ero_bed* cell_area[node_id]
+ # Update deposition rate based on adjusted fluxes
+ Hd = H_loc - H_Before
+ depo_effective[node_id] = (v*qs_out_adj/q[node_id])/(1 - phi)
+ # Deposition should be larger or equal to increase in soil depth
+ depo_effective[node_id] = max(depo_effective[node_id], Hd/dt)
+ ero_sed_effective[node_id] = depo_effective[node_id] - Hd/dt
return vol_SSY_riv
diff --git a/src/landlab/components/space/space_large_scale_eroder.py b/src/landlab/components/space/space_large_scale_eroder.py
index 07f9a75242..618bbd758d 100644
--- a/src/landlab/components/space/space_large_scale_eroder.py
+++ b/src/landlab/components/space/space_large_scale_eroder.py
@@ -7,68 +7,68 @@
from landlab import Component
from landlab import RasterModelGrid
+from landlab.components.depression_finder.lake_mapper import _FLOODED
+from landlab.components.space.ext.calc_sequential_ero_depo import _sequential_ero_depo
from landlab.grid.nodestatus import NodeStatus
from landlab.utils.return_array import return_array_at_node
-from ..depression_finder.lake_mapper import _FLOODED
-from .ext.calc_sequential_ero_depo import _sequential_ero_depo
-
-ROOT2 = np.sqrt(2.0) # syntactic sugar for precalculated square root of 2
+ROOT2 = np.sqrt(2.0)
TIME_STEP_FACTOR = 0.5 # factor used in simple subdivision solver
class SpaceLargeScaleEroder(Component):
- """Stream Power with Alluvium Conservation and Entrainment (SPACE) large scale eroder
+ """Stream Power with Alluvium Conservation and Entrainment large scale eroder.
- The SPACE_large_Scale_eroder is based on the SPACE component and is designed
- to be more robust against large time steps and coded in such a way that mass
- conservation is explicitly conserved during calculation.
+ The :class:`~.SpaceLargeScaleEroder` is based on the SPACE component
+ and is designed to be more robust against large time steps and coded in
+ such a way that mass conservation is explicitly conserved during calculation.
See the publication:
- Shobe, C. M., Tucker, G. E., and Barnhart, K. R.: The SPACE 1.0 model: a
- Landlab component for 2-D calculation of sediment transport, bedrock
- erosion, and landscape evolution, Geosci. Model Dev., 10, 4577-4604,
- `https://doi.org/10.5194/gmd-10-4577-2017 `_, 2017.
+
+ Shobe, C. M., Tucker, G. E., and Barnhart, K. R.: The SPACE 1.0 model: a
+ Landlab component for 2-D calculation of sediment transport, bedrock
+ erosion, and landscape evolution, Geosci. Model Dev., 10, 4577-4604,
+ `https://doi.org/10.5194/gmd-10-4577-2017
+ `_, 2017.
Unlike other some other fluvial erosion componets in Landlab, in this
- component (and :py:class:`~landlab.components.ErosionDeposition`) no
+ component (and class:`~landlab.components.ErosionDeposition`) no
erosion occurs in depressions or in areas with adverse slopes. There is no
ability to pass a keyword argument ``erode_flooded_nodes``.
- If a depressions are handled (as indicated by the presence of the field
- "flood_status_code" at nodes), then deposition occurs throughout the
+ If depressions are handled (as indicated by the presence of the field
+ ``"flood_status_code"`` at nodes), then deposition occurs throughout the
depression and sediment is passed out of the depression. Where pits are
encountered, then all sediment is deposited at that node only.
- Note: In the current version, we do not provide an adaptive time stepper.
- This will be addded in future versions of this component.
+ .. note::
+
+ In the current version, we do not provide an adaptive time stepper.
+ This will be addded in future versions of this component.
- For more explanation and examples,
- check out the correponding notebook of this component
+ For more explanation and examples, check out the correponding notebook of
+ this component.
Examples
---------
- >>> import numpy as np
>>> from landlab import RasterModelGrid
- >>> from landlab.components import PriorityFloodFlowRouter, SpaceLargeScaleEroder
- >>> import matplotlib.pyplot as plt # For plotting results; optional
- >>> from landlab import imshow_grid # For plotting results; optional
-
- >>> num_rows = 20
- >>> num_columns = 20
- >>> node_spacing = 100.0
- >>> mg = RasterModelGrid((num_rows, num_columns), xy_spacing=node_spacing)
- >>> node_next_to_outlet = num_columns + 1
- >>> np.random.seed(seed=5000)
- >>> _ = mg.add_zeros("topographic__elevation", at="node")
- >>> _ = mg.add_zeros("soil__depth", at="node")
- >>> mg.at_node["soil__depth"][mg.core_nodes] = 2.0
- >>> _ = mg.add_zeros("bedrock__elevation", at="node")
- >>> mg.at_node["bedrock__elevation"] += (
- ... mg.node_y / 10.0 + mg.node_x / 10.0 + np.random.rand(len(mg.node_y)) / 10.0
+ >>> from landlab.components import PriorityFloodFlowRouter
+ >>> from landlab.components import SpaceLargeScaleEroder
+
+ >>> mg = RasterModelGrid((5, 4), xy_spacing=100.0)
+
+ >>> mg.at_node["soil__depth"] = [
+ ... [0.0, 0.0, 0.0, 0.0],
+ ... [0.0, 2.0, 2.0, 0.0],
+ ... [0.0, 2.0, 2.0, 0.0],
+ ... [0.0, 2.0, 2.0, 0.0],
+ ... [0.0, 0.0, 0.0, 0.0],
+ ... ]
+ >>> mg.at_node["bedrock__elevation"] = mg.y_of_node / 10.0 + mg.x_of_node / 10.0
+ >>> mg.at_node["topographic__elevation"] = (
+ ... mg.at_node["bedrock__elevation"] + mg.at_node["soil__depth"]
... )
- >>> mg.at_node["bedrock__elevation"][:] = mg.at_node["topographic__elevation"]
- >>> mg.at_node["topographic__elevation"][:] += mg.at_node["soil__depth"]
+
>>> mg.set_closed_boundaries_at_grid_edges(
... bottom_is_closed=True,
... left_is_closed=True,
@@ -93,49 +93,29 @@ class SpaceLargeScaleEroder(Component):
... sp_crit_sed=0,
... sp_crit_br=0,
... )
- >>> timestep = 10.0
- >>> elapsed_time = 0.0
- >>> count = 0
- >>> run_time = 1e4
- >>> sed_flux = np.zeros(int(run_time // timestep))
- >>> while elapsed_time < run_time:
+
+ >>> node_next_to_outlet = mg.shape[1] + 1
+ >>> sed_flux = []
+ >>> for _ in range(500):
... fr.run_one_step()
- ... _ = sp.run_one_step(dt=timestep)
- ... sed_flux[count] = mg.at_node["sediment__flux"][node_next_to_outlet]
- ... elapsed_time += timestep
- ... count += 1
+ ... _ = sp.run_one_step(dt=10.0)
+ ... sed_flux.append(mg.at_node["sediment__flux"][node_next_to_outlet])
...
- Plot the results.
-
- >>> fig = plt.figure()
- >>> plot = plt.subplot()
- >>> _ = imshow_grid(
- ... mg,
- ... "topographic__elevation",
- ... plot_name="Sediment flux",
- ... var_name="Sediment flux",
- ... var_units=r"m$^3$/yr",
- ... grid_units=("m", "m"),
- ... cmap="terrain",
- ... )
- >>> _ = plt.figure()
- >>> _ = imshow_grid(
- ... mg,
- ... "sediment__flux",
- ... plot_name="Sediment flux",
- ... var_name="Sediment flux",
- ... var_units=r"m$^3$/yr",
- ... grid_units=("m", "m"),
- ... cmap="terrain",
- ... )
- >>> fig = plt.figure()
- >>> sedfluxplot = plt.subplot()
- >>> _ = sedfluxplot.plot(
- ... np.arange(len(sed_flux)) * timestep, sed_flux, color="k", linewidth=1.0
- ... )
- >>> _ = sedfluxplot.set_xlabel("Time [yr]")
- >>> _ = sedfluxplot.set_ylabel(r"Sediment flux [m$^3$/yr]")
+ Look at the results.
+
+ >>> mg.at_node["sediment__flux"].reshape(mg.shape) # doctest: +SKIP
+ array([[0. , 0. , 0. , 0. ],
+ [0. , 0.26889843, 0.02719885, 0. ],
+ [0. , 0.09130624, 0.13962599, 0. ],
+ [0. , 0.06421368, 0.09527108, 0. ],
+ [0. , 0. , 0. , 0. ]])
+ >>> sed_flux[::100] # doctest: +SKIP
+ [2053.1526698811363,
+ 601.0591808186842,
+ 405.95978528106235,
+ 272.8232194793999,
+ 176.63159183013272]
References
----------
@@ -150,7 +130,7 @@ class SpaceLargeScaleEroder(Component):
None Listed
- """ # noqa: B950
+ """
_name = "SpaceLargeScaleEroder"
@@ -197,6 +177,39 @@ class SpaceLargeScaleEroder(Component):
"mapping": "node",
"doc": "Sediment flux (volume per unit time of sediment leaving each node)",
},
+ "sediment__erosion_flux": {
+ "dtype": float,
+ "intent": "out",
+ "optional": False,
+ "units": "m/s",
+ "mapping": "node",
+ "doc": (
+ "Sediment erosion flux from bed to water column (depth eroded per"
+ " unit time)"
+ ),
+ },
+ "sediment__deposition_flux": {
+ "dtype": float,
+ "intent": "out",
+ "optional": False,
+ "units": "m/s",
+ "mapping": "node",
+ "doc": (
+ "Sediment deposition flux from water column to bed (depth deposited"
+ " per unit time)"
+ ),
+ },
+ "bedrock__erosion_flux": {
+ "dtype": float,
+ "intent": "out",
+ "optional": False,
+ "units": "m/s",
+ "mapping": "node",
+ "doc": (
+ "Bedrock erosion flux from bedrock to water column (depth eroded per"
+ " unit time)"
+ ),
+ },
"soil__depth": {
"dtype": float,
"intent": "inout",
@@ -547,6 +560,9 @@ def run_one_step_basic(self, dt=10):
K_sed_vector = np.broadcast_to(self._K_sed, self._q.shape)
+ ero_sed_effective = np.zeros_like(K_sed_vector)
+ depo_effective = np.zeros_like(K_sed_vector)
+
vol_SSY_riv = _sequential_ero_depo(
stack_flip_ud_sel,
r,
@@ -563,6 +579,8 @@ def run_one_step_basic(self, dt=10):
self._sed_erosion_term,
self._br_erosion_term,
K_sed_vector,
+ ero_sed_effective,
+ depo_effective,
self._v_s,
self._phi,
self._F_f,
@@ -576,8 +594,18 @@ def run_one_step_basic(self, dt=10):
cores = self._grid.core_nodes
z[cores] = br[cores] + H[cores]
+ self.grid.at_node["sediment__erosion_flux"][:] = ero_sed_effective
+ self.grid.at_node["sediment__deposition_flux"][:] = depo_effective
+ self.grid.at_node["bedrock__erosion_flux"][:] = self._Er
+
return vol_SSY_riv, V_leaving_riv
def run_one_step(self, dt):
- vol_SSY_riv, V_leaving_riv = self.run_one_step_basic(dt)
+ """
+ Returns:
+ - vol_SSY_riv (float): Suspended sediment yield leaving the domain as wash load
+ - V_leaving_riv (float): Volume of bedload sediment leaving the domain.
+ """
+
+ (vol_SSY_riv, V_leaving_riv) = self.run_one_step_basic(dt)
return vol_SSY_riv, V_leaving_riv
diff --git a/tests/components/space/test_space_large_scale_eroder.py b/tests/components/space/test_space_large_scale_eroder.py
index a00941185a..6d7f889721 100644
--- a/tests/components/space/test_space_large_scale_eroder.py
+++ b/tests/components/space/test_space_large_scale_eroder.py
@@ -1317,6 +1317,121 @@ def test_MassBalance():
cores = mg.core_nodes
area = mg.cell_area_at_node
+ # ... and run it to steady state (10000x1-year timesteps).
+ for _ in range(10000):
+ fa.run_one_step()
+ soil_B = cp.deepcopy(H)
+ bed_B = cp.deepcopy(br)
+ vol_SSY_riv, V_leaving_riv = sp.run_one_step(dt=dt)
+ diff_MB = (
+ np.sum((bed_B[cores] - br[cores]) * area[cores])
+ + np.sum((soil_B[cores] - H[cores]) * area[cores]) * (1 - sp._phi)
+ - vol_SSY_riv * dt
+ - V_leaving_riv
+ )
+
+ br[mg.core_nodes] += U * dt # m
+ soil[0] = 0.0 # enforce 0 soil depth at boundary to keep lowering steady
+ z[:] = br[:] + soil[:]
+
+ # Test Every iteration
+ testing.assert_array_almost_equal(
+ z[cores],
+ br[cores] + H[cores],
+ decimal=5,
+ err_msg="Topography does not equal sum of bedrock and soil! Decrease timestep",
+ verbose=True,
+ )
+ testing.assert_array_less(
+ abs(diff_MB),
+ 1e-8 * mg.number_of_nodes,
+ err_msg=(
+ "Mass balance error SpaceLargeScaleEroder! Try to resolve by "
+ "decreasing timestep"
+ ),
+ verbose=True,
+ )
+
+ # Check mass balance on effective erosion and deposition values
+ soil_new_calc = (
+ soil_B
+ + (
+ mg.at_node["sediment__deposition_flux"]
+ - mg.at_node["sediment__erosion_flux"]
+ )
+ * dt
+ )
+ np.testing.assert_almost_equal(
+ H[mg.core_nodes], soil_new_calc[mg.core_nodes], decimal=8
+ )
+
+
+# %%
+@pytest.mark.slow
+def test_MassBalance_lower_pore_density():
+ # %%
+ # set up a 15x15 grid with one open outlet node and low initial elevations.
+ nr = 15
+ nc = 15
+ 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")
+
+ mg["node"]["topographic__elevation"] += (
+ mg.node_y / 100000 + mg.node_x / 100000 + np.random.rand(len(mg.node_y)) / 10000
+ )
+ mg.set_closed_boundaries_at_grid_edges(
+ bottom_is_closed=True,
+ left_is_closed=True,
+ right_is_closed=True,
+ top_is_closed=True,
+ )
+ 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 = FlowAccumulator(
+ mg, flow_director="D8", depression_finder="DepressionFinderAndRouter"
+ )
+
+ # Parameter values for detachment-limited test
+ K_br = 0.002
+ K_sed = 0.002
+ U = 0.0001
+ dt = 10.0
+ F_f = 0.2 # all detached rock disappears; detachment-ltd end-member
+ m_sp = 0.5
+ n_sp = 1.0
+ v_s = 0.25
+ H_star = 0.1
+
+ # Instantiate the Space component...
+ sp = SpaceLargeScaleEroder(
+ mg,
+ K_sed=K_sed,
+ K_br=K_br,
+ F_f=F_f,
+ phi=0.4,
+ H_star=H_star,
+ v_s=v_s,
+ m_sp=m_sp,
+ n_sp=n_sp,
+ sp_crit_sed=0,
+ sp_crit_br=0,
+ )
+ # Get values before run
+ z = mg.at_node["topographic__elevation"]
+ br = mg.at_node["bedrock__elevation"]
+ H = mg.at_node["soil__depth"]
+ cores = mg.core_nodes
+ area = mg.cell_area_at_node
+
# ... and run it to steady state (10000x1-year timesteps).
for _ in range(10000):
fa.run_one_step()
@@ -1351,3 +1466,16 @@ def test_MassBalance():
),
verbose=True,
)
+
+ # Check mass balance on effective erosion and deposition values
+ soil_new_calc = (
+ soil_B
+ + (
+ mg.at_node["sediment__deposition_flux"]
+ - mg.at_node["sediment__erosion_flux"]
+ )
+ * dt
+ )
+ np.testing.assert_almost_equal(
+ H[mg.core_nodes], soil_new_calc[mg.core_nodes], decimal=8
+ )