Spreading of Fire – Beau Christ

Today’s question

A fire starts at one tree in the middle of a forest. If the forest is 50% trees, how much of it burns? What if it is 70% trees? We will build a cellular automaton in which each tree only “sees” its four neighbors, then repeat the simulation to find out how tree density changes the outcome.

Reading and code

Read textbook Module 10.3, “Spreading of Fire.” The book lists Module 10.2, “Diffusion,” as a prerequisite; the two grid ideas we need from it (neighborhoods and boundary conditions) are explained on this page. The randomness is the same runif(1) < p test from Modules 9.2 and 9.5. Download the complete R script after coding along.

A model made of cells

Our earlier models tracked a few numbers over time: a population, a temperature, the sizes of S, I, and R. A cellular automaton instead tracks a grid of cells. Each cell holds a state, and every cell follows the same local rule at each time step.

In the fire model, each cell has one of three states:

Value	Name	Meaning
0	`EMPTY`	bare ground, no tree
1	`TREE`	a tree that is not burning
2	`BURNING`	a tree that is on fire

We give the values names so the code says what it means. forest[i, j] == BURNING is easier to read (and check) than forest[i, j] == 2.

EMPTY <- 0
TREE <- 1
BURNING <- 2

Build a forest

Each cell gets a tree with probability prob_tree, the tree density. That is one uniform draw per cell. Before running it: in a 10 × 10 forest with prob_tree = 0.7, about how many trees do you expect?

init_forest <- function(n, prob_tree) {
  forest <- matrix(EMPTY, nrow = n, ncol = n)
  forest[runif(n * n) < prob_tree] <- TREE  # one draw per cell
  forest
}

set.seed(1031)
small_forest <- init_forest(10, 0.7)
small_forest

      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]    0    0    0    1    0    1    1    1    1     1
 [2,]    1    0    0    1    0    1    1    1    1     1
 [3,]    1    1    0    1    0    1    1    1    1     1
 [4,]    0    1    0    0    1    1    0    1    1     1
 [5,]    1    1    0    0    1    1    1    1    0     1
 [6,]    1    1    1    1    1    1    0    1    1     1
 [7,]    1    1    1    1    0    1    1    1    0     1
 [8,]    1    0    0    0    1    1    1    0    0     1
 [9,]    1    0    1    1    1    1    1    1    1     1
[10,]    0    1    1    1    1    1    1    0    0     1

sum(small_forest == TREE)

[1] 72

To see a grid, we color each cell by its value. image() draws a matrix with row 1 at the bottom, so show_forest() flips it so the picture matches the printed matrix.

forest_colors <- c("wheat", "forestgreen", "orangered")  # EMPTY, TREE, BURNING

show_forest <- function(forest, title = "") {
  n <- nrow(forest)
  image(1:n, 1:n, t(forest[n:1, ]), col = forest_colors,
        zlim = c(EMPTY, BURNING), axes = FALSE, asp = 1,
        xlab = "", ylab = "", main = title)
}

show_forest(small_forest, "A 10 × 10 forest, density 0.7")

The updating rule

At each time step, a cell’s next value depends on the cell itself and its von Neumann neighborhood: the cells directly north, east, south, and west. Diagonal cells don’t count.

The book’s rules (pp. 438–439):

An empty cell stays empty.
A burning tree burns down and becomes empty.
A tree next to at least one burning neighbor catches fire, unless it is immune. It is immune with probability prob_immune (wet wood, for example).
A tree with no burning neighbors can still be struck by lightning with probability prob_lightning. It burns only if it is struck and not immune, which has probability prob_lightning * (1 - prob_immune).

spread <- function(site, N, E, S, W, prob_immune, prob_lightning) {
  if (site == EMPTY) {
    return(EMPTY)
  }
  if (site == BURNING) {
    return(EMPTY)                      # a burning tree burns down
  }
  # site is a TREE
  if (N == BURNING || E == BURNING || S == BURNING || W == BURNING) {
    if (runif(1) < prob_immune) {
      return(TREE)                     # fire nearby, but immune this step
    } else {
      return(BURNING)
    }
  }
  if (runif(1) < prob_lightning * (1 - prob_immune)) {
    return(BURNING)                    # struck by lightning and not immune
  }
  TREE
}

Test the rule on cases where you already know the answer. Set both probabilities to 0 so there’s no randomness:

spread(TREE, BURNING, EMPTY, TREE, TREE, 0, 0)  # burning neighbor to the north

[1] 2

spread(TREE, TREE, TREE, TREE, TREE, 0, 0)      # no fire nearby

[1] 1

spread(BURNING, TREE, TREE, TREE, TREE, 0, 0)   # burns down

[1] 0

What about the edges?

A cell on the top row has no north neighbor. The book uses periodic boundary conditions: the grid wraps around, so the top row’s north neighbor is the bottom row and the left column’s west neighbor is the right column. Think of the forest drawn on a donut.

To make that easy, we copy the last row above the first, the first row below the last, and do the same with columns. The result is a grid one cell bigger on every side:

extend_periodic <- function(forest) {
  n <- nrow(forest)
  forest_ns <- rbind(forest[n, ], forest, forest[1, ])    # north and south
  cbind(forest_ns[, n], forest_ns, forest_ns[, 1])        # east and west
}

tiny <- matrix(1:9, nrow = 3, byrow = TRUE)
tiny

     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6
[3,]    7    8    9

extend_periodic(tiny)

     [,1] [,2] [,3] [,4] [,5]
[1,]    9    7    8    9    7
[2,]    3    1    2    3    1
[3,]    6    4    5    6    4
[4,]    9    7    8    9    7
[5,]    3    1    2    3    1

Predict before running: in the extended grid, what are the north and west neighbors of the original cell holding 1? Check your answer against the output.

One time step, then many

To advance the forest one step, extend it, then apply spread() to each original cell. Cell [i, j] of the forest sits at [i + 1, j + 1] in the extended grid. Notice that we read neighbors from the old grid and write results into a new grid. If we updated in place, a tree that just caught fire could spread fire again in the same step.

step_forest <- function(forest, prob_immune, prob_lightning) {
  n <- nrow(forest)
  extended <- extend_periodic(forest)
  new_forest <- forest

  for (i in 1:n) {
    for (j in 1:n) {
      r <- i + 1
      c <- j + 1
      new_forest[i, j] <- spread(extended[r, c],
                                 extended[r - 1, c],  # north
                                 extended[r, c + 1],  # east
                                 extended[r + 1, c],  # south
                                 extended[r, c - 1],  # west
                                 prob_immune, prob_lightning)
    }
  }
  new_forest
}

Now start a 40 × 40 forest with density 0.7, set the center tree on fire, and take 30 steps. We save every grid in a list so we can look back at any time. There is no lightning and no immunity yet, so the only randomness is where the trees are.

set.seed(1032)
n <- 40
forest <- init_forest(n, 0.7)
center <- n / 2
forest[center, center] <- BURNING

number_of_steps <- 30
grids <- list(forest)              # grids[[1]] is time 0

for (t in 1:number_of_steps) {
  forest <- step_forest(forest, prob_immune = 0, prob_lightning = 0)
  grids[[t + 1]] <- forest
}

par(mfrow = c(2, 3), mar = c(1, 1, 2, 1))
for (t in c(0, 4, 8, 12, 20, 30)) {
  show_forest(grids[[t + 1]], paste("t =", t))
}

The fire moves outward as a ragged front, leaving bare ground behind it. Some trees survive because no burning neighbor ever touched them. A diagonal neighbor isn’t enough.

Ask a question one forest cannot answer

How does tree density change the fraction of trees that burn? One forest at one density is just an example, so we need many forests at each density, the same Monte Carlo idea as the random walk.

fraction_burned() builds a forest, lights its center, runs until nothing is burning, and reports the fraction of the original trees that burned. The fire can only ever burn out (burning trees become empty, and nothing regrows), so the while loop always ends.

fraction_burned <- function(n, prob_tree, prob_immune = 0) {
  forest <- init_forest(n, prob_tree)
  center <- ceiling(n / 2)
  forest[center, center] <- BURNING
  trees_at_start <- sum(forest != EMPTY)

  while (any(forest == BURNING)) {
    forest <- step_forest(forest, prob_immune, prob_lightning = 0)
  }
  (trees_at_start - sum(forest == TREE)) / trees_at_start
}

Predict first: sketch the curve you expect for fraction burned vs. tree density from 0.3 to 0.9. A straight line? Something else?

set.seed(1033)
densities <- seq(0.3, 0.9, by = 0.05)
runs_per_density <- 20
mean_burned <- numeric(length(densities))

for (k in seq_along(densities)) {
  results <- numeric(runs_per_density)
  for (run in 1:runs_per_density) {
    results[run] <- fraction_burned(30, densities[k])
  }
  mean_burned[k] <- mean(results)
}

data.frame(density = densities, mean_fraction_burned = round(mean_burned, 3))

   density mean_fraction_burned
1     0.30                0.015
2     0.35                0.024
3     0.40                0.026
4     0.45                0.043
5     0.50                0.108
6     0.55                0.272
7     0.60                0.733
8     0.65                0.906
9     0.70                0.987
10    0.75                0.994
11    0.80                0.998
12    0.85                1.000
13    0.90                1.000

plot(densities, mean_burned, type = "o", pch = 16, col = "orangered",
     ylim = c(0, 1), xlab = "Tree density (prob_tree)",
     ylab = "Mean fraction of trees burned",
     main = "One fire in a 30 × 30 forest (20 runs per density)")

The curve is not a straight line. Below a density of about 0.5, the fire almost always dies out near where it started. Above about 0.65, it almost always sweeps the whole forest. In between, a small change in density makes a large change in the outcome. That is a tipping point (physicists call it a percolation threshold), and it comes from nothing more than each tree checking four neighbors.

Look at the runs near the threshold: would 20 runs be enough there? Why might the spread of results be largest at those densities?

If we have time: change one assumption

Pick one and compare with the curve above, keeping the same seed:

Immunity. Call fraction_burned(30, densities[k], prob_immune = 0.25). Does the tipping point move left or right? Why?
Lightning. In the 40 × 40 simulation, start with no burning tree and use prob_lightning = 0.0005. Where do fires start, and does the forest ever stop burning?
Boundaries. Replace extend_periodic() with a version that surrounds the forest with EMPTY cells, so fire can’t wrap around the edges. Which densities does that change most?

What the cellular automaton adds

Each cell follows a simple local rule, but the pattern of the whole forest (a spreading front, surviving patches, a sharp threshold) is not written anywhere in the code. It emerges from many cells interacting. As with Monte Carlo, one run shows what can happen; repeated runs tell us what usually happens under the model’s assumptions.