NATL INST. OF STAND & TECH

A111D7 IDhnD

MIST

PUBLICATIONS

NIST Technical Note 1611

Fire-Front Propagation Using the

Level Set Method

Ronald G. Rehm

Randall J. McDermott

1 00 VIm7I National Institute of Standards and Technology • U.S. Department of Commerce

.115-753

NIST Technical Note 1611

Fire-Front Propagation Using the

Level Set Method

Ronald G. Rehm

Randall J. McDermott

U.S. Department ofCommerce

Building and Fire Research Laboratory

National Institute ofStandards and Technology

Gaithersburg, MD 20899-8663

March 2009

Issued April 2009

/ W \

SfATES O*

•"

U.S. Department of Commerce

Gary Locke, Secretary

National Institute of Standards and Technology

Patrick D. Gallagher, Deputy Director
Certain commercial entities, equipment, or materials may be identified in this document in order to describe an experimental procedure or concept adequately. Such

identification is not intended to imply recommendation or endorsement by the National Institute of Standards and Technology, nor is it intended to imply that the

entities, materials, or equipment are necessarily the best available for the purpose.

National Institute of Standards and Technology Technical Note 1611

Natl. Inst. Stand. Technol. Tech. Note 1611, 15 pages (April 2009)

CODEN: NSPUE2
Fire-Front Propagation Using the Level Set Method

Ronald G. Rehm and Randall J. McDermott

Building and Fire Research Laboratory

National Institute of Standards and Technology

Gaithersburg, MD 20899

March 25, 2009

Abstract

Propagation of an outdoor fire front in wildland or in a combination of wildland and structural fuels

(the so-called wildland-urban interface or WUI fire), can be modelled as an initial-value problem using

either a Lagrangian or an Eulerian description. The equations associated with each description are presented, and the methods used to solve the equations are discussed. Some comparisons between the two

methods are also made. The emphasis in this report is on the Eulerian equations and on the level-set numerical method. Earlier studies had presented the Lagrangian formulation, and a method-of-lines solution. Advantages of the Eulerian/level-set method are discussed, and several examples that illustrate these advantages are presented.

1 Introduction

In earlier papers, [l]-[4], mathematical models for the propagation of wind-driven fire fronts in a combination of wildland and structural fuels, i.e., models for Wildland-Urban Interface (WUI) fires, have been

presented and discussed. The fire front in the wildland fuel is often treated as a front of zero thickness

and is tracked in a Lagrangian description. The initial front is spatially discretized into a set of nodes

whose time evolution are determined as the solution of a set of ordinary differential equations (ODEs)

by the method-of-lines (MOL). Results for the behavior of the front were reported for wind-blown fires in the presence of single or multiple burning structures on level or uneven terrain, Rehm [2] and Rehm

and Mell [4]. The model presented in these papers, [l]-[4], is unique in that it attempts to describe the behavior

of wind-blown WUI fires, not just wildland fires. Although Lagrangian models for front propagation in wildland fuels, i.e., wildfires, have been around for a few decades, no one has attempted to extend the methodology to WUI fires. In Section 2 of this report, for comparison, we summarize the Lagrangian

formulation and its solution for wind-blown fire-front propagation presented in references [l]-[4]. Then,

in Section 3, we present an Eulerian formulation of fire- front propagation and, in Section 4, a solution to these equations by the numerical technique known as the level-set method.

2 Lagrangian Equations

For simplicity in this section, we consider only fire-front propagation in the presence of burning struc- tures over level ground (no topography). We follow the formulation presented in the Appendix of Rehm

[2]. The governing equations are the ordinary differential equations (ODEs) describing the propagation

1
of an element of the fire front in the horizontal plane:

dR ;=t _w

(1) — = {U-n)n

The equations are given in vector form R = xix +yiy

, where ix, iy

are unit vectors in the x and y directions. U = Uxlx + Uyly is the spread rate (or the rate of spread (ROS) vector in m/s) ofthe fire front at the location

(x,y), and nx ,ny are the components of the unit normal to the fire front directed toward the unburnt fuel. 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, can depend on several variables including the total wind

speed at that location. Let V = Vx ix + Vy

iy be the total wind velocity at a specified height. Assume that the linear relation for Australian grass, Cheney et al [9], given below, is valid for the normal ROS and

the local normal wind velocity:

U„ = r (l + cfV„)

where U„ = U -n and V„ = V n with Pq — 0. 165 m/s and Cf = 3.24s/m. Then

dR

(2) — =f (l + cf

V-n)n

If the fire front curve at any specified time t is described by the vector function (x(s,t),y(s,t)), where

s is a parameter specifying the curve, then, the unit tangent vector can be written as

(3) X S ' (** /3*7 **' dy,\

{Ts" +

|V4>| \ds FsV x ds

where |V(()| = ^(dx/ds)2 + (dy/ds)2. The unit normal vector is

(4) "=iT7li(-^ 1_ (Jtyt

+

^^T

In component form, the equations can be written

dx

Tt

=Unnx

(5) -± = U„ny dy

dt

where

U„ = r f 1 +Cf

3 Eulerian Equations

(-dy/ds)Vx + (dx/ds)Vy ^(dx/ds)2 + (dy/ds) 2

At this point, 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 ty(x,y,t) = of the front function <f>. 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 ty(x,y, t) = from a specified

point along the front, then

S)s+@)S"

so that

jig ~ — gj and £ ~ g| . Then, the unit tangent vector can be written

(7) *=Ti7ir("5^+^ |V(|)| V ^

h+

d^
while the unit normal vector becomes

( 8) " = 7^It(

|V<()| \dx' &+^Jy x dy

where |V<t>| = ^/(d$/dx)2 + {d$/dy)2. Write out the complete Eulerian equations:

(9) -^+ U^+Uy^=0

where

dt +Ux

dx

+Uy

dy

Ux = r (l + cfV-n)nx = f{\?\,V-n/\V\)nx,

(10) Uy = r (l+cf

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

?o(l +c/V n). Generally, the speed function depends upon the wind speed \V\ and the angle 8 between

the normal to the front and the wind velocity: V • n/\V\ = cos0.

4 Level Set Numerical Method

In this section we present our method for solving the level set equation (9). 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 x Ly with uniform grid spacings Ax and Ay in the x and y directions,

respectively. The scalar §(iAx,jAy,nAt) and the speed function (i.e. rate of spread) U(iAx,jAy,nAt) are node centered: i — 0:NX , j = : Ny

, n is time step index. No flux, or zero gradient, boundary conditions

are applied along each of the boundaries. The initial condition for <j) is arbitrary, except that there must

be some initial curve or curves representing the fire front initially for which ty(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 tn

, the first step in the numerical procedure is to compute the node-centered scalar gradient by a central difference,

(11) ,'*Y Ku-ti-u

&x

i 2Ax

„v ,'&>V _^l-«/-l

SyJij 2Ay

where 8/§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 (8) and

compute the speed function via (10). Step 3: Given the front velocity at each node from Step 2, U"-, 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 [10]. 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 Ax centered at the node location. The gradient is determined from a central difference of the scalar face values for the cell,

n

,.

§<!> _ <tW - tywest U ' bx~ Ax

The scalar face values are determined from the flux limiter scheme. In this example, let us compute the value for the east face (/+ \) and assume that Ux > 0. In the limiter scheme, the first step is to compute

the local and upwind data variations,

(14) &loc —

b+ij'-ij.

(15) KP = $ij ~ to-ljThese values are used to form the local data ratio,

(16) r = Kp_

In practice, we set r = if A/oc = 0. The face value is now simply determined from

( 1 7) tyeast = tyj +

2

5(<

1>I'+ ' J ~ M'

where B is the flux limiter function evaluated from one of the following, Toro [10]:

B(r) — first-order upwinding,

(18) B(r) = max(0,min(l,/")) minmod,

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

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 [11]. This scheme is simply a linear combination of two

Forward Euler steps. Let F(ty) = U V<j) denote the advective terms of the PDE obtained via Steps 2 and

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

tfj = tfj-AtFijip),

(19) 4#! - Ulj + Utfj-AtFiM))

The fire front is obtained from the zero level crossing of the scalar field. In Matlab this level set is easily visualized using the contour utility.

5 Initial 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 [5], Fendell and Wolff [6], and Rehm and Mell

[4]. A general form for the speed function has been suggested by Fendell and Wolff [6]. 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 [2] and Rehm and Mell [4], used a speed function appropriate for Australian

grass fires ro(l +c/V n), since both the form and the parameters were found by Cheney et al. [9] to match reasonably well field measurements. 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 [6], [7].
The two speed functions used in the simulations below are:

(20) f(V,Q)=ro(l+cf

V-n),

and the simplified speed function of Mallet et al [7], expressed in our notation as:

f(V,Q) = r (l + cf

y/Vcosn e) |8| < re/2,

(21) f(V,e) = r (a+(l-a)sin6) k/2<\Q\<k.

Here, n and a are parameters, taken to be n = 1.5 and a = 0.5 for the simulations shown in Figures 1 and

2. The remaining simulations, shown in Figures 3, 4 and 5, used the speed function given in Eq. (20). We now discuss some simulations that 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.(20), with a wind speed of 3 m/s. For the Lagrangian formulation

with the MOL solution procedure, the propagation speed was found to agree with the value calculated from Eq. (20) to within one tenth of one percent, whereas the Eulerian/level-set method was found to agree with the value calculated by Eq. (18) to approximately 3.5 per cent (Nx — Ny — 30) and 3 percent

(Nx =Ny = 40).

Figures 1 and 2 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 ydirection. The initial front is taken to have a radius of 10 m, centered at x = 0,y = 50. The speed function

given in (21) with Pq = 0A65m/s,Cf — 3.24s/m,n = 1. 5, a — 0.5, is used. In Figure 1, the front is shown

at five equal time intervals 6s 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 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 second or 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 to be discretized in both the x- and the ydirections, and boundary conditions along each of the edges of this domain must be specified. For the case shown in Figures 1 and 2, a domain 100 m on a side (Lx = Ly = 100 m with < y < Ly

, —Lx/2 <

x < Lx /2) was considered, with the discretization being 30 nodes (Nx = Ny — 30) in each direction. Alternately, the Lagrangian formulation requires that the initial fire front be discretized periodically along

its circumference. For the case shown in Figures 1 and 2, 120 nodes were used along the circumference.

Variation of any of the parameters noted above can cause significant differences in the simulations,

making detailed comparisons between simulation results, as well as simulation efficiencies, difficult. All

simulations were carried out using Mathematica [8]. The simulations for a given resolution using the

level set method generally required more cpu time than did the MOL solutions with comparable resolution. However, the parameters for each type of simulation were not varied systematically to determine

cpu time dependencies. Also, we did not attempt to determine the expected cpu time differences between

the two methodologies.

In Figure 2, the fire fronts from each simulation are compared directly. In the first frame of Figure 2, 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 fronts 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. 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. Both simulations seem to do a reasonable job quantitatively of

determining the fire head location. In the next two figures, capabilities touted by the advocates of the level-set methodology are illus- trated. The merger of two fire fronts is illustrated in Figure 3, and in Figure 4, fronts are shown spreading

to form an island or pocket of unburned fuel, that is subsequently consumed. Figure 3 shows four frames
from a level-set simulation of the merger of two fire fronts, one a straight-line front progressing in the

positive y-direction, and the other generated by a growing spot fire ahead of the straight-line front. Figure

4 shows 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.

The ability to qualitatively describe each of these situations arises naturally from the LSM and is generally regarded as great advantage of the method. Although the behavior of the resulting fronts seems

physically realistic, there is no assurance that this appealing qualitative behavior is in fact quantitatively

correct. Furthermore, there is little prospect for quantitative evaluation of the 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 mergings of fronts would be

tedious.

Finally, in Figure 5, a line fire front is generated from the merger of three initial circular fires that have been allowed to propagate under the influence of a 3 m/s ambient wind in the y-direction. Each

initially circular front grows preferentially in the y-direction (head direction), but also grows in the x- or

flank direction until the three separate fronts merge into a single, scalloped front. As noted before, the difference between the head and flank ROS is dependent upon the speed function used in the simulation.

6 Discussion*

Equation (9) determines the advance with time of the function §(x,y,t), whose level set § = 0, represents

the fire front. Note that the only quantity that has physical meaning is the location of the curve §(x,y,t) =

0. This fact underlies both the attractiveness of the level set approach and the issues discussed below.

In order to solve this equation on a finite domain, boundary conditions must be specified. Moreover,

the solution obtained must depend at least to some degree on the boundary conditions chosen. However,

since the function <)) has no meaning away from the front, the boundary conditions chosen cannot be justified by appeal to physical principles related to fire dynamics. In sufficiently simple situations, solutions obtained using this technique can be compared with those obtained by directly tracking the location of

the front (the "Lagrangian" technique) as discussed above. Agreement between the two approaches is then offered as evidence of the correctness of the level set technique, as done in the previous section of

this report. However, there is no guarantee that boundary conditions that work for a single isolated fire front are necessarily correct in more complex situations. Another uncertainty arises when the solution obtained describes the merging and creation of "islands" of unburned fuel, as shown in the examples above. The fact that solutions of this type can be readily

obtained using level set techniques is often cited as an indication of the superiority of the method. This

may be true, but it rests on the assumption that the merging process, where the direction normal to the front is essentially undefined, is correctly described. Since front tracking models are by definition undefined in such regions, this conclusion can only be justified if an analysis based on a more complete,

physically based model can be analyzed in some asymptotic limit to produce a similar result. This does

not seem to have been done for fire fronts. Such an analysis would doubtless be difficult to perform,

since it is far from clear what more detailed model should be employed as the starting point for any such

study.

At this point in time, a reasonable assessment would be that front tracking and level set approaches

result in two different models. For sufficiently simple scenarios, they are equivalent. However, the level

set approach clearly offers much greater generality than front tracking. In the absence of definitive analytical or experimental data, this would clearly be the preferred approach, once the role of boundary

conditions has been properly addressed. It may well be that by choosing a sufficiently large computational domain, the issue of boundary conditions can be finessed. If so, then the level set model can stand on its own as a simplified physics based approach for fire spread over larger domains.

* This discussion of the level set approach was prepared by Dr. Howard R. Baum.
References

1. R.G. Rehm and D.D. Evans, "Physics - Based Modeling of Wildland - Urban Interface

Fires," in "Remote Sensing and Modeling Applications to Wildland Fires," a volume

Edited by John Qu in Springer's Geosciences Series. To appear 2009.

2. Rehm, Ronald G., "The Effects of Winds from Burning Structures on Ground-Fire Propagation at the Wildland-Urban Interface," Combustion Theory and Modelling, Volume 12, Issue 3, 2008, Pages 477-496.

3. Mell, William E., Samuel L. Manzello, Alexander Maranghides and Ronald G. Rehm,"WildandUrban-Interface Fires: Overview and Research Needs," International Journal of Wildland

Fire. Invited paper, under review.

4. Rehm, Ronald G. and William (Ruddy) Mell, "A Simple Model for Wind Effects of Burning Structures and Topography on WUI Surface-Fire Propagation," accepted for publication in the International Journal of Wildland Fire.

5. Sethian, J. A., Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Dynamics, Computer Vision, and Material Science, Cambridge Monographs on Applied and Computational Mathematics, Cambridge University

Press, 1999.

6. Fendell, F.E. and M.F. Wolff, (200 l)"Wind-Aided Fire Spread, Forest Fires, Behavior and

Ecological Effects," Edward A. Johnson and Kiyoko Miyanishi (eds) Chapter 6, pp 1 7 1 - 223.

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

Propagation with Level Set Methods," 22 pp.

8. Mathematica Version 6 software: see http://www.wolfram.com/

9. 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.

10. Toro, E.F., "Reimann Solvers and Numerical Methods for Fluid Dynamics: A Practical

Introduction," Second Edition, Springer (1999).

11. Gottlieb, S., C.W. Shu, E. Tadmor, "Strong Stability-Preserving High-order Time Discretization Methods," SIAMReview, Vol. 43, Number 1, pages 89-1 12, (2001).
Y

100

1

'

' i ' ' ' i ' '_

80 If/"-OA

60 *"vMI

llllll

J/JJ/J 40 ^eSzfr

20 -

~

i i , , -40 -20 20 40

Figure 1. 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).
100 100

xo

60

40 e

80

60

40 G 20 - 20 -

-40 -20

Y

20 40

100

1

!

1

80 -

60 -

40 -

20 -

i 1 i

-

-40 -20 20 40

Figure 2. Comparison of the solutions obtained by the two methods, initially and every six seconds.
Figure 3. 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.

10
-40 -20 20 40 -40 -20 20 40

Figure 4. 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.

11
Figure 5. A wind-driven fire front shown at four times during the level-set simulation. The front is initially generated by 3 circular fire fronts placed symmetrically with respect to x. The head of each

circular front propagates in the direction of the wind, while the flanks propagate laterally and merge

into a single front. For this simulation, the ambient wind is 3 m/s in the y-direction.

12