Note

This page was generated from a jupyter notebook.

Wolf-Sheep-Grass Model with Soil Creep

This notebook demonstrates coupling of an ABM implemented in Mesa and a grid-based numerical model written in Landlab. The example is the canonical “wolf-sheep-grass” example of an agent-based model. Here we add an additional twist: when sheep eat grass, the soil beneath becomes more easily mobile. This then influences soil transport: the transport efficiency is higher where the grass is “damaged”. An additional feedback lies in the thickness of the soil: grass will not grow if the soil is too thin.

The rules in this example are deliberately simple. The main goal of this tutorial is to illustrate the mechanics of building an integrated model that combines agent-based elements (via Mesa) with continuum-based elements (via Landlab) on a shared grid.

(Greg Tucker, June 2020; most recent update November 2021)

Running the Mesa Wolf-Sheep-Grass model by itself

To start, here’s an example of how to run a Mesa model in a notebook. First, we’ll run a check to make sure Mesa is installed and available; if it is not, follow the instructions in the message to install it, then re-start the kernel (Kernel => Restart) and continue.

try:
    import mesa

    print("Mesa version", mesa.__version__)
except ModuleNotFoundError:
    print(
        """
Mesa needs to be installed in order to run this notebook.

Normally Mesa should be pre-installed alongside the Landlab notebook collection.
But it appears that Mesa is not already installed on the system on which you are
running this notebook. You can install Mesa from a command prompt using either:

`conda install -c conda-forge mesa`

or

`pip install mesa`
"""
    )
    raise

Next, we’ll import the Mesa example Wolf-Sheep-Grass model from the examples collection (more info here.

from mesa.examples.advanced.wolf_sheep.agents import GrassPatch
from mesa.examples.advanced.wolf_sheep.model import WolfSheep
from mesa.experimental.devs import ABMSimulator

Create an instance of the WolfSheep model, with the grass option set to True:

simulator = ABMSimulator()
ws = WolfSheep(simulator=simulator, grass=True)

Define a function to set up an array representing the growth status of grass on the model grid (in other words, extract the information from the model’s GrassPatch agents), as well as a function to plot the current grass status. This is really a translation of data structures: the Mesa model stores data inside agents, which themselves reside at particular grid cells. Here we want to extract the information pertaining to the status of each cell’s GrassPatch—is it fully grown or “damaged”—and store that information in a simple 2D numpy array.

import copy

import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np

ground_cover_cmap = copy.copy(mpl.colormaps["YlGn"])


def generate_grass_map(model):
    grass_map = np.zeros((model.grid.width, model.grid.height))
    for cell in model.grid:
        (x, y) = cell.coordinate
        cell_content = cell.agents
        for agent in cell_content:
            if type(agent) is GrassPatch:
                if agent.fully_grown:
                    grass_map[x][y] = 2
                else:
                    grass_map[x][y] = 1
    return grass_map


def plot_grass_map(grass_map):
    plt.imshow(grass_map, interpolation="nearest", cmap=ground_cover_cmap)
    plt.colorbar()

Run the model and display the results:

for _ in range(25):
    ws.step()

gm = generate_grass_map(ws)
plot_grass_map(gm)

One-way coupling: using the grass cover in a soil-creep model

Here we initialize and run the W-S-G model for a short duration. We then extract its map of fully grown versus damaged grass, and use that to set the soil creep coefficient in a model of downslope soil creep. The point here is just to show that it’s pretty easy to use a grid from a Mesa model as input to a Landlab-built model.

simulator = ABMSimulator()
ws = WolfSheep(simulator=simulator, grass=True)

for _ in range(3):
    ws.step()

gm = generate_grass_map(ws)
plot_grass_map(gm)

Import from Landlab a RasterModelGrid (which will be Landlab’s version of the model grid), the imshow_grid function (for plotting Landlab grid fields), and the LinearDiffuser component (which will implement down-slope soil creep).

import copy

import matplotlib as mpl

from landlab import RasterModelGrid, imshow_grid
from landlab.components import LinearDiffuser

# Create a grid the same size as the W-S-G model's grid
rmg = RasterModelGrid((ws.grid.width, ws.grid.height))

# Create elevation field and have it slope down to the south at 10% gradient
elev = rmg.add_zeros("topographic__elevation", at="node")
elev[:] = 0.1 * rmg.y_of_node

# Have one open boundary on the south side
rmg.set_closed_boundaries_at_grid_edges(True, True, True, False)

# Remember the starting elevation so we can calculate cumulative erosion/deposition
initial_elev = np.zeros(rmg.number_of_nodes)
initial_elev[:] = elev

# Create a field for the creep coefficient, and set parameters for two
# rates: slow (full grass cover) and fast (partial or "eaten" grass cover)
creep_coef = rmg.add_zeros("creep_coefficient", at="node")
fast_creep = 0.1
slow_creep = 0.001

# Assign the higher creep coefficient to cells where the grass has
# been eaten and not yet recovered; the slower value is assigned to
# "fully grown" grass patches.
creep_coef[gm.flatten() == 1] = fast_creep
creep_coef[gm.flatten() == 2] = slow_creep

# Instantiate a LinearDiffuser (soil creep) Landlab component
diffuser = LinearDiffuser(rmg, linear_diffusivity=creep_coef)

# Set the time step duration
dt = 0.2 * rmg.dx * rmg.dx / fast_creep
print(f"Time step duration is {dt} years.")

# Run the soil creep model
for i in range(50):
    diffuser.run_one_step(dt)

# Calculate and plot the erosion/deposition patterns
ero_dep = elev - initial_elev
maxchange = np.amax(np.abs(ero_dep))
imshow_grid(
    rmg,
    ero_dep,
    vmin=-maxchange,
    vmax=maxchange,
    cmap=copy.copy(mpl.colormaps["coolwarm_r"]),
    colorbar_label="Cumulative deposition (+) or erosion (-), m",
)

# Plot the grass cover again
imshow_grid(
    rmg, gm, cmap=ground_cover_cmap, colorbar_label="Ground cover (1 = bare, 2 = grass)"
)

imshow_grid(
    rmg,
    elev,
    cmap=copy.copy(mpl.colormaps["pink"]),
    colorbar_label="Elevation above base of slope (m)",
)

Interestingly, erosion tends to occur at locations where grass cover upslope captures incoming soil.

So far, however, this is just one-way feedback: the previously damaged grass patches, as calculated in the wolf-sheep-grass ABM, become susceptible to erosion, but this does not (yet) feed back into future grass growth or erosional loss. Let’s turn to that next.

Two-way feedback

Here, we explore two-way feedback by running the two models iteratively. We track soil thickness, and “damage” any grass where the soil is thinner than a given amount. We also limit soil flux according to its thickness, so that absent soil cannot move.

These rules are deliberately simple. One could make the model more realistic by, for example, setting the grass regrowth time (a property of the GrassPatch agents in the ABM) to a value that depends on the thickness of the soil (a Landlab field).

simulator = ABMSimulator()
ws = WolfSheep(
    simulator=simulator,
    initial_sheep=20,
    initial_wolves=10,
    grass=True,
    grass_regrowth_time=15,  # give grass a fighting chance...
)

initial_soil_depth = 0.3
min_depth_for_grass = 0.2
hstar = 0.2
fast_creep = 0.1
slow_creep = 0.001

# Create a grid the same size as the W-S-G model's grid
rmg = RasterModelGrid((ws.grid.width, ws.grid.height))

# Create elevation field and have it slope down to the south at 10% gradient
elev = rmg.add_zeros("topographic__elevation", at="node")
elev[:] = 0.1 * rmg.y_of_node

# Have one open boundary on the south side
rmg.set_closed_boundaries_at_grid_edges(True, True, True, False)

# Remember the starting elevation so we can calculate cumulative erosion/deposition
initial_elev = np.zeros(rmg.number_of_nodes)
initial_elev[:] = elev

# Also remember the elevation of the prior time step, so we can difference
prior_elev = np.zeros(rmg.number_of_nodes)

# Create a field for the creep coefficient, and set parameters for two
# rates: slow (full grass cover) and fast (partial or "eaten" grass cover)
creep_coef = rmg.add_zeros("creep_coefficient", at="node")

# Create a soil-thickness field
soil = rmg.add_zeros("soil__depth", at="node")
soil[:] = initial_soil_depth

# Instantiate a LinearDiffuser (soil creep) Landlab component
diffuser = LinearDiffuser(rmg, linear_diffusivity=creep_coef)

# Set the time step duration
dt = 0.2 * rmg.dx * rmg.dx / fast_creep
print("Time step duration is {dt} years.")

Next we define a new function limit_grass_by_soil that will render any GrassPatches “non-fully-grown” if the soil is thinner than a specified minimum value. In other words, we represent soil limitation with a simple threshold in which the grass in any cell with soil thinner than the threshold can never be fully grown. Again, a more realistic way to do this might be to reduce the regrowth rate, but our simple threshold treatment will serve for the purpose of showing how we can use data from a Landlab field to influence data associated with spatially distributed agents in a Mesa model:

def limit_grass_by_soil(wsg_model, soil, min_soil_depth):
    soilmatrix = soil.reshape((wsg_model.width, wsg_model.height))
    for cell in wsg_model.grid:
        (x, y) = cell.coordinate
        cell_content = cell.agents
        if soilmatrix[x][y] < min_soil_depth:
            for agent in cell_content:
                if type(agent) is GrassPatch:
                    agent.fully_grown = False

Run the integrated model in a time loop. Our algorithm performs the following sequence of calculations in each iteration:

• Get a copy of the current grass status as a 2D array

• Update the soil-creep coefficient Landlab field according to the grass status and the soil thickness

• Run soil creep for one time step and update the soil thickness (we could have used a DepthDependentLinearDiffuser for this, but here a simpler approach will suffice)

• Set grass in any cells with insufficient soil to be non-fully-grown

• Run the wolf-sheep-grass model for one time step

The data exchange happens in two function calls. generate_grass_map translates grass status data from the Mesa model’s data structure to a Landlab field, and limit_grass_by_soil translates Landlab’s soil thickness field into a restriction on grass status in the Mesa model’s GrassPatch agents.

# Main loop
for _ in range(50):
    # Assign the higher creep coefficient to cells where the grass has
    # been eaten and not yet recovered; the slower value is assigned to
    # "fully grown" grass patches.
    gm = generate_grass_map(ws)
    creep_coef[gm.flatten() == 1] = fast_creep
    creep_coef[gm.flatten() == 2] = slow_creep

    # Adjust the creep coefficient to account for soil depth
    creep_coef *= 1.0 - np.exp(-soil / hstar)

    # Run the soil-creep model
    prior_elev[:] = elev
    diffuser.run_one_step(dt)

    # Update the soil cover
    soil += elev - prior_elev

    # Update the grass cover
    limit_grass_by_soil(ws, soil, min_depth_for_grass)

    # Run the W-S-G model
    ws.step()

The next few plots examine the results to illustrate how the interaction of soil creep and grass consumption by mobile agents (sheep) has influenced the landscape:

# Calculate and plot the erosion/deposition patterns
ero_dep = elev - initial_elev
maxchange = np.amax(np.abs(ero_dep))
imshow_grid(
    rmg,
    ero_dep,
    vmin=-maxchange,
    vmax=maxchange,
    cmap="coolwarm_r",
    colorbar_label="Depth of soil accumulation (+) or loss (-), m",
)

# Soil thickness
imshow_grid(rmg, soil, colorbar_label="Soil thickness, m")

# Ground cover
imshow_grid(
    rmg, gm, cmap=ground_cover_cmap, colorbar_label="Ground cover (1 = bare, 2 = grass)"
)

Here soil erosion at the top of the slope inhibits grass cover, while soil accumulation at the base of the slope allows grass to continue to grow.

---

Generated by nbsphinx from a Jupyter notebook.