Escarpment [Py]

A simple example simulating the evolution of an escarpment (on a 2D raster grid) under the action of bedrock channel and hillslope erosion.

This example is pretty similar to Mountain [Py], although solving equation (1) without the uplift term \(U\) on a semi-infinite (periodic) domain. It also starts from a different initial topograhic surface and routes flow using a multiple direction approach.

import fastscapelib as fs
import numpy as np
import matplotlib.pyplot as plt

Setup the Grid, Flow Graph and Eroders

Create a RasterGrid of 101x201 nodes with a total length of 10 km in y (rows) and 20 km in x (columns).

Set fixed-value boundary at the left border, free (core) boundary at the right border and reflective boundaries for the top-down borders.

bs = [
    fs.NodeStatus.FIXED_VALUE,
    fs.NodeStatus.CORE,
    fs.NodeStatus.LOOPED,
    fs.NodeStatus.LOOPED,
]

grid = fs.RasterGrid.from_length([101, 201], [1e4, 2e4], bs)

Create a FlowGraph object with multiple direction flow routing and the resolution of closed depressions on the topographic surface. See Flow Routing Strategies for more examples on possible flow routing strategies.

By default, base level nodes are set from fixed value boundary conditions (all nodes on the left border in this example).

flow_graph = fs.FlowGraph(
    grid,
    [fs.SingleFlowRouter(), fs.MSTSinkResolver(), fs.MultiFlowRouter(1.1)],
)

Setup eroder classes (bedrock channel + hillslope) with a given set of parameter values.

spl_eroder = fs.SPLEroder(
    flow_graph,
    k_coef=1e-4,
    area_exp=0.4,
    slope_exp=1,
    tolerance=1e-5,
)

diffusion_eroder = fs.DiffusionADIEroder(grid, 0.01)

Setup Initial Conditions and External Forcing

Setup initial topography as two nearly flat plateaus (+ small random perturbations) separated by a steep escarpment. Also initialize the array for drainage area.

rng = np.random.Generator(np.random.PCG64(1234))

init_elevation = rng.uniform(0, 1, size=grid.shape)
init_elevation[:, 100:] += 400

elevation = init_elevation
drainage_area = np.empty_like(elevation)

Run the Model

Run the model for a few dozens of time steps (total simulation time: 100k years).

dt = 2e3
nsteps = 50

for step in range(nsteps):
    # flow routing
    flow_graph.update_routes(elevation)
    
    # flow accumulation (drainage area)
    flow_graph.accumulate(drainage_area, 1.0)
    
    # apply channel erosion then hillslope diffusion
    spl_erosion = spl_eroder.erode(elevation, drainage_area, dt)
    diff_erosion = diffusion_eroder.erode(elevation - spl_erosion, dt)
    
    # update topography
    elevation = elevation - spl_erosion - diff_erosion

Plot Outputs and Other Diagnostics

• Topographic elevation

fig, ax = plt.subplots(figsize=(12, 6))
plt.imshow(elevation)
plt.colorbar();

• Drainage area (log)

fig, ax = plt.subplots(figsize=(12, 6))
plt.imshow(np.log(drainage_area), cmap=plt.cm.Blues)
plt.colorbar();