Modeling forest fires and tipping points

Names of group members

// put your names here!

Goals of this assignment

The primary goal of this assignment are to model the spread of a forest fire using an agent-based model. In doing so, we will:

  • Observe how complex emergent phenomena result from simple rules in a model.
  • Examine and quantify the concept of a "tipping point" in a model.

Note: This is a two-day assignment!

Why model forest fires?

While this isn't a huge problem in Michigan, the states in the western United States are currently having a tremendous problem with huge and difficult-to-control forest fires. This comes from a combination of extended drought conditions, dense woodlands, and forest management policies that suppress small fires and thus ensure that large quantities of dead, dry trees and brush are available to burn when a large fire inevitably starts (typically from lighting strikes, but occasionally from negligent campers). In recent years, this has been exacerbated by climate change, which has both caused drought conditions to be more severe and allowed tree-killing diseases and insects to flourish, which produces more dead, easily-burned wood.

These forest fires destroy ecosystems and peoples' homes and other property, and can result in the loss of substantial human and animal life. A key challenge in forest management is to attempt to contain these huge forest fires once they start, in order to protect human lives, settlements, and infrastructure. To that end, it is critical to have models of how fire spreads in various condition - see, for example, the Open Wildland Fire Modeling group.

More generally, the type of model that we're going to create is an example of a "percolation" model, where one substance (in this case, fire) moves through another substance (in this case, forest). This type of problem is interesting in a variety of fields, including geology (oil or water percolating through rock, sand, or soil), human behavior (crowd movement in amusement parks), and in physics (understanding how two materials mix together).

What is a "tipping point"?

This model also demonstrates the concept of a "critical threshold" or a "tipping point". This is a phenomenon that occurs when a small change in an input parameter results in a large change in outcome. This is a phenomenon that shows up in both simple and complex models, and happens in such varied circumstances as forest fires, the spread of disease in populations, and the transfer of information within a population.

Let's try the model!

Before we try to implement it, let's get a feel for how this model works. We're going to do this using the NetLogo Fire model.

Finding the model: Open NetLogo, click on "File" and then "Models Library". In the "Sample Models" directory, click on "Earth Science" and then "Fire"

Running the model: The only parameter that you can vary is "density" - the fraction of cells that contain trees (and thus the density of trees in the forest). You can drag the slider back and forth to change the density. After you do that, click the "setup" button and then "go".

What do you observe? Try setting the "density" value to various numbers between 0 and 99, and see how the wildfire spreads. What happens when the value is low vs. when it is high? And, is there any particular value where the behavior changes very rapidly as you change the density of trees? If so, try to home in on that number and report it below.

Put your observations here!


The rules for our model

In this model, we'll create a two-dimensional, square array with sides that are N cells long that represents the forested region we're interested in. The cells in the array can have three values, 0 (empty), 1 (trees), and 2 (on fire). At the beginning of the model, a user-specified fraction of the cells $\mathrm{f_{trees\_start}}$ (equivalent to the NetLogo density parameter) are randomly filled with trees, and the remaining cells are empty. One edge of the board (say, the entire i=0 column) is set on fire.

Each cell has a "neighborhood" that is composed of its four neighbors to the left, right, above, and below it. (Note: not the diagonal elements!) If a cell is along one of the edges of the array, only consider the neighbors that it has, and don't try to go out of the bounds of the array!

The model takes steps forward in time, where every cell is modified based on the previous step. The model evolves as follows:

  • If the cell was empty last turn, it stays empty this turn.
  • If the cell is a tree and any of its neighbors were on fire last turn, it catches on fire.
  • If the cell was on fire last turn, the fire has consumed all of the trees and it is now empty.

The model evolves forward in time until all of the fires have burned out. After this happens, you can calculate the fraction of the cells that still have trees at the end of the model ($\mathrm{f_{trees\_end}}$) and the fraction of cells that are empty ($\mathrm{f_{empty}}$). The fraction of burned cells, $\mathrm{f_{burned}}$, is just the difference between the fraction of cells that were initially trees and the fraction of cells that are trees at the end of the model - in other words,

$\mathrm{f_{burned}} = \mathrm{f_{trees\_end}} - \mathrm{f_{trees\_start}}$

Your mission!

Your mission is to answer the question: "How does the spread of fire relate to the density of the forest?" More precisely, we're asking "How does $\mathrm{f_{burned}}$ depend on $\mathrm{f_{trees\_start}}$?"

To to this, we want to do the following:

First, as a group create pseudo-code that shows how you plan to implement this model for an arbitrary value of $\mathrm{f_{trees\_start}}$. Make sure that you think about how to set up the initial conditions, how to evolve the model, and how to calculate the fraction of trees and empty cells that remain in the end. Use the whiteboard and make sure to take a picture of it!

Second, implement this model as code. Demonstrate that it works on a model board 40 cells on a side for a few values of $\mathrm{f_{trees\_start}}$ - say, 0.25, 0.5, 0.75, and 1.0. Do you notice any differences in behavior as you change $\mathrm{f_{trees\_start}}$?

Finally, loop over many values of $\mathrm{f_{trees\_start}}$ - say, values from 0.01 to 1.0, in steps of 0.01 - and run the model many times. Keep track of the fraction of cells that are burned, and show how $\mathrm{f_{burned}}$ and $\mathrm{f_{trees\_start}}$ relate to each other.

Note: This is a two-day assignment. Also, feel free to use any code from the Numpy 2D array tutorial, the pre-class assignment, or anything else you've done this semester!


In [ ]:
# put your code here, adding any additional cells as necessary!

# standard includes
import numpy as np
import numpy.random as rand
%matplotlib inline
import matplotlib.pyplot as plt

# Next we are going to import some specific libraries we will 
# use to get the animation to work cleanly
from IPython.display import display, clear_output
import time

In [ ]:
# ADAPTED FROM THE 2D NUMPY TUTORIAL

# function plotgrid() takes in a 2D array and uses pyplot to make a plot.
# this function returns no values!

def plotgrid(myarray):
    
    # first create two vectors based on the x and y sizes of the grid
    x_range = np.linspace(0, myarray.shape[0], myarray.shape[0]) 
    y_range = np.linspace(0, myarray.shape[0], myarray.shape[0])
    
    # use the numpy meshgrid function to create two matrices 
    # of the same size as myarray with x and y indexes
    x_indexes, y_indexes = np.meshgrid(x_range, y_range)
    
    # make a list of all the x and y indexes that are either trees or fire.
    tree_x = x_indexes[myarray == 1];   # trees!
    tree_y = y_indexes[myarray == 1]; 
    fire_x = x_indexes[myarray == 2];   # fire!
    fire_y = y_indexes[myarray == 2]; 
    
    # plot the trees and fire using appropriate colors.
    # make the size of the polygons 
    # larger than the default so they're easy to see!
    plt.plot(tree_x,tree_y, 'gs',markersize=10)   
    plt.plot(fire_x,fire_y, 'rs',markersize=10)  
    
    #Set the x and y limits to include a space of overlap so we don't cut off the shapes
    plt.ylim([-1,myarray.shape[0]+1]) 
    plt.xlim([-1,myarray.shape[0]+1])

In [ ]:
def set_board(board_size=40,prob_tree=0.5):
    '''
    Creates the initial game board.

    Inputs:
        board_size: length of one edge of the board
        prob_tree: probability that a given cell is a tree.

    Outputs a 2D numpy array with values set to either 0, 1, or 2
        (empty, tree, or fire)
    '''
    
    # all cells initialized to 'empty' (0) by default
    game_board = np.zeros((board_size,board_size),dtype='int64')
    
    # loop over board and roll the dice; if the random number is less
    # than prob_tree, make it a tree.
    for i in range(board_size):
        for j in range(board_size):
            if rand.random() <= prob_tree:
                game_board[i,j] = 1

    # set the whole i=0 edge of the board on fire. We're arsonists!
    game_board[0,:] = 2
    
    return game_board


def advance_board(game_board):
    '''
    Advances the game board using the given rules.
    Input: the game board.
    Output: a new game board
    '''
    
    # create a new array that's just like the original one, but initially 
    # set to all zeros (i.e., totally empty)
    new_board = np.zeros_like(game_board)
    
    # loop over each cell in the board and decide what to do.
    for i in range(game_board.shape[0]):
        for j in range(game_board.shape[1]):
    
            # if the cell was empty last turn, it's still empty.
            # if it was on fire last turn, it's now empty.
            if game_board[i,j] == 0 or game_board[i,j] == 2:
                new_board[i,j] = 0
    
            # now, if it's a tree we have to decide what to do.
            if game_board[i,j] == 1:
                
                # initially make it a tree
                new_board[i,j] = 1

                # If one of the neighboring cells was on fire last turn, 
                # this cell is now on fire!
                if i > 0:
                    if game_board[i-1,j] == 2:
                        new_board[i,j] = 2
                    
                if j > 0:
                    if game_board[i,j-1] == 2:
                        new_board[i,j] = 2

                if i < game_board.shape[0]-1:
                    if game_board[i+1,j] == 2:
                        new_board[i,j] = 2
                        
                if j < game_board.shape[1]-1:                    
                    if game_board[i,j+1] == 2:
                        new_board[i,j] = 2

    # return the new board
    return new_board


def calc_stats(game_board):
    '''
    Calculates the fraction of cells on the game board that are 
    a tree or are empty.
    
    Input: the game board
    
    Output: fraction that's empty, fraction that's covered in trees.
    '''
    
    # use numpy to count up the fractions that are empty
    frac_empty = (game_board == 0).sum() / game_board.size

    # do the same for trees
    frac_tree = (game_board == 1).sum() / game_board.size
    
    # return it!
    return frac_empty, frac_tree

In [ ]:
'''
This cell does a quick test run with a user-specified fraction of trees, 
and makes an animation.  We don't necessarily expect the students to do 
something this fancy, but it would certainly be nice.  All of the plotting
stuff has been adapted from the Numpy 2D tutorial animation section.
'''

# some initial parameters
prob_tree=0.7
board_size = 40

# set up a figure
fig = plt.figure(figsize=(6,6))

# set up a game board
game_board = set_board(board_size=board_size, prob_tree=prob_tree)

# plot the initial board
plotgrid(game_board)

# need to decide when to stop - we're just going to ask if something
# is still on fire!
on_fire = True

# while there are fires, keep going
while on_fire == True:

    # advance the game board
    game_board = advance_board(game_board)
    
    # do a bunch of plotting stuff, using code from the Numpy 2D array tutorial.
    plotgrid(game_board)
    time.sleep(0.05)  # set this to be small because I don't want to wait!
    clear_output(wait=True)
    display(fig)
    fig.clear()

    # calculate the fractions that are empty and are trees
    frac_empty, frac_trees = calc_stats(game_board)

    # if 100% of cells are empty or are trees, there is no fire and we can stop.
    if frac_empty + frac_trees == 1.0:
        on_fire = False

# clean it up.
plt.close()                 # Close dynamic display

In [ ]:
board_size = 40

f_tree = []
f_burned = []

for tree_fraction in np.arange(0.01,1.01,0.01):

    print("tree fraction is:", tree_fraction)
    
    game_board = set_board(board_size=board_size, prob_tree=tree_fraction)

    on_fire = True
    while on_fire == True:
        game_board = advance_board(game_board)
        frac_empty, frac_trees = calc_stats(game_board)
        if frac_empty + frac_trees == 1.0:
            on_fire = False

    #frac_empty, frac_trees = calc_stats(game_board)
    f_tree.append(tree_fraction)
    f_burned.append(frac_empty - (1.0-tree_fraction))
    

plt.plot(f_tree, f_burned)
plt.xlabel("tree fraction")
plt.ylabel("burned fraction (normalized)")

A question

Describe how you see the model's behavior. Do you observe a "tipping point" in $\mathrm{f_{burned}}$ - that is to say, an abrupt change in behavior in $\mathrm{f_{burned}}$ - as you increase $\mathrm{f_{trees\_start}}$? Does it agree with what you saw when you were experimenting with the NetLogo model?

put your answer here!

ANSWER: As the fraction of trees goes up, the fire burns much more steadily across the forest. At very low fractions, it basically burns itself out, while at large fractions it burns all the way across the board. There is a clear "tipping point" around f_trees_start of 0.6. (The precise value is 0.59)

Assignment wrapup

Please fill out the form that appears when you run the code below. You must completely fill this out in order to receive credit for the assignment!


In [ ]:
from IPython.display import HTML
HTML(
"""
<iframe 
	src="https://goo.gl/forms/gXClWJLmneatfgew2?embedded=true" 
	width="80%" 
	height="1200px" 
	frameborder="0" 
	marginheight="0" 
	marginwidth="0">
	Loading...
</iframe>
"""
)

Turn it in!

Turn this assignment in to the Day 23 dropbox in the "in-class activities" folder.