wolf_sheep/wolf_sheep_with_soil_creep.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Wolf-Sheep-Grass Model with Soil Creep\n",
    "\n",
    "This notebook is meant to demonstrate coupling of an ABM implemented in Mesa and a grid-based numerical model written in Landlab.\n",
    "\n",
    "Greg Tucker, June 2020"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running the Mesa Wolf-Sheep-Grass model by itself\n",
    "\n",
    "Here we just show that you can run a Mesa model in a notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# To install dependencies for running on mybinder.org\n",
    "import sys\n",
    "!{sys.executable} -m pip install -r requirements.txt"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from wolf_sheep.model import WolfSheep \n",
    "from wolf_sheep.agents import GrassPatch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ws = WolfSheep(grass=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def generate_grass_map(model):\n",
    "    grass_map = np.zeros((model.grid.width, model.grid.height))\n",
    "    for cell in model.grid.coord_iter():\n",
    "        cell_content, x, y = cell\n",
    "        for agent in cell_content:\n",
    "            if type(agent) is GrassPatch:\n",
    "                if agent.fully_grown:\n",
    "                    grass_map[x][y] = 2\n",
    "                else:\n",
    "                    grass_map[x][y] = 1\n",
    "    return grass_map\n",
    "\n",
    "def plot_grass_map(grass_map):  \n",
    "    plt.imshow(grass_map, interpolation='nearest')\n",
    "    plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ws.run_model(step_count=25)\n",
    "gm = generate_grass_map(ws)\n",
    "plot_grass_map(gm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## One-way coupling: using the grass cover in a soil-creep model\n",
    "\n",
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ws = WolfSheep(grass=True)\n",
    "ws.run_model(step_count=3)\n",
    "gm = generate_grass_map(ws)\n",
    "plot_grass_map(gm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from landlab import RasterModelGrid, imshow_grid\n",
    "from landlab.components import LinearDiffuser"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a grid the same size as the W-S-G model's grid\n",
    "rmg = RasterModelGrid((ws.grid.width, ws.grid.height))\n",
    "\n",
    "# Create elevation field and have it slope down to the south at 10% gradient\n",
    "elev = rmg.add_zeros('topographic__elevation', at='node')\n",
    "elev[:] = 0.1 * rmg.y_of_node\n",
    "\n",
    "# Have one open boundary on the south side\n",
    "rmg.set_closed_boundaries_at_grid_edges(True, True, True, False)\n",
    "\n",
    "# Remember the starting elevation so we can calculate cumulative erosion/deposition\n",
    "initial_elev = np.zeros(rmg.number_of_nodes)\n",
    "initial_elev[:] = elev\n",
    "\n",
    "# Create a field for the creep coefficient, and set parameters for two\n",
    "# rates: slow (full grass cover) and fast (partial or \"eaten\" grass cover)\n",
    "creep_coef = rmg.add_zeros('creep_coefficient', at='node')\n",
    "fast_creep = 0.1\n",
    "slow_creep = 0.001\n",
    "\n",
    "# Assign the higher creep coefficient to cells where the grass has\n",
    "# been eaten and not yet recovered; the slower value is assigned to\n",
    "# \"fully grown\" grass patches.\n",
    "creep_coef[gm.flatten()==1] = fast_creep\n",
    "creep_coef[gm.flatten()==2] = slow_creep\n",
    "\n",
    "# Instantiate a LinearDiffuser (soil creep) Landlab component\n",
    "diffuser = LinearDiffuser(rmg, linear_diffusivity=creep_coef)\n",
    "\n",
    "# Set the time step duration\n",
    "dt = 0.2 * rmg.dx * rmg.dx / fast_creep\n",
    "print('Time step duration is ' + str(dt) + ' years.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run the soil creep model\n",
    "for i in range(50):\n",
    "    diffuser.run_one_step(dt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate and plot the erosion/deposition patterns\n",
    "ero_dep = elev - initial_elev\n",
    "maxchange = np.amax(np.abs(ero_dep))\n",
    "imshow_grid(rmg, ero_dep, vmin=-maxchange, vmax=maxchange, cmap='coolwarm_r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the grass cover again\n",
    "imshow_grid(rmg, gm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imshow_grid(rmg, elev)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Interestingly, erosion tends to occur at locations where grass cover upslope captures incoming soil."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Two-way feedback\n",
    "\n",
    "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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ws = WolfSheep(grass=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "initial_soil_depth = 0.2\n",
    "min_depth_for_grass = 0.2\n",
    "hstar = 0.2\n",
    "fast_creep = 0.1\n",
    "slow_creep = 0.001"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a grid the same size as the W-S-G model's grid\n",
    "rmg = RasterModelGrid((ws.grid.width, ws.grid.height))\n",
    "\n",
    "# Create elevation field and have it slope down to the south at 10% gradient\n",
    "elev = rmg.add_zeros('topographic__elevation', at='node')\n",
    "elev[:] = 0.1 * rmg.y_of_node\n",
    "\n",
    "# Have one open boundary on the south side\n",
    "rmg.set_closed_boundaries_at_grid_edges(True, True, True, False)\n",
    "\n",
    "# Remember the starting elevation so we can calculate cumulative erosion/deposition\n",
    "initial_elev = np.zeros(rmg.number_of_nodes)\n",
    "initial_elev[:] = elev\n",
    "\n",
    "# Also remember the elevation of the prior time step, so we can difference\n",
    "prior_elev = np.zeros(rmg.number_of_nodes)\n",
    "\n",
    "# Create a field for the creep coefficient, and set parameters for two\n",
    "# rates: slow (full grass cover) and fast (partial or \"eaten\" grass cover)\n",
    "creep_coef = rmg.add_zeros('creep_coefficient', at='node')\n",
    "\n",
    "# Create a soil-thickness field\n",
    "soil = rmg.add_zeros('soil__depth', at='node')\n",
    "soil[:] = initial_soil_depth\n",
    "\n",
    "# Instantiate a LinearDiffuser (soil creep) Landlab component\n",
    "diffuser = LinearDiffuser(rmg, linear_diffusivity=creep_coef)\n",
    "\n",
    "# Set the time step duration\n",
    "dt = 0.2 * rmg.dx * rmg.dx / fast_creep\n",
    "print('Time step duration is ' + str(dt) + ' years.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def limit_grass_by_soil(wsg_model, soil, min_soil_depth):\n",
    "    soilmatrix = soil.reshape((wsg_model.width, wsg_model.height))\n",
    "    for cell in wsg_model.grid.coord_iter():\n",
    "        cell_content, x, y = cell\n",
    "        if soilmatrix[x][y] < min_soil_depth:\n",
    "            for agent in cell_content:\n",
    "                if type(agent) is GrassPatch:\n",
    "                    agent.fully_grown = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Main loop\n",
    "for i in range(50):\n",
    "    \n",
    "    # Assign the higher creep coefficient to cells where the grass has\n",
    "    # been eaten and not yet recovered; the slower value is assigned to\n",
    "    # \"fully grown\" grass patches.\n",
    "    gm = generate_grass_map(ws)\n",
    "    creep_coef[gm.flatten()==1] = fast_creep\n",
    "    creep_coef[gm.flatten()==2] = slow_creep\n",
    "\n",
    "    # Adjust the creep coefficient to account for soil depth\n",
    "    creep_coef *= 1.0 - np.exp( -soil / hstar)\n",
    "\n",
    "    # Run the soil-creep model\n",
    "    prior_elev[:] = elev\n",
    "    diffuser.run_one_step(dt)\n",
    "    \n",
    "    # Update the soil cover\n",
    "    soil += elev - prior_elev\n",
    "\n",
    "    # Update the grass cover\n",
    "    limit_grass_by_soil(ws, soil, min_depth_for_grass)\n",
    "    \n",
    "    # Run the W-S-G model\n",
    "    ws.step()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate and plot the erosion/deposition patterns\n",
    "ero_dep = elev - initial_elev\n",
    "maxchange = np.amax(np.abs(ero_dep))\n",
    "imshow_grid(rmg, ero_dep, vmin=-maxchange, vmax=maxchange, cmap='coolwarm_r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imshow_grid(rmg, soil)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imshow_grid(rmg, gm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}