Skip to content

Commit

Permalink
Merge pull request landlab#1931 from BCampforts/bc/space_return_eff_e…
Browse files Browse the repository at this point in the history 
…ro_depo

Provide effective sediment erosion and deposition rates from SpaceLargeScaleEroder
  • Loading branch information
@mcflugen
mcflugen authored Sep 18, 2024
2 parents 60b1f1a + a52ee1b commit fcbfe6e
Show file tree
Hide file tree
Showing 5 changed files with 300 additions and 173 deletions.
2 changes: 2 additions & 0 deletions news/1931.misc
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@

Provide effective sediment erosion and deposition rates for 'SpaceLargeScaleEroder'
134 changes: 48 additions & 86 deletions ...utorials/landscape_evolution/space/SPACE_large_scale_eroder_user_guide_and_examples.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
]
},
{
Expand All @@ -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."
]
Expand Down Expand Up @@ -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",
")"
]
},
{
Expand All @@ -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": {},
Expand Down Expand Up @@ -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])"
]
},
{
Expand All @@ -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",
Expand All @@ -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\")"
]
},
{
Expand All @@ -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]\")"
]
},
{
Expand All @@ -426,7 +388,7 @@
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
Expand All @@ -440,7 +402,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.1"
"version": "3.12.5"
}
},
"nbformat": 4,
Expand Down
21 changes: 14 additions & 7 deletions src/landlab/components/space/ext/calc_sequential_ero_depo.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Loading

0 comments on commit fcbfe6e

Please sign in to comment.