Climate changes, people change and land also changes with time. We can`t believe that the river networks will remain the same over the next 1000 years or that mountains and depressions will have the same elevation in the next 10000 years. But changes are not related to huge time frames, they can occur in decades or years at lower rates as well. In order to evaluate those changes we need some formulation of the key components of land evolution: fluvial, hillslope and uplift. We have developed a tutorial with Python and the Landlab library to simulate the land evolution at basin scale for 100 thousand years; inputs come from geospatial rasters and output data is exported as Ascii raster files.

You will need a conda environment with geospatial tools for the tutorial. Create the environment following this link:

hatarilabs.com/ih-en/how-to-install-python-geopandas-in-windows-on-a-conda-environment-tutorial

Tutorial

Code

Modeling Land Evolution at Basin Scale with Python and Landlab

# import generic packages
import numpy as np
from matplotlib import pyplot as plt

# import geospatial packages
import rasterio
from rasterio.plot import show

# import landlab components
from landlab import HexModelGrid, RasterModelGrid, imshow_grid
from landlab.components import (
    DepressionFinderAndRouter,
    FlowAccumulator,
    LinearDiffuser,
    StreamPowerEroder,
)

# Open raster image 
rasterObj = rasterio.open('../data/DEM_18S_clip.tif')
show(rasterObj)

#extract array from raster
elevArray = rasterObj.read(1)
plt.imshow(elevArray)

# define parameters for fluvial, hillslope and uplift components

uplift_rate = 0.001  # [m/yr], initially set at 0.001
K_sp = 1.0e-5  # units vary depending on m_sp and n_sp, initially set at 1e-5
m_sp = 0.5  # exponent on drainage area in stream power equation, initially 0.5
n_sp = 1.0  # exponent on slope in stream power equation, initially 1.
K_hs = 0.05  # [m^2/yr], initially 0.05

#create grid from raster attributes

nrows = rasterObj.height  # number of raster cells on each side, initially 150
ncols = rasterObj.width
dxy = rasterObj.transform.a  # side length of a raster model cell, or resolution [m], initially 50

# define a landlab raster

mg = RasterModelGrid(shape=(nrows, ncols), 
                     xy_spacing=dxy,
                     xy_of_lower_left=(rasterObj.bounds[0],rasterObj.bounds[1]))

# show number of rows, cols and resolution

print(nrows, ncols, dxy)

167 179 90.3109901477272

# define temporal distribution 

dt = 1000  # time step [yr]
total_time = 0  # amount of time the landscape has evolved [yr]
tmax = 100000  # time for the model loop to run [yr]

t = np.arange(0, tmax, dt)

# create a dataset of zero values

zr = mg.add_zeros("topographic__elevation", at="node")
imshow_grid(mg, "topographic__elevation")

# apply cell elevation to defined arrray
zr += elevArray[::-1,:].ravel()
imshow_grid(mg, "topographic__elevation")

# initialize the process components of the equation

frr = FlowAccumulator(mg)  # intializing flow routing
spr = StreamPowerEroder(
    mg, K_sp=K_sp, m_sp=m_sp, n_sp=n_sp, threshold_sp=0.0
)  # initializing stream power incision
dfn = LinearDiffuser(
    mg, linear_diffusivity=K_hs, deposit=False
)  # initializing linear diffusion
df = DepressionFinderAndRouter(mg) # Initializing the pit finder

# run the process for every time step 

for ti in t:
    zr[mg.core_nodes] += uplift_rate * dt  # uplift the landscape
    dfn.run_one_step(dt)  # diffuse the landscape
    frr.run_one_step()  # route flow
    # df.map_depressions()
    spr.run_one_step(dt)  # fluvial incision
    total_time += dt  # update time keeper
    print("\r"+str(total_time), end=' ')

# plot final surface elevation

fig = plt.figure(figsize=(18,12))

imshow_grid(
    mg, "topographic__elevation", grid_units=("m", "m"), var_name="Elevation (m)", shrink=0.5
)

title_text = f"$K_{{sp}}$={K_sp}; $K_{{hs}}$={K_hs}; $time$={total_time}; $dx$={dxy}"

plt.title(title_text)

plt.show()

# export data as geospatial rasters

from landlab.io.esri_ascii import write_esri_ascii

files = write_esri_ascii('../Out/basinLandEvolution.asc', mg)