Mathematical Modelling of Running Crown Forest Fires

. Adapted mathematical model for simulation of running crown forest ﬁre propagation is considered. Simplifying assumptions, equations of the model, initial and boundary conditions, ﬁnite diﬀerence approximations are introduced. The results of computer modelling and the peculiarities of forest ﬁre behaviour in heterogeneous forests are discussed


Introduction
Interest in mathematical models of forest fires is caused by a large number of scientific difficulties.It is noted in [3,5,8,11] that forest fire models have been developed since 1940 to the present, but a lot of chemical and thermodynamic questions related to fire behaviour are still to be resolved.Forest fires are divided into underground (peatbog) fires, surface fires, active crown fires, running crown fires (also called independent crown fires), and mass fires [5,8].A number of researchers notice, that running crown fires possess the greatest speed of propagation.They are extremely dangerous and very difficult to fight, thus mathematical modelling is represented as an important problem [5,8,11].The state-of-the-art of forest fire modelling is reviewed in detail in [11].It's indicated that there are only a few ( [11], table 7) successful computer software for crown forest fire calculation, based on the crown fire spread empirical model of Rothermel and the Forestry Canada Fire Danger Group's semi-empirical model.It is possible to note several weaknesses of existing forest fire models: a lack of theoretical justification of accepted mathematical descriptions of physicochemical reactions and transformations; unavailability of analytical reference solutions of boundary value problems; the need of efficient numerical approximations and high speed algorithms.The current state of IT opens the opportunity to develop and formulate forest fire mathematical models, which take into account processes of heat and mass transfer, radiation, chemical reactions and transformations.Such models require significant computing power, they involve quite lengthy calculation time even on high speed computers, i.e. direct application of such codes in real-time systems is impossible.The practical use of forest fire mathematical models in decision support systems can be realized in two directions: filling of a knowledge base [1] of computed forest fire propagation forecasts or creation of simplified (empirical) models for real-time calculations.
The most important theoretical model of the crown fire spread is developed by Grishin ([11]).It is based on the fundamental laws of physics, conservation and theoretical justifications are provided.In the present article we are going to present the implementation of this model in the software complex.Theoretical and methodological foundations of running crown forest fire model realization and of according knowledge base formation are described.We suggest Grishin's mathematical model extension and adaptation, refined descriptions of some model's components.New method of finite difference approximation for the boundary value problem is used and results of the numerical experiments are discussed.The strongest change in environmental condition parameters occurs in the zone called fire front, which propagates at some speed across the territory, covered with forest, and is visually observed as a zone captured by a flame.As a rule the effects of smoke generation over a big territory and clouds formation over a fire zone, because of condensation of the water vapor released during combustible forest materials (CFM) burning, take place.
Describing the geometry of forest fire propagation the following properties are usually mentioned -in a horizontal direction: a region of the burnt area, a burning zone, not burnt area; in a vertical cut: a bottom forest layer, a flame, a forest canopy (tree crowns), boundary layer of the atmosphere ( [7,8,9,10]).Among physical and chemical processes of CFM burning, the following processes are defined: CFM heating, drying, pyrolysis and burning, oxidation of the gaseous, condensed and disperse products of pyrolysis, aeration of soot particles from coke, aeration of ashes particles and a smoke formation.The general scheme of change of an aggregate state and a chemical composition in a forest fire zone is described, for example in [7], it includes the heat supply due to convection, heat conductivity and radiation causes change in condition and composition of CFM.During combustion are formed and present in gaseous and dispersed conditions: carbon monoxide (CO), hydrogen (H 2 ), methane (CH 4 ), carbon dioxide (CO 2 ), water vapor (H 2 O), soot particles (C).In a solid phase during generation of coke, which includes carbon (C) and a mineral part of CFM, are taking place: formation of ashes because of oxidation of coke and smoke particles in gaseous and dispersed phases, generation of soot particles due to pyrolysis and aeration.reactionary properties of the forest canopy, temperature and the sizes of the ignition centre are known.
The following assumptions are usually made: environment is considered to be a five-phase porous medium, consisting of dry organic substance, dispersed water, solid pyrolysis product (coke), ash and a gas phase; the gas phase consists of oxygen, combustible pyrolysis product components, inert air components and also water vapor and inert products of burning; the gradient of temperature across the forest canopy is small in a comparison with a temperature gradient in a longitudinal direction; influence of Coriolis force and centrifugal force is small in comparison with gravity ( [7]).The mathematical model taking into account the listed assumptions is described by the system of equations, concerning which the author notices its complexity and that practical application of such a model is not always possible.
In this paper, the model of running crown forest fires propagation is studied with following additional assumptions: • pressure is considered constant; • wind velocity in the forest canopy basically depends on forest phytocenosis structure parameters and strongly depends on coordinates and characteristics of the forest fire front itself, it is accepted to be equal to the equilibrium velocity, calculated by the formula (4.2.7) [7]; • differences in thermal (and diffusive) streams between the top and bottom borders of a forest canopy are approximated by the Newton's formulas [7,12]; • heat inflow into a forest canopy due to radiation from a flame torch is negligible.
Then the one-temperature mathematical model of a running crown fire can be written in the form of the following system of PDEs: ) (2.10) Here t is the time, x is a coordinate in the system of coordinates connected with the center of an initial fire, the x-axis is directed towards unperturbed wind velocity parallel to a horizontal spreading surface; the coordinate z is counted upwards from a spreading surface; V = (u, v, w) is equilibrium wind velocity vector; T is temperature (in Kelvins), T ∞ is unperturbed ambient temperature; λ T is turbulent thermal conductivity; ϕ j , j = 1, 2, 3, 4 defines volume fractions of the multiphase reactive medium, where ϕ 1 corresponds to dry organic substance, ϕ 2 is water in liquid-drop state combined with CFM, ϕ 3 is coke (condensed pyrolysis product), ϕ 4 is mineral part of CFM (ash); ρ j , j = 1, 2, 3, 4 is the j-th phase density; ρ 5 is the density of a gas phase (a mix of gases), ρ ∞ is unperturbed density of a mix of gases (air density); c ν , ν = 1, 2, 3 is mass concentration of components of a gas phase, where c 1 corresponds to oxygen (O 2 ), c 2 to combustible gases (combustible pyrolysis product components), c 3 to mixes of other gases (inert components of air, water vapor, inert products of reactions of pyrolysis, coke burning and of combustible gases oxidation); c 1∞ and c 2∞ are mass concentrations of oxygen and combustible gases in unperturbed atmosphere; M ν , ν = 1, 2, 3 are molecular masses of gas phase components; M C is molecular mass of carbon, M ∞ is molecular mass of air; R 1 , R 2 , R 3 are mass rates of reactions of dry CFM pyrolysis (chemical decomposition of substance by heating with allocation of combustible gases), of moisture evaporation from CFM (drying), coke burning; R 51 , R 52 , Q are mass rates of generation (disappearance) of oxygen, combustible gases, gas phase; R 5 is mass rate of reaction of burning (oxidation) of combustible gases; c pj , j = 1, 2, 3, 4 is the j-th phase thermal capacity; c p5 is thermal capacity of a gas phase; q 2 , q 3 and q 5 are heat effects of processes of evaporation, of burning of the condensed fuel and of gaseous combustible pyrolysis products accordingly; D T is the diffusion coefficient; ∆h is crown height (∆h = h 3 − h 2 , where h 3 and h 2 are heights of the top and bottom borders of forest canopy accordingly); α is the coefficient of heat exchange between the atmosphere and a forest canopy; α c is the coke number of CFM; κ R defines integrated absorptance; σ is the Stefan-Boltzmann constant, Equations (2.1)-(2.10)are written in a form which differs from the traditional one (see, [7]).In our model some changes are done and some members of the equations are grouped together to identify the physical processes described.The presented form of equations (derivative in time, convection, diffusion, and the right-hand side) simplifies the understanding of finite difference approximations which are used below.The numerical algorithm is implemented in the computer program using Wolfram Mathematica system.The accepted form of the equations is also convenient because it makes more understandable where and what nonlinearities take place.

Initial and Boundary Conditions
Boundary conditions are set as follows ( [7]): (3.1) ) Here P is the forest fire zone, −∞ and +∞ indicate burnt and unburned areas removed from fire zone to adequate distance (unperturbed values c 1∞ , c 2∞ and T ∞ ).The processes of reactions, phase transitions, physical and chemical transformations in a part of forest untouched by fire processes require the following conditions ( [7]): where T * is the given value of temperature.In other words, processes described in the model do not start, until the temperature reaches a certain level.For a given simulated forest region it is possible to define the density of dry organic substance ρ 1 , the bulk density ρ 0 of typical layer of CFM and its moisture content W = (m 0 − m 1 )/m 1 , where m 0 and m 1 are mass of CFM in natural and absolutely dry conditions accordingly.Then the initial values of volume fractions are calculated by the formulas: where ζ defines ash content of combustible forest materials.According to [7], the values of ash content are in the range 0.001 < ζ < 0.01.It is noticed, that the influence of ash content on modelled processes is weak, therefore further we will accept ϕ 4 = 0. Also, we assume that ϕ 3H = 0, since in an unburnt zone coke is not formed yet.Using the assumption that CFM are completely burnt down, relations for final volume fractions can be described by the formulas: For the case of one-dimensional process, when x-axis coincides with the wind direction, the initial distributions T 0 , c 10 , c 20 , ϕ 10 , ϕ 20 , ϕ 30 , which are used in Math.Model.Anal., 15(2):161-174, 2010.computational experiments, are presented in Fig. 1.Such initial distributions of temperature, volume fractions and component concentrations are confirmed by the results of physical modelling of crown forest fires propagation, and also by the results of calculations on the "self-coordinated mathematical model" given in [7].The formulated initial conditions correspond to characteristics of forest fires in case when ignition has occurred in the zone of the finite size, and further burning extends by different possible scenarios: a switch to the established regime or to an extinction of a fire.

Corrected Mathematical Model
To complete the system of equations it is necessary to write down the dependences describing the speeds of CFM pyrolysis, drying, coke burning and chemical reactions in a gas phase.For example, according to [7,9]: where k 01 , k 02 , k 03 are pre-exponential factors of chemical reactions ( [5,17]), E 1 , E 2 , E 3 is the activation energy of chemical reactions, R is universal gas constant, s σ defines specific surface of the condensed product of pyrolysis (of coke), ν Γ is the proportion of gaseous combustible pyrolysis products.An important component of the formulated mathematical model is the expression for gaseous pyrolysis products burning speed R 5 .There is no uniform approach on this question in the scientific literature.Often it is assumed that it is sufficient to consider only one reaction [7,9,17]: since it brings the greatest energy release during the forest fire.
In [12] the following features of the mentioned reaction are marked: "Dry (without water vapor) CO + O 2 reacts with a low speed.However, speed of reaction becomes fast when a small amount of water vapor or hydrogen is added to a mixture.Then radical OH, atomic hydrogen H and oxygen O serve as the primary active centers.Reaction of continuation of a chain with simultaneous CO 2 generation is the reaction CO + OH → CO 2 + H" and "Practically during combustion of natural fuels water vapor and hydrogen are always presented at mixes and provide high speed of CO burning reaction".Thus the formulas describing speed of reaction (4.5) require a specification.To describe the speed of reaction (4.5) in [7] two different expressions are offered in various papers.They essentially differ.A possibility of their practical use raises the doubts, moreover the analysis of formula (5.4.19) in [7] shows the discrepancy of measuring system units in the parts describing the same reaction at various concentrations of oxygen.
In [9] the kind of the formula and the value pre-exhibitors of reaction (4.5) are taken from [12], but the value of activation energy is taken from [6,7].The formula is given as: Our mathematical model consists of equations (2.1)-(2.10),(4.1)-(4.4),boundary conditions (3.1)-(3.8),and also it includes the formula for R 5 .If there is not enough oxygen for the total burning out of all carbon monoxide, using formula (4.6) for the description of reaction R 5 gives physically incorrect result, since the concentration of oxygen can reduce below zero level.At any (including very small) values c 1 , speed of burning of combustible gases and destruction of the oxygen, described by expression (4.6) can be very high.
The computational experiments made for the fitting of mathematical model enabled us to derive a formula for reaction (4.5).Using a method of stationary concentrations [17] for chemical reaction of carbon monoxide burning (4.5), concentration c 2 can be expressed through stationary concentration c 1 as follows: The given formula means that for the total combustion of carbon monoxide (CO), mass concentration of which equals c 2 in some volume, the required mass concentration in the same volume of oxygen (O 2 ) is defined by the expression (4.7).Using (4.7) the adaptation of formula (4.6) is proposed in which process of carbon monoxide burning goes according the "excess/deficiency" principle: It should be noted that such a model is not used in any known published work on the theory of combustion; the delimitation of its applicability requires further studies.Performed computational experiments indicate that the use of formula (4.8) allows us to avoid reduction of oxygen concentration below zero, and at the same time it adequately models the running crown fire propagation limiting conditions by the wind velocity, CFM density and moisture content.

Finite Difference Approximation of the Model
The boundary value problem formulated above is solved numerically using the explicit finite difference approximations: Equations for concentration and temperature are approximated as follows: The given approximations describe the one-dimensional problem; the generalization for a multi-dimensional case is obvious because explicit difference schemes are applied.The values of functions at nodes i + 1/2, i − 1/2 are calculated by the formula of the two-point upstream weighting ( [14,15]).
It is well-known that explicit schemes are only conditionally stable.In order to get unconditionally stable schemes the implicit approximations should be used, see e.g.[2] where porous media mathematical models are applied to simulate wood drying processes.
Methodical calculations were conducted to justify and optimize the proposed method of solving the boundary value problem.The impact of independent members of quasi-linearization methods and grid selection on the accuracy of the discrete solution were examined.Two-point upstream weighting formula was used to improve the accuracy of computed solutions near fronts.Comparison of numerical solutions obtained with two-point, central-and backwarddifference approximation formulas showed that it is preferable to use the twopoint upstream weighting scheme given in [14] for the approximation of forest fire mathematical model.We note that such schemes are most efficient for solution of problems of physicochemical methods for enhancing oil recovery, see also [13], where the analysis is done of upwind and high-resolution schemes for numerical solution of convection dominated problems in porous media.
In order to fit the mathematical model, results of computational experiments of crown forest fire propagation were compared with the experimental results of the Institute of Forest and Wood of the USSR Academy of Sciences and of Tomsk State University published in scientific and technical literature, and by comparing simulation results with Van Wagner's experiments, conducted in natural conditions at the plantations of "red pine" [7,16].In particular, the calculations reproduced conditions in a zone of crown fire, including the characteristics of fire front dynamics during establishing to a regime of steady propagation.Computer modelling also reproduced with high accuracy the speed of fire front for all cases examined, corresponding to different wind velocities.
The use of explicit approximations imposes severe restrictions on the time step size ( [4,14]).Control of numerical instability and misconvergence (strong nonlinearities of right hand sides in the equations for temperature and oxygen concentration) was performed by analyzing the graphics of solutions.Application of numerical experiments conduction technology using toolkit developed for knowledge base creation and support [1] simplified this analysis greatly.

Modes of Running Crown Fire Propagation and Conditions of its Extinction
The analysis of the numerical results calculated on proposed model (2.1)-(4.8)shows that two various scenarios of forest fire behaviour are possible: selfextinction and steady propagation.Let's present and comment the solutions of 1D problem of the dynamics of fire propagation in the direction of the wind (towards the x-axis).It is assumed that the parameters of forest and therefore the parameters of fire front in the direction of y-axis (perpendicular to the wind velocity) are homogeneous.At the first scenario, fire extinction occurs,  Fig. 1-3 illustrate the dynamics of calculated solutions in the case of extinction scenario.After forest fire occurrence it propagates some distance in the wind direction, the maximum temperature decreases due to heat loss and after reaching the critical value the burning stops (see, Fig. 2).Areas of the low concentration of oxygen and nonzero concentration of combustible pyrolysis components "separate" from the stopped fronts of pyrolysis and drying, these areas are extending in a direction of a flow, concentrations influenced by convection are gradually approaching the corresponding values in the unperturbed atmosphere.
Another scenario is also possible, when fire propagates in a direction of a wind without extinction.Calculated solutions show the areas of typical behaviour, corresponding to arguments stated above about the geometry of fire, i.e. in the wind direction presence of different zones can be traced: burnt out, burning, unburnt.In a burning zone it is possible to allocate a site of intensive increase of temperature, also a relatively narrow zone where the temperature is maximum and a site of intensive decrease in the temperature, caused by heat loss for CFM heating and drying.To the left of the zone of maximum of temperature there is coke, to the right there is CFM.On the site where temperature increases coke burns down, dry CFM is pyrolysised, released gases burnt out.In the vicinity of the maximum value of temperature the combustible pyrolysis product components concentration c 2 intensively increases (at the same site volume fraction of dry CFM considerably decreases), oxygen concentration c 1 due to oxidation of combustible pyrolysis products falls down close to zero, there's not enough oxygen for carbon monoxide (CO) total combustion.Further on the site where temperature decreases the concentration c 2 relatively slowly goes down to initial level.Profiles of solutions (temperature T , multiphase reacting medium volume fractions ϕ j , j = 1, 2, 3, 4, mass concentration of gas phase components c ν , ν = 1, 2, 3) quite rapidly takes the asymptotic form.Further in the case of uniform distribution of CFM layer, profiles of the density and moisture content in the unburnt area are transferred with minor curvature alterations in a wind direction with constant rate of spread ω n .Solutions for the different moments of time are shown in Figs.4-5.Let's notice, that similar distributions (u ∞ = 8 m/s, W = 66.66 %, ρ 0 = 0.2 kg/m 3 , ∆h = 5 m, T = 300 K, ρ 1 = 500 kg/m 3 , ρ 2 = 1000 kg/m 3 , ) are received at calculations on "the self-coordinated model" [7].The calculated profiles, the fire rate of spread ω n = 0.625u ∞ are consistent with the results given in [7].7 The Boundaries of Scenarios in the {ρ, W, u} Space Using the described computer model the question on boundaries, switch points of scenarios of fire process-steady propagation or extinction were studied.Whenever possible it is necessary to localize the switch line, when one scenario belongs to parameters on one side, another on the opposite.The study involves the calculation of a big number of variants in a wide range of determining parameters and the analysis of the obtained results in hyperspace (CFM density, moisture content, the equilibrium wind velocity) of sequence of points.
Using the toolkit, developed by authors on the base of Wolfram Mathematica kernel, for creation and support of knowledge bases of mathematical models [1] allows to do series of numerical experiments in the automated mode and to analyze the obtained results.

Conclusions
Let us present some conclusions, interpretations of results confirmed by presented graphics.In Fig. 6, regions of steady fire propagation are schematically illustrated on a plane {wind velocity, moisture content} at three various CFM layer density values.We note, that at any fixed equilibrium wind velocity there is some critical moisture content value W > W * at which fire propagation stops due to the big losses of thermal energy on CFM heating and drying.In practice, when fighting or localizing the combustion by spraying water, the conditions for fire extinction by moisture content parameter are artificially created.The obtained results of computations show, that to create fire extinction conditions it is enough to increase the moisture content value above 0.7.The artificial creation of limiting conditions of propagation are also possible by taking special values of two other parameters (forest fuel breaks and buffer zones, CFM burning before combustion front, directed explosions).The graphics presented, in particular, show that, for example, at equilibrium wind velocity 5 m/s and CFM density less than 0.1 kg/m 3 the steady propagation of running crown fire is impossible (shown in Fig. 6 as vertical lines).Definition in hyper-plane {ρ, u} of approximate position of propagation or extinction scenarios boundary is possible for any moisture content value.In order to predict the critical value of moisture content it is necessary, for instance, to do series of calculations, changing the CFM initial density, and fixing equilibrium wind velocity, then the change of the temperature distribution at the successive time moments is controlled.
Let us separately notice, that the data can be used to predict the dangerous combustible situations in a drought season on a concrete blown site of forest.
Of particular interest is the study of the dynamics of fire spread in piecewise homogeneous medium.Corresponding conditions describe the practical challenges of forest fire fighting, such as the availability of forest buffer strips, i.e. strips free of forest vegetation (fuel breaks, rails, roads, watercourses), and fire barriers, i.e. strips with artificially increased water content, as well as strips of trees of non-combustible sorts [3,5].In particular, the numerical experiments allow to get the width of a fuel break or fire barrier (and its moisture content) required for fire extinction.For example, in a typical spruce forest (ρ 0 = 0.2 kg/m 3 , ρ 1 = 700 kg/m 3 , ∆h = 15 m) in case of wind velocity u ∞ = 10 m/s and moisture content W = 40 % the availability of fuel break with width less than 4 m leads only to a temporary decrease in temperature of a crown fire, after passing the fuel break temperature of a fire front increases again, the fire goes back to a regime of steady propagation.When the fuelbreak's width is 4 m (or more) the fire extincts, the effect of partial combustion of branches (the generation of small amount of pyrolysis products ϕ 3 ) on the other side of a fuel break takes place.
Let us note, that using the described model critical values that characterize the width of a strip with increased water content (for instance, as a result of water dropped from plane) can be easily calculated.

Figure 1 .
Figure 1.Initial distributions of volume fractions of the multiphase environment and concentration of gas phase components.The temperature is shown by a continuous thick line; numbers 1, 2, 3 denote the curves corresponding to dry organic substance of combustible forest materials, water in an liquid-drop condition, coke; 4, 5 illustrate concentration of oxygen, carbon monoxide.

Figure 2 .
Figure 2. Dynamics of temperature decrease at moments t = 0, 0.5, 1, 2, 2.5, 3 and 4 s.distributions of mass concentrations c 1 and c 2 of gas phase components and of temperature T are approaching to corresponding unperturbed values c 1∞ , c 2∞ and T ∞ .

Figure 3 .
Figure 3.Typical distributions of characteristics of forest fire propagation at extinction scenario at moments of time t = 1, 2, 3 s.

Figure 4 .
Figure 4. Typical profiles of forest fire characteristics at steady propagation scenario at moments of time t = 1, 3 s.

Figure 5 .
Figure 5. Profiles of temperature and oxygen concentration at steady forest fire propagation at moments of time t = 1, 2, 3, 5, 7.5 and 10 s.