Citation
Reaction diffusion equations and numerical wildland fire models

Material Information

Title:
Reaction diffusion equations and numerical wildland fire models
Creator:
Kim, Minjeong
Publication Date:
Language:
English
Physical Description:
xvii, 139 leaves : ; 28 cm

Subjects

Subjects / Keywords:
Reaction-diffusion equations ( lcsh )
Wildfires -- Mathematical models ( lcsh )
Differential equations, Partial ( lcsh )
Differential equations, Partial ( fast )
Reaction-diffusion equations ( fast )
Wildfires -- Mathematical models ( fast )
Genre:
bibliography ( marcgt )
theses ( marcgt )
non-fiction ( marcgt )

Notes

Bibliography:
Includes bibliographical references (leaves 131-139).
General Note:
Department of Mathematical and Statistical Sciences
Statement of Responsibility:
by Minjeong Kim.

Record Information

Source Institution:
|University of Colorado Denver
Holding Location:
|Auraria Library
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
761029679 ( OCLC )
ocn761029679
Classification:
LD1193.L622 2011d K55 ( lcc )

Full Text
REACTION DIFFUSION EQUATIONS AND NUMERICAL WILDLAND FIRE
MODELS
by
Minjeong Kim
M.S., Kyungpook National University, South Korea, 2001
A thesis submitted to the
University of Colorado Denver
in partial fulfillment
of the requirements for the degree of
Doctor of Philosophy
Applied Mathematics
2011


This thesis for the Doctor of Philosophy
degree by
Minjeong Kim
has been approved
by
Jan Mandel
Loren Cobb
Jiangguo (James) Liu
Weldon Lodwick
Date


Kim, Minjeong (Ph.D., Applied Mathematics)
Reaction Diffusion Equations and Numerical Wildland Fire Models
Thesis directed by Professor Jan Mandel
ABSTRACT
In this thesis, we formulate a mathematical model based upon partial differential
equations that describes wildland fires. We first investigate the reaction rate and then
the PDE-based wildland fire model is discussed. The traveling wave solutions of the
wildland fire model and measured temperature profile are presented. We also discuss
the identification of coefficients, and a dimensionless form of wildland fire model.
Standard numerical schemes are used for the numerical simulations, and various
numerical tests are performed in order to examine the properties of the numerical
model and numerical discretization. An efficient implementation of a finite differ-
ence scheme in 2D is described and discussed. Furthermore, we present a fireline
propagation model by using the level set method, and observed problems in the level
set method are discussed.
Lastly, we analyze a reaction-diffusion model in ID in phase space. We present
some concepts of dynamic systems and prove the existence of a traveling wave in
phase space under some assumptions. We obtain rigorous inequalities which show the
existence of a traveling wave solution in phase space, and we finally present numerical
solutions of the traveling wave in phase space.


This abstract accurately represents the content of the candidates thesis. I recommend
its publication.
Signed
Jan Mandel


ACKNOWLEDGMENT
I would like to give my special thank to Professor Jan Mandel as my advisor, and as
a mentor. Without his support, guidance and endless patience, my dissertation would
have been impossible. I am also indebted to Professor Lynn Bennethum and Professor
Loren Cobb for their advices. I thank to Professor Sangdong Kim for motivating me
to study abroad.
I would also like thank the faculty and my fellow students at the Department
of Mathematical and Statistical Sciences for providing great working conditions. In
particular, thank to Dr. Bedrich Sousedik and Dr. Christopher Harder for studying
and discussions.
Lastly, I thank to my family, and in particular, my husband and my parents.


Contents
1 Introduction 2
2 PDE-based model 7
2.1 Introduction.............................................. 7
2.2 Formulation of the model................................. 10
2.3 Derivation of the model ................................. 13
2.3.1 Heat and fuel supply balance equations............. 13
2.3.2 Reaction rate...................................... 17
2.3.3 Heat balance with full fuel supply..................20
2.4 Identification of coefficients .....................23
2.4.1 Reaction rate coefficients..........................23
2.4.2 Cooling coefficient ................................26
2.4.3 Scales and non-dimensional coefficients.............29
2.4.4 Related work....................................... 36
2.5 Calibration of coefficients in one dimension..............37
vi


3 Numerical methods for the PDE model
41
3.1 Efficient numerical implementation............................41
3.2 Computational tests ..........................................48
3.2.1 Line ignition with oscillatory wind ...................50
3.2.2 Spot fire without wind.................................50
3.2.3 Effect of the wind speed...............................52
3.2.4 Some observations .................................... 56
3.3 Numerical convergence tests.................................. 58
3.4 Numerical difficulties........................................61
4 The fireline propagation model 68
4.1 Level set based wildland fire model.......................... 70
4.1.1 Some basics about the level set method................ 71
4.2 Numerical schemes for fireline propagation model............. 74
4.2.1 Numerical schemes for the Level Set equation...........74
4.2.2 Observations and future work...........................77
5 Phase-space analysis of the traveling wave 83
5.1 Traveling wave in phase space.................................84
5.1.1 Parameterized system of ODEs in phase space............85
5.2 Properties of dynamical systems...............................93
5.3 Traveling wave in the fire model as a dynamical system .... 100
5.3.1 Invariant manifolds...................................100
5.3.2 Nullclines............................................103
vii


5.4 Inequalities for the speed of the traveling wave..........Ill
5.5 Existence of a traveling wave speed.......................118
5.6 Numerical solution in phase space.........................123
6 Conclusion 129
viii


List of Figures
2.1 Sample reaction heat balance function f(T) from equation
(2.21)............................................................21
2.2 Reaction heat balance potential U.................................22
2.3 Solution with the Arrhenius reaction rate. Due to nonzero
reaction rate at ambient temperature, fuel starts disappearing
and thus a propagating combustion wave does not develop.
The coefficients are k = 2.1360 x 10_1m2s_1 A*3, A = 1.8793 x
102Ks~\ B = 5.5849 x 102A, C = 4.8372 x lO^A-1, and
Cs = 1-6250 x lO^s"1. Tc = 1200A and T{ = 670A...............24
2.4 Solution with Arrhenius reaction rate modified by temperature
offset Ta to force zero reaction rate at ambient temperature.
A propagating combustion wave develops. The coefficients
are k = 2.1360 x 10-1m2s-1 A-3, A = 1.8793 x 102As-\
B = 5.5849 xlO2 A, C = 4.8372x10- A"1, and Cs = 1.6250 x
lO-V1. Tc = 1200A and Tj = 670A.............................. 27
IX


2.5 Temperature profile of a traveling wave. The wave moved
about 398m from the initial position in 2300s...............28
2.6 The fuel fraction remaining after the combustion wave as a
function of dimensionless parameters A and (3...............32
2.7 Ratio of the leading and the trailing edge of combustion wave
with respect to A and (3 from (2.31). The width is measured
as the length of the part of the wave above temperature 500A7 33
2.8 Width of the trailing slope as a function of A and /?. The rela-
tion between the wave speed and (diffusion coefficient)0,5, i.e. \fk.
We see approximately linear relations.........................35
2.9 Periodic patterns of the solutions by varied (3. (a) (3 = 0.4172,
(b) (3 = 0.4636, (c) (3 = 0.4945, (d) /3 = 0.5718. Each graph
shows the temperature profile by timeframe...........38
2.10 Time-temperature profile (dotted line) measured in a grass
wildland fire at a fixed sensor location, digitized from Kremens
et al. [2003], and a computed profile (solid line) from simulation. 39
3.1 Five-point stencil........................................... 44
3.2 Pentadiagonal matrix Jtt by using standard five-point finite
difference discretization.....................................45
3.3 Diagonal matrix Jts, Jst, and Jss............................46
x


3.4 Column storage of the sparse Jacobian matrix. Five columns,
labelled as i m, i 1, i, i + 1, and i + m, are from the
pentadiagonal matrix Jtt The last three columns store the
diagonal matrices Jts> Jst, and Jss........................47
3.5 North to south straight line ignition and a sinusoidal wind
with speed 3m/s after 400s. The dark areas are burnt regions. 49
3.6 Spot fire without wind simulations by Section 2.2 equations
(2.2-2.3). The mesh size is 4m and time step is Is after 600s.
(a) is the numerical fire evolution and (b) is a contour of (a).
(c) is fuel consumption in (a).............................51
3.7 Effect of wind speed 0.005m/s fire simulation of Section 2.2
equations (2.2-2.3) after 800s. Time step is Is and mesh size is
4m. With a small wind speed the fire also moves more slowly
and it is closer to circular shape, (b) is a contour plots of (a).
(c) is a fuel consumption in (a)...........................53
3.8 Effect of wind speed 0.05m/s fire simulation of Section 2.2
equations (2.2-2.3) after 400s. The numerical solution moves
from southwest to northeast. Time step is Is and mesh size is
8m. (b) is a contour plots of (a), (c) is a fuel consumption in
(a)........................................................ 54
xi


3.9 Effect of wind speed 3m/s fire simulation of Section 2.2 equa-
tions (2.2-2.3). Fire moves from southwest to northeast. Time
step is Is and mesh size is 16m. (b) is contour plots of (a),
and (c) is a fuel consumption of (a). Note that the sharp front
edges in (a).................................................. 55
3.10 Wind effect of fuel consumption, (a) Fuel consumption with
domain size 256m x 256m and wind speed 3m/s after 20s.
(b) Fuel consumption with wind speed 5m/s. The dark area
means fuel consumption and light yellow is unburnt region.
We see that the fire spreads faster with a higher wind speed. 57
3.11 Errors from Table 3.1 in the logarithmic scale(left), and the
convergence norm ratios(right). Errors from Table 3.2 in the
logarithmic scale(left), and the convergence norm ratio(right). 60
3.12 Errors from Table 3.3 in the logarithmic scale(left), and the
convergence norm ratios (right). Errors from Table 3.4 in the
logarithmic scale(left), and the convergence norm ratio(right). 62
3.13 Fire simulation with large mesh size of 80m and wind speed
of 0.005m/s. The time step is 2s, and ignition area is 400m
by 320m. (a) is the initial plot of (b).......................65
3.14 Fire simulation with mesh size of 2.5m after 600s. All other
conditions are the same in Figure 3.13........................66
xii


3.15 Fire simulation with large mesh size of 80m and wind speed of
1 m/s. While the case of small wind speed (3.13), a traveling
fire develops............................................................66
4.1 the normal to the curve(4.6)...................................... 72
4.2 A circular fire propagation with constant wind. The contour
is a fireline...................................................... 77
4.3 (a) A trace of circle movement along wind direction. Initial
value circle, (b) The circle moved to the south-east. As time is
passing, the circle keep shrinking. The re-initialization avoids
the shrinking problem in (c)..................................... 78
4.4 A trace of motion in a rotation velocity field. Re-initialization
method is added in (a) and (b)................................... 79
4.5 (a) Points on the zero level set on the mesh grid lines are stored
as a set A and points on the zero level set function, which is
not on the mesh grid lines, are stored as a set B. (b) Closest
points to the grid points are stored as a set C..................80
4.6 Motion of a circle in a rotational velocity field. The Narrow
Band scheme with re-initialization is added in (a) and (b).
There still exists some mass loss while the circle is rotating
and eventually the circle does not match at the end.............. 82
5.1 Intersection of the solution of equations (5.35) (5.36) with
the plane V = 0 at (zq, Cq)
xiii
98


5.2 The unstable, stable, and center subspaces, Eu, Es, and Ec,
and manifolds, WU(Xo), VT5(Xo), and WC(Xo) respectively. 102
5.3 Nullclines Vz = 0 of (5.32) with different wave speed for 5 = 0,
5 = 0.5, and 5=1 and the traveling wave speed c = 0.45 in
(a), c = 0.4 in (b), c = 0.35 in (c). The parameters are
A = 0.25 and (3 = 0.2...........................................104
5.4 Nullclines V = £(Se~7 XT) and V = 0, for a fixed 5 with
direction fields. Scaled to unit length in the (V,T) projection. 105
5.5 Direction field of equation (5.35) in the plane 5=1 with the
nullcline Vz = 0 and tangent line to the nullcline at the origin
in (a), in the plane 5 = 0.8 in (b), and in the plane 5 = 0.5 in
(c) also shown. In (a) and (b), Tig is the ignition temperature.
The parameters are c = 0.4, A = 0.25, and (3 = 0.2..............106
5.6 The solutions of the systems (5.35)-(5.36), projected on the
(T, V) plane. The solution is not a traveling solution because
the solution does not crosses the nullcline Vz 0..............107
5.7 The solution of equations (5.35)-(5.36). A traveling wave so-
lution is shown. The view is rotated around the V-axis. The
parameters are c = 0.7992, A = 0.06, and (3 0.15671... 109
xiv


5.8 Nullcline Vz = 0 and the solution of equation (5.35). The
view is rotated around the V-axis. The maximum value of the
nullcline is negative and the solution keeps going upward from
the critical point. Thus, no traveling wave can exist. The
5.9 Solution of equations (5.35)-(5.36) in the plane 5=1 and
S = 0.9. We see that A) in the region where
V > ^(e_1 A)T. The solution of equations (5.35)-(5.36)
keeps going upward, and this solution is not the traveling wave
solution. The parameters are c = 0.4, A = 0.25, and (3 0.4. 113
5.10 The solution of equation (5.35)-(5.36) in 3D. In (a), c is too
large, and the solution with the plane T = 0. In (b), c is
too small and the solution with the plane V = 0. In (c), the
end point of the solution goes to the limit point (5.37), and
this solution is the traveling wave solution. The parameters
parameters are c = 0.075, A = e *, and (3 = 0.15671
110
A = 0.06 and (3 = 0.15671
114
5.11 Function e t AT in Theorem 8, which is denoted by 0(T)
here. The parameters are c = 0.6 and (3 = 0.4................
115
5.12 (a) p with varied A in equation (5.50). (b) with varied A.
The speed of the traveling wave c = 0.07992 and c <
(c) The traveling wave solution of equations (5.35)-(5.36) with
the plane (T, V)
. 116
xv


5.13 (a) hfi- with varied A and the parameter c = 0.6. We see
y 1
that the parameter c > (b) The solution of equations
(5.35) -(5.36) with the plane (T, V) and the parameter c = 0.6.
The solution of equations (5.35)-(5.36) keeps increasing, and
clearly it is not the traveling wave solution of the equation
(5.35) -(5.36). The parameters are 0 = 0.4 and A = 0.254.
5.14 Function <&(c) is continuous.............................
117
121
xvi


List of Tables
3.1 Convergence ratio for decreasing time step, wind speed 2m/s,
the mesh size h = 2m....................................59
3.2 Convergence rate for decreasing time step, wind speed 4m/s,
the mesh size h = 4m....................................61
3.3 Numerical convergence behavior for both the time step dt and
the mesh size h, wind speed 0.5m/s......................61
3.4 Numerical convergence behavior for both the time step dt and
the mesh size h. Wind speed is 8m/s.....................63
xvn


Chapter 1
Introduction
Reaction-diffusion equations arise in a large area of mathematical models
involving chemical reactions. As a practical application, there is a wide
interest in such equations, with numerous recent studies being performed
on, in particular, biological models such as the epidemic model and bacterial
colonies.
We present a highly simplified wildland fire model as an application of
reaction diffusion equations. This model has two physical variables, the fuel
temperature and the fuel mass. The speed of propagation of the fire is a
property of the solution of a system of differential equations, and it is thus
determined by the coefficients of the equations.
One reason for considering a PDE-based model is that even simple reaction-
diffusion equations are capable of complex nonlinear, unsteady behavior such
as pulsation and bifurcation. We note that the physical behavior can be
2


achieved by very simple models that can reproduce qualitative behavior of
some features of fire.
On the other hand, it is not easy to determine the coefficients of the
PDEs. In this thesis, properties of the combustion wave are used to calibrate
the parameters of the model.
Exact solutions to the PDE wildland fire model are not known. Thus,
computational simulation is important. In addition these PDE models can-
not exactly replicate the physical phenomena. Furthermore, numerical so-
lutions can exhibit errors related to the implementation itself. Therefore,
numerical tests and observations on the physical correctness of numerical so-
lutions are necessary. Efficient implementations can dramatically reduce the
costs of the simulation.
Although the reaction-diffusion-advection equations with enriched finite
element methods are designed and error estimates are given in Asensio et al.
[2004], Codina [1998], Franca et al. [2005], and Franca et al. [2006], these
works do not account for nonlinear effects at the leading edge. Error esti-
mates are provided in Asensio and Ferragut [2000] and Ferragut and Asensio
[2002] in the context of mixed finite elements, which is applied to a single
species combustion equation and more general reaction-diffusion problems,
but these works do not consider a fuel-balance equation. In this thesis, we
present a nonlinear reaction-diffusion equation with fuel consumption, and
we obtain traveling wave solutions.
Although it is known that coupling with the atmosphere is essential for
3


wildland fire behavior, such interaction is not considered in this thesis. We
instead limit our attention to investigating the wildland fire models by them-
selves and finding numerical solutions.
Because a common feature of the solutions of reaction-diffusion equations
is the development of a sharp wave, with the solution being almost constant
elsewhere, interface tracking techniques such as level set methods (Sethian
[1999] and Sussman et al. [1994]) are relevant here as well. This method is
a basis of the implementation in WRF-FIRE (Mandel et al. [2009], Mandel
et al. [2011]), and this empirical model developed by NCAR (Clark et al.
[2004], Coen [2005a]). In this empirical model, the fire propagation speed
is normal to the fireline and is a function of wind and terrain slope, and an
exponential decay of fuel from the time of ignition is assumed. Thus, the fire
propagation speed can be input directly from known fire behavior for various
fuels.
We observe that the numerical solution of PDEs is a traveling wave so-
lution. This observation gives us a motivation to study the traveling wave
solution of the equations. Thus, we are interested in the traveling wave solu-
tions of reaction-diffusion equations. The traveling wave solution has a form
in connection with a time and a speed of wave, and this form of solution
allows us to have ODEs instead of PDEs. In general, ODEs are easy to
analyze, and it often allows us to find an explicit solution of the equation.
The parameters of the equations determine a behavior of the solution
in phase space. The exact solution of reaction-diffusion equations should
4


represent the behavior of the model or physical properties. The analysis
of the direction field in phase space clearly indicates the properties of the
solution in phase space which depend on parameters. These dependencies
provide rigorous properties of traveling wave solutions in phase space.
This thesis is organized as follows. In Chapter 2, we formulate a simple
PDE-based wildland fire model and identify the coefficients instead of taking
these from physical and chemical properties of fuel.
In Chapter 3, we first introduce a familiar discretization of unsteady state
reaction-diffusion-convection equations, and then provide an efficient imple-
mentation. Various computational tests and numerical convergence tests will
be performed. Finally, we will discuss observed numerical difficulties.
In Chapter 4, we introduce the level set method and apply it to the wild-
land fire problem. An observed numerical drawback, the mass loss problem,
is introduced from some numerical experiments. We then apply simple res-
olutions to cure this difficulty, and then it turns out that preventing this
difficulty is a interesting area in studying the level set method.
In Chapter 5, we turn our attention to the phase-space method for the
reaction-diffusion equations. The phase-space method is often used to reduce
second order equations to first order equations. Then the first order equations
become autonomous dynamic systems. We introduce the basic concepts of
dynamics and investigate a phase-space direction field, which is 3D here. We
prove the existence of a traveling wave under suitable assumptions.
The authors personal contributions in this work are as follows. Chapter
5


2: Derivation of the simple PDE-based model (2.2)-(2.3) in (2.9)-(2.12) and
identification of coefficients. Analysis of dimensionless parameters. Chapter
3: Matlab and Fortran implementation of equations (2.2)-(2.3). Compu-
tational tests of numerical convergence and observations of numerical diffi-
culties. Furthermore, the efficient implementation of equations (2.2)-(2.3)
allowed for faster data assimilation in the DDDAS(Dynamic Data Driven
Application Systems) project. Chapter 4: Survey of level set method and
the investigation of problems of level set method in experiments. Contri-
butions to the implementation of the method and fuel burning in the fire
code in WRF-Fire. Chapter 5: Rigorous results in equations (5.35)-(5.37).
Phase-space analysis of the traveling wave using invariant manifolds. Joint
work with advisor: Proofs of Lemmas and Theorems showing that a traveling
wave does not exist under certain conditions, and thus proving inequalities
for the parameters of the traveling wave. Proof of the existence of traveling
wave solutions in phase space under some assumptions.
6


Chapter 2
PDE-based model
2.1 Introduction
Modeling wildland fire requires us to understand the properties of fire be-
haviors. This modeling contains complexities of fuel characteristics and heat
transfer mechanism with the effect of wind. We formulate a highly simple
mathematical wildland fire model. It is a plausible mathematical model of
wildland fire which helps to understand the physical meanings of a compli-
cated real wildland fire.
Various aspects of special cases of the reaction-diffusion equations (2.2)-
(2.3) in Section 2.2 below have been studied in a number of papers. A formal
expansion in an Arrhenius reaction model to get the wave speed and a pre-
diction of whether a small fire will or will not spread was used by Weber
[1991]. The ignition wave speed and the extinction wave speed were com-
7


puted numerically by Mercer and Weber [1995]. The speed of combustion
waves were analyzed by asymptotic expansion and the stability of the travel-
ing wave front was obtained by linearized system of ODEs and its eigenvalues
by Gubernov et al. [2003].
The speed of the traveling wave in a simplified model with the reaction
started by ignition at a given temperature was derived by Campos et al.
[2004] under the assumption that a temperature-independent reaction rate
is proportional to the fuel remaining, similar to the model by Balbi et al.
[1999]. This temperature-independent reaction rate is often used in a PDE-
based model, and the reaction rate of the model has a proportionality to the
fuel remaining.
The speed of a combustion traveling wave fronts was investigated by using
relaxation and shooting methods by Gubernov et al. [2003] and Gubernov
et al. [2004] from the solid to gaseous fuel (i.e., also with fuel diffusion instead
of just temperature diffusion) for different values of parameter, j3 (the ratio
of activation to heat release). They found that the large /? and gaseous fuel
give a good correspondence between analytical and numerical results.
A nonlinear eigenvalue problem for a traveling wave in a different combus-
tion problem, with fuel reaction, was derived by Norbury and Stuart [1988a]
and Norbury and Stuart [1988b], The equation was solved numerically by
the shooting method, and the existence and stability of the traveling wave
solution is studied.
For an analytical study of the evolution of waves to a traveling waveform,
8


see Sherratt [1998], and for a numerical study, see Gazdag and Canosa [1974]
and Zhao and Wei [2003]. Proofs of the existence of a solution and attractors
for reaction-diffusion equations are given by [Robinson, 2001, Ch. 8 and 11],
although traveling waves are not mentioned. The existence of a solution for a
finite time is proved by Asensio and Ferragut [2002] for the reaction-diffusion
equation of (2.2) without convection, but the nonlinear diffusion term was
replaced by V (kT3VT). Again, traveling waves and fuel depletion were not
considered.
Equation (2.2) without fuel depletion (i.e., with constant S) and without
wind (i.e., v = 0) is a special case of the nonlinear reaction-diffusion equation,
V (VT) + / (T). (2.1)
Reaction-diffusion equations of the form of equation (2.1) are known to pos-
sess traveling wave solutions, which switch values close to stationary states
given by / (T) = 0 (Gilding and Kersner [2004] and Infeld and Rowlands
[2000]). The simplest model problem is Fishers equation with f (T) =
T(l T), for which the existence of a traveling wave solution and a for-
mula for its speed were proven by Kolmogorov et al. [1937].
Two reaction models (solid and pyrolysis gas) are considered by Seron
et al. [2005], who argue that the modeling of pyrolysis by a separate reaction
is essential for capturing realistic fire behavior. For more complicated theo-
retical models of this type that include other components, e.g., water vapor,
9


see Grishin [1996].
For a similar system, but with fuel diffusion, the existence and the speed
of traveling waves were obtained by asymptotic methods already in classi-
cal works, summarized by the monograph Zeldovich et al. [1985]. For the
system (2.2)-(2.3), an approximate combustion wave speed by asymptotic
methods was obtained by Weber et al. [1997]. However, the heat loss term,
C (T Ta), was not included in the approximation. No heat loss case is
equivalent to our setting C = 0 in equations (2.2)-(2.4).
The proposed model in this chapter is based on a reaction-diffusion-
convection equation (2.2)-(2.4) in Section 2.2. A more detailed derivation
of the model from physical principles such as the conservation of mass and
energy is contained by Section 2.3.
In Section 2.4, we identify coefficients, using a dimensionless form of
reaction-diffusion-convection equation (2.29)-(2.30). Non-dimensional pa-
rameters presented in this section are used to identify coefficients.
This chapter is partly based on Mandel et al. [2008].
2.2 Formulation of the model
We consider a model of a fire in a layer just above the ground. First, we
define the following terms:
T (K) is the temperature of the fire layer (K is degrees Kelvin),
10


S [0,1] is the fuel supply mass fraction (the relative amount of fuel re-
maining: the mass of remaining fuel divided by the total mass of fuel),
k (m2s 1) is the thermal diffusivity,
A (Ks 1) is the temperature rise per second at the maximum burning rate
with full initial fuel load and no cooling present,
B (K) is the proportionality coefficient in the modified Arrhenius law below,
C (K-1) is the scaled coefficient of the heat transfer to the environment,
Cs (s-1) is the fuel relative disappearance rate,
Ta (K) is the ambient temperature, and
v (ms-1) is the wind speed given by atmospheric data or model.
The model equations consist of the conservation of energy and balance of
fuel:
^ = V (kVT) v VT + A (5e-B/(T-T) C (T Ta)), (2.2)
(2.3)
with the initial values
S (t0) = 1 and T (t0) = Tinit.
(2.4)
11


The diffusion term V (&VT) models short-range heat transfer by radia-
tion in a semi-permeable medium, v VT models heat advected by the wind,
Se-B/(T~Ta) jg ra^e fuej jg consumed due to burning, and the fuel deple-
tion gains heat when the fire is ignited, i.e., T > Ta. AC (T Ta) models
the convective heat lost to the atmosphere. The reaction rate e~B^T~Ta'> is
obtained by modifying the reaction rate e~B!T from the Arrhenius law by
an offset to force zero reaction at ambient temperature, with the resulting
reaction rate smoothly depending on temperature. The reaction rate from
the Arrhenius law will be discussed in Section 2.3.2.
It is known that systems of a form similar to (2.2 2.3) admit traveling
wave solutions. An approximation to the temperature reaction equation
gives the size of the reaction zone and the slope of the temperature curve.
The temperature in the traveling wave has a sharp leading edge, followed
by an exponentially decaying cool-down trailing edge. This was observed
numerically but we were not able to find a rigorous proof in the literature
for exactly this (2.2 2.3) case.
Models that do not guarantee zero combustion at ambient temperature
suffer from the cold boundary difficulty: by the time a combustion wave
gets to a given location, the fuel at that location is depleted by the ongoing
reaction at ambient temperature. So, no perpetual traveling combustion
waves can exist, and there are only pseudo-waves that travel only for a
finite time and then vanish (Berestycki et al. [1991], Mercer et al. [1996], and
Zeldovich et al. [1985]).
12


2.3 Derivation of the model
We consider fire in a ground layer of some unspecified finite small thickness
h. The fire layer consists of the fuel and air just above the fuel. All mod-
eled quantities are treated as two dimensional, homogenized in the vertical
direction over the ground layer. We will not attempt to substitute the coef-
ficients from material properties because of the degree of simplification and
the uncertainty present in the homogenization. Instead, physical laws will be
used to derive the form of the equations and the coefficients will be identified
later from the dynamic behavior of the solution. We first derive the system
of PDEs based on conservation of energy and fuel reaction rate in Section
2.3.1 and then discuss the choice of the reaction term in Section 2.3.2
2.3.1 Heat and fuel supply balance equations
In this section, we write the equations in a form suitable for identification of
the coefficients, which will be formally done in Section 2.4 The chemical
reactions are a heat source and the chemical reaction of burning is generating
the heat. We model the burning as a reaction where the reaction rate only
depends on the temperature. Then the reaction rate with a coefficient of
proportionality, Cs (1 /s), is Csr(T), and r is a dimensionless reaction rate.
Let Fc > 0 (kgm~2) be the concentration of fuel remaining. Then the
fuel consumption rate is proportional to the reaction rate where the amount
13


of fuel available,
-^ = -F'Csr(T)
We introduce the mass fraction of fuel by
(2.5)
S
Fc
Fo
where F0 is the initial fuel quantity.
The heat generated per unit surface area is then proportional to the fuel
lost,
qg = AiFcCsr (T) (Wm~2), (2.6)
where Ay (J kg~x) is the heat released per unit mass of fuel.
Heat transfer is due to radiation and convection to the atmosphere. The
short-range heat transfer due to radiation and turbulence is modeled by dif-
fusion. The two dimensional heat flux vector then is
qr = AqVT (Wm"1). (2.7)
Equation (2.7) models the heat transfer from the area of high temperature
to the area of low temperature through the unit length, and equation (2.7)
guarantees this heat transfer. The constant Aq (WK~l) will be identified
later.
Heat per unit area lost due to natural convection to the atmosphere and
14


to buoyancy is given by Newtons law of cooling,
qc = Ca(T Ta) (Wm~2) , (2.8)
where Ta is the ambient temperature (K) and Ca (Wm~2K~x) is the heat
transfer coefficient. In this model, it is assumed that the convective heat
transfer is dominant, and so the effect of radiation into the atmosphere is
included in (2.8).
The material time derivative, which is a derivative taken along a path
following the flow with velocity field, v, of the temperature is given by
DT
~Dt
dT
= Ti+'"vt-
(2.9)
Here dT/dt is the Eulerian (or spatial) time derivative of temperature, v is
the (homogenized) velocity of the air, p is the homogenized surface density
of the fire layer (kg m-3), and cp is the homogenized specific heat of the fire
layer (J kg~1K~l). Again, none of the coefficients h, p, or Cp can be assumed
to be actually known. The units of the product hpcp are JK~lm~2. The
velocity vector v is obtained from the state of the atmosphere as data, or in
future work by coupling with an atmospheric model.
From the divergence theorem, we now obtain the conservation of energy
in the fire layer as
DT
hpCp~Dt = V Qr + Q9 ~ Qc ^21^
15


The goal is to obtain a system of equations with a minimal number of
coefficients and in as simple form as possible. In addition, we wish to relate
the behavior of the solution to the non-dimensional parameters in Section
2.4.3 below, rather than to material and physical properties of the medium,
which are in general unknown. Substituting in the appropriate values for
the heat sources and fluxes, (2.7), (2.6), and (2.8), into (2.9) and (2.10), and
some simple algebra, we obtain the energy balance and fuel reaction rate
equations
= V-(fcVT)-v VT + A(Sr(T)-C0(T-Ta)), (2.11)
^ = -CsSr(T), (2.12)
with
k = ki/(hpcp), A = AiCs/{hcpp), and C0-Ca/Ai.
Alternatively, we could have started with the disappearance of fuel on the
left hand side of the heat balance equation (2.10) (fuel that has burned does
not need to be heated), which would lead to an equation of the form
(1 + CjS) ^ = V (kVT) v VT + A (Sr (T) -C0(T- Ta))
at
instead of (2.11). We have chosen not to do so since our goal is to work with
the simplest possible model, whose coefficients can be identified. The fuel dis-
16


appearance in the heat equation affects the temperature profile significantly
only in the reaction zone, which is the highest part of the temperature curve.
Before the ignition, there is a full fuel load, 5=1, and after a fairly short
reaction time, well before most of the cooling takes place, the remaining fuel
settles to some residual value which then remains constant. The effect of
the decreased heat capacity of the remaining fuel is then absorbed into the
cooling term ACq.
2.3.2 Reaction rate
In this section we discuss the reaction rate with an Arrhenius type of equa-
tion, and we derive a modified equation by introducing dimensionless tem-
perature and parameter. The Arrhenius equation from physical chemistry is
given by
r{T) = Ae(~-tr\ (2.13)
where r is a rate constant of reaction and T is a temperature, and the reaction
rate falls exponentially with decreasing in the temperature. This equation
gives the reaction rate as a function of temperature, where A is the pre-
exponential factor, Ra is the gas constant, Ea is an activation energy, and
Ea/Ra has units K. Equation (2.13) is valid only for a single reaction and
lower temperature regime, but it is used here as an approximation form.
The activation energy Ea is mostly obtained empirically. One consequence
of (2.13) is that the reaction has a nonzero rate at some temperature above
17


absolute zero.
As already explained in Section 2.2, the cold boundary problem is caused
by the fact that the time scale of burning is much smaller than the oxidation
rate at ambient temperature (Williams [1985]). A dimensionless parameter e*
of equation (2.17) below is used in Mercer et al. [1996], Asensio and Ferragut
[2002], Mercer and Weber [1995], and the numerically investigated results in
the dependency of dimensionless parameters with respect to a wave speed
can be found in Mercer et al. [1996] and Mercer and Weber [1995]. However,
none of these papers explains why they used equation (2.16) as the reaction
rate, and they did not mention why the dimensionless parameter e* needs to
be introduced.
We summarize the assumptions and a process of deriving the dimension-
less parameter e* (Zeldovich et al. [1985]). The idea is that the parameter
e* is approximating the reaction rate in a narrow range of temperature. The
narrow range of temperature, AT = T T, where the temperature T* is
taken from a neighborhood where the reaction occurs. This temperature
T* can be Ta in a study of spontaneous ignition, or it can be the combus-
tion temperature in a study of propagating combustion. In the problem, T,
is determined during problem solving since T* is the unknown. Here, the
temperature T is the ambient temperature Ta.
By the method of expanding the exponent, which is originally suggested
18


in Williams [1985], we introduce the dimensionless temperature difference u,
u =
Eg
RaT?
(T T*).
(2.14)
From the identity
we obtain
1 +at
1 T.
= 1 -
AT
T.

Ea 1 Ea
RgT ~ T(T + AT) TaT, 1 + ^ ~ RaTt
(
)
(2.15)
Equation (2.15) allows us to rewrite the reaction rate as
___ Eg __ Eg U
c =6 KaT* gHTeTtTj
(2.16)
where u is given by (2.14) and the parameter
e.
RaTt
Eg '
(2.17)
We simply rewrite the Arrhenius reaction rate (2.16) as
r(T) = e-B'T,
(2.18)
where the coefficient B has units K. This equation is valid only for gas fuel
premixed with a sufficient supply of oxygen. This approximation ignores fuel
19


surface effects but it is widely used nonetheless. We establish the reaction
rate by neglecting the initial temperature, i.e., by assuming that the reaction
rate is equal to zero at the ambient temperature and modify (2.18) so that
no oxidation occurs below some fixed temperature, To (Section 2.4), and take
instead
Note that the fuel consumption rate is a smooth function of T, which is
favorable for a numerical solution, unlike the rate in Asensio and Ferragut
[2002], where a cutoff function was used.
2.3.3 Heat balance with full fuel supply
Consider first the hypothetical case in which T is constant in space, so that
only heat due to the reaction (burning) and natural convection contribute
non-zero terms in the heat equation, (2.2), and essentially the full initial fuel
supply Fo is present at all times (the rate of fuel consumption is negligible,
Cs ~ 0, so S ~ 1). Then, (2.11) reduces to
Constant values of temperature which are solutions to (2.20) are called equi-
librium points, and at these points the heat produced by the reaction equals
(2.19)
= A (e-B/(T-To) -C(T- Ta)).
(2.20)
20


T =300 T =1200 T=670
a c i
Figure 2.1: Sample reaction heat balance function f(T) from equation (2.21).
the heat lost to the environment,
f(T) = e-BI(T-To) -C{T-Ta) = 0. (2.21)
Equation (2.21) has at most three roots (Frank-Kamenetskii [1955]), see
for example, Figure 2.1.
The first zero, denoted as Tp, called the lower temperature regime by
Frank-Kamenetskii [1955], is a stable equilibrium temperature. If the tem-
perature goes below this temperature then the heat generated from the re-
action dominates and the temperature rises, since some reaction is present
even at ambient. If the temperature goes above this temperature then cool-
ing dominates and the temperature decreases. The middle zero, 7), is an
21


10
T =300 T =1200 T.=670
a c i
Figure 2.2: Reaction heat balance potential U
unstable equilibrium point. If the temperature goes above Ta and below
then cooling dominates, and the temperature decreases. Above 7), the heat
due to chemical reactions dominate and the temperature increases. We refer
to Ti as the auto-ignition temperature, the temperature above which the re-
action is self-sustaining (Quintiere [1998]). The stable equilibrium Tc is the
maximum stable combustion temperature, assuming replenishing of the sup-
ply of fuel and oxygen. The temperature Tc is called the high temperature
regime by Frank-Kamenetskii [1955]. The stability properties of the equi-
librium points are also clear from the graph of the potential U(T), defined
by U'(T) = f(T). The stable equilibrium points are local minima of the
potential, while the auto-ignition temperature is a local maximum and thus
an unstable equilibrium (Figure 2.2).
22


2.4 Identification of coefficients
We wish to calibrate the coefficients from physically observable quantities of
the fire rather than physical material properties. It is not simple to obtain
a reasonable behavior of the solution from substituting physical coefficients
values into the equations. Further, as explained in Section 2.3, because of
a number of simplifying assumptions employed and because of the homoge-
nization of coefficients over a fire layer of unspecified thickness, it is not quite
clear what the material properties should be anyway.
We first use basic reaction dynamics and a reduced model to find rough
approximate values of the coefficients that produce a reasonable solution.
After that, we transform the equations to a non-dimensional form, which
allows us to separate the coefficients into those that determine the qualitative
behavior of the solution and those that determine the scales. We use the
approximate coefficients obtained from the reduced models as initial values
for the identification of the coefficients by the nondimensionalization method
to match observed temperature profiles.
2.4.1 Reaction rate coefficients
While the coefficients B and C/F0 in (2.21) are generally unknown, they
can be found from the equilibrium temperature points. Suppose that two
roots Ti and Tc of f(T) from (2.21) are given such that T0 Then simple algebra gives a unique B and C in (2.21) from the conditions
23


Fuel supply mass fraction
Figure 2.3: Solution with the Arrhenius reaction rate. Due to nonzero
reaction rate at ambient temperature, fuel starts disappearing and thus
a propagating combustion wave does not develop. The coefficients are
k = 2.1360 x 10-1m2s-1 A-3, A = 1.8793 x 102Ks~\ B = 5.5849 x 102A,
C = 4.8372 x 10-5K~\ and Cs 1.6250 x 10_1s_1. Tc = 1200K and
Ti = 670K.
24


/ (Ti) = / (Tc) = 0, thus
e~B/Ti (7 (T< Ta) = 0,
e-B/Tc-C(Tc-Ta) = 0.
Eliminating C gives
0 = e~BtTc (T{ Ta)
e-B/T
(Tc-Ta),
which is equivalent to
Tj-Ta
Tc-Ta
_ e-B(l/Ti-l/Tc) = Q5
hence
Tc-T0 Ti-To
-B/(Ti-T0)
and C = ------
Ti-Ta
(2.22)
When T{ is too close to T, the resulting energy balance equation f(T) = 0
has only two roots. This, however, does not occur for the values of Ta, Tt,
and Tc of interest.
First consider the solution of (2.11-2.12) and the reaction rate (2.19)
with T0 = 0. Then the reaction is the Arrhenius rate known from chemistry
and there is a nonzero reaction rate at T = Ta. This results in fuel loss
everywhere, and, in our computational experiments, no traveling combustion
25


wave developed in Figure 2.3 since, after a relatively short time, there was
not enough fuel to sustain combustion due to the cold boundary difficulty
(see Section 2.2).
Since the values of B and C are obtained from realistic Ti and Tc, the fuel
disappearance happens quickly and the traveling wave is sustained (Figure
2.3), we force the reaction rate, r(T), to be zero at ambient temperature by
choosing the offset T0 = Ta. Using the offset by Ta is essentially the same
as assuming that the ambient temperature is absolute zero, as commonly
done in combustion literature Weber et al. [1997]. In this case, a propagating
combustion wave develops on Figures 2.4 and 2.5. Therefore, we use T0 = Ta.
ft should be noted that traveling combustion waves, such as in Figure
2.5, are caused by the combined effect of reaction and diffusion; convection
does not play a role in this section. The reaction heat diffuses forward on
the leading edge, heating the fuel ahead of the wave, until the reaction ahead
of the wave can sustain itself, thus causing the combustion to spread. On
the trailing edge, the reaction drops off due to fuel depletion, after which
temperature decays due to cooling.
2.4.2 Cooling coefficient
The coefficients B and C in the modified reaction form (2.21) have been
determined from reasonable values of Tc and Ti by (2.22). Now we want to
determine the remaining coefficient, A, which can be done if we know the
characteristic cooling time. Consider the trailing edge of a traveling combus-
26


reaction
Figure 2.4: Solution with Arrhenius reaction rate modified by temperature
offset Ta to force zero reaction rate at ambient temperature. A propagating
combustion wave develops. The coefficients are k 2.1360 x 10_1m2s_1 K~3,
A = 1.8793 x 102Ks~\ B = 5.5849 x 102A, C = 4.8372 x lO"5^-1, and
Cs = 1.6250 x lO^s-1. Tc = 1200K and T* = 670A.
27


Traveling combustion waves
Figure 2.5: Temperature profile of a traveling wave. The wave moved about
398m from the initial position in 2300s.
tion wave, after all or most of the fuel has been depleted, temperature drops,
and heat generated by the reaction and diffusion drop to an insignificant
level. From that point on, the temperature approximately satisfies
dT
- = -AC(T-Ta).
Thus, at the trailing edge, given T at some time to, we have
T(t) = Ta + (T(t0) Ta) e~AC^
and we can define the characteristic cooling time, tc, (see similar discussions
in Balbi et al. [1999]), to be the time which the fire layer takes to cool by a
28


factor of 1/e, i.e.,
T(t0 + tc) Ta = (T(t0) Ta),
which gives that ACtc 1, or
A
(2.23)
2.4.3 Scales and non-dimensional coefficients
We now write the model in terms of nondimensional variables, which control
the qualitative behavior of the system (2.2-2.3). Again, we do not consider
the wind here yet, and so
= V (kVT) + A (^Se~T^ C (T Ta)j , (2.24)
fjQ B
= -SCse~T=K, T > Ta. (2.25)
at
Given time temperature scale Tj, space scale xx, and time scale t,x, con-
sider the substitution
T=1 x = t = -. (2.26)
Tj xi 11
29


The equation (2.24)-(2.25) becomes
^2 = (kVf) + A (Se 5% CTiT) ,
*i dt x\ \ ) \ )'
^ = -SCse~7n, T > Ta.
at
Multiplying through by tx/Tx and tx respectively, we have the equations in
the dimensionless form,
dT t\k ~ /~ ~\ txA -
= -Vv VT J + ASe
dt xf \ ) Ti
dS -JL.
-^ = -hSCse TT^,
dt
T > 0,
ACtiT,
which depends on the six dimensionless coefficients,
h Cs,
(2.27)
2 = d'V (Vf) + c3Se~c*/f cbT,
2 = -c&Se~Cilf, f > 0.
dt
_ tik _Ta txA B
^1 2 ^2 rp i Qj rji ) ^4 rj~) 1 ^5 Cg
^1 Tl Tj 7l
so (2.24)-(2.25) becomes
We choose the scales aq, fi, and 7\ from what we expect the dominant
mechanisms to be and to simplify the equations. It is natural to expect the
temperature scale to be governed by the reaction rate, thus 7\ = B, and so
30


C4 = 1. Next we scale the time so that c3 = 1, hence t\ T\/A = B/A.
Finally, we expect that the spatial scale is governed by the diffusion term
and scale to make c\ = 1, which gives aq = (kt\)^2 = kxl2Bxl2 A~1!2. Then
c5 = ACt\ ACB/A and c$ tiCs = BCs/A.
In conclusion, the scaling (2.26) with
Tj = B, Xl = kl/2B^2A~^2, h = BA~\ (2.28)
transforms (2.24-2.25) into the non-dimensional form
= V (vf) + Se~1/f AT,
dt V /
S = -(3Se-Vf, T> 0,
dt
with two dimensionless coefficients
A = CB,
(2.29)
(2.30)
(2.31)
Therefore, the qualitative behavior of the solution is determined only by
the nondimensional coefficients A and (3, which can be varied independently.
The nondimensional form (2.29-2.30) suggests a strategy for identification
of the coefficients k, A, B, C, Cs first match nondimensional properties of
the traveling combustion wave, such as the remaining fuel fraction, or the
ratio of the width of the leading edge to the trailing edge. As we see from
Figures 2.6 and 2.7, the relationship between nondimensional parameters, A
31


x 10'
0.06
x 10'
Figure 2.6: The fuel fraction remaining after the combustion wave as a func-
tion of dimensionless parameters A and 0.
32


0.05
1 0.01
0.02
0.06
Figure 2.7: Ratio of the leading and the trailing edge of combustion wave
with respect to A and /3 from (2.31). The width is measured as the length of
the part of the wave above temperature 500/C.
33


and ft, and nondimensional properties (e.g. remaining fuel fraction, the ratio
of the width of the leading edge to the trailing edge) of the combustion wave
is approximately linear.
This relation helps to get an initial temperature curve. Once we have
found the initial temperature curve, we can optimize the two curves by con-
structing a distance function between two curves. We used FMINSEARCH in
MATLAB, which minimizes the distance function starting at an initial value
with parameters. The optimized measured temperature curve can then be
found by this process.
In Figure 2.6, the narrow range, which was perturbed by very small
amounts from the optimized nondimensional parameters, of A and /?, was
used in order to have an approximately linear relation. Thus, we see that
the fuel remaining is very small. The width of the wave can be measured
e.g. as the distance of the points where the temperature equals 50% of the
maximum with respect to A and /3, and this measured width is used in Figure
2.7.
The nondimensional traveling wave solution T(t,x), S(t,x) has some
(nondimensional) maximal temperature Tmax, width w, and speed v, while
the data (a measured temperature profile) has the maximal temperature
Tmax, the width ui, and the speed v of the traveling wave. This determines
the scales
T, =
w vw
X\ = , and 11 = ^3.
w vw
34


3.8-,
Figure 2.8: Width of the trailing slope as a function of A and j3. The rela-
tion between the wave speed and (diffusion coefficient)0'5, i.e. \fk. We see
approximately linear relations.
35


The substitution
~ x
x =
and
into the system (2.24-2.25) with the coefficients
A = Ti/ti, B = Tlt C = X/TU Cs P/h, and k = x\/ (Tfh)
which has the desired nondimensional properties as well as the correct max-
imal temperature, width, and speed of a traveling combustion wave.
2.4.4 Related work
A dimensionless system similar to (2.29-2.30) was studied in Weber et al.
[1997] when A = 0, i.e., combustion insulated against heat loss. The speed
of the traveling wave as a function of /? was determined numerically and
by an asymptotic expansion in Weber et al. [1997]. The observed existence
of a traveling combustion wave only for small values of (1 is also shown in
Weber et al. [1997]. When increasing /?, the solution becomes periodic, then
the period doubles in Figure 2.9. We have observed that increasing /? has a
(2.32)
admits the scaled solution
and
36


similar effect for equations (2.29-2.30) when A > 0. Also, we have observed
that a sustained combustion wave is possible only when A is small enough.
The dimensionless parameters e, A, and (3 were introduced in Mercer
and Weber [1995], Mercer et al. [1996], Weber et al. [1997], and Asensio and
Ferragut [2002], as we discussed in Section (2.3.2). The derived dimension-
less equation can still have the ubiquitous cold boundary problem due to a
natural cooling term, depending on how the dimensionless equation was de-
rived (Weber et al. [1997]). By setting the ambient temperature to absolute
zero in the natural cooling term, the cold boundary difficulty can be avoided.
However, none of the paper analyzed with the natural cooling term at their
work. In a different way, the derived dimensionless equation of (Weber et al.
[1997] and Mercer and Weber [1997]) avoid the cold boundary problem in the
heat loss term (Mercer and Weber [1995], Mercer et al. [1996], and Asensio
and Ferragut [2002]), but it still requires extra work in dealing with new
dimensionless parameter e, besides A and (3. Our model only has two di-
mensionless parameters, A and (3. Having fewer parameters in the equations
is a huge benefit for understanding the model.
2.5 Calibration of coefficients in one dimen-
sion
We have found the initial values B = 5.5849 x 104K, and C = 5.9739 x
10-4A-1 by using the values Tj = 670A and Tc = 1200A' in equation (2.22),
37


(a) (b)
(c) (d)
Figure 2.9: Periodic patterns of the solutions by varied (3. (a) (3 = 0.4172,
(b) (3 = 0.4636, (c) (3 = 0.4945, (d) (3 = 0.5718. Each graph shows the
temperature profile by timeframe.
38


1000
Figure 2.10: Time-temperature profile (dotted line) measured in a grass wild-
land fire at a fixed sensor location, digitized from Kremens et al. [2003], and
a computed profile (solid line) from simulation.
and then the value A = 1.5217 x lO1/^-1 from (2.23), using the value
tc 110s from Kremens et al. [2003]. Not every initial condition gives rise
to a traveling combustion wave (Mercer et al. [1996]).
Inspired by Mercer and Weber [1997], we use an initial condition of the
form T(x,t0) = Tce~^x~x^+ Ta, where xq is in the center of the interval
and a = 10\/2m.. This initial condition is smooth. Thus, it does not excite
possible numerical artifacts. It has numerically local support and for a mod-
est a provides ignition sufficient to develop into two sustained combustion
waves traveling from the center.
We then found empirically suitable values of k and Cg that result in
traveling combustion waves, computed the nondimensional coefficients A and
/3, adjusted them using Figure 2.10 (dotted line), and scaled using (2.32) to
39


match the maximal temperature and the width of the wave in Figure 2.10
(dotted line), and the speed of the traveling combustion wave, 0.17m/s,
from Kremens et al. [2003]. There was a small amount of wind in Kremens
et al. [2003], however we have not considered the wind here. The resulting
coefficients are k = 2.1360 x 10~1m2s~1 K~3, A = 1.8793 x 102/fs-1, B
5.5849 x 102K, C = 4.8372 x lO"5^"1, and Cs = 1.6250 x lO-1*-1. The
corresponding nondimensional coefficients were A = 2.7000 x 10-2 and 0 =
4.829 x 10_1.
The computed traveling combustion wave (Figures 2.4, 2.5 and 2.10, solid
line) is a reasonable match with the observation (Figure 2.10, dotted line).
The trailing edge of the computed temperature profile (Figure 2.10, solid line)
was not so well matched but this model is quite limited and other matches
reported in the literature are similar (Balbi et al. [1999]). The real data looks
like the superposition of two exponential decay modes, possibly the fast one
from cooling and the long one from the heat stored in water in the ground.
We have also noted that when the ratio A/Cs increases, the temperature
in the traveling combustion wave increases, increasing the thermal diffusiv-
ity coefficient k increases the width and the speed of the combustion wave
and that the maximum temperature in the traveling wave decreases if Cs
increases. Sufficiently small value of Cs is needed for sustained combustion.
We have noted that the numerical solution by finite differences becomes un-
stable when the ratio k/h, where h is the mesh size, is too small.
40


Chapter 3
Numerical methods for the
PDE model
3.1 Efficient numerical implementation
There are several approaches for solving systems of nonlinear PDEs. We
used a standard finite difference method, with an upwind scheme in the
space variables for stability purposes. The trapezoidal method in time with
Newtons method was implemented. An explicit method, such as Runge-
Kutta, would be much more expensive than the trapezoidal method because
a very small time step is required due to the diffusion and reaction terms. In
Newtons method, the GMRES method with FFT as a preconditioner was
also used.
The discretized systems of equations (2.2-2.3) in space is ^
41


where
H~1----S+Ul + ^) (max(t;x,0)uM + min(t;x, Oju^)
B
(max(vy, 0)ui j + min^, 0)rt^) + Au2je U1~T AC{u\j Ta),
-Csu2e
(3.1)
where u = [*], u\ = T and u2 = S, i = 1, n, j = 1, n at time t.
The quantities u[i and are defined by
uU =
u\ u\ 1
and ut, =
+ u\+1 ~ ui
h
and the quantities ux and u^-
UU =
u{ u{ 1
-------and it,, =
h M
+ ui+1 ~ u\
h
The quantities of u2 , u2i, u2p and u2 j are defined similarly.
Applying the trapezoidal method in time, we get
^n+l 2 (V^(^n+1, ^n+l) T ^pi^ni ^n))
(3.2)
Equation (3.2) is solved for un+, by minimizing the residual r, given by
T ^n+1 (^n T 2 (V^(^n+1) ^n+1) T by Newtons method.
42


Write a equation (3.3) as Fr(un+i) = 0. Then in every step of Newtons
Method, we solve the linear system
Fr(ui) + = o,
(3.4)
for k is iteration number of Newtons Method. Equation (3.4) can be
written as
and
(/-
At difn+l
2 du
)A u = r,
(3.5)
where / ^ = E;+1), rn = -Fr{ukn), Au uk. Then the linear system (3.5) will be solved for Au by using
the GMRES method. Once we have the system of linear equations (3.5), we
use a reduced system by using the Schur complement.
Define
K(nk) =
^ Jtt Jts ^
JsT Jss
(3.6)
where the solution is uk [j]. In case we develop the simple fire model
by adding more variables, the matrix (3.6) will be extended easily, since this
way of implementation is capable of accommodating variables by extending
the block matrix. Because Jss is the diagonal inverse matrix of a diagonal
43


Node numbering
Figure 3.1: Five-point stencil
matrix with the reciprocals of the diagonal entries, we use the reduced matrix
ft1 = Jtt ~ Jts Jss Jst-
In every Newtons iteration, the Jacobian (3.6) should be calculated, and
this process expends many CPU cycles and requires a large memory space.
Thus, we need to find an efficient way to implement (3.6) by taking advan-
tage of the sparsity of the Jacobians in (3.6). Suppose we have a uniformly
discretized square domain so that we have m2 points to each directions. Then
the Jacobian matrix (3.6) has the size 2m2 x 2m2.
Using the standard finite difference discretization to the diffusion term
(Figure 3.1), the block matrix Jtt has a five-point stencil, so it is a pen-
tadiagonal matrix (Figure 3.2). The block matrix Jts is m2 x m2 diagonal
44


Matrix Jn
1 2 m l i-m i-i i ) i + m m-
mm
:








.*



Figure 3.2: Pentadiagonal matrix Jtt by using standard five-point finite
difference discretization.
matrix (Figure 3.3), and Jst, and Jss are also m2 x m2 diagonal matrices
(Figure 3.3).
In order to reduce the memory space, which is occupied by zeros, we store
the Jacobian matrix as Figure 3.4. The matrix in Figure 3.4 has a column
data structure of size m2 x 8. The first five columns of Figure 3.4 are a block
matrix Jtt and last three columns are Jts, Jst, and Jss in order.
The first five columns, except the third column, of Figure 3.4 are off-
diagonals of matrix Jtt The first column of Figure 3.4 has m zeros in the
beginning but the fifth column has m zeros in the end. The second column
of Figure 3.4 has one zero in the beginning but the fourth column has one
zero in the end. The third column of Figure 3.4 is a diagonal of matrix Jtt-
Once we have this column storage of data structure, we perform the matrix-
45


Matrix JTS t J^ nruj Jss
\ 2 ? mi i-m i -1 i + 1 i* m m
1
2
3
m 1
i m
i -1
/
Figure 3.3: Diagonal matrix JTs, Jst, and J55.
vector multiplication and take the inverse of the Jacobian. For example, Jgg
is obtained by taking the inverse of each element of Jss column in Figure
3.4.
This sparse data structure allows us to execute faster simulations. Faster
simulations provide great benefits for getting numerical solutions, and help
us to reproduce the advance of wildfires in real time.
By eliminating the fuel variable, the reduced system (3.5) becomes
R^Aui = rx JTS Jgs r2, (3.7)
Ati2 Jss r2 Jss Jst ri, (3-8)
where r = .
To compute (3.7) and (3.8) by GMRES, we follow the algorithm below:
46


Jn. Jn Jsr Jss
i-rtt 1-1 i+i 1 + m
1 I 0 0 \
2 0
3 0
m 0
m* t


m:-m 0
0
0
m \ 0 0 /
Figure 3.4: Column storage of the sparse Jacobian matrix. Five columns,
labelled as i m, i 1, i, i + 1, and i + m, are from the pentadiagonal matrix
Jtt The last three columns store the diagonal matrices Jts, Jsr, and Jss-
Define Jmaxtrix column storage of data structure in Figure 3.4
b2 = Jmaxtrix(:, 6) [Jmaxtrix(:, 8)]-1 Jmaxtrix(\, 7)
Ri = Jmaxtrix(:, 1 : 5)
Ri = Ri(:, 3) b2
r = ri Jmaxtrix(:,6) [Jmaxtrix(:, 8)]_1 r2 in equation (3.7)
Call GMRES, V\ gmres(Ri,r)
Get V2 = [Jmaxtrix(:, 8)]_1 r2
[Jmaxtrix(\, 8)]_1 Jmaxtrix(:,7) v\ in equation (3.8)
Au= [vi\v2]
47


Update u u Au
The GMRES minimizes a residual of linear equations (3.7) with the FFT
preconditioner. The GMRES algorithm calls a subroutine to multiply a vec-
tor by a matrix, and here we use column storage matrix. After we have
minimized the residual ri, which belongs to iq, we get the residual r2, which
belongs to u2-
Remark 1 Adding more variables, (T, S,Y\t Yn), are easily extended by
getting more block matrices,
Jtt Jts JtYi JTYn
JST Jss JsYi JsYn
JYiT JyiS JyiYi Jyxy
\ -WnT JYnS JYnYi JYnYn
from the Jacobian matrix of (3.6). In this sense, using sparsity of (3.6)
gives us great advantages in the simulation when we deal with large system
of nonlinear partial differential equations with several variables.
3.2 Computational tests
We obtained a simple model, (2.2-2.3), which gave us the well matching
temperature curve to experimental temperature curve, presented in Figure
48


IX I#
Figure 3.5: North to south straight line ignition and a sinusoidal wind with
speed 3m/s after 400s. The dark areas are burnt regions.
49


2.10. In this section, we perform various model tests. First, we investigate the
numerical model in order to see if the discretization model has grid directional
preferences or not. We have tested with similar examples as in Clark et al.
[2004], such as a spot fire without a wind speed and a sinusoidal pattern line
ignition. Second, we observe numerical solutions with variable wind speeds.
Since we have identified the coefficients without the wind term (Section 2.5
and Figure 2.10), it is interesting to explore the effect of varying the wind
speed.
3.2.1 Line ignition with oscillatory wind
In this section, we carry test for robustness. A sinusoidal pattern with line
ignition was tested in Figure 3.5. We see that our discretized systems do not
have any directional preference. Similar numerical experiments are shown in
Clark et al. [2004], and they also do not have any directional preference.
3.2.2 Spot fire without wind
The purpose of this test is to validate the discretized numerical model. A
line ignition with oscillatory wind was compared with analytical solutions in
Clark et al. [2004] to validate the numerical scheme. Since we do not have
an exact solution for equations (2.2)-(2.3), we also performed similar tests as
in Clark et al. [2004].
In this test, we used domain size 256 meters by 256 meters, time step 1
50




m}h'.
^ ;nri
r &'
(a)
6tC ieit-nds
Figure 3.6: Spot fire without wind simulations by Section 2.2 equations (2.2-
2.3). The mesh size is 4m and time step is Is after 600s. (a) is the numerical
fire evolution and (b) is a contour of (a), (c) is fuel consumption in (a).
51


second, and mesh size 4 meters. The fire is ignited in the middle of the do-
main. The region of fuel consumption after 600 seconds is in Figure 3.6c. As
we see, the fire spreads out symmetrically in all directions, and it is symmetric
in all directions. Thus the code does not have a directional preference since
the circular motion of numerical solution in all directions is well performed.
Fuel consumption for a spot fire is presented in Figure 3.6c,. As time
advances, the fuel consumption spreads out in all directions, so that the fuel
is consumed mostly in the middle. From the middle to near the boundary of
the circle is a burnt region, which belongs to the front of the traveling wave.
Near the boundary of the circle, we see a lighter region than in the middle,
which belongs to the tail of the traveling wave. The fire is cooling down due
to the Newtonian cooling at the end.
3.2.3 Effect of the wind speed
In this section, we examine the effect of changing the wind speed on the
numerical solutions. Without wind, we get the circular shape of the fire in
Figure 3.6a and 3.6b. The temperature of the fire goes up with the fuel
consumption and the fire propagates in all direction due to the diffusion.
Finally, the fire cools down because of the cooling. With a small wind, the
wind speed is 0.05m/s to the north east, we see a moon shape propagation
as time advances (Figure 3.7a, Figure 3.7b, and Figure 3.7c). Higher wind
speeds, 0.5m/s and 3m/s, with the same wind direction, give us the half
moon shapes in Figure 3.8a and Figure 3.9a, respectively.
52


(a)
(b)
800 seconds
x(m)
(c)
Figure 3.7: Effect of wind speed 0.005m/s fire simulation of Section 2.2
equations (2.2-2.3) after 800s. Time step is Is and mesh size is 4m. With a
small wind speed the fire also moves more slowly and it is closer to circular
shape, (b) is a contour plots of (a), (c) is a fuel consumption in (a).
53


II .
* II hi
(a)
(b)
400 seconds
x(m)
(c)
Figure 3.8: Effect of wind speed 0.05m/s fire simulation of Section 2.2 equa-
tions (2.2-2.3) after 400s. The numerical solution moves from southwest to
northeast. Time step is Is and mesh size is 8m. (b) is a contour plots of (a),
(c) is a fuel consumption in (a).
54


(a)
(b)
JUG 4*tar Jt
Figure 3.9: Effect of wind speed 3m/s fire simulation of Section 2.2 equations
(2.2-2.3). Fire moves from southwest to northeast. Time step is Is and mesh
size is 16m. (b) is contour plots of (a), and (c) is a fuel consumption of (a).
Note that the sharp front edges in (a).
55


Remark 2 We have also observed different fire shape in Figure 3.9c com-
pared to the results in Clark et al. [2004] The coefficients we obtained are
without a convection term. This means that the calibrated coefficients are
dominated by diffusion and reaction terms in order to match the temperature
profile. These problems are reflected in Figure 3.9c and we easily see the fire
spread out to the whole domain. With a proper spread rate, which means
different coefficients, we may get results closer to Clark et al. [2004]-
3.2.4 Some observations
We observed some interesting experimental results while we were testing our
model. Although the simple PDE-based model was set up and the coefficients
are calibrated without considering the wind speed, we have found that the
model performs reasonably with nonzero wind speed, and we have obtained
different plausible shapes of fire by varying the wind speed in Section 3.2.3.
In reality, a large wind speed provides more oxygen and increases fuel
consumption, while the effect of heat carried away due to wind speed. Thus,
the fuel depletion continues, while the fire advances. However, we have ob-
served that the large wind speed caused less fuel consumption so that the fire
propagated without completely consuming all the fuel (Figure 3.10a-3.10b).
This result is because our model did not implement the interactions with
the reaction rate and it did not count this. Of course, the fire will be extin-
guished if the wind speed is too large, because a large wind speed can also
take away the heat.
56


(a) (b)
Figure 3.10: Wind effect of fuel consumption, (a) Fuel consumption with
domain size 256m x 256m and wind speed 3m/s after 20s. (b) Fuel con-
sumption with wind speed 5m/s. The dark area means fuel consumption
and light yellow is unburnt region. We see that the fire spreads faster with a
higher wind speed.
57


3.3 Numerical convergence tests
In this section, we describe several numerical convergence tests. Because we
do not have an exact solution, we investigate convergence rates by decreasing
both the mesh size and the time step. These tests show us the effect of time
size and mesh size on the numerical solution under the given CFL number,
which implies that given stability region.
In this section, we also observe convergence rate in connection with CFL
and mesh Peclet numbers,
CFL =
v\ At
h
and Pec =
\v\h
k '
The mesh Peclet indicates the ratio between diffusion to advection. Clearly,
the mesh Peclet numbers becomes large or small not only because of the ratio
of advection to diffusion but also because of the mesh size h.
Every test has the domain size 1024m x 1024m and the end time is 32s.
We note that relatively small time steps lead to more accurate numerical
solutions, as reflected by the H-H^ norm in Table (3.1)-(3.2). Relatively small
errors are shown in Table (3.1)-(3.2) under the different CFL numbers but
the equal mesh Peclet number, while large errors are shown in Table (3.3)-
(3.4) under the same CFL number but the different mesh Peclet numbers.
We note that the mesh Peclet in Table (3.1)-(3.2) is relatively smaller than
the mesh Peclets in Table (3.3)- (3.4).
First, we test convergence rates for decreasing time from 3.2s to 0.2s.
58


Mesh Pec = 13.3333 U l2 L(X>
dt CFL error ratio error ratio error ratio
3.2 3.2
1.6 1.6 0.360 4.232 114.746
0.8 0.8 0.090 0.251 1.097 0.259 34.933 0.304
0.4 0.4 0.022 0.252 0.276 0.252 9.131 0.261
0.2 0.2 0.005 0.243 0.066 0.241 2.231 0.244
Table 3.1: Convergence ratio for decreasing time step, wind speed 2m/s, the
mesh size h = 2m.
Each test has wind speed 2m/s, 4m/s and mesh size 2m, 4m. Although
tests have the same CFL number, they have different mesh Peclet since we
want to perform a numerical convergence test in time under the same stability
condition.
Table 3.1 and 3.2 show about 0(h2) convergence. The error was estimated
by
\\U2At-UAt\\a,
where U is a numerical solution of equation (2.2)-(2.3) in Section 2.2 and a
is 1,2, and oo.
Second, we test the numerical convergence behavior, with varying the
mesh size and the time step. Each test has wind speed 0.5m/s, 0.8/ms, and
8m/s. The observed convergence is 0(h1+P) due to a stabilized convection
term, upwind scheme, and trapezoidal method in time. The error was ob-
tained by
\\U2At,2h U&t,h\\a ,
59


Figure 3.11: Errors from Table 3.1 in the logarithmic scale(left), and the
convergence norm ratios(right). Errors from Table 3.2 in the logarithmic
scale(left), and the convergence norm ratio(right).
60


Mesh Pec = 53.3333 Lx l2 Loo
dt CFL error ratio error ratio error ratio
3.2 3.2
1.6 1.6 0.928 7.809 144.669
0.8 0.8 0.236 0.255 2.060 0.263 35.008 0.241
0.4 0.4 0.059 0.252 0.520 0.252 9.0724 0.259
0.2 0.2 0.014 0.246 0.127 0.244 2.229 0.245
Table 3.2: Convergence rate for decreasing time step, wind speed 4m/s, the
mesh size h 4m.
CFL = 0.05 Lx l2 Loo
h dt error ratio error ratio error ratio
32 3.2 5.036 69.780 1357.84
16 1.6 2.389 0.474 39.273 0.562 1015.33 0.747
8 0.8 1.360 0.569 23.596 0.600 586.914 0.578
4 0.4 0.647 0.476 12.677 0.537 438.005 0.746
2 0.2 0.236 0.365 4.944 0.389 184.736 0.421
Table 3.3: Numerical convergence behavior for both the time step dt and the
mesh size h, wind speed 0.5m/s.
where U is a numerical solution of equation (2.2)-(2.3) in Section 2.2 and a
is 1,2, and oo.
3.4 Numerical difficulties
A realistic domain size of 90m x 90m and with different mesh sizes and time
steps were used by Asensio and Ferragut [2002], who had actual numerical
results only with a non-dimensionalized model. For comparison, the Hayman
fire burned 60,000 acres. The experimental design using 6 nested domains
61



I
s
c
Figure 3.12: Errors from Table 3.3 in the logarithmic scale(left), and the
convergence norm ratios (right). Errors from Table 3.4 in the logarithmic
scale(left), and the convergence norm ratio(right).
62


CFL = 0.8 Li L>2
h dt error ratio error ratio error ratio
32 3.2 9.005 64.594 1532.53
16 1.6 6.266 0.695 53.694 0.831 1266.04 0.826
8 0.8 4.217 0.673 40.843 0.760 933.48 0.737
4 0.4 2.424 0.574 26.392 0.646 563.045 0.603
2 0.2 1.027 0.423 12.294 0.465 251.448 0.446
Table 3.4: Numerical convergence behavior for both the time step dt and the
mesh size h. Wind speed is 8m/s.
had 10km on the horizontal gird spacing in Coen [2005b]. We see huge gaps
of domain size between the real fire and experimented numerical solutions
Section 3.2.
Numerical difficulties are caused by the small mesh size and time step
needed to advance the fire. For example, in the thin reaction region, the
resolution needs to be small: the spread rate is known to be approximately
0.3m/s. However, a large domain is necessary at the same time in the sim-
ulation of wildland fires. Thus a refined grid mesh in the combustion area
was used (Dupuy and Morvan [2005]).
The complexities of fire modeling also cause difficulties not only for for-
mulating the fire model, but also in getting numerical solutions. That is the
reason we derived a simple PDE-based on model.
In the numerical experiments, the area of ignition also affects a proper
resolution of the effect of wind speed on the fire. Hence we must choose
an appropriate mesh size and time step for the given domain size and wind
63


speed.
We now explore numerical difficulties related to solving systems of non-
linear PDEs with a realistic domain size. Since such realistic domains are
large, we need mesh sizes to be as large as possible while still maintaining a
good resolution of the fire.
At the experimental level, we first set the domain size as 4000m x 4000m,
the wind speed, 0.005m/s, and the time step, 2s. Then we choose mesh sizes
of 80m, 10m, and 2.5m and the end time is 600s. In these tests, the same
ignition area and fuel gap were applied. We found that the mesh size of 80m
is not enough to propagate the fire (Figure 3.13), and that we require a grid
size at least of 2.5m (Figure 3.14).
When we increase the wind speed to lm/s and the grid size to 80m
(Figure 3.15), a simulated numerical solution does not have numerical prob-
lems. The increased wind speed accelerated the propagation of fire (Fig-
ure 3.15). These results are evident if we consider a basic formula such as
distance = velocity x time, which explains why we cannot avoid numerical
difficulties in experiments.
In the first example in Figure 3.13, we have a mesh size of 80m with wind
speed 0.005m/s. In order to transfer the heat from one node to the next
node, which is far from 80m, the required time is 16000s. Of course, 600s is
not enough to propagate the fire, but this does not mean that the fire should
extinguish. This result gives us a motivation to choose a small mesh size like
2.5m, and 600s is enough to cause the interactions between nodes with 2.5m
64


0 seconds
1500-
g ' |
0)
2 1000^
co
(a)
300 seconds
1400
g 1200
| 1000
CO
800.
V(m) 0 0 x(m)
(b)
Figure 3.13: Fire simulation with large mesh size of 80m and wind speed of
0.005m/s. The time step is 2s, and ignition area is 400m by 320m. (a) is
the initial plot of (b)
65


600 seconds
Figure 3.14: Fire simulation with mesh size of 2.5m after 600s. All other
conditions are the same in Figure 3.13.
300 second*
Figure 3.15: Fire simulation with large mesh size of 80m and wind speed of
1 m/s. While the case of small wind speed (3.13), a traveling fire develops.
66


mesh size. However, we have to pay extra costs with this small mesh size in
computation.
In second example in Figure 3.15, we consider the mesh size of 80m, and
wind speed is 1 m/s. From the formula, we know that the required time is
just 80s and it is clear that there are no numerical problems in Figure 3.15.
67


Chapter 4
The fireline propagation model
In this chapter, we introduce a semi-empirical model. This model postulates
that the fireline evolves with a given spread rate in the normal direction. The
semi-empirical model assumes that the fire propagation speed is normal to
a fireline as a function of wind and terrain slope. Fuel decays exponentially
from the time of ignition, and this model imposes the fuel depletion rate. This
semi-empirical model can be implemented by using the level set method in
the Section 4.1.
In the semi-empirical model, we consider fire burning in the area Du(t)
in the plane with boundary Dr{t), called the fireline, with the outside nor-
mal n = n(rr, y, <), (x,y) Dr(t). The time of ignition ti(x,y) at a point
(x,y) Dii(t) is defined as the time when the point is at the fireline,
(x,y) Dv (U (x,y)).
The spread rate is given by the form (Rothermel [1972], Clark et al.
68


[2004])
0, if R < 0,
< max? if R ^ ^max 1 (4.1)
R, V otherwise,
mm{co,R0 +R4) +Rs}, R# = a (v n)b, Rs = dVz n,(4.2)
where Rq is the spread rate in the absence of wind, Rmax is a maximal speed
of propagation, Co is the minimal spread rate, R is the wind correction, v n
is the component of the wind v in the spread direction n, Rs is the terrain
correction, and Vz n is the slope of the terrain 2 in the spread direction n.
Let Fo (x, y) be the initial fuel supply and W (x, y) be the time constant
of the fuel. Then the heat flux from the fire to the atmosphere is given by
the amount of fuel burned
H =- A (x, y) ^SF (x, y, t),
(4.3)
where
SFo (*, y) if (x> y) e £)()
SF{x,y,t)=l , (4.4)
I SFo (x, y), otherwise
by assuming that the fuel decreases exponentially with time from ignition
while the fuel is burning. The coefficients Co, Ro, Rmax, a, b, d, W, and
A are determined from laboratory experiment, and A characterizes the fuel
69


types and properties. This model was developed in Clark et al. [2004], where
a tracer scheme was used to advance the fireline. Tracer schemes do well
when advecting the shape of the fire in wind, but they require a complicated
code for numerous special cases when the topology of the fireline changes,
such as when fire fronts merge. Also, tracer schemes are not well suited for
data assimilation, because changing the location of the fireline represented by
tracers would require a special, complicated code, and the state of the model
cannot be adjusted by making corrections to gridded arrays, as is usual in
data assimilation methods.
4.1 Level set based wildland fire model
Another possible implementation of the propagation of the fireline is the level
set method. The level set method evolves a function ip = p(x,y,t), called
the level set function, such that the burning area is the level set
Du(t) = {(x,y) :p(x,y,t) < 0}
and the fireline is the level set
Dr{t) = {(x, y) : p (x, y, t) = 0} .
70


The level set function satisfies the differential equation
-jg + fl l|V*>|| = o, (4.5)
where the spread rate R = R(x,y,t) is given by (4.1). Equation (4.5) is
solved numerically.
4.1.1 Some basics about the level set method
A curve Dr = Dr(t) in R2 is represented by the level set
Dr = {{x,y) :tp{x,y,t) = 0} .
The function

speed R(x,y,t), the spread rate, in the normal direction. The speed may
depend on local properties such as curvature and global properties of the
front. The level set method converts the evolution of the curve into a differ-
ential equation for the level set function, which is then solved numerically.
In applications, the curve Dp is closed, and we choose the inside of the curve
to be where p < 0 and it is the area of fire. The gradient gives the direction
normal to the level curves. If Vp ^ 0, we have for the normal direction
(Figure 4.1) that
_ V = IWr
(4.6)
where |*| is the Euclidean norm.
71


Figure 4.1: the normal to the curve(4.6)
The propagation of the curve Dp gives the position of the curve at time
t, which is a fireline at time t since the curve evolves with speed R in the
normal direction. It is known that p is governed by the level set equation
(Osher and Sethian [1988], [Sethian, 1999, Ch. 1]):
^ + R(x,y,t) |V To derive equation (4.7), we consider a point x(t) = (x(t),y(t)) G Dp(t).
Then p(x(t),t) = 0, and from the chain rule,
0
d . .
ji^ dp
w + V,.V)x,
(4.8)
72


where Vtx = ^). Because Vtp = |V dip
+ |'Vy?| n V(x = 0. (4.9)
The speed R in the normal direction is R(x, y, t) = n Vtx. We see that the
level set function p is now governed by the partial differential equation (4.9).
Note that the level set equation requires the propagation speed R(x, y, t)
to be determined at all points (x, y), not just on the curve Dr(t), even if the
curve propagation itself depends only on the speed R given on Dr(t). One
typical example of level set functions is a signed distance function,
p (x) = dist(x, L),
with inside and + outside the curve. The signed distance function satisfies
|V<^| = 1. For stability reasons, the ideal case of

distance function, but this property is not preserved in general by advancing
the level set function in time.
The curve Dp will be also called the interface and the ultimate goal is to
track the motion of this interface as it propagates. An extension velocity is
the velocity field at the fireline in the computational domain. A good choice
of extension velocity field is to take the value at the closest point on the level
set curve, which tends to preserve p close to the signed distance [Osher and
Fedkiw, 2003, p. 27]. There are various approaches in choosing the extension
velocity field, such as extrapolating the velocity to the closest grid point on
73


the front and choosing the fluid velocity in the fluid simulations. However, it
is known that these techniques can cause smoothing of the level set function
in detailed simulations so that rebuilding the signed distance function, called
re-initialization, needs to be performed at each time step in order to avoid
these problems [Sethian, 1999, p. 129].
4.2 Numerical schemes for fireline propaga-
tion model
4.2.1 Numerical schemes for the Level Set equation
The solution of the Hamilton-Jacobi equation (4.7) can be approximated by
several numerical schemes. A simple choice of a numerical scheme for dis-
cretization of G(V method, which is a central difference approximation in space, [Osher and
Fedkiw, 2003, p. 50],
V-
G = G(


7yV +

)-**(

'z )-ty(
V>- VvV
),
where and Vx are numerical approximations of the partial derivative Vx,

Ax
J rr+ and VT Ax
74


and similarly for Vy coefficients, defined as
e1 = max |Gi(Vxv?, V Here G\ and G2 are partial derivatives of G with respect to ipx and

spectively. Then we have the advection equation in two dimensions, G\ = VX
and G2 = Vy, where V(x,y) (Vx, Vy). Although the Lax-Friedrichs method
is a less diffusive scheme, it still creates spurious oscillations at the bound-
ary and causes considerable diffusion into the solution. Another conservative
method is artificial diffusion with upwinding, but it also causes smoothing of
sharp corners [Sethian, 1999, p. 56].
The Hamiltonian of equation (4.7) in one dimension is H(ipx) R (x) \^px\-
By letting u ut + [H(u)]x = 0 (4.10)
is a conservation law in one dimension. Then the Hamiltonian H(u) =
R(x)\u\ is convex since > 0. This Hamiltonian is approximated by
a numerical flux function and the system of conservation laws with the flux
function requires a positive definite matrix for well-posedness (Friedrichs and
Lax [1971]), and thus requires a convex flux function. Construction of ap-
propriate flux functions which preserves the conservation equation (4.10) can
give us more accurate solutions.
75


In the numerical solution of the fireline propagation model (4.7), we used
the Godunov scheme, which is an upwinded finite difference approximation
in space, and Heuns method, which is a Runge-Kutta method of order 2, in
time. Heuns method is given by
n+1
(pn + AtG( G( = tpn +
where G = R(x,y) |Vy?| in equation (4.7), and the approximation of the
gradient is computed by the Godunov method [Osher and Fedkiw, 2003, p.
58],
VxV?
V*y? if V+y> < o and Vxip < 0,
if V+y> > 0 and V~y> > 0,
< 0 if V~

0
otherwise if |V~y?| > |V*y>,
if i v-^i < i i,
(4.11)
where Vy? = [Vxy>, Vyy?] and similarly for Vj,y>. Furthermore,
|V.
76


X(m)
Figure 4.2: A circular fire propagation with constant wind. The contour is a
fireline.
4.2.2 Observations and future work
A numerically stable scheme that includes upwinding, such as the Godunov
method (4.11), is required to properly compute the norm of the gradient
ip. When the gradient in the level set equation is approximated by central
differences, the numerical method fails quickly.
Heuns method is a second-order Runge Kutta method, with second-order
accuracy. In the time discretization, Heuns method is required, because of
the conservation of law in equation (4.7). Heuns method is known to avoid
the overestimation or underestimation of (p by averaging the two.
From two simple experiments, such as a circle movement along a wind
velocity field and a rotating velocity field, we found a shrinking of the cir-
cle. In these simple experiments, we considered the spread rate R (x, y, t) of
77


t=0
(a)
E.
>
t=3
\\\\\\\\\\
VV\\\
v.\\\\
\\\\
jk w
w
W wWWWWW v x\
\\\\\\\\\\\\\\ x x\
\\\
X(m)
t=3
x\\\\>W\\\\
. \ v v v x'x 'x
\\\W x\)C\^ Vxx x xxx
\\\\\\\\\^XX xVx
WVWxx X X
E.
x\\V'
-x\\\'-
wwwwwww
wwwwwww-
X(m)
(b)
(c)
Figure 4.3: (a)
circle, (b) The
keep shrinking.
A trace of circle movement along wind direction. Initial value
circle moved to the south-east. As time is passing, the circle
The re-initialization avoids the shrinking problem in (c).
78


t=5
t=5
X(m)
X(m)
(a) (b)
Figure 4.4: A trace of motion in a rotation velocity field. Re-initialization
method is added in (a) and (b).
the equation (4.7) to be only the rotating wind velocity field, unlike equa-
tion (4.1). Thus, the shrinking of the circle should not appear. Shrinking
of the circle did not appear in short time experiments, but eventually this
problem appeared. Thus, the initial circle (Figure 4.3a) went away (Figure
4.3b). While we are taking a rotating velocity field with the initial circle, we
observed that the circle could not move without loss of mass (Figure 4.4a).
It turns out that preventing mass loss is quite an interesting subject in the
study of level set methods (Enright et al. [2002], Enright et al. [2005]). Several
approaches have been developed, such as hybrid particle level set methods,
the quadrature-free discontinuous Galerkin method, and flux-based level set
methods (Enright et al. [2005], Frolkovic and Mikula [2007], and Marchandise
et al. [2006]).
In order to resolve the mass loss problem, a quick remedy is a re-initialization
79



-"
- r'

-A1


A B A
B / c / '
aT a
(a) (b)
Figure 4.5: (a) Points on the zero level set on the mesh grid lines are stored
as a set A and points on the zero level set function, which is not on the mesh
grid lines, are stored as a set B. (b) Closest points to the grid points are
stored as a set C.
scheme which simply reconstructs the signed distance function from a new
front every few steps (Figure 4.3c and 4.4b), and re-initializes the entire do-
main with the new front. Re-initializing the level set function can maintain
the signed distance function property every few steps, and it also helps to
avoid the mass loss problem in a diagonal movement (Figure 4.3c). The
algorithm of re-initialized level set method is listed below (Figure 4.5).
Find mesh grid points in the zero level set and store those points as set
A.
Find grid points where the zero level set and grid lines intersect and
store those points as set B.
Find the closest points to the grid points on the segments, which con-
nect points in set B and A in the mesh cells, and store those points as
set C.
80


Find the minimum distances from the points in the set A, B, and C to
each grid point and store the minimum distances to each grid point.
Construct the signed distance function for the grid points with stored
minimum distances.
With a large domain, this re-initialization method is expensive computa-
tionally since every point has a distance to each grid point, and this distance
needs to be calculated with every point in the set A, B and C. The Narrow
Band scheme, which is one of the most efficient schemes, can be used to re-
duce computation time. The level set function depends only on points of
the narrow band. The Narrow Band implementation is done on points of the
narrow band along the curve (see Adalsteinsson and Sethian [1995] and Peng
et al. [1999]). The points of the narrow band are found from the zero level
set at each time step, and then we only update the values on those points.
Lastly, we re-initialize the level set function. The choice of the proper width
of the narrow band is still an issue, but it is known that choosing six grid
points from the zero level set is reasonable based on experiments ([Sethian,
1999, p. 56]). This scheme results in a much faster method and also makes
extending the velocity of the fireline cheaper and easier since it is extended
to fewer points.
Note that the re-initialization method is also not enough to resolve the
mass loss problem (Figure 4.6) in a rotating flow.
81


!=7.8
X(m)
(a)
(b)
:'r;s t Th:,rBMd
^t,,e circie is ^ -itua;bLeTotr:ie
82


Chapter 5
Phase-space analysis of the
traveling wave
In nonlinear systems, it is hard to find an analytical solution. Even with
an analytic solution, such a solution is often far too complicated to gain
insight into its properties. The phase-space method allows us to determine
the qualitative behavior of some solutions.
Analysis of the traveling wave in phase space yields a simplification of
PDEs, and the solution is a geometrical interpretation of the state variables
in phase space.
In this chapter, we first show that a traveling wave solution of PDEs is
equivalent to a solution of ODEs in phase space, and some basics of dynamics
and general properties of dynamical systems will be discussed. Lastly, we give
inequalities for the speed of the traveling wave and the existence of a traveling
83


wave.
5.1 Traveling wave in phase space
In this section, we briefly introduce the phase-space method and we obtain
the parameterized system of ODEs in phase space in order to investigate the
traveling wave solutions of the combustion PDEs from Chapter 2.
The phase-space method is one approach to find solutions of an au-
tonomous system, where the system does not explicitly depend on the in-
dependent variable. The method introduces additional variables, that allow
us to reduce the equations to nonlinear ODEs. This approach allows us to
analyze the solutions of the ODEs, which is a dynamic system, instead of the
PDEs.
Phase space consists of vectors formed by variables. The solution in the
phase-space is a curve, which is called a trajectory or phase curve.
A two-dimensional phase-space analysis of a reaction-diffusion equation
was introduced by Grindrod [1991], and an application to a fire model similar
to ours was discussed by Weber et al. [1997], assuming no heat loss and zero
absolute zero ambient temperature, making the equations easier to analyze.
They found the lower bound of the traveling wave speed numerically, and
the existence of the traveling wave speed was also investigated numerically
in the phase-plane. Here we do a three-dimensional phase-space analysis with
heat loss and no diffusion fuel properties, which leads to a more complicated
84