# Data-Driven Modeling of Wildfire Spread with Stochastic Cellular Automata
and Latent Spatio-Temporal Dynamics

Nicholas Grieshop111University of Missouri<EMAIL_ADDRESS>Christopher K. Wikle 222University of Missouri

###### Abstract

We propose a Bayesian stochastic cellular automata modeling approach to model
the spread of wildfires with uncertainty quantification. The model considers a
dynamic neighborhood structure that allows neighbor states to inform
transition probabilities in a multistate categorical model. Additional spatial
information is captured by the use of a temporally evolving latent spatio-
temporal dynamic process linked to the original spatial domain by spatial
basis functions. The Bayesian construction allows for uncertainty
quantification associated with each of the predicted fire states. The approach
is applied to a heavily instrumented controlled burn.

Keywords: Cellular automata, spatio-temporal statistics, wildfire modeling

##  1 Introduction

In 2021 there were 59,000 wildfires recorded in the USA, covering an area of
over 7 million acres (Hoover and Hanson,, 2022). Anthropogenic climate change
is lengthening the wildfire season and leading to an increase in wildfire
frequency. The damage caused by these fires and the expense of mitigation
strategies costs billions of dollars annually. Improving models that can
predict the spread of wildfires is important for developing and managing fire
fighting resources and in being able to provide timely warning to those who
may be in danger from an advancing fire.

There are numerous approaches to modeling the spread of wildfires (see the
overview in Sullivan,, 2007). For example, in semi-empirical physical models
such as Rothermel, (1972), the rate of fire spread is parameterized as a
function of the heat flux into a given area and the amount of energy necessary
to cause combustion. The heat flux into a given area is a function of the fuel
source and local factors such as wind speed and elevation. These types of
thermodynamic models are incorporated into wildfire simulators such as FARSITE
(Finney,, 1998) where modifications to the approach of Rothermel are
implemented, including consideration of different treatments of crown versus
surface level fires and the addition of a spotting component. Spotting occurs
when a wildfire jumps and spreads to a nonadjacent area. Although there has
been discussion of how fire spotting can be modeled from physical principles
(e.g., Koo et al.,, 2010), incorporating realistic spotting in wildfire models
remains an important topic of research.

From a mathematical perspective, wildfire propagation can be modeled using the
so-called “level set approach.” In this method, the fire front is the object
of interest, and the propagation of this front over space and time is modeled
by considering the movement dynamics via an implicit level-set function.
Information such as wind speed and elevation, as motivated by a physical
model, can be incorporated into the front spread, as done in Mallet et al.,
(2009) or Alessandri et al., (2020), where the parameters involved in the rate
of spread are estimated. In Muñoz-Esparza et al., (2018), complex fuel sources
are considered in the level set approach. Historically, level-set methods rely
on these types of empirical parameterizations and do not directly use data
from the current fire to inform parameters. A recent exception is the work of
Dabrowski et al., (2022), who assimilate data in real-time into a level-set
model forced by a Rothermel relationship using a Bayesian filtering approach.
In addition, Yoo and Wikle, (2022) demonstrate a hierarchical Bayesian data-
driven approach to learn the spatially-varying propagation speed in the
direction normal to the fire front in a level-set framework.

An alternative approach for fire spread modeling is based on cellular automata
(CA), where the spread is controlled by simple rules. In a CA, the temporal
domain consists of equally spaced time points, and the spatial domain is
divided into discrete cells. CA models have a long history in applied
mathematics and computer science, at least going back to von Neumann, (1966).
In the traditional CA approach, the evolution of the states of each spatial
cell from time $t-1$ to time $t$ is governed by a set of simple rules. A
famous deterministic demonstration is given by Conway’s Game of Life, Gardner,
(1970), where at each time point, the transition rule for a cell considers its
eight neighbors, and a state transition occurs depending on the sum of its
neighbors’ states. From this simple set of rules, complex time evolution
behavior can be observed without additional input from the users - aside from
setting the initial state of the cells. Two challenges in developing of a CA
model are to specify the number of neighbors of a cell at time $t-1$ and to
discover or specify a set of rules to update the state of the cell for time
$t$. In addition, one must decide if the transition rules are best represented
as deterministic or stochastic.

Specific representations of the CA method have been incorporated into a
variety of statistical models. For example, agent-based models in which
individual agents are considered “cells” and operate according to a specific
set of probabilistic rules have been used to model the spread of rabies in
raccoon populations (Hooten and Wikle,, 2010; Wikle and Hooten,, 2015).
Simmonds et al., (2020) present a review of this approach applied to
hydrological data and Banks and Hooten, (2021) present a review with
discussion of formal uncertainty quantification for such models.

The traditional CA approach has been applied to wildfire modeling with
differing definitions of the neighborhood structure and transition rules. For
example, a traditional Moore’s neighborhood (cardinal and ordinal neighbors)
has often been used (e.g., Karafyllidis and Thanailakis,, 1997; Banks and
Hooten,, 2021). Alternatively, Zhang and Wang, (2021) considered a two-tiered
neighborhood structure where local information was captured using a small-
scale neighborhood, and an additional larger-scale neighborhood was used for
long-range effects. The spatial cells need not be defined as regular squares,
as shown in Hernández Encinas et al., (2007) and Johnston et al., (2008). This
approach is further discussed in Liu et al., (2018).

Once the neighborhood structure has been defined, the rules of state
propagation must be specified or learned. In many models, such as those
discussed in Currie et al., (2019), Liu et al., (2018), Karafyllidis and
Thanailakis, (1997), and Lautenberger, (2013), the rules of fire propagation
are specified as simple physical relationships and are not learned.

The method outlined in this manuscript builds upon the work of Currie et al.,
(2019) by utilizing physical principles from semi-empirical models and
environmental covariates when constructing a neighborhood structure.
Importantly, rather than specify these relationships directly, our approach is
data-driven in that both the neighborhood and transition rules are learned. In
particular, we incorporate a novel dynamic neighborhood structure motivated by
the Rothermel, (1972) relationship, in addition to transition rules informed
by a latent low-rank spatio-temporal dynamic process. The model is implemented
in a spatio-temporal Bayesian hierarchical framework (e.g., Cressie and
Wikle,, 2015) to provide formal uncertainty quantification of the predicted
fire spread and covariate effects, analogous to the agent-based methods
discussed in Hooten and Wikle, (2010) and reviewed in Banks and Hooten,
(2021).

Section 2 presents our approach using the CA model with physically motivated
covariates and a stochastic latent process, which is then implemented in a
Bayesian hierarchical model to reflect the uncertainty of data, model, and
parameters. Section 3 gives the results of a simulation, demonstrating the CA
method’s ability to learn the rules of cellular evolution, and in Section 4
the method is applied to a real-world fire. Section 5 discusses further
extensions to this framework and other possibilities for determining
neighborhood structure.

##  2 Bayesian Cellular Automata Model Methodology

This section presents the Bayesian hierarchical cellular automata model for
fire spread as well as a description of the model implementation and model
evaluation metrics.

###  2.1 Stochastic Cellular Automata Model

Consider a fire that can occur within a rectangular spatial domain, $D_{s}$,
that is partitioned with $n$ spatial locations (cells),
$\\{\mathbf{s}_{i}:i=1,\ldots,n\\}$. We assume the fire evolves in discrete
time indexed by $t\in\\{1,2,3,\ldots,T\\}$, and that each cell can take one of
$J$ states; here we restrict our analysis to the case where $J=3$ as described
below, corresponding to unburned, burning, and burned wildfire states. We
denote the state of cell $\mathbf{s}_{i}$ at time $t$ by
$\\{S_{t}(\mathbf{s}_{i})=j:j=1,\ldots,J\\}$. We then assume that these
potential states at each discrete time and space location follow a multinomial
distribution with probabilities $p_{jt}(\mathbf{s}_{i})$. That is, for each
spatial location $\mathbf{s}_{i}$ and time $t$, the state probabilities must
satisfy $\sum_{j=1}^{J}p_{jt}(\mathbf{s}_{i})=1$. The challenge is then to
estimate these spatio-temporal probabilities.

Following the method of ordered categorical responses in Albert and Chib,
(1993), we relate the observed discrete ordered cateogries,
$S_{t}(\mathbf{s}_{i})=j:j=1,\ldots,J$, to an unobserved latent spatio-
temporal process, $Z_{t}(\mathbf{s}_{i})$. Specifically, we relate the
probability of an observation being in state $j$, $p_{jt}(\mathbf{s}_{i})$, to
the latent process $Z_{t}(i)$ with a normal cumulative density function (CDF)
link, $\Phi$, and cutpoints $\lambda_{j}$;
$\sum_{j=1}^{j}p_{jt}(\mathbf{s}_{i})=\Phi(\lambda_{j}-Z_{t}(i))$. For
identifiability reasons, the first cutpoint, $\lambda_{1}$, is set equal to 0,
and the last cutpoint, $\lambda_{J}$ set equal to infinity, with the remaining
cutpoints learned.

In the original Albert and Chib, (1993) framework, the latent process, $Z$, is
assumed to vary only due to covariates. Here, we consider spatio-temporal
covariates that can affect the probability of state transition, but also allow
the state transition probabilities to depend on a latent low-rank spatio-
temporal dynamic process to capture dynamics that are not well-specified by
the neighbor-based covariates (e.g., see Cressie and Wikle,, 2015). In
particular, we specify the latent process by

$\mathbf{Z}_{t}=\mathbf{X}_{t}\bm{\beta}+\mathbf{H}\mathbf{Y}_{t}+\bm{\epsilon}_{t},$
(2.1)

where $\mathbf{Z}_{t}\equiv[Z_{t}(1),\ldots,Z_{t}(n)]^{\prime}$,
$\mathbf{X}_{t}$ is an $n\times p$ matrix of p covariates corresponding to
each spatial location at time $t$ (see below), $\bm{\beta}$ are associated
regression coefficients, $\mathbf{Y}_{t}$ is an $r$-dimensional dynamic
process (where $r<<n$), $\mathbf{H}$ is an $n\times r$ spatial basis function
matrix that maps this low-rank process to the $n$ spatial locations, and
$\bm{\epsilon}_{t}$ is a zero-mean error term constrained to be uncorrelated
with variance 1, $\bm{\epsilon}_{t}\sim N(\bm{0},\mathbf{I})$, as in Albert
and Chib, (1993).

Note that the data augmentation approach considered here is just one method to
model ordinal data (e.g., see Schliep and Hoeting, (2015) for an overview of
additional methodologies for spatial data). For example, the Pólya Gamma
approach of Polson et al., (2013) has proven to be an efficient approach for
Bayesian modeling of categorical data. In our case, the use of the data
augmentation scheme from Albert and Chib, (1993) was motivated by our
consideration of a single unobserved latent process, $\mathbf{Z}_{t}$. As
shown in Appendix D, an implementation of our model from the Pólya Gamma
perspective provides essentially equivalent inference and prediction as the
Albert and Chib, (1993) approach presented here.

####  2.1.1 Latent Spatio-Temporal Dynamics

At the next level of the model hierarchy, we specify the evolution of the
dynamic process $\mathbf{Y}_{t}$ in (2.1) as a vector autoregression

$\mathbf{Y}_{t}=\mathbf{M}\mathbf{Y}_{t-1}+\bm{\eta}_{t},$ (2.2)

where $\bm{\eta}_{t}\sim N(\bm{0},\mathbf{Q})$ and $\mathbf{M}$ is an $r\times
r$ transition operator. As discussed in Wikle et al., (2019), it is important
to allow the transition operator to be unstructured to capture realistic
spatio-temporal dynamics (e.g., transient growth). This is facilitated
computationally by the low-rank basis representation, which is often
considered in spatio-temporal dynamic models given that the underlying
dynamics typically exist on a lower-dimensional manifold than that on which it
is modeled (e.g., Wikle et al.,, 2019).

We discuss implementation issues associated with the choice of covariates and
basis functions in 2.2. In addition, the next level of the model hierarchy
requires specification of prior distributions for {$\gamma_{j}:j=2,\dots J-1$}
parameters, $\bm{\beta}$, $\mathbf{M}$, $\mathbf{Q}$, and $\mathbf{Y}_{0}$.
These are presented in Section 2.2 as well.

###  2.2 Model Implementation

This section discusses implementation choices associated with the local
covariates in (2.1), the choice of spatial basis functions in (2.1), the prior
distributions for model parameters, the MCMC algorithm, and our approach to
model evaluation.

####  2.2.1 Specification of Local Covariates

The latent $\mathbf{Z}_{t}$ process in (2.1) is a function of local time-
varying covariates, $\mathbf{X}_{t}$, whose specification is motivated by
empirical thermodynamic relationships that are sometimes used to characterize
fire spread. Specifically, Rothermel, (1972) proposed perhaps the most-used
thermodynamic relationship for the spread of fires (see also Andrews,, 2018,
for further discussion):

$R=\frac{I_{R}\xi(1+\theta_{w}+\theta_{s})}{\rho_{b}\zeta Q_{ig}},$ (2.3)

where the rate of fire spread, $R$ (ft/min), is defined as a function of
reaction intensity ($I_{R}$), larger-scale propagating flux ratio (i.e., a
long-range wind factor $\theta_{w}$), a slope factor ($\theta_{s}$), and
scaled by properties of the fuel ($\rho_{b},\zeta,Q_{ig}$). Wind speed
influences the size and shape of the neighborhood to match the physical
properties of the Rothermel equation. Rather than use this equation directly,
as it is usually used in fire modeling, we utilize this relationship to modify
the neighborhood structure of our CA model dynamically in time as a function
of wind speed. That is, when the wind speed in either the east-west or north-
south direction is greater than $2m/s$, then the wind expands the neighborhood
in the corresponding direction, as shown in Figure 1. This novel dynamic
neighborhood structure, where the neighborhood can change at each time point,
allows for the model to capture larger-scale dynamics, as in Zhang and Wang,
(2021), but is motivated by well-known physical and thermodynamic
relationships.

Figure 1: Neighborhood structure used for the S5 fire considered in Section 4;
if the wind is from the southwest (sw) direction with x and y wind vectors
$>2$m/s then the neighborhood includes the three cells from the sw. Likewise,
if the wind is from the south (s) with a wind speed in the y direction $>2$
m/s and the x wind vector $\leq 2$m/s, then only an additional cell from the
south is included in the structure.

Thus, the matrix $\mathbf{X}_{t}$, an $n$ by $p=3$ matrix, consists of the
state of the cells in their dynamic neighborhood. The use of three states is
again motivated by the Rothermel equation because the heat flux into a cell is
a function of the fuel property and an unburnt cell is fundamentally different
from a burnt cell due to the change in fuel properties of the cell. Denoting
$\mathcal{N}_{k,t}$ to be the neighborhood structure of cell k at time $t$, we
can then write $X_{k,t,1}=\sum_{i}I(S_{t-1}(s_{i})=1),i\in\mathcal{N}_{k,t}$,
$X_{k,t,2}=\sum_{i}I(S_{t-1}(s_{i})=2),i\in\mathcal{N}_{k,t}$, and
$X_{k,t,3}=\sum_{i}I(S_{t-1}(s_{i})=3),i\in\mathcal{N}_{k,t}$, where
$S_{t-1}(s_{i})$ is the state at location i (which is a neighbor of cell k) at
the previous time.

####  2.2.2 Specification of Spatial Basis Functions

The reduced-rank dynamic process $\mathbf{Y}_{t}$ in (2.1) relies on $r$
spatial basis functions in the $n\times r$ matrix $\mathbf{H}$. There are many
choices for such basis functions (e.g., Wikle et al.,, 2019). Here, we
consider a constructed empirical orthogonal function (EOF) spatial basis. As
discussed in Cressie and Wikle, (2015), EOFs can be good spatial basis
functions for dynamical processes because they provide an optimal (from the
Karhuen-Loève decomposition) low-rank spatial dimension reduction with
realistic covariance structure (note, sensitivity to the choice of basis
functions is briefly discussed in Section 4). The calculation of EOFs requires
time replicates. Although we have time replicates for fire data, if one is
predicting the spread outside the range of the fire front, the basis may not
provide adequate coverage of unburnt areas. We address this with a novel
constructed EOF algorithm.

The constructed EOF approach requires that we initially simulate a plausible
realization of fire spread. Specifically, the final state of the fire at the
last observation time $T$ is used as an input to a simple model. The simple
model, $f_{sim}$, can then be run for the desired number of time steps,
$\tau$, and the state of the fire can be predicted using the simulation. The
simulation model predicts the probability of each state at time $T+1$, $T+2$,
…, $T+\tau$ and the most probable state used to evolve the fire. Once a
collection of predicted states is computed, each cell can randomly be assigned
a temperature value, $\mathbf{Z}_{t}:t=T+1,\ldots$, from that same state using
$g_{sim}$. This approach for constructing EOFs is given in Algorithm 1 and the
details of the specification of $f_{sim}$ and $g_{sim}$ is further discussed
in Appendix B.

Importantly, the simulation model, $f_{sim}$, is not meant to be an ideal
simulator for a particular fire. The use of the simplified simulator is simply
to evolve the fire front sufficiently to allow the EOF basis to cover unburnt
regions in the actual CA model as presented in Section 2.

The choice of EOFs as a basis set is subjective, but it does have some
advantages compared to other bases, such as spatial bisquare bases (Wikle et
al.,, 2019). The primary advantage of using the EOF basis functions is that
they are guaranteed to capture the spatial dependence most efficiently for a
given number of basis functions. Having fewer spatial basis function
coefficients is essential in our application due to the estimation of the
latent dynamical process parameters (i.e., $\mathbf{M},\mathbf{Q}$) given a
relatively few number of time points. See Appendix E for a demonstration of
the sensitivity of our model to the specification of a low-rank bisquare basis
set.

Algorithm 1 Construction of the EOFs

$\mathbf{S}:\text{n spatial locations at T time points}$

$\mathbf{Z}_{t}(s=1),\mathbf{Z}_{t}(s=2),\mathbf{Z}_{t}(s=3)$

while $t\leq\tau$ do

$\mathbf{S}_{T+t}\leftarrow f_{sim}(\mathbf{S}_{T+t-1})$

$\mathbf{Z}_{T+t}\leftarrow g_{sim}(\mathbf{S}_{T+t})$

end while

$\mathbf{H}\leftarrow EOF(\mathbf{Z}_{1,2,...,T+\tau})$

####  2.2.3 Prior Specification and MCMC Algorithm

To complete the Bayesian hierarchical model, non-informative priors are
assigned to the model parameters. The local covariate effects, $\bm{\beta}$,
are assigned independent $N(0,5)$ priors. Elements of the evolution matrix,
$\mathbf{M}$, are assigned independent $N(0,2)$ priors, and the inverse
innovation covariance matrix is assigned an inverse Wishart prior,
$\mathbf{Q}^{-1}\sim W((\mathbf{I}_{r}r)^{-1},\mathbf{I}_{r})$. The cutoffs,
$\lambda_{j}$, are given improper flat priors, as in Albert and Chib, (1993).

Once the local covariates, spatial basis functions, and prior distributions
are chosen, inference and prediction is accomplished through Markov chain
Monte Carlo (MCMC) methods as outlined in Algorithm 2 and presented in detail
in Appendix A.

Algorithm 2 General Sampling algorithm, see Appendix A for full details.

$\mathbf{H}$, $\mathbf{X}_{t}$, $\mathbf{Z}_{t}$

procedure Sampling

while sampling == TRUE do

$\mathbf{Y}_{t}|\cdot\sim N(\cdot)$

$\mathbf{Q}^{-1}|\cdot\sim Wis(\cdot)$

$\bm{\beta}|\cdot\sim N(\cdot)$

$\mathbf{Z}_{t}|\cdot\sim
N(\mathbf{X}_{t}\bm{\beta}+\mathbf{H}\mathbf{Y}_{t},\mathbf{I})$

$\lambda|\cdot\sim
Unif(\max[\max(Z_{i}:S_{i}=j,\lambda_{j-1}],\min[\min(Z_{i}:S_{i}=j+1,\lambda_{j+1}])$

end while

end procedure

####  2.2.4 Model Evaluation

The model is evaluated using two different forecast verification metrics that
are specifically designed for categorical data – the ranked probability score
(RPS) and Gilbert Skill Score (GSS), as discussed in Jolliffe and Stephenson,
(2012).

GSS is a measure based on the $2\times 2$ contingency table of the forecasted
states. Our model has 3 predicted states, but we can compress these into two
classes (“burning” and “not burning”) so that the contingency table can be
rewritten as:

| Reference |
---|---|---
Prediction | Burning | Not Burning
Burning | a | b
Not Burning | c | d
Table 1: $2\times 2$ Contingency table for computing GSS

GSS is then calculated by:

GSS $\displaystyle=\frac{a-\frac{(a+c)(a+b)}{n}}{a+b+c-\frac{(a+c)(a+b)}{n}},$

where the term $(a+c)(a+b)/n$ can be considered to be the expected number of
true predictions due to chance. A larger GSS indicates a better model, with 1
being a perfect prediction. GSS has been used for forecast evaluation in
meteorology, climatology, and related fields (e.g., Schaefer,, 1990;
Stephenson,, 2000).

In addition to being able to predict a future state, our model is able to
forecast the probability of each state. Thus, it is appropriate to consider an
evaluation approach that considers the probability. The RPS (e.g., Epstein,,
1969; Mason,, 2004; Ferro et al.,, 2008) is a classic scoring rule that
considers the probability of each cell being state $j$. The RPS for a single
prediction location is calculated by:

RPS
$\displaystyle=\frac{1}{J-1}\sum_{j=1}^{J}[(\sum_{k=1}^{j}p_{k}-\sum_{k=1}^{j}I_{k})]^{2},$

where J is the number of categories, in this case 3. The RPS for each forecast
location is calculated and the average value across all locations is reported.
A lower RPS indicates a better model with 0 indicating a perfect prediction.

##  3 Simulation Experiment

The purpose of this simulation is to ensure that the ordinal categorical CA
model can retrieve the true probabilities of transitioning. Thus, we consider
a simulated fire evolving in time that is categorized into three states;
unburnt, burning, and burnt. The simulated fire considers a spatial domain of
$20\times 30$ cells with 45-time points - a similar size to the real-world
application presented below. The rule for the fire propagation (i.e., from
state unburnt $\to$ burning) was a function of the number of burning Moore
neighbors, see Table 2 for the transition probabilities based on the number of
burning neighbors. This rule roughly approximates the physical propagation of
a fire front as given by the Rothermel equation (Eq. 2.3) as discussed in
Section 2.1. The Bayesian hierarchical model described in Section 2 was fitted
to this dataset. There was an intermittent wind term in the positive y
direction, so the neighborhood structure was a function of this wind term,
using the neighbors as shown in Figure 1. All state transitions were a
function of this dynamic neighborhood structure, so the latent
$\mathbf{Y}_{t}$ process was omitted. The MCMC was run for 10,000 iterations,
with the first half discarded as burn-in.

The first 40 observations from the simulated dataset were used to train the
model, and the in-sample prediction, as seen in Table 2, shows that the method
can retrieve the true probabilities of state transitions.

Metric | Unburnt | Burning | Burnt
---|---|---|---
Prediction Correct | 15193 (81) | 4035 (226) | 3796 (69)
GSS | 0.975 | 0.897 | 0.923
Table 2: Simulation Results - Prediction Correct is the number of correct
(incorrect) by count (e.g., the model correctly predicted state “unburnt”
15,193 times and predicted state “unburn” 81 times, but the truth was
otherwise). The simulation had an RPS of 0.022 (versus the naïve RPS of 0.747
from equal probabilities for each state)

The overarching goal of this CA model is to learn the probabilities of the
three-state fire spread model and use these learned rules to forecast the fire
spread. Table 3 shows that the CA approach presented here can correctly
recover the probability of state transition based on the local neighborhood.
The model can also be used for forecasting. The next five time points beyond
the training period were used to test forecast performance, and the RPS was
0.100 versus the naïve RPS of 0.579 based on equal probabilities.
Incorporating the CA modeling framework within a Bayesian inference paradigm
allows us to quantify the uncertainty of the transitions between states. The
HPDs in Table3 demonstrate the uncertainty associated with the transition
probabilities; thus, forecasted and predicted states from this model include
uncertainty measures that are unavailable in traditional deterministic CA
models.

# unburnt | # burning | Prob Burning | HPD
---|---|---|---
0 | 8 | 1.0000 | (1,1)
1 | 7 | 1.0000 | (1, 1)
2 | 6 | 0.9999 | (0.9999, 1)
3 | 5 | 0.9994 | (0.9989, 0.9999)
4 | 4 | 0.9641 | 0.9533, 0.98163)
5 | 3 | 0.6368 | (0.6127, 0.7188)
6 | 2 | 0.1357 | (0.1295, 0.1889)
7 | 1 | 0.0054 | (0.0047, 0.0109)
8 | 0 | 3.167e-05 | (1.974e-05, 1.165e-04)
Table 3: The true probability of advancing from state unburnt to burning based
solely on the number of neighbors burning and unburnt. The truth from the
simulation and the recovered learned probabilities from the model are
presented as 95% HPD intervals, with bold intervals indicating those which
contain the true parameter.

##  4 S5 Fire Example

The real-world data for this application comes from the RXCadre series of
experimental burns from Florida, USA. An infrared camera recorded the area’s
temperature, and local weather conditions were recorded (Ottmar et al.,,
2015). The data were originally $240\times 320$ resolution with 105-time
points. For the analysis presented here, the data were averaged onto a 24x32
grid with the mean temperature of the pixels within each larger cell used as
data. These temperatures were then categorized into three states: burnt,
burning, and unburnt. The criteria for the classification were based on the
temperature of the cell, with progression from unburnt to burning after the
temperature crossed a threshold above the background average of 300K. After a
burning cell returned to the average background temperature of 300K, the cell
was considered to have transitioned to a burnt state. Using discrete states
instead of modeling the temperature allows the model considered here to be
applied to a broader range of data sources (e.g., those that just present the
fire front boundary at a given time).

The first 39 observations were used to train the model, and then the trained
model was used to forecast 5 time points forward. The MCMC sampler was run for
10,000 iterations, with the first 5,000 discarded as burn-in. First, a model
using only the covariates $\mathbf{X}_{t}$, with no latent temporal process
was fitted, and the results are shown in Figure 2. The model can capture the
fire’s growth in the positive x and y direction, but the local covariates are
insufficient to model realistic spread. For example, the model predicts the
transition from the burning state to the burnt state to be more rapid than the
observations, as seen in the right side of Figure 2.

The limitation illustrated by only considering the covariates to model fire
spread can be addressed by the addition of the latent dynamical process
$\mathbf{Y}_{t}$ linked by the spatial basis functions. The first 5 EOFs were
constructed following the procedure in Section 2.2.2; Figure 3 shows the
general structure.

Figure 2: S5 fire forecasting using covariates only with no latent temporal
process ($\mathbf{Y}_{t}$). The use of only local covariates is insufficient
to explain the spread of the fire. The left side is the truth at the beginning
and end of the forecast period. The middle panel shows the mean predicted
probability of the cell being in the burning state. The right side of the
figure shows the mean predicted probability for the burnt state. The top of
the panel shows the first forecast step at timepoint 39, and the bottom of the
panel is the five timepoint forecast at timepoint 43. Figure 3: The first 5
EOFs after forecasting the fire using the simulation model. These are the
basis functions used to link the latent temporary evolving $Y$ process back to
the original spatial domain. Figure 4: S5 fire forecasting using the proposed
method with the addition of the temporally evolving latent $\mathbf{Y}_{t}$
dynamic process. The left panel shows the true states, and the middle panel
shows the mean prediction probability of the burning states. The right panel
is the mean probability of burnt state, which is significantly reduced when
compared to Figure 2. Figure 5: HPD intervals for the first 3 time forecasts
for the three possible states. The top row is for the first forecast time
period. The first two panels in each row show the lower and upper values of
the HPD for the unburnt state, the middle two panels are the bounds for the
burning state, and the right two panels are the bounds for the burnt state.
The subsequent rows are for the next two forcast time periods. As the length
of the forecast increases, all of the HPDs become wider, as expected.

As shown in Table 4 and visually in Figure 4, adding the latent dynamical
process $\mathbf{Y}_{t}$ on the expansion coefficients from the first 5 EOFs
leads to an improvement in prediction accuracy. The fire in Figure 4 is
forecast to progress more rapidly in the positive x and y directions,
providing a better match to the truth than the forecast with only local
covariate effects - this is further supported by the model scores in Table 4.
The GSS increases for each category with the largest improvement in capturing
the growth of the fire. The RPS decreased, showing an improved model fit with
the addition of the latent dynamical process. Figure 5 demonstrates the
strength of using a Bayesian hierarchical model. The model is able to provide
bounds on the state estimation for each forecast time. The results in Appendix
C investigate an additional fire with different fire spread characteristics,
but still from the RXCadre experiment. This example also demonstrates that the
model with the latent spatio-temporal dynamics performs better than the
covariate-only model.

Metric | Unburnt | Burning | Burnt
---|---|---|---
Covariates Only |  |  |
Prediction Correct | 1332 (76) | 1074 (309) | 260 (789)
GSS | 0.778 | 0.248 | 0.127
Simulated EOFs |  |  |
Prediction Correct | 1341 (77) | 1821 (480) | 84 (37)
GSS | 0.8111 | 0.70 | 0.185
Table 4: S5 fire results - Prediction Correct is the number of correct
(incorrect) forecast states. The covariate-only model had an RPS of 0.354
(versus the naïve RPS of 0.712), and the latent dynamic model had an RPS of
0.212

##  5 Conclusion

The Bayesian hierarchical cellular automata model presented here offers a few
key benefits. A three-state model is more justifiable from a physical
perspective than a two-state model with only burning and unburned states. This
is because burnt cells would have smaller addition to the heat flux in an
adjacent cell than a model that only allows for the transition from unburnt to
burning. The CA model prsented here could incorporate local covariates such as
(but not limited to) wind speed, by dynamically varying the definition of
neighbors. In addition, the model matrix, $\mathbf{X}_{t}$, can easily be
expanded to include terms measuring the cell’s fuel if the fire spreads over a
non-homogeneous region or include slope properties - as both are terms in the
Rothermel equation. Importantly, adding the unobserved latent dynamic process
$\mathbf{Y}_{t}$ explains some of the dynamics not captured by the local
neighborhood structure, and this latent process is linked using a novel
construction of low-rank EOF basis functions. Finally, the model can easily be
used to forecast a fire’s spread and quantify prediction uncertainty, and this
uncertainty of the prediction can be allowed to propagate throughout the
length of the forecast.

Future extensions of the method could explore different definitions of
neighbors. In the model presented here, the neighborhood structure was novel
in that it was defined based on the Rothermel equation, but data-driven
methods could be used to learn this structure. In addition, adding a fire
spotting process to the model could aid in capturing some of the dynamic fire
front shifts that happen in large-scale real-world fires. Mechanisms to allow
for missing data need to be considered in cases where the data is not reliably
recorded on a regular spatial and temporal scale.

##  Appendix A MCMC Sampler

The full-conditional distributions for our MCMC algorithm are presented below.

###  A.1 Priors

$\displaystyle[\bm{\beta}]$ $\displaystyle\sim
N(\bm{\mu}_{b},\bm{\Sigma}_{b})$ $\displaystyle[\mathbf{Y}_{0}]$
$\displaystyle\sim N(\bm{\mu}_{0},\bm{\Sigma}_{0})$
$\displaystyle[\mathbf{Q}^{-1}]$ $\displaystyle\sim
W((\bm{\nu}_{Q}C_{Q})^{-1},\bm{\nu}_{Q})$ $\displaystyle[\mathbf{m}]$
$\displaystyle\sim N(\bm{\mu}_{m},\bm{\Sigma}_{m})$

Hyperparameters were fixed with $\bm{\mu}_{b}=\bm{\mu}_{m}=\bm{\mu}_{0}=0$,
$\bm{\Sigma}_{b}=\bm{\Sigma}_{m}=diag(2)$, $\nu_{Q}=1$, and
$\bm{\Sigma}_{0}=diag(5)$

###  A.2 Sampler

$\mathbf{Z}_{t}=\mathbf{X}_{t}\bm{\beta}+\mathbf{H}\mathbf{Y}_{t}+\bm{\epsilon}$

  1. 1.

$\mathbf{Y}_{0}\sim N(V_{0}a_{0},V_{0})$

     * •

$V_{0}=(\mathbf{M}^{\prime}\mathbf{Q}^{-1}\mathbf{M}+\Sigma_{0}^{-1})^{-1}$

     * •

$a_{0}=\mathbf{M}^{\prime}\mathbf{Q}^{-1}\mathbf{Y}_{1}+\Sigma_{0}^{-1}\mu_{0}$

  2. 2.

$\mathbf{Y}_{t}\sim N(V_{t}a_{t},V_{t})$

     * •

$V_{t}=(\mathbf{H}^{\prime}\mathbf{R}^{-1}\mathbf{H}+\mathbf{Q}^{-1}+\mathbf{M}^{\prime}\mathbf{Q}^{-1}\mathbf{M})^{-1}$

       * –

Note this is where the computation savings come in. We can invert an r x r
matrix faster than an n x n.

     * •

$a_{t}=\mathbf{H}^{\prime}\mathbf{R}^{-1}(\mathbf{Z}_{t}-\mathbf{X}_{t}\bm{\beta})+\mathbf{Q}^{-1}\mathbf{M}\mathbf{Y}_{t-1}+\mathbf{M}^{\prime}\mathbf{Q}^{-1}\mathbf{Y}_{t+1}$

  3. 3.

$Y_{t=T}\sim N(V_{T}a_{T},V_{T})$

     * •

$V_{T}=(\mathbf{H}^{\prime}\mathbf{R}^{-1}\mathbf{H}+\mathbf{Q}^{-1})^{-1}$

     * •

$a_{T}=\mathbf{H}^{\prime}\mathbf{R}^{-1}(\mathbf{Z}_{t}-\mathbf{X}_{t}\bm{\beta})+\mathbf{Q}^{-1}\mathbf{M}\mathbf{Y}_{T-1}$

  4. 4.

$\mathbf{Q}\sim W(Q_{1},Q_{2})$

     * •

$Q_{1}=(\sum(\mathbf{Y}_{t}-\mathbf{M}\mathbf{Y}_{t-1})(\mathbf{Y}_{t}-\mathbf{M}\mathbf{Y}_{t-1})^{\prime}+\nu_{Q}C_{Q})^{-1}$

     * •

$Q_{2}=\nu_{Q}+T$

  5. 5.

$vec(\mathbf{M})\sim N(V_{m}a_{m},V_{m})$

     * •

$V_{m}=((\mathbf{Y}^{\prime}_{0:T-1}\otimes\mathbf{I}_{n})^{\prime}(\mathbf{I}_{T}\otimes\mathbf{Q})^{-1}(\mathbf{Y}^{\prime}_{0:T-1}\otimes\mathbf{I}_{n})+\Sigma_{m}^{-1})^{-1}$

     * •

$a_{m}=(\mathbf{Y}^{\prime}_{0:T-1}\otimes\mathbf{I}_{n})^{\prime}(\mathbf{I}_{T}\otimes\mathbf{Q})^{-1}vec(\mathbf{Y}_{1:T})+\Sigma_{m}^{-1}\mu_{m}$

  6. 6.

$\beta^{i}\sim N(V_{b}a_{b},V_{b})$

     * •

$V_{b}=(\mathbf{X}^{\prime}\mathbf{I}\mathbf{X}+\Sigma_{b}^{-1})^{-1}$

     * •

$a_{b}=(vec(\mathbf{Z}-\mathbf{H}\mathbf{Y})^{\prime}\mathbf{I}_{nT,nT}\mathbf{X}+\mu_{b}\Sigma_{b}^{-1})^{\prime}$

  7. 7.

$\mathbf{Z}|S\sim N(\mathbf{X}\bm{\beta}+\mathbf{H}\mathbf{Y},\mathbf{I})$

     * •

Truncated on the left and right by $\lambda_{j-1},\lambda_{j}$ for $S_{i}=j$

     * •

Note for identifiability reasons $\lambda_{1}=0$. $\lambda_{0}=-\infty$ and
$\lambda_{J+1}=\infty$

  8. 8.

$\lambda_{i}\sim\text{Unif}(\max[\max(Z_{i}:S_{i}=j,\lambda_{j-1}],\min[\min(Z_{i}:S_{i}=j+1,\lambda_{j+1}])$

##  Appendix B EOF Construction

The method to construct the EOF is as follows. From time $t=1$ to the end of
the training period, time $t=T$, a multinomial model, $f_{sim}$, was fit with
only the cell’s neighbors states as covariates. This model was then used to
forecast the most probable state of the fire up to time $T+\tau$, where $\tau$
is the forecast length. To finalize the construction of the EOF, the
temperature of each cell needs to be known. From time $t=1$ to $t=T$ this is
given from the data, but for time $T+1$ to $T+\tau$ these need to be
simulated. Using the temperature at the training time points, the simulated
temperature was randomly sampled from this distribution, $g_{sim}$. For
example, if a cell was forecasted to be state “burning”, the temperature was
randomly sampled (imputed) from the distribution of all cells that were state
“burning” from time $t=1$ to $t=T$.

##  Appendix C S7 fire

In addition to the S5 fire discussed in Section 4, we consider another fire
from the RXCadre experiment that demonstrated different characteristics than
the S5 fire This fire, S7, had a fire front that was more steadily burning
with the fire lacking some of the large jumps in fire boundary that were
experienced in the S5 fire. In addition, the cells transitioned to the burnt
state more rapidly than in the S5 fire. As can be seen in Figure 6, there was
some exceptionally rapid growth that occurred that the model was unable to
capture. However, as demonstrated in Table 5 the addition of the simulated
EOFs coupled with the dynamic $\mathbf{Y}_{t}$ process was able to outperform
the model with only local covariates. The GSS was higher in all categories
which indicates a better fit and the RPS was lower which supports the same
conclusion.

Figure 6: S7 simulated EOF forecasts Metric | Unburnt | Burning | Burnt
---|---|---|---
Covariates Only |  |  |
Prediction Correct | 2459 (16) | 0 (3) | 751 (311)
GSS | 0.709 | 0.000 | 0.442
Simulated EOF |  |  |
Prediction Correct | 2698 (205) | 60 (34) | 710 (133)
GSS | 0.729 | 0.133 | 0.758
Table 5: S7 fire results - Prediction Correct is the number of correct
(incorrect) by count. The covariate model had an RPS of 0.282 and the
simulated EOF model an RPS of 0.151 (versus the naive RPS of 0.733).

##  Appendix D Pólya Gamma Augmentation

The use of the data augmentation scheme from Albert and Chib, (1993) was
motivated by the use of a single unobserved latent process, $\mathbf{Z}_{t}$.
Other data augmentation strategies for categorical data have been proposed,
such as Polson et al., (2013), which utilize a latent Pólya Gamma process for
each category. Both augmentation schemes produce similar results as
demonstrated by considering the simulation from Section 3. The HPD intervals
for the Pólya Gamma augmentation scheme are shown in Table 6 and are similar
to the results in Table 3.

# unburnt | # burning | Prob Burning | HPD
---|---|---|---
0 | 8 | 1.0000 | (0.958,0.998)
1 | 7 | 1.0000 | (0.975, 0.998)
2 | 6 | 0.9999 | (0.985, 0.998)
3 | 5 | 0.9994 | (0.989, 0.997)
4 | 4 | 0.9641 | 0.957, 0.980)
5 | 3 | 0.6368 | (0.621, 0.753)
6 | 2 | 0.1357 | (0.093, 0.152)
7 | 1 | 0.0054 | (0.006, 0.012)
8 | 0 | 3.167e-05 | (0.000, 0.001)
Table 6: The true probability of advancing from state unburnt to burning based
solely on the number of neighbors burning and unburnt. The truth from the
simulation and the recovered learned probabilities for the Polya Gamma data
augmentation approach presented as 95% HPD intervals, with bold intervals
indicating those which contain the true parameter.

The Pólya Gamma data augmentation scheme produced similar results to the
Albert and Chib method while forecasting for a the same length. The Pólya
Gamma scheme had an RPS of 0.101 compared to the naïve RPS of 0.579 and the
RPS of the proposed method of 0.100 as presented in Section 3.

##  Appendix E Bisquare Basis Comparison

As a sensitivity comparison, we considered bisquare basis functions with
equally spaced knot points in the domain. Due to the equal spacing imposed
constraint, 6-knot locations were chosen, which is close to the number of EOFs
but is a conservative comparison. In this case, the model does not capture the
evolving fire front as well, and the corresponding forecast predictions are
worse when compared to the simulated EOF approach, with an RPS of 0.2502
compared to the method with simulated EOFs that had an RPS of 0.212.

## References

  * Albert and Chib, (1993) Albert, J. H. and Chib, S. (1993).  Bayesian Analysis of Binary and Polychotomous Response Data.  Journal of the American Statistical Association, 88(422):669–679.
  * Alessandri et al., (2020) Alessandri, A., Bagnerini, P., Gaggero, M., and Mantelli, L. (2020).  Parameter Estimation of Fire Propagation Models Using Level Set Methods.  Applied Mathematical Modelling.
  * Andrews, (2018) Andrews, P. L. (2018).  The Rothermel Surface Fire Spread Model and Associated Developments: A Comprehensive Explanation.  Technical report, U.S. Department of Agriculture.
  * Banks and Hooten, (2021) Banks, D. L. and Hooten, M. B. (2021).  Statistical Challenges in Agent-Based Modeling.  American Statistician, 75(3):235–242.
  * Cressie and Wikle, (2015) Cressie, N. and Wikle, C. K. (2015).  Statistics for spatio-temporal data.  John Wiley & Sons.
  * Currie et al., (2019) Currie, M., Speer, K., Hiers, J. K., O’brien, J. J., Goodrick, S., and Quaife, B. (2019).  Pixel-level statistical analyses of prescribed fire spread.  Canadian Journal of Forest Research, 49(1):18–26.
  * Dabrowski et al., (2022) Dabrowski, J. J., Huston, C., Hilton, J., Mangeon, S., and Kuhnert, P. (2022).  Towards data assimilation in level-set wildfire models using bayesian filtering.  arXiv preprint arXiv:2206.08501.
  * Epstein, (1969) Epstein, E. S. (1969).  A scoring system for probability forecasts of ranked categories.  Journal of Applied Meteorology (1962-1982), 8(6):985–987.
  * Ferro et al., (2008) Ferro, C., Richardson, D., and Weigel, A. (2008).  On the effect of ensemble size on the discrete and continuous ranked probability scores.  Meteorological Applications, 15:19 – 24.
  * Finney, (1998) Finney, M. A. (1998).  FARSITE: Fire area simulator-model development and evaluation.  Technical report.
  * Gardner, (1970) Gardner, M. (1970).  Mathematical games.  Scientific American, 223(4):120–123.
  * Hernández Encinas et al., (2007) Hernández Encinas, L., Hoya White, S., Martín del Rey, A., and Rodríguez Sánchez, G. (2007).  Modelling forest fire spread using hexagonal cellular automata.  Applied Mathematical Modelling, 31(6):1213–1227.
  * Hooten and Wikle, (2010) Hooten, M. B. and Wikle, C. K. (2010).  Statistical agent-based models for discrete spatio-temporal systems.  Journal of the American Statistical Association, 105(489):236–248.
  * Hoover and Hanson, (2022) Hoover, K. and Hanson, L. A. (2022).  Wildfire statistics.  Technical report, Congressional Research Service.
  * Johnston et al., (2008) Johnston, P., Kelso, J., and Milne, G. J. (2008).  Efficient simulation of wildfire spread on an irregular grid.  International Journal of Wildland Fire, 17(5):614–627.
  * Jolliffe and Stephenson, (2012) Jolliffe, I. T. and Stephenson, D. B. (2012).  Forecast verification: a practitioner’s guide in atmospheric science.  John Wiley & Sons.
  * Karafyllidis and Thanailakis, (1997) Karafyllidis, I. and Thanailakis, A. (1997).  A model for predicting forest fire spreading using cellular automata.  Ecological Modelling, 99(1):87–97.
  * Koo et al., (2010) Koo, E., Pagni, P. J., Weise, D. R., and Woycheese, J. P. (2010).  Firebrands and spotting ignition in large-scale fires.
  * Lautenberger, (2013) Lautenberger, C. (2013).  Wildland fire modeling with an Eulerian level set method and automated calibration.  Fire Safety Journal, 62(PART C):289–298.
  * Liu et al., (2018) Liu, Y., Liu, H., Zhou, Y., and Sun, C. (2018).  Spread vector induced cellular automata model for real-time crown fire behavior simulation.  Environmental Modelling and Software, 108:14–39.
  * Mallet et al., (2009) Mallet, V., Keyes, D. E., and Fendell, F. E. (2009).  Modeling wildland fire propagation with level set methods.  Computers and Mathematics with Applications, 57(7):1089–1101.
  * Mason, (2004) Mason, S. J. (2004).  On using “climatology” as a reference strategy in the brier and ranked probability skill scores.  Monthly Weather Review, 132(7):1891 – 1895.
  * Muñoz-Esparza et al., (2018) Muñoz-Esparza, D., Kosović, B., Jiménez, P. A., and Coen, J. L. (2018).  An Accurate Fire-Spread Algorithm in the Weather Research and Forecasting Model Using the Level-Set Method.  Journal of Advances in Modeling Earth Systems, 10(4):908–926.
  * Ottmar et al., (2015) Ottmar, R., Hiers, J., Butler, B., Clements, C., Dickinson, M., Hudak, A., O’Brien, J., Potter, B., Rowell, E., Strand, T., and Zajkowski, T. (2015).  Measurements, datasets and preliminary results from the rxcadre project ? 2008, 2011 and 2012.  International Journal of Wildland Fire, 25.
  * Polson et al., (2013) Polson, N. G., Scott, J. G., and Windle, J. (2013).  Bayesian inference for logistic models using pólya–gamma latent variables.  Journal of the American statistical Association, 108(504):1339–1349.
  * Rothermel, (1972) Rothermel, R. C. (1972).  A mathematical model for predicting fire spread in wildland fuels.
  * Schaefer, (1990) Schaefer, J. T. (1990).  The critical success index as an indicator of warning skill.  Weather and Forecasting, 5(4):570 – 575.
  * Schliep and Hoeting, (2015) Schliep, E. M. and Hoeting, J. A. (2015).  Data augmentation and parameter expansion for independent or spatially correlated ordinal data.  Computational Statistics & Data Analysis, 90:1–14.
  * Simmonds et al., (2020) Simmonds, J., Gómez, J. A., and Ledezma, A. (2020).  The role of agent-based modeling and multi-agent systems in flood-based hydrological problems: A brief review.  Journal of Water and Climate Change, 11(4):1580–1602.
  * Stephenson, (2000) Stephenson, D. B. (2000).  Use of the “odds ratio” for diagnosing forecast skill.  Weather and Forecasting, 15(2):221 – 232.
  * Sullivan, (2007) Sullivan, A. L. (2007).  A review of wildland fire spread modelling, 1990-present 2: Empirical and quasi-empirical models.  International Journal of Wildland Fire.
  * von Neumann, (1966) von Neumann, J. (1966).  Theory of self-reproducing automata.  Edited by Arthur W. Burks.
  * Wikle and Hooten, (2015) Wikle, C. K. and Hooten, M. B. (2015).  Hierarchical agent-based spatio-temporal dynamic models for discrete-valued data.  Handbook of Discrete-Valued Time Series, pages 349–366.
  * Wikle et al., (2019) Wikle, C. K., Zammit-Mangion, A., and Cressie, N. (2019).  Spatio-temporal Statistics with R.  Chapman and Hall/CRC.
  * Yoo and Wikle, (2022) Yoo, M. and Wikle, C. K. (2022).  A bayesian spatio-temporal level set dynamic model and application to fire front propagation.
  * Zhang and Wang, (2021) Zhang, B. and Wang, H. (2021).  A new type of dual-scale neighborhood based on vectorization for cellular automata models.  GIScience and Remote Sensing, 58(3):386–404.