NISTIR 7803 

Mathematical Modeling of 

Wildland-Urban Interface 

Fires 

Ronald G. Rehm 

Randall J. McDermott 

. 
 
NISTIR 7803 

Mathematical Modeling of 

Wildland-Urban Interface 

Fires 

Ronald G. Rehm

RGR Consulting, LLC  

Ronald G. Rehm

Randall J. McDermott

Engineering Laboratory

Fire Research Division

June 2011 

U.S. Department of Commerce 

Rebecca M. Blank, Acting Secretary 

National Institute of Standards and Technology 

Patrick D. Gallagher, Under Secretary for Standards and Technology and Director 
 
Disclaimer 

Certain trade names or company products are mentioned in the text to specify adequately 

the experimental procedure and equipment used. In no case does such identification 

imply recommendation or endorsement by the National Institute of Standards and 

Technology, nor does it imply that the equipment is the best available for the purpose. 

Regarding Non-Metric Units: The policy of the National Institute of Standards and 

Technology is to use metric units in all its published materials. To aid the understanding 

of this report, in most cases, measurements are reported in both metric and U.S. 

customary units.


Abstract

This paper presents a simple model to demonstrate the effect on grass-fire propagation of the winds induced by structural fires in a Wildland-Urban Interface (WUI) setting. The model combines an empirical formula for wind-driven

grass-fire spread and a physics-based analytical solution to the Euler equations to determine the ground-level wind

produced by the burning structure. The scaling of the wind is based on the heat release rate (HRR) of the structural

fire as well as other parameters. Also considered are an ambient wind and a topographical wind, assumed to be proportional to the ground slope. Data on grass and structure fires required by the model are discussed. The model can be

presented using either a Lagrangian or an Eulerian description of the fire front. Methods used to solve each description

are presented and compared, with the so-called Level-Set Method (LSM) used in the latter case. The LSM has the

distinct advantage that it can follow front propagation during pinching off and disappearance of a portion of the front

as well as mergers with other fronts, and examples are presented of each type of behavior. Prediction of the fire front

propagation in the Lagrangian description is illustrated by several examples: a front passing a single burning structure

on flat terrain, a front passing a burning structure on a hill, and a front passing several burning structures. The model

predicts that a fire front will be accelerated toward the burning structure upon approach and decelerated after passing

the structure, thereby spending more time near the burning structure than if the structure were absent. The model also

shows, as expected, that the slope of a hill will accelerate or decelerate the front, depending upon the slope magnitude

and sign. Finally, several burning structures multiply the effects of an individual burning structure, causing the front

to linger longer in the vicinity of the structures.

2
Introduction

Over the past few years, Wildland-Urban Interface (WUI) fires have become of much greater concern.

The WUI includes people and property, and, therefore, the costs of damage in the built environment, especially when measured in terms of injury to people, death and property damage, are far higher than in

wildland settings, National Interagency Fire Center (2006). The WUI has been found to be large and growing rapidly in the United States as population expands and housing development in wildland areas continues.

People are building more houses in a wildland setting every year, Environmental Literacy Council (2007).

For example, data from year 2000 show that about 36 % or 42 million homes in the U.S. are in the WUI

and the numbers are growing rapidly, Radeloff et al (2005) and Gustafson et al (2005). Also, as the average

temperature, both within the United States and globally, has increased, earlier snowmelts have taken place,

extending the number of weeks every year that forests are exposed to high temperatures and dangerously

dry conditions, Andrews et al. (2007). Increasingly hot and dry conditions in the United States have led,

in turn, to more acres of wildland fuel consumed by wildfire with greater threat to the WUI, Andrews et al

(2007), National Interagency Fire Center (2006).

The philosophy of management of wildlands in the United States has undergone a major change over the

past few decades, Andrews et al (2007), Berry (2007). Early in the twentieth century, management practice

was to suppress wildfire completely. The Smokey-the-Bear campaign emphasized reporting and suppression

of all fire in wildland settings. This practice, in turn, allowed the buildup of near-ground vegetative material

in forests so that, when a fire did occur, it burned with a much higher intensity than if this material had been

expunged by more frequent, but less intense fires. Today, the philosophy is that fire is a naturally occurring

process in forests and that periodic, low-intensity fires help to keep the forest floor clear of vegetative debris

that can turn low-intensity fires into high-intensity ones. The forests that have been cleared periodically by

small fires are now believed to be healthier and much more able to resist the negative effects of wildfire.

Observations show that WUI fires can behave differently than wildland fires, Murphy et al (2007).

Wildfires are found to spread as a fire line or fire front with distinct features, such as intensity, rate of spread

(ROS), flame height, etc. that depend upon properties of the wildland fuel, the meteorological properties

and topography. WUI fires are produced when wildfires invade a community. For example, Figure 1 shows

a photograph, taken by John Gibbins of the San Diego Post Tribune, of a wildland fire front approaching the

Scripps Ranch residential community during the Cedar Fire in October 2003.

In both wildland and WUI fires, spot fires ignited by brands are a primary fire-spread mechanism. In

WUI fires, brands can ignite a house, producing a vigorous structure fire and substantial additional brand production while the surrounding trees remain untouched. Examples of such behavior, including photographs,

are presented in the USDA report on the Angora Fire that occurred southwest of Lake Tahoe during 24-26

June, 2007, Murphy et al (2007). The brands were estimated to travel up to 1/4 mile downwind in this fire

and were suspected of igniting additional houses.

While the danger and costs of WUI fires have increased dramatically, National Interagency Fire Center

(2006), the tools to address these fires have not been developed, Mell et al (2007). The length and time

scales associated with wildland and WUI fires vary greatly, ranging from the millimeter length scales and

sub-second time scales associated with combustion of fuel elements to tens of kilometers and days associated

with smoke transport and total burned area. This wide variation in scales, coupled with the highly nonlinear

behavior of many of the physical processes, represents the major challenge faced in attempting to model

such fires.

Models that hope to predict WUI rather than wildland fire spread are likely to be more complex because

of the heterogeneity of the fuel. The usual conceptual models for the interaction of a wildfire with structures

3
regard the structures as isolated and surrounded by wildland fuels and imply a density of houses that is so

low that the burning of these houses has no effect on the progression of the wildfire, GAO Report (2005).

In reality, often there are many burning structures in addition to the wildfire itself that can contribute to

the fire spread. For example, in the photograph of Gibbins, Figure 1, imagine that some unburned structures

were to be ignited by brands. Also assume that the fire front were to invade the community. (The Cedar Fire

actually did not invade the community because a road acted as a fire break.) Such a scenario is an example

of the type of WUI fire envisioned for the analysis presented here.

Recently, an important study by Maranghides and Mell (2009) has appeared that emphasizes for the

first time both the discrete nature of the fuels in WUI fires and the utilization of a systematic approach

to the collection of post-fire data in and around a community with the intention of estimating the spatial

and temporal characteristics of a WUI fire. The fire, known as the 2007 California Firestorm, occurred

during October 2007, and was one of the top four fire incidents due to the number of structures destroyed

and acres burned. The community studied was the Trails development of Rancho Bernardo, north of San

Diego, California, which included 274 homes, of which 245 were within the fire perimeter. The two fires

that invaded the community a little over two hours apart were the Witch and Guejito fires, named after the

location of the initiation of each fire, the Witch and the Guejito Creeks respectively. The weather conditions

during this time showed a rapid wind shift and increased wind (from the west to the east with a new sustained

westward wind of 38 km/h and gusts up to 69 km/h) accompanied by reductions in humidity (from 30 to 40

% before the wind shift to 16 % after). The winds contributed significantly to the rapid fire spread.

Many important findings have arisen from this work. The study established a timeline for the fire spread

through the wildland fuels surrounding the community first and then through the community. The Guejito

fire was estimated to approach the community through the wildland fuel at a rate of 9 km/h and then to drop

to a ROS of 0.35 km/h through the community itself. Embers from the approaching wildland fire began to

hit the community an hour before the main fire front, traveling a distance of 9.0 km and igniting three homes

and several patches of ornamental vegetation. The fire spread up to 500 m into the interior of the community.

Although the embers from the wildland fire initiated several ignitions within the community, the majority of

the damage was caused by the fire front arrival. Interventions by firefighters and homeowners were found

to be very significant, saving 60 % of the defended structures. Direct ember ignitions accounted for one

out of every three homes destroyed, while two out of every three destroyed homes were ignited directly

or indirectly by embers. 40 % of the homes on the perimeter of the community were destroyed compared

to 20 % in the interior. The peak ignition rate for structures occurred when the wildland fire reached the

community and was 21 per hour. There were 74 structures destroyed with 29 (or 40 %) estimated to be

burning at the same time.

Currently operational models for wildland fires might be satisfactory for some predictions of large wildfire behavior, where average data on wildland fuels is available from Geographical Information Systems

(GIS) over elements measuring 30 m per side. Likewise, predictions of wildland or WUI fire behavior over

areas measuring a few kilometers per side using field models might soon be possible, provided the detailed

data necessary to make such predictions can be obtained. However, it is very unlikely that the capability

for useful predictions of large-scale WUI fire spread by any model will be available within the next several

years.

The objective of this paper is to present a simple physics-based model that demonstrates the effect on

ground-fire propagation of the winds induced by structural fires. In this proposed model, the fuel system is

regarded as having two components, a surface-fuel portion consisting of grass alone and the structural fuel

portion. Propagation of fires in the surface fuel portion is treated as in the wildfire models discussed below.

The structures on the other hand are treated as discrete fuel elements, which generate three-dimensional

4
buoyant plumes and entrain fresh air to sustain the combustion. These entrainment winds influence the

surface-fire propagation. The approximate model utilizes an empirical formula for wind-driven grass-fire

spread and a physics-based analytical solution to the Euler equations to determine the ground-level wind

produced by the burning structure. This simple model does not resort to a costly Computational Fluid

Dynamics (CFD) simulation. The model is based on research presented in studies reported by Rehm (2006).

In the following section of the paper, we briefly discuss physics-based models of fire behavior, focussing

on those related to structure fires. Next, empirical relations and data needed for the model are given. In particular, an empirical relationship is given which relates the rate of spread (ROS) and the local wind velocity

for Australian grass fires; also estimates are made of the heat release rate (HRR) for structural fires. Then,

the coupled model is presented for the propagation of a grass fire driven by winds produced by the local meteorological conditions, entrainment flow from one or more burning houses, and from topographic features

such as a hill. For the interested reader, details of the mathematical formulation of the model are given in

the Appendix. Fire front propagation for this model is illustrated by examples in the next section. Finally,

discussion and conclusions are presented.

Physics-Based Models of Fire Behavior

Physics-based mathematical models of fire behavior have undergone significant development since the early

70s, Rehm (2006). These models can be divided broadly into two categories: indoor, or enclosure fires, and

outdoor fires. There are significant differences between these two types of fires. For example, indoor fires

are generally limited by the oxygen that can reach the fire through openings in the enclosure. Outdoor fires,

on the other hand, are limited by the supply of wildland fuel. Also, fuel moisture content is a critical factor

controlling ignition and spread of wildfires, whereas fuel moisture content is much less important for the

description of indoor fire behavior.

Here we discuss primarily models for structural fires. This limitation is introduced for two reasons.

First, there are several recent papers that review models of wildfire. For example, a recent review of wildfire

predictive capability by Andrews et al (2007) has appeared in Scientific American. More technical reviews

of wildfire models are given by Mell et al (2007), Pastor et al (2003)and by Perry (1998) among others.

Second, models for structural fires have been developed to deal with fires inside single or multiple

enclosures, as noted above, and not to predict behavior of fully involved structural fires. We are aware

of no data on burning rate or heat release rate (HRR) for any scenario including ignition through burnout

for a whole structure. Yet such information is required to develop a model for WUI fires. Acquisition

of such data is very expensive since the complete burning of buildings is not generally sanctioned by fire

authorities, and if sanctioned, not adequately instrumented to obtain the data required. The ignition location,

the building materials, the interior fire-spread scenario, the prevailing wind conditions, and the response of

the structure to severe burning conditions (such as wall, floor or ceiling penetration and/or collapse) would

all be necessary data, for example, to adequately calibrate and test a model for a house burning down. A

database to test the robustness of any outdoor structural-fire model would require many such full-scale burns.

No such database exists.

Indoor fire models are further subdivided into two categories, zone models and field models. The formulation of both zone and field models start from the basic conservation laws of mass, momentum, energy and

species, together with radiative transport. Zone models have been utilized to study fire dynamics in structures for about four decades. They take advantage of approximate mathematical submodels of the physical

processes that occur in enclosure fires to simplify the conservation laws. As a result, they end up with

nonlinear ordinary differential equations together with complex, nonlinear algebraic relations connecting

5
dependent variables. These simplifications reduce both the data and the computational resources required

to predict the progression of a fire in a structure. The model CFAST, Jones et al (2004), is a recent example

from NIST of this class of models.

Field models begin with a partial differential equations (PDEs) description of the conservation laws and

radiative transport, and attempt to integrate these equations directly using techniques from Computational

Fluid Dynamics (CFD) and other disciplines. During the past three decades, field models have been developed and applied very successfully to fires in structures. Field models require considerably more data and

computational resources to make predictions and can provide substantially more detail about fire behavior

than zone models. There are several examples of field models; a recent example from NIST of this class of

models is the Fire Dynamics Simulator (FDS), McGrattan et al (2004).

Physics-based models of wildland or outdoor fires have generally followed a similar, but somewhat delayed, evolution to that of indoor fires. The simpler models have generally been based on the Rothermel

model, Rothermel (1972), which was developed over three decades ago. Field models for outdoor fires, on

the other hand, have only begun serious development over the past several years, McGrattan et al (1996),

Linn (1997), and Mell et al (2007). In general, field models describing both indoor and outdoor fires have

lagged in development over the simpler models due to the previously prohibitive computer-resource requirements and the corresponding data requirements for such models.

Generally, the simple models of wildland fires treat the fuel as a locally homogeneous surface distribution of mass that varies slowly in the horizontal directions. The three-dimensional elevation of the ground

surface is taken into account in determining the fire progression. However, usually, there is little consideration of the large-scale vertical distribution of the fuel, and the small-scale fuel structure is parameterized.

Heat Release Rates of Grass and Structure Fires

Rehm, et al (2001) reviewed the literature on the potential energy content of various wildland fuels and

compared these numbers with the potential energy content of structures. The purpose of that comparison

was to estimate the density of structures required for the potential energy content to be equal to that of a

particular wildland fuel. That work emphasized the importance of the potential energy content in the burning of structures as part of the overall energy available for combustion. However, that study did not consider

any dynamical processes such as ignition of the fuels, or HRRs required to sustain and propagate the fire.

Furthermore, there was no estimate of the duration and completeness of the combustion processes during a

WUI fire.

This section extends that previous work by considering time-dependent processes for two WUI fuels.

An empirical relation for the propagation speed for grass fires is used, and the burning times or the HRR for

ventilation-controlled structure fires is estimated.

Grass Fires

We utilize an empirical relationship presented in the paper of Cheney et al. (1998) for the rate of spread

(ROS) of the fire line in Australian grass. This ROS, rw, in m/s, is determined as a function of the ambient

wind speed, Va in m/s, the effective length W of the head portion of the fireline in m and the dead fuel

moisture percentage Mfin %:

(1) rw = 0.165(1+3.24 Va) exp[(−0.859−2.04 Va)/W] · exp(−0.108 Mf) [m/s]

For simplicity in the examples below, we have utilized the limit of very small moisture content and very

long fire-lines, so that rw ≈ 0.165(1+3.24 Va).

6
Structure Fires

An estimate of the energy release rate during a house fire in the 1991 Oakland Hills fire was reported

by Trelles and Pagni (1997). According to that estimate, a house burns at a peak rate of 45 MW for 1 hour

(yielding about 160 GJ), and then dies down over another 6 hour period. The die-down of the fire is approximated as two steps, one 10 MW for 3 hours and the last as 5 MW for 3 more hours. The total burn time

is 7 hours, and the total energy released by the house is about 320 GJ. Assuming brush around each house

another 5 MW is released for one hour, yielding an additional 18 GJ. (The additional HRR and total energy

produced by the burning brush around the structure are included in the calculations of Trelles and Pagni

(1997). However, these estimates are unnecessary for our purposes here.) If we include this additional HRR

and total energy, and if the house is assumed to be 15 m by 15 m by 5 m, then we estimate the total potential

fuel loading per unit area to be of order 1.4 GJ/m2, the peak HRR per unit area to be of order 0.2 MW/m2.

For comparison, oil yields a heat release rate per unit area of approximately 2 MW/m2, see for example,

McGrattan et al. (1996) and Baum et al. (1994).

Confirmation of this estimate for the magnitude of the peak HRR for a burning structure can be found in

the chapter on compartment fires in the book on fire behavior by Quintiere (1997). Here, Quintiere describes

the stages of a fire in a compartment and estimates the peak heating rates possible during the latter stages of

a compartment fire when the fire is fully developed and the burning rate is ventilation limited, or restricted

by the amount of air entering the enclosure through the vents. During a compartment fire, the flow in and

around the enclosure is driven by buoyancy, which is generated by the burning taking place both inside and

outside the compartment.

In one example, Quintiere estimates a total HRR of 9 MW for a compartment in which the fuel load

is taken to be proportional to the floor area, that in this case is 12 m2. He points out that, this peak HRR

could increase to over 60 MW if the fuel was proportional not only to the floor area, but to the whole inside

area of the compartment. Furthermore, a multi-room structure, with a fuel loading of the more modest

type, producing 9 MW for each room, could also easily exceed the roughly 50 MW peak HRR estimated

by Trelles and Pagni (1997). Therefore, based on compartment fire analysis, it seems very plausible that

structure fires could have peak HRR reaching several times the estimate of 50 MW, and the duration of these

peak HRR would be measured in tens of minutes to hours.

These estimates say nothing about the fact that the roof of the structure might develop a hole or even

collapse under prolonged vigorous burning. In that case, the fire might then resemble more a burning crib

than an enclosure fire. They also do not address the issue of the effects of winds on peak heat release rate or

burn duration. It seems likely that winds would increase the peak HRR and reduce the burn duration, but the

magnitude of these changes is not known. For our purposes here, the estimates above will be used without

trying to assess these other effects.

Two comparisons can be made between structure fuels and grass fuels using the information given above

and cited in the references. First, the energy density of grasses can be compared with the housing energy

density. Rehm et al. (2002) quotes Albini (1984) who gives a value of 5.5 [m/s] for the ROS and a value

of 1 [MW/m] for the fireline intensity I of a grass fire. Using the burnable energy density of the grass as

eg = I/ROS, we find eg ≈ 0.2 MJ/m2. For Australian grass, Mell (2007) uses ROS = 5.5 m/s and I = 5

MW/m, so that eg ≈ 1 MJ/m2. These estimates can be compared to the corresponding energy density for a

rectangular array of same-size houses of specified housing density, see Rehm (2006). A housing density of

2.5 houses per hectare (one acre lot per house), which can represent a suburban housing density in an affluent

suburb, with an energy of 300 GJ per structure, yields a potential energy density es ≈ 75 MJ/m2. A housing

7
density of 10 houses per hectare (1/4 acre lots), which is a much more typical suburban neighborhood, yields

an energy density es ≈ 300 MJ/m2. For these housing densities, the potential energy of the structures is much

higher than that of the grass.

Next, the heat release rate (HRR) or Q0 for grass line fires can be compared to the HRR for an array

of structures, again see Rehm (2006). If we take the distance between houses Lˆ as a characteristic fireline

length and take 1 acre lots, then, for the structure, Q0 ≈ 50 MW, while for the corresponding grass fire

front, HRR ≈ I ·Lˆ = 64 MW. For 1/4 acre lots, the fireline HRR is about 32 MW. The dynamics of the wind

fields generated by burning fuels under low-wind conditions is determined by the HRR, and this argument

demonstrates that the wildland and structure contributions to the dynamics can be comparable under common conditions.

A Simple Model for WUI Surface-Fire Spread

Plume Model of Baum and McCaffrey

In the paper of Baum and McCaffrey (1989), there is a fundamental analysis of the structure of a plume

and its associated flow field produced by a pool fire in a quiescent atmosphere. The quiescent-atmosphere

assumption implies that the plume remains upright and axially symmetric and requires the heat release rate

of the burning house produces buoyant velocities large compared with the local wind around the structure.

In this paper, an empirical correlation for centerline temperature and velocity is determined from the compilation of data obtained from a large number of pool-fire experiments carried out by many investigators over

a wide range of pool-fire diameters. Based on the buoyant, inviscid equations of motion and this correlation, the analysis obtains the scaling relations for the characteristic length and velocity scales for a pool-fire

plume. Furthermore, a detailed velocity profile is determined from a solution to these equations.

The model of Baum and McCaffrey (1989) is for a single buoyancy-driven plume in an inviscid, quiescent fluid of density ρ0, temperature T0 and pressure P0 at ground level. The magnitude of the heat release

rate of the source is designated as Q0, and the specific heat of the air is denoted as Cp. The model starts with

the equations for mass, momentum and energy, assuming axial symmetry. The velocity field is then decomposed into two components, one arising from the divergence and other from the curl. The divergence of the

velocity results from thermal expansion of the gas, and the curl is the vorticity, and these components can

be related to the plume centerline temperature and velocity correlations. From this analysis, the following

set of scaling relations arise:

D

∗ =



Q0

ρ0CpT0

√g

2/5

V

∗ =

p

gD∗(2)

where D

∗=length scale [m], Q0=heat source [W], ρ0=ambient density [kg/m3

], Cp=specific heat at constant

pressure [J/kg], T0=ambient temperature [degrees K], g=acceleration of gravity [m/s2] and V

∗=velocity scale

[m/s].

Finally, a detailed solution for the velocity field, which is valid both inside and outside the plume, is

found by Baum and McCaffrey. This velocity field at ground level is shown in Figure 2.

Baum and McCaffrey applied their model to the study of winds generated by mass fires. The plume

model was also used by Ohlemiller and Corley (1994) to estimate the thermally-induced winds generated

8
during large-scale mass-fire experiments carried out by Forestry Canada. In both cases, the estimated winds

were found to be consistent with the measured winds.

Similarly, Trelles and Pagni (1997) used the model to estimate the winds generated by multiple burning

houses at several times during the Oakland Hills fire of 1991. These predicted winds were then compared

with measured wind data at the same times, and it was found that significant wind changes occurred, consistent with the model predictions, at nearly the same times. Specifically during the Oakland Hills fire, over

a 15 minute interval, from 11:45 AM to 12 noon on 20 October 1991, the number of houses burning was

found to increase from 38 to 259, producing dramatic changes in the winds consistent with the increased

burning.

The reader should note that the fire front generated by a surface fire itself also induces entrainment

winds, see Dold et al (2006), and Mell et al (2007). However, as in all of the simpler wildfire models to date,

except for that of Dold, the fireline entrainment winds have not been taken into account. This approximation

represents a distinct limitation of the simpler wildfire models and this proposed model. However, the proposed model represents an attempt to couple wildfire and structural fire dynamically in a model that could

make predictions fast enough to be included in current operational models, such as BehavePlus, Carlton et

al (2004) and FARSITE, Finney (2004).

Fire-front Propagation on a Three-dimensional Surface

For the spread of a wildfire, it is usual to consider a fire front of arbitrary shape on a horizontal surface

propagating normal to itself into unburned fuel. Behind the front, the fuel is assumed to be burned, and the

front is taken to be thin relative to other dimensions of the problem. The model for the front propagation can

then be formulated mathematically in two related but different descriptions. One is the so-called Lagrangian

description and the other is an Eulerian description, see for example, Sethian (1999) and Fendell and Wolff

(2001). In the former formulation, the advance of each Lagrangian particle on the front is related to the

empirically determined normal rate of spread (ROS) of a fire at the locally determined wind speed. This

normal ROS is the local front velocity perpendicular to the front. It is the most straightforward description

and requires following only a one-dimensional, time-dependent array of these Lagrangian particles. The latter formulates the problem as a time-dependent, convection-diffusion partial differential equation, for which

the fire front at any time is a curve representing a constant value of a dependent variable of the problem.

According to Sethian (1999), this formulation offers some advantages for following the front progression.

However, a distinct disadvantage of this formulation is that it requires solution of a partial differential equation (PDE) in multiple spatial dimensions and time. In this section, we utilize the Lagrangian description.

The governing equations are the ordinary differential equations (ODEs) describing the propagation of

an element of the fire front along the surface:

d~R

dt

(3) = (U~ ·~n)~n

The equations are given in vector form ~R = x(s,t)~i+y(s,t)~j +z(s,t)~k, where~i,~j,~k are unit vectors in the x,

y, z directions. U~ = Ux

~i+Uy~j +Uz

~k is the ROS vector of the fire front at the location (x, y,z), and nx,ny,nz

are the components of the unit normal to the fire front directed toward the unburnt fuel. s is the arc length

along the curve. The vector (U~ ·~n)~n is the normal ROS. The model predicts the fire front propagation.

Let the height of the surface be specified as a function of the horizontal location. Then, only the first

two of the three vector-component equations need to be solved since the fire-front curve is constrained to the

surface. Therefore, the surface function is used to eliminate the height in the component equations for x and

9
y, yielding two ODEs for x(s,t) and y(s,t). These are solved as outlined below, and the height is determined

at any horizontal location from the equation for the surface.

At each point, the fire front is advanced in the direction normal to the front at a speed determined by the

local ROS for the fire. This ROS depends on the local wind speed normal to the front, which depends on the

ambient velocity, an equivalent topographically induced velocity and the entrainment velocities produced

by all burning structures. The topographically induced wind accounts approximately for the observation

that the ROS increases uphill and decreases downhill due to buoyancy on uneven terrain, Rothermel (1972).

Any local velocity induced by the fire front on itself is not taken into account. The vertical component

of the velocity field is made consistent with the topography, as described in detail in the Appendix. For

computational purposes, the fire front is discretized and then moved incrementally to its new location. We

start with an approximation to the normal ROS, and then numerically solve the governing equations. We use

the Method of Lines (MOL) and a centered difference scheme for the spatial discretization of the fire line.

In the examples, we assume that the fire line initially is a straight line along the x-axis, running between -L

and L (see the schematic diagram in Figure 3). Details of the mathematical equations for following the fire

front are presented in the Appendix.

The model presented here and in earlier papers by the authors, is unique in that it attempts to describe the

behavior of wind-blown WUI fires, not just wildland fires. Although Lagrangian models describing front

propagation in wildland fuels, i.e., wildfires, have been around for a few decades, no one has attempted

to extend the methodology to study WUI fires. In the next section, results from both the Lagrangian and

Eulerian formulation of wind-blown WUI fires are presented and compared, and various observations are

made about these results. The equations for both formulations are presented in Appendices A and B. The

Level-Set methodology (LSM) is described in detail in Appendix C.

Model Results

Fire Front Propagation Past a Burning Structure

Two results using the methodology outlined above for computing a fire front propagating past a single

burning structure are shown in Figure 4. These results were reported earlier by Rehm (2006) and are given

here for comparison. In the plot at the left, the fire front progression is shown for a structure burning at 200

MW intensity with an ambient wind speed of 2 m/s blowing toward the top of the diagram. The structure is

shown as a square 12 m on a side, and the fire front is shown every 25 s starting as a straight line 30 m below

the center of the burning structure. The length of the initial fire line is 2L = 60m. For this case, the entrainment velocity produced by the burning structure accelerates the fire front as it approaches the structure and

retards it after it passes. The fire front effectively spends more time in the vicinity of the burning structure.

In the plot on the right, the conditions are the same except that the fire front starts initially 10 m behind

the structure; i.e., it is assumed that the peak heat release rate for the burning structure is achieved only after

the fire front has passed it. In this case, the fire front is simply retarded by the entrainment flow as it tries to

escape from the vicinity of the burning structure.

Fire Front Propagation over a 2D Hill with a Burning Structure

The next example using the methodology outlined above and presented in detail in the Appendix shows

the computation of a fire front propagating along an elliptical hill in the presence of an ambient wind and

the entrainment wind from a burning structure. The elevation of the hill is assumed to vary in the x- and

10
y-directions with the functional form

Z(x, y) = H

1+ (x/Wx)

2 + (y/Wy)2

(4)

where H is the vertical amplitude of the hill in meters, here taken to be 50 meters, Wx is the length scale

measured in meters in the x-direction and Wy is the length scale measured in meters in the y-direction; Wx

and Wy give the horizontal dimensions of the hill, each assumed to be 100 meters (so the example hill is

circular).

The equations of motion for the fire front then include an ambient wind, here assumed to be 2 m/s in the

y-direction only and an entrainment wind generated by a burning structure having a heat release rate (HRR)

of 100 MW located at x = 0, y = −70, and the effective wind generated by the topography. Figure 5 shows

the result of such a computation for which the initial fire line is located at y = −90, meters from the peak of

the hill and extends initially between x = −50 and x = 50 meters in the cross-stream direction.

Figure 5 shows the progression of the fire at equal intervals of time as it moves up the hill under the

influence of three wind components, the ambient wind, the entrainment wind produced by the burning structure and the effective slope-generated wind. Using this model, the fire front progresses at a uniform speed,

dependent on the ambient wind in the y-direction, in the absence of the hill and the burning structure. With

no ambient wind and no burning structure, the model predicts that the topography-induced wind increases

the ROS of the front up the hill and decreases it down the hill nonuniformly in the x-direction. Also for this

example, comparison with and without the ambient wind shows that the ambient wind increases the ROS

of the front as it propagates over the hill. The entrainment wind from the burning structure accelerates the

front ahead of the structure and decelerates it behind the structure nonuniformly in the x−direction.

Fire Front Propagation over Several Burning Structures

Finally, we turn our attention to a larger area in which there are multiple burning structures at arbitrary,

but specified locations. We wish to see how these burning structures can influence grass fire-front propagation as predicted by this model. We have simulated a fire front propagating on a 250 m by 250 m portion

of Worley, Idaho (USA) located within the Couer d’Alene Indian Tribe Reservation Boundaries. Data sets

that determine the location of non-fuel surfaces, vegetative fuels and structural fuels to 1 meter or 2 meter

resolution for testing NIST modeling efforts have been developed under NIST sponsorship by the Couer

d’Alene Indian Tribe GIS Program, Couer d’Alene (2007). The portion of the reservation used in this study

is shown on the left of Figure 6. The 250 m by 250 m layout with eight structures is shown on the right

of Figure 6. For simplicity, we have assumed that the slightly sloped land is flat and that it is covered with

Australian grass, rather than trees and roads, as shown.

The results of this simulation are shown in Figure 7. We assumed a fire line initially 200 m long in the

x-direction with an ambient wind of 2 m/s in the y-direction and that each house burned with a HRR of 100

MW. We allowed the simulation to run for 700 s of real time, requiring 423 CPU seconds. The time between

each fire front in the simulation plot shown in Figure 7 is 70 s.

In Figure 7, note the rapid acceleration of the fire front toward the burning structures initially and the

very slow progression of the front late in the simulation when the front tries to escape the influence of the

eight burning structures. The net effect is for the fire front to become trapped in the vicinity of the burning

structures. Therefore,it is likely that any unignited fuel in this vicinity would be exposed longer to the fire

front and have an increased likelihood of ignition. As noted above, however, the model does not address

ignition or burnout of the fuels.

11
Results using the Level-Set Method

The speed function  ̃f(V,θ) and the assumption that the front moves normal to itself into the unburned fuel determine how the fire front will propagate, Sethian (1999), Fendell and Wolff (2001),

and Rehm and Mell (2009). A general form for the speed function has been suggested by Fendell

and Wolff (2001). It appears that the speed function can be specified by three measurements: (1)

the head-fire speed (measured in the direction of the wind), (2) the flank speed (measured in the

cross-wind direction), and (3) the tail speed (measured opposite to the wind speed). However, neither the form nor the parameters used in this form are established generally for wildland fires, in

large part because the field measurements are difficult to make. Earlier studies, Rehm (2008) and

Rehm and Mell (2009), used a speed function appropriate for Australian grass fires  ̃r0(1+c f

~V ·~n),

since both the form and the parameters were found by Cheney et al. (1998) to match field measurements reasonably well. Some results reported here use this speed function for the simulations,

while others use a simplified version of the more general speed function of Fendell and colleagues

(2001), (2007).

The two speed functions used in the simulations below are:

 ̃f(V,θ) = r ̃0(1+c f (5) ~V ·~n),

and the simplified speed function of Mallet et al (2007), expressed in our notation as:

 ̃f(V,θ) = r ̃0(1+c f

√

V cosmθ) |θ| ≤ π/2,

 ̃ (6) f(V,θ) = r ̃0(α+ (1−α)sinθ) π/2 < |θ| ≤ π.

Here, m and α are parameters, taken to be m = 1.5 and α = 0.5 for the simulations shown in Figures 8 and 9, presented in the next section comparing the Lagrangian and Eulerian solutions. The

remaining simulations, shown in Figures 10, 11 and 12, and also presented in the next section,

used the speed function given in Eq. (5).

Lagrangian and Eulerian Solutions Compared

Simulations were carried out to compare the Lagrangian and Eulerian methodologies. The initial comparison was for the trivial case of a straight-line front in a wind blowing perpendicular to

the front with periodic boundary conditions in the cross-wind direction. The spread function was

for Australian grass, Eq. (1), with a wind speed of 3 m/s. For both the Lagrangian formulation

using the Method of Lines (MOL) solution technique and for the LSM procedure, the front progressed with a propagation speed found to be within a few percent of that calculated from Eq. (1),

even with coarse meshes such as JJ = 10, or 21 nodes for the MOL solution, and Nx = Ny = 30 for

the LSM. Next, a more challenging example was chosen.

Figures 8 and 9 below show comparisons between the two methods for the simulation of fire

front propagation for a nontrivial case, an initially circular fire front blown by a uniform wind

of 3 m/s in the y-direction. The initial front is taken to have a radius of 10 m, centered at x =

0, y = 50. The speed function given in Eq. (6) with  ̃r0 = 0.165 m/s, c f = 3.24 s/m, m=1.5, α =

0.5, is used. In Figure 8, the front is shown at five equal time intervals 6 s apart. The Eulerian

formulation, using the Level-Set solution method, is shown on the left; at each time, the zero

contour is obtained by interpolation from the Eulerian solution. For comparison, the Lagrangian

12
formulation using the method of lines solution of the same problem is shown on the right. For this

method, successive Lagrangian nodes have simply been joined by straight lines; therefore, the fire

front in the windward direction gradually loses its smooth appearance as time increases. (A higher

order interpolation of the front between Lagrangian nodes would provide a smoother contour with

better visual agreement.)

Each simulation requires the selection of several parameters to carry out the simulation. For

example, the Eulerian formulation requires the simulation domain be discretized in both the x- and

y-directions, and boundary conditions along each of the edges of this domain must be specified.

Therefore, the LSM formulation is two-dimensional and has Nx × Ny computational nodes. For

the case shown in Figures 8 and 9, a domain 100 m on a side (Lx = Ly = 100 m with 0 ≤ y ≤

Ly,−Lx/2 ≤ x ≤ Lx/2) was considered, with the discretization being 30 nodes (Nx = Ny = 30) in

each direction. Therefore the total number of computational nodes for the LSM is 900. On the

other hand, the Lagrangian formulation requires that the initial fire front be discretized along its

circumference, which is one-dimensional and periodic. For the case shown in Figures 8 and 9, 120

nodes were used along the circumference.

In Figure 9, the fire fronts from each simulation are compared directly. In the first frame of

Figure 9, the initial fire fronts are shown together, while the second compares the fronts after 6

s, the third frame after 12 s, etc. (Note that the fire-front contours obtained from each simulation

have not been smoothed.) As a reference, during the 30 s of simulation time, the head velocity is

expected to be 1.09 m/s, for a wind speed of 3 m/s, so that after 30 s, the head of the fire should be at

x = 0, y = 92.7 m. The level-set simulation determines the head location to be at y = 93.9 m while

the Lagrangian, MOL simulation, determines the head location to be at y = 92.0 m, exhibiting that

both simulations seem to do a reasonable job quantitatively of determining the fire head location.

Variation of any of the parameters noted above can cause significant differences in the simulations, making detailed comparisons between simulation results and between simulation efficiencies

difficult. However, several computations with different values of the parameters were performed

for the case of the initially circular front shown in Figures 8 and 9, and generally it was found that

the Lagrangian formulation with the MOL solution required less computational time than the LSM

for similar accuracy. This conclusion is in part due to the fact that the Lagrangian formulation of the

fire front is described by the temporal evolution of a one-dimensional set of nodes along the initial

front, whereas the Level-Set Method requires calculating the time evolution for a scalar function

on a two-dimensional set of nodes. Furthermore, the computational advantage of the Lagrangian

formulation is enhanced because all of the simulations were carried out using the software package

Mathematica, Wolfram (1999), in which the Ordinary Differential Equation (ODE) solver used by

the MOL technique is carefully programmed and compiled and therefore very computationally efficient. However, as shown in the examples below, the LSM can handle a much wider variety of

problems, ones in which fronts merge or disappear, and these more challenging problems cannot

easily be solved using Lagrangian methods.

In the next two figures, these capabilities of the Level-Set methodology are illustrated. The

merger of fire fronts is illustrated in Figures 10 and 11, where fronts are shown spreading and

joining to form an extended front or to form an island or pocket of unburned fuel that is subsequently consumed. Figure 10 shows four frames from a Level-Set simulation of the merger of two

fire fronts, one a straight line progressing in the positive y-direction, and the other generated by a

growing spot fire ahead of the straight-line front. Figure 11 shows four frames from a Level-Set

13
simulation of the merger of three fire fronts, producing an island or a pocket of unburned fuel,

which is subsequently consumed.

Finally, in Figure 12, six frames from a level-set simulation of a straight-line fire front merging

with three point-ignition fronts. In the first two frames, the straight-line front and a single pointignition front are shown, with the line moving upward at a uniform speed and the point ignition

growing with time. In the third frame, a second point ignition appears to the right and down-wind

of the original point ignition. In the fourth frame, the front from the original point ignition is seen

to merge with the line front while the second point ignition front grows with time. In the fifth

frame, a third point ignition is shown to the left and down-stream of the other two ignitions. In

the final frame, both of the first two ignition fronts are seen to merge with the line front, while the

third ignition front grows with time.

Discussion and Conclusions

Our analysis utilizes the plume model of Baum and McCaffrey (1989) to describe the flow field

generated by a single burning house and to estimate the effects of this flow field on the progression of a surface fire. The import of the analysis, we believe, is that it demonstrates with a simple

physics-based model and an inexpensive computational scheme that a house, once ignited, becomes part of the fuel system and affects fire-line progression. It also allows us to investigate the

changes in the fire-line spread as various parameters are changed, such as the number and location

of burning structures.

The time scale associated with burning wildland fuels is measured in tens of seconds to a few

minutes. In contrast, the time scale for a burning structure is measured in tens of minutes to hours,

This disparity in time scale is due, in part, to the fact that the fuel distributions for a structure

and for wildland fuel differ substantially, Rehm et al (2002). For a structure, for example, the

fuel loading per unit area is usually much higher and the footprint much smaller than for a parcel

of wildland fuel. Because of this difference, burning structures and burning wildland fuels will

not couple well in general. For example, it is shown here that burning houses can substantially

influence grass fire propagation through entrainment winds. In contrast, burning grass, excluding

the potential for ignition of the structure (which is not addressed in this model), is not expected to

alter structure burning. (Fires propagating in wildland fuels will also generate entrainment winds

that might influence structure fires; but this model cannot yet address this more difficult issue.)

Because of the disparity in time scales, however, it will be difficult to utilize a field model

to compute multiple burning houses and vegetative (WUI) fires over large areas in any detail due

to constraints on computational resources. An advantage of this model is that it is very fast and

resolves the fire front with a minimum number of nodes, even when several burning structures are

assumed involved. All examples shown have been computed in real time or less in a computational

development environment (using Mathematica, see Wolfram (1999)), orders of magnitude faster

than current field models. Therefore, this methodology for calculating WUI fire spread could

potentially be incorporated into current operational models. Using field models for operational

guidance on large-area WUI fires in the near term is unlikely.

The heat release rate of a structure fire determines the strength of its plume and defines a

characteristic length scale and a characteristic velocity scale for the entrainment of the plume.

This model requires that the plume stand upright, and, therefore, that the ambient wind be much

14
less than the characteristic plume velocity.

The surface level entrainment velocity resulting from the plume model of Baum and McCaffrey

decays slowly at large distances from the burning structure. The mass-fire results of Baum and

McCaffrey illustrate this fact. In the present model, the slow decay can be characterized by an

amplification factor, which is the ratio of the entrainment wind due to multiple burning structures

to the entrainment wind produced by a single burning structure, Rehm (2006). The entrainment

factor is found to grow as the HRR of each house increases, as the number of houses increases,

and as the housing density increases, as discussed by Rehm (2006).

The progression of a grass-fire front (and probably other vegetation fire fronts also) can be

altered substantially by burning structures as shown by the examples presented here. These examples also show that, within the limitations of this simple model, these changes can be tracked in

a computationally inexpensive and robust fashion. Even though the entrainment induced by the

burning grass is not taken into account in this model, several parameters determine the details of

the fire-front propagation, and only a few of these variations have been examined here.

A burning structure on a hill will increase the ROS of the fire front up the hill more quickly as

it approaches the structure. It will also accelerate the front nonuniformly, with greater acceleration

of the front closest to the structure and less so the further the front is from the structure. Likewise,

as the front recedes from the burning structure, it retards the front nonuniformly.

The entrainment caused by multiple burning structures can cause the fire front to spend significantly more time in the vicinity of these structures. The fire front is first accelerated into the

vicinity of the burning structures as it approaches and then is retarded as it tries to escape the vicinity. Because of this additional time lingering near these burning structures, the fire front may tend

to ignite other fuels such as nearby structures which had not been ignited previously. As noted

before, however, this simple model does not include mechanisms for ignition or burnout of the

fuels. Therefore, this observation is suggested, but cannot be confirmed by the model.

Two methods for formulating and computing fire front propagation in a continuous fuel are presented. The Lagrangian formulation with the method of lines solution to the equations is one, and

the Level-Set method using an Eulerian solution procedure is the other. Details of the mathematical

formulations and solution procedures are presented in the Appendices. Examples of solutions obtained by both methods are presented, discussed and compared. While the Lagrangian formulation

and solution seems to be computationally more efficient, the Eulerian method has the powerful

advantage that it allows fronts to merge and to pinch off. Examples of this latter behavior are

presented and discussed.

The ability to qualitatively describe mergers or pinching off of the fire fronts arises naturally

from the Level-Set method and is a great advantage of the method. Although the behavior of

the resulting fronts seems physically realistic, there is no assurance that the behavior is in fact

quantitatively correct, and given the current lack of data on wildland fire-front mergers, there is

little prospect for quantitative evaluation of this model behavior.

On the other hand, in the Lagrangian formulation, an algorithm that specifies either a front

merger or the development of a fuel pocket behind a front under a variety of conditions is a daunting conceptual and programming task. To our knowledge, no physical description of a fire-front

merger exists, nor do measurements of the phenomenon, for example. Furthermore, even if there

were a valid theoretical description of the phenomena, programming the behavior of complex

15
mergings of fronts would be tedious.

Acknowledgements

This research could not have been carried out without our many colleagues in the Fire Research Division of NIST. They have provided several of the good ideas presented here. The long term collaborations

of the senior author with Drs. Howard R. Baum and Kevin B. McGrattan have influenced all of this work.

More recent collaborators, Dr. Ruddy Mell and Alex Maranghides, have guided the direction of the research

relative to the WUI fire studies. Dr. Anthony Hamins has continued to provide encouragement and support

to pursue WUI fire issues. Finally, this research of the senior author was supported in part by NIST Grant

Number 60NANB7D6187, Dr. Ruddy Mell, Grant Monitor.

References

Albini, Frank A. 1984. Wildland Fires, American Scientist 590-597.

Andrews, Patrica, Mark Finney and Mark Fischetti. 2007. Predicting Wildfires, Scientific American

47-55.

Baum, H.R. and B. McCaffrey (1989) Fire Induced Flow Field - Theory and Experiment, Fire Safety

Science - Proceedings of the Second International Symposium, pp. 129-148, Hemisphere Publ.

Baum, H.R., K.B. McGrattan and R.G. Rehm (1994) Simulation of Smoke Plumes from Large Oil

Fires. Twenty-Fifth Symposium (International) on Combustion, The Combustion Institute 25:1463-

1469.

Berry, Alison (2007) Forest Policy Up in Smoke: Fire Suppression in the United States, A PERC

Report, Property and Environmental Research Center, Bozeman, Montana 59718, 1-21.

http://www.perc.org/pdf/Forest/

Carlton, Don, Pat Andrews and Collin Bevins (2004) BehavePlus Fire Modeling System, User’s

Guide, Quick Start Tutorial, USDA Forest Service, Rocky Mountain Research Station, Systems for

Environmental Management, Online Help development, http://fire.org/

Cheney, N.P., J.S. Gould and W.R. Catchpole (1998) Prediction of Fire Spread in Grasslands, International Journal of Wildland Fire, 8, 1-13.

Dold, John, Anna Zinoviev, and Rodney Weber, Nonlocal Flow Effects in Bushfire Spread Rates, V

International Conference on Forest Fire Research, D.X. Viegas (Ed.) 2006, 8 pages.

Environmental Literacy Council, Forest Fires, Environmental Literacy Council, 1625 K Street, NW,

Suite 1020 Washington, DC 20006-3868 http://www.enviroliteracy.org/article.php/46.html.

Fendell, F.E. and M.F. Wolff (2001) Wind-Aided Fire Spread, Chapter 6, pp 171-223, Forest Fires,

Behavior and Ecological Effects Edward A. Johnson and Kiyoko Miyanishi (eds), Academic Press,

San Diego.

Finney, Mark A., “FARSITE: Fire Area Simulator – Model Development and Evaluation,” United

States Department of Agriculture, U.S. Forest Service, Rocky Mountain Research Station, Research

Paper RMRS-RP-4, March 1998, Revised February 2004. http://fire.org/

16
Forney, Glenn P. and Kevin B. McGrattan (2004) Users Guide for Smokeview Version 4 - A Tool for

Visualizing Fire Dynamics Simulation Data, NIST Special Publication 1017.

Gottlieb, S., C.W. Shu, E. Tadmor, “Strong Stability-Preserving High-order Time Discretization Methods,” SIAM Review, Vol. 43, Number 1, pages 89-112, (2001).

Gustafson, Eric. J., Roger B. Hammer, Volker C. Radeloff and Robert S. Potts, “The relationshiip

between environmental amenities and changing human settlement patterns between 1980 and 2000 in

the Midwestern USA,” Landscape Ecology, Vol. 20, 2005, pp 773-789.

GAO Report (2005) TECHNOLOGY ASSESSMENT; Protecting Structures and Improving Communications during Wildland Fires, General Accountability Office Report GAO-05-380.

Jones, Walter W., Richard D. Peacock, Glenn P. Forney, Paul A. Reneke (2004) CFAST Consolidated

Model of Fire Growth and Smoke Transport (Version 6) Technical Reference Guide, NIST Special

Publication 1030. http://fire.nist.gov/bfrlpubs/.

Linn, R.R. (1997) Transport model for prediction of wildfire behavior. Los Alamos National Laboratory Scientific Report LA13334-T.

Mallet, V., D.E. Keyes and F.E. Fendell, (14 October 2007) “Modeling Wildland Fire Propagation

with Level Set Methods,” 22 pp.

Maranghides, Alexander and William Mell, (April 2009) “A Case Study of a Community Affected

by the Witch and Guejito Fires,” NIST Technical Note 1635, National Institute of Standards and

Technology, Gaithersburg, MD, 59 pp.

McGrattan, K.B., H.R. Baum and R.G. Rehm (1996) Numerical Simulation of Smoke Plumes from

Large Oil Fires. Atmospheric Environment 30 4125-4136.

McGrattan, K.B.; S. Hostikka; J.E. Floyd; H.R. Baum; R.G. Rehm (2007) Fire Dynamics Simulator

(Version 5), Technical Reference Guide, NIST Special Publication 1018-5, 100 p., National Institute

of Standards and Technology, Gaithersburg, MD. http://fire.nist.gov/bfrlpubs/.

Mell, William, Mary Ann Jenkins, Jim Gould and Phil Cheney (2007) A Physics Based Approach to

Modeling Grassland Fires, International Journal of Wildland Fire 16,1-22. http://fire.nist.gov/wui.

Mell, William E., Samuel L. Manzello, Alexander Maranghides. David Butry and Ronald G. Rehm

(2010) “The Wildand-Urban-Interface Fire Problem - Current Approaches and Research Needs,” International Journal of Wildland Fire, Vol. 19, No. 2, 238-251.

Murnane, Richard J. (2006) Castrophe Risk Models for Wildfires in the Wildland-Urban Interface:

What Insurers Need Natural Hazards Review 150-156.

Murphy, Kathy, Tim Rich and Tim Sexton (2007) An Assessment of Fuel Treatment Effects on Fire

Behavior, Suppression Effectiveness, Structure Ignition on the Angors Fire USDA Report R5-TP-025

32 pp.

T. Ohlemiller and D. Corley (1994) Heat Release Rate and Induced Wind Field in a Large Scale Fire,

Combustion Science and Technology, 94 315-330.

G.L.W. Perry, (1998) Current approaches to modelling the spread of wildland fire: a review, Progress

in Physical Geography, 22 222-245.

17
E. Pastor, L. Zarate, E. Planas and J. Arnaldos (2003) Mathematical models and calculational systems

for the study of wildland fire behaviour, Progress in Energy and Combustion Science, 29 139-153.

Quintiere James (1997) Principles of Fire Behavior, (Delmar Publishers, an International Thomson

Company, Albany, New York).

Radeloff, Volker C., Roger B. Hammer and Susan I. Stewart, “Rural and Suburban Sprawl in the U.S.

Midwest from 1940 to 2000 and Its Relation to Forest Fragmentation,” Conservation Biology, Vol.

19, No. 3, June 2005, pp 793-805.

Rehm, R.G., A. Hamins, H.R. Baum, K.B. Mcgrattan and D.D. Evans (2002) COMMUNITY-SCALE

FIRE SPREAD, Proceedings of the California’s 2001 Wildfire Conference:10 Years After the 1991

East Bay Hills Fire, 10-12 October 2001, Blonski, K.S., M.E. Morales and T.J. Morales, (editors),

Oakland California Technical Report 35.01.462. Richmond CA; University of California Forest Products Laboratory, pp 126-139; Also, NISTIR 6891, July, 2002.

Rehm, Ronald G. (2006) The Effects of Winds from Burning Structures on Ground-Fire Propagation

at the Wildland-Urban Interface, NIST Report GCR 06-892, pp 1-31.

Rehm, Ronald G. (2008) “The Effects of Winds from Burning Structures on Ground-Fire Propagation

at the Wildland-Urban Interface,” Combustion Theory and Modelling, Volume 12, Issue 3, Pages

477-496.

Rehm, Ronald G. and William (Ruddy) Mell (2009) “A simple model for wind effects of burning

structures and topography on wildland-urban interface surface-fire propagation,” International Journal

of Wildland Fire 18, No. 3, 290-301.

Rehm. Ronald G. and Randall J. McDermott (April 2009), “Fire-Front Propagation Using the Level

Set Method,” NIST Technical Note 1611, National Institute of Standards and Technology, U.S. Department of Commerce,

Rothermel, R.C. (1972) A Mathematical Model for Predicting Fire Spread in Wildland Fuels, USDA

Forest Service Research Paper INT-115.

Rothermel, R.C. (1983) How to Predict the Spread and Intensity of Forest and Range Fires, USDA

Forest Service Technical Report INT-143.

Sethian, J.A. (1999) Level Set Methods and Fast Marching Methods, Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science, (Cambridge University

Press, Cambridge, England).

Trelles, Javier, and Patrick J. Pagni (1997) Fire-Induced Winds in the 20 October 1991 Oakland Hills

Fire, Fire Safety Science – Proceedings of the Fifth International Symposium, Melbourne, Australia,

Yuji Hasemi, Editor. 911-922.

Toro, E.F., “Reimann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction,”

Second Edition, Springer (1999).

Wildland Fire Statistics, Historically Significant Wildland Fires, and Total Wildland Fires and Acres

(1960-2006) (2006), National Interagency Fire Center, 3833 South Development Avenue, Boise,

Idaho 83705.

http://www.nifc.gov/fire info/historical stats.htm

http://www.nifc.gov/fire info/fires acres.htm

18
Wolfram, Stephen, The Mathematica Book, Fourth Edition, Wolfram Media and Cambridge University Press, 1999.

19
Appendix A: Lagrangian Equations

In this Appendix, Lagrangian equations are presented for determining the fire front propagation on a specified surface or topography in the presence of burning structures. The governing equations are the ordinary

differential equations (ODEs) describing the propagation of an element of the fire front along the surface:

d~R

dt

(7) = (U~ ·~n)~n

The equations are given in vector form ~R = x(s,t)~i+y(s,t)~j +z(s,t)~k, where~i,~j,~k are unit vectors in the x,

y, z directions. U~ =Ux

~i+Uy~j+Uz

~k is the rate of spread (ROS) vector of the fire front at the location (x, y,z),

and nx,ny,nz are the components of the unit normal to the fire front directed toward the unburnt fuel. s is the

arc length of the curve.

The linear relation given in Eq.(1) between the ROS and the local wind velocity is assumed. For a line

of sufficient head width W and for very low moisture content of the grass, this reduces to:

Un = ROS0(1+cfVT n)

where Un = U~ ·~n and VT n =~VT ·~n with ROS0 = 0.165 [m/s] and cf = 3.24. Then

d~R

dt

= ROS0(1+cf

~VT ·~n)~n

For notational simplicity, we will not explicitly show the time dependency of these functions as we

obtain the tangent and the normal vectors to the fire front. The tangent to this curve is determined by~τ ≡

d~R

ds

and can be written as:

~τ ≡

d~R

ds

=

dx

ds

~i+

dy

ds

~j +



∂Z

∂x

dx

ds

+

∂Z

∂y

dy

ds

(8) ~k

Define the normal vector to the fire front by the parameter n as follows:

~N ≡

d~R

dn

=

dx

dn

~i+

dy

dn

~j +



∂Z

∂x

dx

dn

+

∂Z

∂y

dy

dn

(9) ~k.

~N is determined by the requirement that the dot product of the normal and tangential vectors be zero:

~τ ·~N =

d~R

ds ·

d~R

dn = 0.

Define the following quantities:

dZ

ds

≡



∂Z

∂x

dx

ds

+

∂Z

∂y

dy

ds

(10)

and



dA

ds 2

≡



dx

ds2

"

1+



∂Z

∂y

2

#

+



dy

ds2

"

1+



∂Z

∂x

2

#

+



dZ

ds 2

"

2+



∂Z

∂x

2

+



∂Z

∂y

2

#

−2

∂Z

∂x

∂Z

∂y

dx

ds

dy

ds

(11)

Then, we can write the unit normal vector~n to the fire front on the prescribed surface as follows:

~n ≡

~N

|~N|

=

1

|

dA

ds |

dy

ds

+

dZ

ds

∂Z

∂y



~i−



dx

ds

+

dZ

ds

∂Z

∂x



~j +



∂Z

∂x

dy

ds

−

∂Z

∂y

dx

ds

~k



(12)

20
At the wildland urban interface, houses as well as wildland fuels can be burning. In this model, we

account for three elementary wind fields, the ambient wind, a topographically-induced wind and the entrainment wind produced by all burning structures. We assume that the total velocity at any location can

then be determined by adding linearly the contributions of all the elementary wind fields. The total velocity

locally is taken to be the sum of these velocity contributions:

VT,x(x, y,z) = Va,x +Vt,x(x, y) +Ve,x(x, y)

VT,y(x, y,z) = Va,y +Vt,y(x, y) +Ve,y(x, y)

(13) VT,z(x, y,z) = Va,z(x, y) +Vt,z(x, y) +Ve,z(x, y)

We take for simplicity the horizontal components of the local ambient wind (Va,x,Va,y) to be uniform in

space. Then, the z-component of the ambient wind is given by Va,z(x, y) = Va,x

∂Z

∂x +Va,y

∂Z

∂y

.

In the absence of any ambient wind and over flat terrain, a fire front is found to propagate with the

uniform rate of spread ROS0, as discussed above. When there are topographical features, Z(x, y), but no

ambient wind, it is found that the ROS of the fire front increases uphill and decreases downhill because of

buoyancy effects. This observed behavior can be treated by defining a topographically-induced horizontal

velocity that is proportional to the gradient of the hill, Rothermel (1972). Here, we take this equivalent

horizontal velocity to be given by the relations Vt,x(x, y) ≡ α

∂Z

∂x

and Vt,y(x, y) ≡ α

∂Z

∂y

, where α is a proportionality constant, Rothermel (1972). The z-component of the topographically-induced velocity is given by

Vt,z(x, y) = Vt,x

∂Z

∂x +Vt,y

∂Z

∂y

.

For a fire front exposed to the velocity field generated by a single burning structure of HRR Q0, the

characteristic length and velocity scales are D

∗

and V

∗

as discussed earlier. Let ~r

0 denote the vector distance

from the center of the structure to the element of the fire front. The velocity at this point will be Ve(r

0

) =

V

∗G(r0/D∗

), where G(r

0/D∗

) is the dimensionless velocity and D

∗

is the length scale defined above, and

the dimensionless vector distance,~r, is~r =~r

0/D

∗

.

The detailed solution for the dimensionless velocity function at ground level, G(r) was obtained analytically by Baum and McCaffrey in terms of special functions. For computational purposes, however, this

solution was replaced in the example calculations presented here by the functional form given below, which

closely approximates the analytical solution:

G(r) = ar +br2/2+cr3/3+dr4/4 0 ≤ r ≤ r0;

G(r) = ar0 +br2

0/2+cr3

0/3+dr4

0/4−

f0(r −r1)

2

2(r1 −r0)

+

f1(r −r0)

2

2(r1 −r0)

+

f0(r1 −r0)

2

r0 ≤ r ≤ r1;

G(r) = ar0 +br2

0/2+cr3

0/3+dr4

0/4+

(f0 + f1)

2(r1 −r0)

r1

r

(1/3)

r1 ≤ r;

where r0 = 0.8,r1 = 1.0, f0 = 0.407199, f1 = 0.045029,a = −2.39441,b = 11.2283, c = −13.6154,d =

4.9468.

Therefore, for a single burning structure at (x = h, y = H), the induced horizontal entrainment velocity

components at any point (x, y) are

Ve,x =

x−h

p

(x−h)

2 + (y−H)2

V

∗G(

q

(x−h)

2 + (y−H)2/D

∗

)

Ve,y =

y−H

p

(x−h)

2 + (y−H)2

V

∗G(

q

(x−h)

2 + (y−H)2/D

∗

(14) )

21
The vertical component of the entrainment velocity is taken to be Ve,z(x, y) = α(

∂Z

∂x

Ve,x(x, y) + ∂Z

∂y

Ve,y(x, y)).

The total entrainment from all of the burning structures is obtained by simple summation, Rehm (2006).

If, for example, the jth structure has an entrainment velocity V

∗

j

and a characteristic length scale D

∗

j

determined by the heat release rate Qj of the burning structure (as described in Rehm (2006)), and if the location

of the structure is given by x = hj, y = Hj, then the total entrainment velocity is given by

Ve,x(x, y) =

J

∑

j=1

x−hj

p

(x−hj)

2 + (y−Hj)2

V

∗

j G

q

(x−hj)

2 + (y−Hj)2/D∗

j



Ve,y(x, y) =

J

∑

j=1

(y−Hj)

p

(x−hj)

2 + (y−Hj)2

V

∗

j G

q

(x−hj)

2 + (y−Hj)2/D∗

j



(15)

where we have assumed that the entrainment velocity from each burning structure is only dependent upon

the vector distance in a horizontal plane between the observational location (x, y) and the location of the

burning structure (hj,Hj).

Only the first two of the three vector-component equations for the fire-front spread need to be solved,

since the front is constrained to the surface z = Z(x, y) (assuming that we can describe the surface explicitly

in this form). The functional form for the surface is used to eliminate z in the component equations for x

and y, yielding two ODEs for x(s,t) and y(s,t). These are solved as described below, and z is determined at

(x, y) from the equation for the surface.

At each point, the fire front is advanced in the direction normal to the front at a speed determined by

the local ROS for the fire. This ROS, in turn, depends on the wind speed at that location. For computational

purposes, the fire front is discretized and then moved incrementally to its new location. We consider the fire

line initially to be a straight line along the x-axis, running between -L and L and divide this interval into

2 I panels each of length δ, where δ = 1/I. We start with an approximation to the normal ROS, and then

numerically solve the governing equations. We use the Method of Lines (MOL) and a centered difference

scheme for the spatial discretization at all interior nodes of the fire front. For the end nodes, we use a

one-sided difference scheme with the neighboring interior node.

22
Appendix B: Eulerian Equations

Here we turn to the development of the equations for the front propagation in Eulerian form. As in the book

of Sethian [5], assume that the front is defined by the contour φ(x, y,t) = 0 of the front function φ. Let Ux,Uy

be the components of the spread velocity normal to the front, U~ = Ux

~ix +Uy

~iy.

If, as above, we take the parameter s to be the distance along the surface φ(x, y,t) = 0 from a specified

point along the front, then



∂φ

∂x



dx

ds

+



∂φ

∂y



dy

ds

= 0(16)

so that dx

ds ≈ −∂φ

∂y

and dy

ds ≈

∂φ

∂x

. Then, the unit tangent vector can be written

~τ =

1

|∇φ|



−

∂φ

∂y

~ix +

∂φ

∂x

~iy



(17)

while the unit normal vector becomes

~n =

1

|∇φ|



∂φ

∂x

~ix +

∂φ

∂y

~iy



(18)

where |∇φ| =

p

(∂φ/∂x)

2 + (∂φ/∂y)2

.

Write out the complete Eulerian equations:

∂φ

∂t

+Ux

∂φ

∂x

+Uy

∂φ

∂y

= 0(19)

where

Ux = r ̃0(1+cf

~V ·~n)nx ≡ f(|~V|,~V ·~n/|~V|)nx,

Uy = r ̃0(1+cf (20) ~V ·~n)ny ≡ f(|~V|,~V ·~n/|~V|)ny.

The spread vector, (Ux,Uy), is normal to the fire front and has a magnitude equal to the “speed function”

defined by Wolff and Fendell [6], For the examples shown below, the speed function is taken to be  ̃r0(1+

cf

~V ·~n). Generally, the speed function depends upon the wind speed |~V| and the angle θ between the normal

to the front and the wind velocity: ~V ·~n/|~V| ≡ cosθ.

23
Appendix C: Level Set Numerical Method

In this Appendix we present our method for solving the level set equation (19). The method is generally

second-order in space and time but uses flux limiting schemes to preserve monotonicity of the scalar field.

The 2D spatial domain is Lx × Ly with uniform grid spacings ∆x and ∆y in the x and y directions,

respectively. The scalar φ(i∆x, j∆y,n∆t) and the speed function (i.e. rate of spread) U~ (i∆x, j∆y,n∆t) are

node centered: i = 0:Nx, j = 0:Ny, n is time step index. No flux, or zero gradient, boundary conditions are

applied along each of the boundaries. The initial condition for φ is arbitrary, except that there must be some

initial curve or curves representing the fire front initially for which φ(x, y,t = 0) = 0.

In general, the speed function depends on the scalar gradient at the node location. However, to preserve

monotonicity of the scalar field, the scalar gradient is obtained from a flux limiter which is based on the speed

function. We avoid this circular problem by computing the speed function using a scalar gradient obtained

from a central difference (no limiter). This speed function is then used in the PDE and also in determining

the limited scalar gradient for the PDE. The basic steps of the numerical procedure are as follows:

Step 1: Given the scalar field at time t

n

, the first step in the numerical procedure is to compute the

node-centered scalar gradient by a central difference,



δφ

δx

n

i, j

=

φ

n

i+1, j −φ

n

i−1, j

2∆x

,(21)



δφ

δy

n

i, j

=

φ

n

i, j+1 −φ

n

i, j−1

2∆y

,(22)

where δ/δx, for example, represents the numerical approximation to the partial derivative. In what follows

we will drop the time stamp until discussion of the time integration scheme.

Step 2: Use the differences obtained in Step 1 to determine the fire front unit normal vector (18) and

compute the speed function via (20).

Step 3: Given the front velocity at each node from Step 2, U~ n

i, j

, determine the monotonicity preserving

scalar gradient for the PDE. The limiter scheme choices in our level set code (in order of increasing accuracy)

are: (1) first-order upwinding, (2) minmod, and (3) Superbee, Toro (1999). As an example of how the flux

limiter is implemented, consider the computation of the limited scalar gradient in the x direction. We imagine

a cell of width ∆x centered at the node location. The gradient is determined from a central difference of the

scalar face values for the cell,

δφ

δx

=

φeast −φwest

∆x

.(23)

The scalar face values are determined from the flux limiter scheme. In this example, let us compute the

value for the east face (i+

1

2

) and assume that Ux > 0. In the limiter scheme, the first step is to compute the

local and upwind data variations,

∆loc = φi+1, j −φi, j,(24)

∆up = φi, j −φi−1, j.(25)

These values are used to form the local data ratio,

r =

∆up

∆loc

.(26)

24
In practice, we set r = 0 if ∆loc = 0. The face value is now simply determined from

φeast = φi, j +

1

2

B(φi+1, j −φi, j),(27)

where B is the flux limiter function evaluated from one of the following, Toro (1999):

B(r) = 0 first-order upwinding,

B(r) = max(0,min(1,r)) minmod,

B(r) = max(0,min(2r,1),min(r,2)) Superbee.

(28)

Note that when B = 1, the scheme recovers second-order central differencing.

Step 4: Time integration. We have now discussed all the details associated with computing the right

hand side (RHS) of the scalar partial differential equation (PDE). To advance the field in time, we use a

second-order Runge-Kutta scheme, Gottlieb et al (2001). This scheme is simply a linear combination of two

Forward Euler steps. Let F(φ) = U~ · ∇φ denote the advective terms of the PDE obtained via Steps 2 and 3

above; ∇ denotes the limited discrete gradient operator. The time integration proceeds as follows:

φ

∗

i, j = φ

n

i, j −∆tFi, j(φ

n

),

φ

n+1

i, j =

1

2

φ

n

i, j +

1

2



φ

∗

i, j −∆tFi, j(φ

∗

)



.(29)

The fire front is obtained from the zero level crossing of the scalar field.

25
Figure Captions

Figure 1. Photograph of a fire front approaching the Scripps Ranch residential community during the Cedar

Fire in October 2003. (Photograph by John Gibbins of the San Diego Post Tribune.)

Figure 2. TOP: Dimensionless entrainment velocity at ground level induced by a structure fire. The velocity

is plotted as a function of the dimensionless radial distance from the structure. BOTTOM: Dimensionless

entrainment velocity vectors at ground level induced by the burning structure. See Equation (2) in the text

for the length scale used for the nondimensionalization.

Figure 3. A schematic diagram showing the initial location, the ambient velocity (Va) and the extent (between the vertical bars in the x-direction) of the horizontal fire line, The burning structure is initially a

vertical (y) distance h from the front.

Figure 4. The fire front at several equal time intervals for a burning structure with an ambient wind of

2 m/s, a 200 MW fire and L = 30 m. Distances are shown in meters. TOP: House is burning before arrival

of fire front. BOTTOM: House ignites and becomes fully involved after passage of fire front.

Figure 5. The fire front at several equal time intervals propagating up a hill under a an ambient wind,

an entrainment wind from the burning structure and the topography-induced wind caused by the hill. Distances are measured in meters.

Figure 6. Diagram showing the location of an area under study in Worley, Idaho, at the left and a small

section of the study area that includes the positions of eight structures. This small section is used as an

example for model fire predictions in the WUI.

Figure 7. The propagation of a fire front through the small area shown in the previous figure, assuming

all of the structures are burning. Distances are measured in meters.

Figure 8. LEFT: The progression of a wind-blown fire front from an initially circular fire is shown. The

model is formulated as a two-dimensional Eulerian partial differential equation (PDE) and solved with a

level-set methodology. For initial conditions, the front is taken to be circular with a radius of 10 m, centered at x = 0, y = 50. The scenario shows what might be expected during the growth of a spot fire when a

steady wind is blowing in the positive y-direction. RIGHT: The same scenario as shown on the left, except

that these results were obtained using the Lagrangian formulation for the front progression and solving the

resulting ordinary differential equations (ODEs) by the method of lines (MOL).

Figure 9. Comparison of the solutions obtained by the two methods, initially and every six seconds.

Figure 10. Four frames from a level-set simulation of the merger of two fire fronts, one a straight-line

front progressing normal to itself in the positive y-direction, and the other generated from a spot fire ahead

of the line fire. The order of the frames is upper left first, upper right second, lower left third and lower right

last.

Figure 11. Four frames from a level-set simulation of the merger of three fire fronts, producing an island or a pocket of unburned fuel, which is subsequently consumed.

Figure 12. Six frames from a level-set simulation of a straight-line fire front merging with three point26
ignition fronts. In the first two frames the straight-line front and a single point-ignition front are shown with

the line moving upward at a uniform speed and the point ignition growing with time. In the third frame, a

second point ignition appears to the right and down-wind of the original point ignition. In the fourth frame,

the front from the original point ignition is seen to merge with the line front while the second point ignition

front grows with time. In the fifth frame, a third point ignition is shown to the left and down-stream of the

other two ignitions. In the final frame, both of the first two ignition fronts are seen to merge with the line

front, while the third ignition front grows with time.

27
Figure 1. Rehm and McDermott

28
Figure 2a, 2b. Rehm and McDermott

0.5 1 1.5 2 2.5 3

r

!0.3

!0.25

!0.2

!0.15

!0.1

!0.05

V

!4 !3 !2 !1 1 2 3 4

X

!4

!3

!2

!1

1

2

3

4

Y

29
Figure 3. Rehm and McDermott

X

Y

Va

h

30
Figure 4a, 4b. Rehm and McDermott

!30 !20 !10 10 20 30

X

10

20

30

40

50

60

Y

200MW

250s

2m!s

!30 !20 !10 10 20 30

X

!10

10

20

30

Y

200MW

250s

2m!s

31
Figure 5. Rehm and McDermott

32
Figure 6. Rehm and McDermott

33
Figure 7. Rehm and McDermott

700s 100MW

2m!s

!100 !50 50 100

X

50

100

150

200

250

Y

34
Figure 8. Rehm and McDermott

-40 -20 0 20 40

0

20

40

60

80

100

X

Y

-40 -20 20 40

X

0

20

40

60

80

100

Y

35
Figure 9. Rehm and McDermott

-40 -20 0 20 40

0

20

40

60

80

100

X

Y

-40 -20 0 20 40

0

20

40

60

80

100

X

Y

-40 -20 0 20 40

0

20

40

60

80

100

X

Y

-40 -20 0 20 40

0

20

40

60

80

100

X

Y

-40 -20 0 20 40

0

20

40

60

80

100

X

Y

-40 -20 0 20 40

0

20

40

60

80

100

X

Y

36
Figure 10. Rehm and McDermott

-40 -20 0 20 40

0

20

40

60

80

100

-40 -20 0 20 40

0

20

40

60

80

100

-40 -20 0 20 40

0

20

40

60

80

100

-40 -20 0 20 40

0

20

40

60

80

100

37
Figure 11. Rehm and McDermott

-40 -20 0 20 40

0

20

40

60

80

100

-40 -20 0 20 40

0

20

40

60

80

100

-40 -20 0 20 40

0

20

40

60

80

100

-40 -20 0 20 40

0

20

40

60

80

100

38
Figure 12. Rehm and McDermott

-40 -20 0 20 40

0

20

40

60

80

100

-40 -20 0 20 40

0

20

40

60

80

100

-40 -20 0 20 40

0

20

40

60

80

100

-40 -20 0 20 40

0

20

40

60

80

100

-40 -20 0 20 40

0

20

40

60

80

100

-40 -20 0 20 40

0

20

40

60

80

100

39