Note

This page was generated from a jupyter notebook.

WARNING: This tutorial has not been updated to work with Landlab 2.0 and is thus not tested to verify that it will run.

Tutorial For Cellular Automaton Vegetation Model Coupled With Ecohydrologic Model

This tutorial demonstrates implementation of the Cellular Automaton Tree-GRass-Shrub Simulator (CATGRaSS) [Zhou et al., 2013] on a digital elevation model (DEM). This model is built using components from the Landlab component library. CATGRaSS is a spatially explicit model of plant coexistence. It simulates local ecohydrologic dynamics (soil moisture, transpiration, biomass) and spatial evolution of tree, grass, and shrub Plant Functional Types (PFT) driven by rainfall and solar radiation.

Each cell in the model grid can hold a single PFT or remain empty. Tree and shrub plants disperse seeds to their neighbors. Grass seeds are assumed to be available at each cell. Establishment of plants in empty cells is determined probabilistically based on water stress for each PFT. Plants with lower water stress have higher probability of establishment. Plant mortality is simulated probabilistically as a result of aging and drought stress. Fires and grazing will be added to this model soon.

This model (driver) contains: - A local vegetation dynamics model that simulates storm and inter-storm water balance and ecohydrologic fluxes (ET, runoff), and plant biomass dynamics by coupling the following components: - PrecipitationDistribution - Radiation - PotentialEvapotranspiration - SoilMoisture - Vegetation

  • A spatially explicit probabilistic cellular automaton component that simulates plant competition by tracking establishment and mortality of plants based on soil moisture stress:

    • VegCA

To run this Jupyter notebook, please make sure that the following files are in the same folder: - cellular_automaton_vegetation_DEM.ipynb (this notebook) - Inputs_Vegetation_CA.txt (Input parameters for the model) - Ecohyd_functions_DEM.py (Utility functions)

[Ref: Zhou, X, E. Istanbulluoglu, and E.R. Vivoni. “Modeling the ecohydrological role of aspect-controlled radiation on tree-grass-shrub coexistence in a semiarid climate.” Water Resources Research 49.5 (2013): 2872-2895]

In this tutorial, we are going to work with a landscape in central New Mexico, USA, where aspect controls the organization of PFTs. The climate in this area is semi-arid with Mean Annual Precipitation (MAP) of 254 mm [Zhou et. al 2013]. We will do the following: - Import a landscape - Initialize the landscape with random distribution of PFTs - Run the coupled Ecohydrology and cellular automata plant competition model for 50 years - Visualize and examine outputs

Let’s walk through the code:

Import the required libraries:

[ ]:
%matplotlib inline
[ ]:
import time

import numpy as np
from Ecohyd_functions_DEM import Create_PET_lookup, Empty_arrays, Initialize_, Plot_

from landlab import RasterModelGrid as rmg, load_params
from landlab.io import read_esri_ascii

Note: ‘Ecohyd_functions_DEM.py’ is a utility script that contains ‘functions’, which instantiates components and manages inputs and outputs, and help keep this driver concise. Contents of ‘Ecohyd_functions_DEM.py’ can be a part of this driver (current file), however left out to keep driver concise.

We will use two grids in this driver. One grid will represent the actual landscape or domain (e.g., created from a DEM). Another grid, with one cell for each of the plant functional types (PFTs), will be used to create Potential Evapotranspiration (PET) lookup tables.

  • grid: This grid represents the actual landscape. Each cell can be occupied by a single PFT such as tree, shrub, grass, or can be empty (bare). In this example we assume that the elevation field and the vegetation field has the same resolution.

  • grid1: This grid will be used to compute plant-specific PET at a point. Spatially distributed PET Lookup arrays (for all days of the year) will be created for each PFT based on these point values.

Note: In this tutorial, the physical ecohydrological components and cellular automata plant competition will be run on grids with same resolution. To develop differential spatial resolutions for the two models, see the tutorial ‘cellular_automaton_vegetation_flat.ipynb’.

[ ]:
(grid, elevation) = read_esri_ascii("DEM_10m.asc")  # Read the DEM
grid1 = rmg((5, 4), xy_spacing=(5.0, 5.0))  # Representative grid

Include the input file that contains all input parameters needed for all components. This file can either be a Python dictionary or a text file that can be converted into a Python dictionary. If a text file is provided, it will be converted to a Python dictionary. Here we use an existing text file prepared for this exercise.

[ ]:
InputFile = "Inputs_Vegetation_CA_DEM.txt"
data = load_params(InputFile)  # Creates dictionary that holds the inputs

Instantiate Landlab components to simulate corresponding attributes. In this example, we shall demonstrate the use of seasonal rainfall and PFT-specific potential evapotranspiration. The instantiated objects are: - PD_D: object for dry season rainfall, - PD_W: object for wet season rainfall, - Rad: Radiation object computes radiation factor defined as the ratio of total shortwave radiation incident on a sloped surface to total shortwave radiation incident on a flat surface. - Rad_PET: Representative radiation object which is used only as an input for PET. - PET_PFT: Plant specific PET objects (we use a cosine function, fitted to calculated PET, as a function of Day Of the Year (DOY) to reduce computation overhead). This value is spatially distributed by using a radiation factor. - SM: Soil Moisture object simulates root-zone average soil moisture at each cell using inputs of potential evapotranspiration, live leaf area index, and vegetation cover. - VEG: Vegetation dynamics object simulates net primary productivity and biomass and thus leaf area index at each cell based on inputs of root-zone average soil moisture. - vegca: Cellular Automaton plant competition object. This object simulates the spatial dynamics of PFTs. It is run once every year at the end of the growing season. This object is initialized with a random cellular field of PFT. Each year, this object updates the field of PFTs based on probabilistic establishment and mortality rules employed at each cell of the modeled DEM.

Note: Almost every component in Landlab is coded as a ‘class’ (to harness the advantages of object oriented programming). An ‘object’ is the instantiation of the ‘class’ (for more information, please refer any object oriented programming book). A ‘field’ refers to a Landlab field (please refer to the Landlab documentation to learn more about Landlab fields).

Now let’s instantiate all Landlab components that we are going to use for this tutorial:

[ ]:
PD_D, PD_W, Rad, Rad_PET, PET_Tree, PET_Shrub, PET_Grass, SM, VEG, vegca = Initialize_(
    data, grid, grid1, elevation
)

Lets look at the initial organization of PFTs

[ ]:
import matplotlib as mpl
import matplotlib.pyplot as plt

cmap = mpl.colors.ListedColormap(["green", "red", "black", "white", "red", "black"])
bounds = [-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5]
norm = mpl.colors.BoundaryNorm(bounds, cmap.N)
description = "green: grass; red: shrub; black: tree; white: bare"
plt.figure(101)
grid.imshow(
    "vegetation__plant_functional_type",
    at="cell",
    cmap=cmap,
    grid_units=("m", "m"),
    norm=norm,
    limits=[0, 5],
    allow_colorbar=False,
)
plt.figtext(0.2, 0.0, description, weight="bold", fontsize=10)

Specify an approximate number of years for the model to run.

IMPORTANT: This code in numerically extensive. It might take an hour or more to run this simulation for 300 years. It is suggested to run the simulation for 50 years which might take less than 7 minutes to execute.

[ ]:
n_years = 50  # Approx number of years for model to run
# Calculate approximate number of storms per year
fraction_wet = (data["doy__end_of_monsoon"] - data["doy__start_of_monsoon"]) / 365.0
fraction_dry = 1 - fraction_wet
no_of_storms_wet = (
    8760 * (fraction_wet) / (data["mean_interstorm_wet"] + data["mean_storm_wet"])
)
no_of_storms_dry = (
    8760 * (fraction_dry) / (data["mean_interstorm_dry"] + data["mean_storm_dry"])
)
n = int(n_years * (no_of_storms_wet + no_of_storms_dry))

Create empty arrays to store spatio-temporal data over multiple iterations. The captured data can be used for plotting model outputs.

[ ]:
P, Tb, Tr, Time, VegType, PET_, Rad_Factor, EP30, PET_threshold = Empty_arrays(
    n, n_years, grid, grid1
)

To reduce computational overhead, we shall create a lookup array for plant-specific PET values for each day of the year, and slope and aspect grid.

[ ]:
Create_PET_lookup(
    Rad, PET_Tree, PET_Shrub, PET_Grass, PET_, Rad_Factor, EP30, Rad_PET, grid
)

Specify current_time (in years). current_time is the current time in the simulation.

[ ]:
# # Represent current time in years
current_time = 0  # Start from first day of Jan

# Keep track of run time for simulation—optional
Start_time = time.clock()  # Recording time taken for simulation

# declaring few variables that will be used in storm loop
time_check = 0.0  # Buffer to store current_time at previous storm
yrs = 0  # Keep track of number of years passed
WS = 0.0  # Buffer for Water Stress
Tg = 365  # Growing season in days

The loop below couples the components introduced above in a for loop until all “n” number of storms are generated. Time is advanced by the soil moisture object based on storm and interstorm durations that are estimated by the strom generator object. The ecohydrologic model is run on each storm whereas cellular automaton vegetation component is run once every year.

Note: This loop might take around 10 minutes (depending on your computer) to run for a 50 year simulation. Ignore any warnings you might see.

[ ]:
# # Run storm Loop
for i in range(0, n):
    # # Update objects
    # Calculate Day of Year (DOY)
    Julian = int(np.floor((current_time - np.floor(current_time)) * 365.0))
    # Generate seasonal storms
    # for Dry season
    if Julian < data["doy__start_of_monsoon"] or Julian > data["doy__end_of_monsoon"]:
        PD_D.update()
        P[i] = PD_D.get_storm_depth()
        Tr[i] = PD_D.get_precipitation_event_duration()
        Tb[i] = PD_D.get_interstorm_event_duration()
    # Wet Season—Jul to Sep—NA Monsoon
    else:
        PD_W.update()
        P[i] = PD_W.get_storm_depth()
        Tr[i] = PD_W.get_precipitation_event_duration()
        Tb[i] = PD_W.get_interstorm_event_duration()

    # Spatially distribute PET and its 30-day-mean (analogous to degree day)
    grid["cell"]["surface__potential_evapotranspiration_rate"] = (
        np.choose(grid["cell"]["vegetation__plant_functional_type"], PET_[Julian])
    ) * Rad_Factor[Julian]
    grid["cell"]["surface__potential_evapotranspiration_30day_mean"] = (
        np.choose(grid["cell"]["vegetation__plant_functional_type"], EP30[Julian])
    ) * Rad_Factor[Julian]

    # Assign spatial rainfall data
    grid["cell"]["rainfall__daily_depth"] = P[i] * np.ones(grid.number_of_cells)

    # Update soil moisture component
    current_time = SM.update(current_time, Tr=Tr[i], Tb=Tb[i])

    # Decide whether its growing season or not
    if Julian != 364:
        if EP30[Julian + 1, 0] > EP30[Julian, 0]:
            PET_threshold = 1
            # 1 corresponds to ETThresholdup (begin growing season)
        else:
            PET_threshold = 0
            # 0 corresponds to ETThresholddown (end growing season)

    # Update vegetation component
    VEG.update(PETthreshold_switch=PET_threshold, Tb=Tb[i], Tr=Tr[i])

    # Update yearly cumulative water stress data
    WS += (grid["cell"]["vegetation__water_stress"]) * Tb[i] / 24.0

    # Record time (optional)
    Time[i] = current_time

    # Cellular Automata
    if (current_time - time_check) >= 1.0:
        if yrs % 5 == 0:
            print(f"Elapsed time = {yrs} years")
        VegType[yrs] = grid["cell"]["vegetation__plant_functional_type"]
        grid["cell"]["vegetation__cumulative_water_stress"] = WS / Tg
        vegca.update()
        SM.initialize()
        VEG.initialize()
        time_check = current_time
        WS = 0
        yrs += 1
VegType[yrs] = grid["cell"]["vegetation__plant_functional_type"]

Time_Consumed is an optional variable that gives information about computer running time

[ ]:
Final_time = time.clock()
Time_Consumed = (Final_time - Start_time) / 60.0  # in minutes
print(f"Time_consumed = {Time_Consumed} minutes")

Save the outputs using numpy.save(). These files have ‘.nc’ extension, which can be loaded using numpy.load().

[ ]:
# # Saving
sim = "VegCA_DEM_26Jul16_"
# Save_(sim, Tb, Tr, P, VegType, yrs, Time_Consumed, Time)

Lets look at outputs.

Plots of the cellular field of PFT at specified year step can be found below where:

GRASS = Green; SHRUB = Red; TREE = Black; BARE = White;

At the end, percentage cover for each PFT is plotted with respect to time.

[ ]:
Plot_(grid, VegType, yrs, yr_step=10)

If you run this model for around 900 years, you will observe patterns of PFTs. For example, you will find more trees on north facing slopes and mostly shrubs and grass on south facing slopes, as shown below:

[ ]:
from IPython.display import Image

Image(filename="presentation.png")

If you want to explore this model further, open ‘Inputs_Vegetation_CA.txt’ and change the input parameters (e.g., initial PFT distribution percentages, storm characteristics, etc..).


Generated by nbsphinx from a Jupyter notebook.