Simulating Experiences: Excursions in Programming

``Forest Fires and Reforestation''

Mathematica in Education
Vol.3 No.2
Spring 1994
(c) TELOS/Springer-Verlag

by
 Richard J. Gaylord
Department of Materials Science
University of Illinois at Urbana-Champaign
gaylord@ux1.cso.uiuc.edu

and

Kazume Nishidate
Department of Electrical and Electronic Engineering
Iwate University, Faculty of Engineering
Morioka 020  JAPAN
nisidate@wriron1.s.kanazawa-u.ac.jp

edited by Richard J. Gaylord

Introduction

Turbulent cascading has been suggested as the underlying cause of a wide variety of phenomena, including geological upheavals (eg., volcanic eruptions and earthquakes), species extinction during evolution, and fluid turbulence. A simple one-dimensional probabilistic cellular automaton, known as the forest fire model, displays this behavior.

The Forest Fire CA

System

The forest fire CA employs a one-dimensional lattice of length n, with periodic boundary conditions. Lattice sites may have a value of 0, 1 or 2, where 0 represents an empty site (or hole), 1 represents a healthy tree (or tree), and 2 represents a burning tree. A forest is a set of contiguous sites (ie., a connected segment) with value 1 or 2. A forest preserve consists of forests separated by gaps (or deserts) consisting of clusters of connected holes. The system evolves in a specified number of time steps, in each of which entire forests of trees can catch fire and burn down and trees can sprout on individual empty sites.

(Note: having all of the trees in a forest burn down in a single time step while trees grows independently of one another provides a separation between the time scales for the processes of deforestation and reforestation).

Algorithm

Statement and Implementation

(1) A forest preserve of length n, consisting of  emtpy sites and tree sites is created using

forestPreserve = Table[Random[Integer], {n}]

All of the sites in the forest preserve are updated in each time step, according to the following sequence of steps:

(2a) Trees catch fire with probability f and empty sites sprout trees with probability, p, using a set of transformation rules to turn 0's into 1's with probability p, and 1's into 2's with probability, f

treeGrowIgnite = forestPreserve /. {0 :> Floor[1 + p - Random[]], 1 :> Floor[2 + f - Random[]]}

(2b) All tree sites adjacent to an ignited tree site (ie., trees in the same forest ) ignite. This accomplished using the repeated application of a set of transformation rules

forestIgnite = treeGrowIgnite //. {{a___, 2, 1, b___} -> {a, 2, 2, b}, 
                                                     {a___, 1, 2, b___} -> {a, 2, 2, b},
                                                     {2, c___, 1} -> {2, c, 2}, 
                                                     {1, c___, 2} -> {2, c, 2}}

(note: The first two transformation rules change sequences in the list having the form .., 1, 2,.. and .., 2, 1, .. to .., 2, 2, .. . and the last two transformation rules implement periodic boundary conditions by treating trees at each end of the system as belonging to the same forest.

(2c) All ignited trees burn down.

forestNew = forestIgnite /. 2 -> 0

The sequence of steps 2a-c can be combined in an anonymous function which can be applied to any forest preserve configuration

pyro = 
             (treeGrowIgnite = # /. {0 :> Floor[1 + p - Random[]], 1 :> Floor[2 + f - Random[]]};
              forestIgnite = treeGrowIgnite //. {{2, c___, 1} -> {2, c, 2},
                                                                     {1, c___, 2} -> {2, c, 2},
                                                                      {a___, 2, 1, b___} -> {a, 2, 2, b},
                                                                      {a___, 1, 2, b___} -> {a, 2, 2, b}};
              forestNew = forestIgnite /. 2 -> 0)&

where # represents the forest preserve configuration.

(3) Step 2 is repeated m times using

NestList[pyro, forestPreserve, m]

These pieces of code can be combined into a program.

Program

smokeyTheBear[n_, p_, f_, m_] :=
   Module[{},
       
       forestPreserve = Table[Random[Integer], {n}];

       pyro = 
               (treeGrowIgnite = #/.{0 :> Floor[1 + p - Random[]], 1 :> Floor[2 + f - Random[]]};
                forestIgnite = TreeGrowIgnite //. {{2, c___, 1} -> {2, c, 2},
                                                                        {1, c___, 2} -> {2, c, 2},
                                                                        {a___, 2, 1, b___} -> {a, 2, 2, b},
                                                                        {a___, 1, 2, b___} -> {a, 2, 2, b}};
                 forestNew = forestIgnite /. 2 -> 0)&;
 
      trees = NestList[pyro, forestPreserve, m] 
        ]

Distribution of Forests

It is useful to see the 'forest through the trees'. To identify the various forests in the forest preserve, a program which labels clusters on a one-dimensional lattice with reflecting boundary conditions can be written consisting of the following steps:

(1) The lattice sites in the lattice, lat, are updated using

clusterID[#, RotateLeft[#]]&[lat]

where

clusterID[1, 0] := i++
clusterID[1, 1] := i
clusterID[0, b_] := 0 
Attributes[clusterID] = Listable;

where the two arguments of clusterID are the value of a cell and the value of its right nn. clusterID places an identifying label, i,  on  the sites in each cluster of contiguous sites. The initial value of i is 1 and increases by 1 for each successive cluster.

(2) After clusterID has been applied to the lattice, all of the clusters have been correctly labelled with one exception: clusters occurring at the extreme ends of the lattice have been identified as different clusters when they should be treated as a single cluster. To correctly identify the sites in these two clusters as belonging to the same cluster, we can use the anonymous function

If[MatchQ[0, First[#] | Last[#]], #, # /. Last[#] -> 1]&

These code fragments can be combined into a program.

clusterLabel1D[lis_List] :=
    Module[{i = 1, clusterID, relabel},

     clusterID[1, 0] := i++;
     clusterID[1, 1] := i;
     clusterID[0, b_] := 0;
     Attributes[clusterID] = Listable;
   
     result = clusterID[#, RotateLeft[#]]&[lis];
     
     If[MatchQ[0, First[#] | Last[#]], #,  # /. Last[#] -> 1]&[result]
   ]

The use of the clusterLabel1D function can be be illustrated with a simple forest preserve.

config = {1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 1};
labelledForests = clusterLabel1D[config]

Having identified the various forests in config, their size distribution can be determined as follows:

The sizes of the various forests in config is given by

forestSizes = Map[Count[labelledForests,#]&,
                  Rest[Union[labelledForests]] ]

The result indicate that the forests in config have 3, 1, and 3 trees, respectively.

(note: While forestSizes gives the sizes of forests in the order in which they occur in the forest preserve, it does not give the sizes of the deserts lying between the forests. This can be determined using

Map[({#, 1})&, config] //.
{u___, {v_, r_}, {v_, s_}, w___} -> {u, {v, r + s}, w} /. {{y_, t_}, z___,{y_, u_}}->{{y, t + u}, z}

which indicates a sequence consisting of a three-tree forest, a two-site desert, a lone tree, a single empty site, a three-tree forest and an empty site).

The distribution of forest sizes is computed using forestSites with the frequency function

frequency[x_List] := Map[{#, Count[x, #]}&, Union[x]]

frequency[forestSizes]

which indicates that config contains one lone tree and two 3-tree forests.

A graph of the forest size distribution can be made using ListPlot on the final list.

Computer Simulation Projects

1. A two-dimensional  version of the forest fire algorithm has recently been presented (see Physica A, vol. 204, pp.212-229 (1994)) which appears to show a variety of interesting behaviors. The model is described as follows:

The sites of a square lattice with periodic boundary conditions are either: occupied by trees (ie., have  a value 1), burning trees (ie., have a value 2) or are empty (ie., have a  value 0). during a time step, all of the sites are updated acccording to the following rules:

(1) a tree becomes a burning tree with probability (1 - g) if at least one nn is burning.
(2) a tree becomes a burning tree with probability  f(1 - g) if no nn is burning.
(3) a burning tree becomes an empty site.
(4) an empty site becomes a tree with probability p.

where there are three probability parameters: f is the lightning probability, g is the immunity, and p is the tree growth probability.

It has been found that for this model:
 (a) when f and g are zero, spiral-shaped fire fronts, reminiscent of the excitable media (see chapter 10), form for small vaules of p.  
 (b) when g = 0, a self-organized critical state (see chapter 8) occurs. 
 (c)  when  f = 0, a percolation-like phase transition (see chapter 5) takes place at a critical value of g which depends on p.

The spirals form when empty areas grow at the expense of forests as fire advances while forests grow into empty areas in the absense of fire, so the fire fronts and tree fronts wind around one another. The SOC state results when tree growth occurs more often than lightning and forests burn down faster than trees grow so that fires of all sizes occur. The percolation-like phase transition happens when there is a zero fire density.

Implement this algorithm in a CA and use it to create snapshots of the forest preserve showing these behaviors.

2. The average size of forests  in the two-dimensional forest preserve is of interest. To determine this quantity, we can first identify the forests in the preserve and then measue their sizes.

The labeling procedure for the two-dimensional lattice having reflecting boundaries can be explained using a simple lattice as an example.

(lat = {{1,0,1,0,0},{0,1,0,1,0},{1,1,0,1,0},{1,0,0,0,1}})//MatrixForm

(1) The lattice sites at the 'northwest' (ie. top left) corners of the clusters are first labelled using

clusterID[RotateRight[lat, {1, 0}], RotateRight[lat, {0, 1}], lat]

where

clusterCornerID[1, 0, 0] := i++;
clusterCornerID[a_ , __] := a;
Attributes[clusterCornerID] = Listable

The three arguments of clusterCornerID are the value of a site, the nn above the site, and the value of the nn to the left of the site. 

Using clusterCornerID with lat, we get

i = 2;
clusterCornerID[1,0, 0]:= i++;
clusterCornerID[a_, __]:= a;
Attributes[clusterCornerID] = Listable;
lat = {{1,0,1,0,0},{0,1,0,1,0},{1,1,0,1,0},{1,0,0,0,1}};
(cornerLabels = clusterCornerID[#, 
                                RotateRight[#,{1,0}],
                                RotateRight[#,{0,1}]
                                ]&[lat])//MatrixForm

note: In performing the above labeling, the initial value of i  was taken to be 2 (and as a result,  cluster corner sites were numbered 2, 3, ...) because the value 1 was used to identify non-corner cluster sites. However, i = 1 can be used in the final clusterLabel2D program).

(2) Having labelled the sites at the corner of clusters, we next label the sites (those with value 1) that lie within clusters and we merge contiguous clusters . Both of these are accomplished with the following rules

reLabel[0, ___] := 0;
reLabel[a_, b_, c_, d_, e_] := Max[a, b, c, d, e];
Attributes[reLabel] = Listable;

where the five arguments of reLabel are: the value of a site, the value of the nn above the site,  the value of the nn to the left of the site, the value of the nn beneath the site and the value of the nn to the right of the site.Using FixedPoint to apply reLabel repeatedly to cornerLabels (which is the result of applying clusterCornerID to lat), until the cluster labels no longer change, we get

reLabel[0, ___] := 0;
reLabel[a_, b_, c_, d_, e_] := Max[a, b, c, d, e];
Attributes[reLabel] = Listable;
FixedPoint[reLabel[#,
                   RotateRight[#,{1, 0}],
                   RotateRight[#,{0, 1}],
                   RotateRight[#,{-1,0}],
                   RotateRight[#,{0,-1}]]&, 
           cornerLabels]//MatrixForm

We can now combine these code fragments into a program.

clusterLabel2D[lat_List]:=
   Module[{i=2},

       clusterCornerID[1,0, 0]:= i++;
       clusterCornerID[a_, __]:= a;
       Attributes[clusterCornerID] = Listable;
    
       cornerLabels = clusterCornerID[#, RotateRight[#,{1,0}], RotateRight[#,{0,1}]]&[lat];

       reLabel[0, ___] := 0;
       reLabel[a_, b_, c_, d_, e_] := Max[a, b, c, d, e];
       Attributes[reLabel] = Listable;
 
       FixedPoint[reLabel[#,  
                                       RotateRight[#,{1, 0}], 
                                       RotateRight[#,{0, 1}], 
                                       RotateRight[#,{-1,0}], 
                                       RotateRight[#,{0,-1}]]&, cornerLabels]
               ]

If desired, it is straightforward (see the relabelrules2 function in the  clusterLabel program in Chapter 5) to relabel the identified clusters so that their numbering is sequential with no gaps.

While the clusterLabel2D program is specifically designed for use with a lattice system having reflecting boundary conditions, it can be used with a lattice having absorbing boundary conditions by 'decorating' the lattice with rows and columns of 0's (the sandpile model program in Chapter  8  shows how to do this using the absorbBC function ).

Compare the relative speeds of the clusterLabel2D program and the clusterLabel program in chapter 5 as function of the lattice size.

General References

Paczuski, Maya and Bak, Per. 1993 "Theory of the one-dimensional forest-fire model", Phys. Rev. E, 48, pp. R3214-R3216.

Bak, Per and Paczuski, Maya. 1993 "Why Nature is Complex", Physics World, pp. 396-403.