Exploring a Convection–Diffusion–Reaction Model of the Propagation of Forest Fires: Computation of Risk Maps for Heterogeneous Environments

: The propagation of a forest ﬁre can be described by a convection–diffusion–reaction problem in two spatial dimensions, where the unknowns are the local temperature and the portion of fuel consumed as functions of spatial position and time. This model can be solved numerically in an efﬁcient way by a linearly implicit–explicit (IMEX) method to discretize the convection and nonlinear diffusion terms combined with a Strang-type operator splitting to handle the reaction term. This method is applied to several variants of the model with variable, nonlinear diffusion functions, where it turns out that increasing diffusivity (with respect to a given base case) signiﬁcantly enlarges the portion of fuel burnt within a given time while choosing an equivalent constant diffusivity or a degenerate one produces comparable results for that quantity. In addition, the effect of spatial heterogeneity as described by a variable topography is studied. The variability of topography inﬂuences the local velocity and direction of wind. It is demonstrated how this variability affects the direction and speed of propagation of the wildﬁre and the location and size of the area of fuel consumed. The possibility to solve the base model efﬁciently is utilized for the computation of so-called risk maps. Here the risk associated with a given position in a sub-area of the computational domain is quantiﬁed by the rapidity of consumption of a given amount of fuel by a ﬁre starting in that position. As a result, we obtain that, in comparison with the planar case and under the same wind conditions, the model predicts a higher risk for those areas where both the variability of topography (as expressed by the gradient of its height function) and the wind velocity are inﬂuential. In general, numerical simulations show that in all cases the risk map with for a non-planar topography includes areas with a reduced risk as well as such with an enhanced risk as compared to the planar case.


Introduction
It is the purpose of this contribution to apply a recently-developed efficient numerical method [1] for the solution of a wildland fire model [2] to simulate the propagation of a wildfire in various spatially heterogeneous environments. In particular we demonstrate the use of the method for the computation of so-called risk maps that are based on solving the model under systematic variations of the initial focus of wildfire. The governing model is a particular system of convection-diffusion-reaction partial differential equations (PDEs) that is specified in Section 2.
General principles for the design of models of wildland fire are outlined, for instance, in [3,4]. The model employed herein, due to Asensio and Ferragut [2], is described in detail in [1]. Further references to this model and its variants include [5][6][7]. Applications of this model to real-world wildfire scenarios in several specific geographic regions are documented in [4,8,9]. It is also the basis of the spectral algorithm advanced by San Martin and Torres [10,11], but the particular nature of that numerical method is applicable to a constant diffusion coefficient only. An alternative approach to the description of wildfires through partial differential equations, and their numerical solution by explicit methods is provided in [12,13]. These models and simulators are all based on the principles of combustion theory (see, for instance, [14,15]). Under simplifying assumptions one then obtains a fully three-dimensional reaction-diffusion-convection system governed by thermodynamics [16]. Such systems are usually reduced to two horizontal space dimensions based on the consideration that the ratio of height of vegetation with respect to the characteristic size of a forest is small [17]. Thus, although the combustion process is fully turbulent and three-dimensional, most models used for the simulation of wildfire propagation are formulated and eventually solved in two horizontal space dimensions, although some three-dimensional effects, for instance due to variation of topography, are considered. These include, in particular, [1,2,[4][5][6][7][8][9][10][11][12][13]. This simplifaction, especially for large terrains, also reflects limitations in computational resources available. In this context we refer to Tables I and II in [10] for an overview on models, their assumptions, and the corresponding computational treatment. Explicit treatments of the vertical propagation of a forest fire (but on domains of limited horizontal extension) include [18,19]; see also the references cited in these works.
The simulation of a wildfire by numerical solution of an initial-boundary value problem for the governing convection-diffusion-reaction model is a well-known challenge [20] due to the presence of a diffusion term and stiff reaction terms. As a consequence of this property, explicit finite difference schemes are usually associated with a stability condition (CFL condition) that enforces very small time steps for moderate to fine spatial discretizations. This restriction can be avoided by implicit-explicit (IMEX) time discretizations that were introduced for the present wildfire model in our previous paper [1]. These methods impose a less restrictive limitation of the time step than explicit methods. This is achieved by carefully distinguishing between stiff and non-stiff occurrence of the vector of unknowns in the spatially-discretized equations, and choosing the time discretization of the respective terms in a corresponding either implicit or explicit fashion. However, IMEX methods are more general than semi-implicit methods and are based on interlacing evaluations of the stages of particular explicit and implicit Runge-Kutta (RK) schemes [21,22]. We refer to [23,24] as the earliest references to IMEX methods in the context of time-dependent partial differential equations and [25][26][27][28][29][30][31][32][33][34][35] for various applications of IMEX-RK schemes.
IMEX schemes can be understood as sophisticated time discretizations to be applied to the spatial discretization of a partial differential equation, that is, to the governing equation semi-discretized (discrete in space, continuous in time) in a so-called method-of-lines formulation. The spatial discretization itself needs to be done separately. To this end, we herein employ a fifth-order weighted essentially non-oscillatory (WENO) finite-difference discretization [36][37][38] of convective terms on Cartesian grids that can also be regarded as a second-order fully conservative finite-volume discretization on these grids. Our preference for this technique with respect to others, such as finite volume schemes with flux limiters (see [36,39,40] for these and many other alternative schemes) stems from our familiarity with it, in particular in the context of convection-diffusion models related to eco-epidemiological problems [30,31] and non-local models of pedestrian flow [32]. In this work, discontinuous solutions may arise due to degenerate diffusion, as well as sharp gradients, that are dealt with in a robust manner by the WENO reconstructions. Articles that use some mathematics, to make the techniques used more plausible, also include [41].
From a scientific point of view, the contribution and novelty of the present paper is threefold. Firstly, we complement the description of an efficient new numerical method (provided in detail in [1]), namely a linearly implicit IMEX time discretization involving Strang-type operator splitting and WENO spatial discretization by additional examples that demonstrate that the method is flexible enough to handle nonlinear or degenerating diffusion functions and variable topographies. Secondly, we obtain some qualitative insight into the predictions made by the underlying wildfire model. Thirdly, it is demonstrated how the availability of an efficient numerical simulation method can be employed to construct wildfire risk maps for portions of a given terrain.
The remainder of the paper is organized as follows. In Section 2 the governing mathematical model is summarized. We only introduce the model in its final, dimensionless form. Detailed derivations of the wildfire model and its ingredients are provided in [1,2,4]. Next, in Section 3 we outline the numerical method. Numerical results are presented in Section 4. Specifically, after introducing preliminaries in Section 4.1, we simulate in Section 4.2 four scenarios (Scenarios 1.1 to 1.4) with different definitions of K(u). Then, in Section 4.3 we fix the diffusion coefficient K(u) and evaluate the influence of the topography of the terrain that affects the wind velocity w through its gradient. To this end we simulate Scenarios 2.1 to 2.4 that differ in the topography of the terrain, in particular the steepness of mountains. Then, based on these topographies, we proceed in Section 4.4 with the computation of risk maps, where the specific risk associated with a point or small patch of the computational domain is the smallness of the time needed to burn a determined fixed amount of fuel (here measured as a percentage of the fuel initially available in Ω) by a wildfire having its initial focus at that position. Some conclusions are collected in Section 5.

Summary of the Mathematical Model
The governing model is the system of convection-diffusion-reaction partial differential equations (PDEs) where t is time, x ∈ Ω is the spatial variable where the domain Ω ⊂ R 2 represents the forest in which the fire may propagate, and u = u(x, t) and v = v(x, t) are the scalar unknowns. Here u is the non-dimensionalized temperature and v the non-dimensionalized mass fraction of solid fuel. Moreover, K = K(u) is a given diffusion coefficient, and w is an advection velocity that represents wind speed. The functions f (u, v, x) and g(u, v) are the reactive part of the model. (These ingredients will be specified further below.) The complete wildfire model is described by the system (1) along with zero-flux boundary conditions where n is the unit normal vector to ∂Ω, and the initial conditions We here state the model ingredients in final form, based on the assumption that the dimensional rate r of the chemical reaction of fuel and oxidants into products is given by the Arrhenius equation r = A exp(−E A /(RU)), where A is a constant, E A denotes the activation energy, R is the universal gas constant, and U is absolute temperature (values of parameters are specified in Section 4.1). We mention that precise statements of the assumptions involved in the model development are provided in [1,2,4].
In particular, it is assumed that there is always sufficient amount of oxidant provided by air, so that the conservation law for the oxidant can be ignored (as is stated in the Section 2 of [2]).
If a reference ambient temperature U ref = U ∞ and the non-dimensional inverse of the activation energy ε = RU ∞ /E A are given, then the non-dimensional temperature is u = (U − U ∞ )/(εU ∞ ). Furthermore, time t and the spatial coordinates x and y within x = (x, y) are understood as dimensionless variables. Here the original dimensional quantities are non-dimensionalized by the time scale t 0 = (ε/(q react A)) exp(1/ε), where q react is a non-dimensional reaction heat, and the length scale l 0 = (t 0 k/(ρC)) 1/2 , where k is thermal conductivity, C is specific heat, and ρ is the density of the fuel.
Next, we assume that the wind velocity w is given by where the vector field w 0 (x, t) is a given wind velocity that may depend on spatial position and time (but is considered constant in our numerical experiments), and T (x) is the topography of the domain [10,11]. This choice differs from that of [1], where the terrain was assumed flat, and allows us to include spatial heterogeneity. Furthermore, we assume that spatial propagation of heat occurs through radiation (as described through the Stefan-Boltzmann law) as well as through natural convection. Then the diffusion function K(u) is given by where κ = 4σδU 3 ∞ /k. Here σ is the Stefan-Boltzmann constant and δ is the length of the optical path for radiation through the substance. (The particular functional form (5) is the one derived from first principles in [1,2,4]. In particular, it is based on writing the heat flux by radiation (as defined by the Stefan-Boltzmann law) in dimensionless form as However, for illustrative purposes and since our numerical methods are able to handle a wide range of nonlinear and even degenerate diffusion coefficients, we will also consider alternative algebraic expressions, see Section 4.2). The reaction term f (u, v, x) is given by Here the function c(u) indicates the phase change between the endothermic (solid) and the exothermic (gaseous) phase, where u pc is a given non-dimensional phase change temperature, and α denotes a natural convection coefficient. Finally, the function g(u, v) that accounts for fuel consumption takes the form Summarizing, we may collect the model in its final form and relations necessary in the remainder of the paper as follows. The governing partial differential equations are (1) that are solved along with the boundary conditions (2) and the initial condition (3). As for the coefficient functions in (1), the wind velocity w(x, t) is specified by (4) (specific choices of the topography function T (x) will be introduced in Section 4.3), the diffusion function K(u) by (5) (alternative formulas will be discussed in Section 4.2), and the reaction term f (u, v, x) by (6) and (7). The fuel consumption function g(u, v) is given by (8). Values of the constants involved in K(u), f (u, v, x) and g(u, v) are provided in Section 4.1.

Numerical Method
We assume (for simplicity) a quadratic domain Ω = [0, L] 2 and denote by u = u(x, t) and v = v(x, t) the solution of (1). We define a Cartesian grid with nodes x i = (x i , y j ), i, j = 1, . . . , M, with x i = y i = (i − 1/2)∆x, ∆x = L/M, and we define the index vector i = (i, j) ∈ M := {1, . . . , M} 2 . The unit vectors e 1 = (1, 0) and e 2 = (0, 1) are used to refer to neighboring grid points, for instance x i+e 1 = (x i+1 , y j ) and x i+e 2 = (x i , y j+1 ). We define u as a solution computed at an instant t in the grid points, where u i (t) = u(x i , t), and analogous notation is used for v. For simplicity we denote K i = K(u i ). We may then approximate (1), (5)-(8) in semi-discrete form (that is, in discrete in space but continuous in time form) by the system of ordinary differential equations (ODEs) Here the terms C(u) and D(u)u represent the spatial discretizations of the convective and diffusive terms arising in the first PDE in (1). Their entries are given by C is the numerical flux corresponding to the fifth-order WENO spatial discretization of the convective term ∇ · (wu) [35]. For the boundary conditions (2), we can apply absorbing or reflective boundary conditions in each dimension by imposing ghost nodes.
We discretize (10) implicitly in time for v and explicitly for u. The resulting scheme can be written as In order to evaluate (11) in each iteration step, we first compute the value of v n+1 , then this value is used to compute u n+1 . Now, we associate with (11) the solution operator ψ ∆t defined by The remaining part of the system (9), namely the system of ordinary differential equations is discretized by a linearly implicit IMEX-RK method. Full technical details of this procedure are provided in [1], and are similar to treatments of related convection-diffusion problems [25][26][27][28].
To explain the main idea, let us recall that in principle Runge-Kutta (RK) ODE solvers could be applied to solve (9) numerically. For instance, strong stability preserving (SSP) explicit RK schemes are a popular class of time integrators associated with a favorable stability constraint on the time step ∆t [25][26][27][28][29][33][34][35]. Alternatively, once could employ implicit-explicit Runge-Kutta (IMEX-RK) methods (see [33][34][35]), for which only the diffusion term is treated implicitly. In this case the stability condition on ∆t is less restrictive than for explicit discretizations, but a large system of nonlinear algebraic equations must be solved in each time step. This shortcoming of IMEX-RK methods is avoided by the methodology proposed in [25,26] that is based on linearly implicit-explicit RK schemes, and which is applied here to the discretization of (1). The idea is to distinguish the terms arising in (9) between stiff and non-stiff dependence on the solution vectors u and v. One then chooses the time discretization by an implicit and an explicit RK scheme for the terms involving stiff and non-stiff dependence, respectively. In the product D(u)u, the dependence on u within D(u) is considered non-stiff, while that of the factor u is considered stiff. On the other hand, in the product vζ(u), the term ζ(u) is considered non-stiff, while the term v is considered stiff. Based on this distinction, the resulting linearly implicit IMEX-RK method for (12) is defined by a pair of RK schemes, namely an explicit one (ERK scheme) and a diagonally implicit one (DIRK scheme) that handle the non-stiff and stiff dependencies, respectively, and one alternates between evaluations of the stages of the ERK scheme and solving linear systems to to evaluate those of DIRK scheme. The details of this procedure are outlined in [1], and we employ here the same coefficients, namely the IMEX-RK scheme H-LDIRK3(2,2,2) defined by a pair of particular two-stage ERK and DIRK schemes [25,27]. The resulting discretization of (12) over a time step of length ∆t can be written as where u n+1 is the approximation of (12) by the IMEX-RK method. Then the complete Strang splitting method (cf. [39,42,43]) to solve (9) is formulated as Once again we refer to [1] for details on the numerical method.

Preliminaries
The values of parameters used in all examples are the universal gas constant R = 1.987207 cal/(K mol), an ambient temperature of U ∞ = 303 K, an activation energy E A = 20 kcal/mol = 83.68 kJ/mol which yields ε = 0.02980905 ≈ 0.03, and A = 10 9 s −1 [15]. The mean magnitudes of other constants, for the numerical examples in [2], and which we adopt for our numerical experiments in [1] as well as herein, are Furthermore, Asensio and Ferragut [2] chose the heat of combustion H in such a way that one can assume q react = 1. Utilizing the parameters (13), we then get At the moment we do not have access to the specific value of U pc . However, we may estimate the maximal temperature u max that is attainable. To this end, we examine the ODE system (10), meaning du/dt = −(q react /ε)dv/dt.
On the other hand, the second equation in (10) For instance, assume that at some point in the spatial domain we impose the initial temperature such that u(0) = 30, which corresponds to U = (1 + 30ε)U ∞ = 1.9U ∞ = 570 K. This assumption implies that u max = (1/ε) + 30 = 63.3, or equivalently, a maximum temperature (in absolute value) of U max = (1 + εu max )U ∞ = 870 K. In light of this calculation we will choose either u pc = 0 or 0 < u pc ≤ u(0) < u max , where u(0) is the dimensionless temperature at the point where the fire starts.
In all numerical examples (with the exception of the accuracy test, Scenario 2.0) we employed the method introduced as S-LIMEX in [1] and assume a spatial domain Ω = [0, L] 2 with L = 200. The uniform spatial discretization is chosen as ∆x = 1, so that the computational domain Ω is subdivided into 200 × 200 cells.

Scenarios 1.1 to 1.4: Numerical Experiments with Various Diffusion Coefficients
In this group of scenarios, we explore the influence of the dynamics of the combustion model. To this end we fix the initial distribution of temperature and that of fuel, and assume that the terrain is flat and homogeneous, but vary the diffusion function K(u) to illustrate the effect of different model assumptions. We assume that the initial distribution of fuel is constant, as well as that the topography is flat. Specifically, the initial datum for temperature corresponds to an initial focus on a small square subdomain, i.e., we set With these initial data we simulate a base case and three variants of it. The base case, Scenario 1.1, is based on using the same remaining parameters and model functions as in [1], namely K(u) defined by (5) with κ = 0.1 and ε = 0.035, f (u, v, x) given by (6) with c(u) defined through (7) with u pc = 20 and α ≡ 10 −3 , and g(u, v) given by (8) with q react = 1.
We are interested in comparing results with those obtained for alternative definitions of the diffusion coefficient K(u). To this end, we repeat the simulation in Scenario 1.2 with the same parameters with the exception that κ is five times larger, i.e., κ = 0.5. Next, to motivate two alternative choices of the function u → K(u), we recall first that the maximum temperature u max can be estimated from the ODE version of (1), (5)- (8), that is (10).
The result is that when we solve (10) for t > 0 starting from u(0) and v(0), then u(t) ≤ u max := (q react /ε)v(0) + u(0). If we assume that in our case the relevant values are v(0) = 1 and u(0) = 40, then u max = 1/0.035 + 40 ≈ 68.57 =: u * max . Thus, K(u) assumes values between K(0) = 1.1 and K(u * max ) = 4.9304. If we wish to compare results with those obtained for a constant diffusivity (as stipulated in [10]) D∆u instead of ∇ · (K(u)∇u), then a reasonable value of the diffusion coefficient D for comparison should be Finally, we introduce the possibility of a degenerating diffusion coefficient, that is we allow that K(u) = 0 for isolated values of u or even a u-interval of positive length. For instance, we may assume that the mechanisms of heat transfer are active only wherever u exceeds a critical value u crit . Applying this idea to the diffusion function K(u) of Scenario 1.1, that is (5), we obtain For Scenario 1.4, we utilize (14) with u crit = u pc = 20. Figure 1 shows a plot of the diffusion functions for Scenarios 1.1 to 1.4, and Figure 2 displays the corresponding numerical results. Roughly speaking, Scenario 1.2 predicts a wildfire that has consumed a larger portion of fuel, and has travelled a larger distance, than that of Scenario 1.1. Far from the reaction front the maximum temperature in the combustion zone is slightly smaller. With a constant diffusion coefficient (Scenario 1.3) we obtain a less wide shape of the burnt area (in comparison with Scenario 1.1), and a smaller value of the burnt portion of fuel. Finally, as is to be expected with a strongly degenerating diffusion coefficient, the lateral flanks of the temperature surface of Scenario 1.4 are steeper than for all other scenarios.

Scenarios 2.0 to 2.4: Effect of the Variability of Topography
In Scenario 2.0 we provide a comparison with a reference solution to show the numerical convergence of S-LIMEX scheme and in Scenarios 2.1 to 2.4 we study the net effect of terrain topography on the dynamics of wildfire propagation. To this end, we fix the diffusion function K(u) as given by (5), and choose all other parameters as in Scenario 1.1 as well with the exception of ε = 0.02, with constant vector w 0 = (2.5, 2.5) of velocity of the wind which blows in the north-east direction (that is, in the direction of increasing xand y-coordinates). However, we assume various cases for the topography function T (x) = T (x, y). To this end, and following [10], we define a shape function s(x, y; γ) := exp −(x 2 + y 2 )/γ with a parameter γ > 0 that describes a peak of height 1, centered at the origin. We assume that there are two peaks in the domain of variable height, as expressed by For Scenario 2.0 we chose the parameters (h 1 , h 2 , x 1 , x 2 , y 1 , y 2 , γ 1 , γ 2 ) = (60, 60, 55, 25, 25, 55, 400, 60) and three alternative choices of these parameters, in Scenarios 2.1 to 2.4, respectively, which are specified in the caption of Figure 3 that illustrates the variants of topography functions.  In Table 1 we compute the approximate L 1 -errors e M (u) and e M (v) at three simulated time T = 0.2, T = 0.3 and T = 0.4 with different discretizations by using M = 40, 80, 160 and 320, with respect to reference solution computed with S-LIMEX scheme with M = 640. We observe that the approximate L 1 -errors decrease as the grid is refined, which can be also observed in Figure 4 where we display the numerical approximations for M = 160 and M = 640 at simulated times T = 0.2 and T = 0.4. We observe that the approximation obtained with M = 160 is a already a good approximation of the reference solution.  Table 1).
We display the numerical approximation obtained at simulated times T = 5, 10, 15, 20, and 25. We display the evolutions of the temperature and fuel over a flat domain and topography T 1 in Figure 5, and show the analogous results for topographies T 2 and T 3 in Figure 6. In the cases T 1 , T 2 and T 3 the topography acts as repulsion, which moves away the propagation fires ( Here and in Figure 6, the first and third row show temperature and the second and fourth row show fuel, and small white square indicates the location of the initial wildfire.

Scenarios 3.1 to 3.4: Risk Maps
An analysis that is interesting when studying the effects of topography on the source of fire in a forest fire is the elaboration of the so-called risk maps, which we present in this section. To do this we consider four cases here, the first corresponds to the absence of the variable T where we are in the presence of a perfectly flat terrain and three other cases which correspond to each of the topographies T 1 , T 2 and T 3 of the previous section. Each of the risk maps constructed corresponds to square initial ignitions of length 12 placed on D ij := (x, y) ∈ Ω | |x − (55 + 10i)| ≤ 6, |y − (55 + 10j)| ≤ 6 , i, j = 1, . . . , 5, which yields a total number of 25 initial foci. The risk maps indicate the different values obtained for the time T risk when the focus of the fire takes to consume 5% of the total fuel available in the domain, so a short time corresponds to a high risk and a long time corresponds to a low risk. In each of the four cases (of a flat domain and a domain with one of the topographies T 1 , T 2 and T 3 ) the ignitions are located in the south-west of the original domain. Specifically, to assign a value of T risk to each of the patches of size 10 × 10 D 0 ij := (x, y) ∈ Ω | |x − (55 + 10i)| ≤ 5, |y − (55 + 10j)| ≤ 5 , i, j = 1, . . . , 5, we simulate the wildfire model starting from the initial condition where D ij is specified (16) and the velocity is given by w = w 0 in the flat case, w = w 0 + ∇T i , i = 1, 2, 3 for the non-flat topographies, and w 0 = (2.5, 2.5). The result is in each case a risk map for the domain [60, 110] × [60, 110], with a (fairly coarse) resolution of 25 patches of size 10 × 10.
In the case of flat terrain the propagation of fire and burnt fuel always have the same spatio-temporal evolution but shifted according to the location of the initial focus ( Figure 7) and 5% of the total available fuel is consumed, the time in the 25 foci is identical and it corresponds to T risk = 6.45 and the risk map obtained which can be seen in the top left plots of Figure 8 corresponds to a high-risk scenario (but the risk is the same for each 10 × 10 patch).  In the case of topography T 2 , the minimum value T risk = 6.04 is reached on [60, 70] × [100, 110] and the maximum value T risk = 6.77 that is reached on [80, 90] 2 . The evolution of temperature and fuel can be seen in Figures 11 and 12

Conclusions
In the present work we analyzed the wildland fire model in [5] to simulate the propagation of a wildfire in various spatially heterogeneous environments. In particular, we focus on variable, nonlinear diffusion functions to explore the dynamics of the combustion model and a variable topography function to analyze the effects of the spatial heterogeneity. This model was solved numerically in an efficient way by the numerical method proposed in [1].
Several diffusion functions are explored (Scenarios 1.1 to 1.4) to arrive at the conclusion that these functions are strongly related to the shape of the burnt area. The main result is that increasing diffusivity (with respect to a given base case) significantly enlarges the portion of fuel burnt within a given time while choosing an equivalent constant diffusivity or a degenerate one produces comparable results for that quantity. Considering three different topographies, we were able to show how the terrain affects the spread of the forest fire, acting as an element that manages to increase fuel consumption (Scenarios 2.2 and 2.4) but also as an element that manages to vary the direction of fire spread acting as a repulsion (Scenario 2.3) which pushes it toward areas far from its initial trajectory (Scenario 2.2). On the other hand, due to the shapes of the landscape, in two of these cases it is possible to capture how the elevation of the ground manages to accelerate the spread of fire (see Appendix B in [44]).
Considering that a determining factor in the spread of fire is the value of the slope associated with the terrain, in the sense that the higher the slope, said spread occurs more quickly, it is essential to obtain the times in which a certain constant portion of the fuel has been consumed through the risk maps (see Figure 8). Finally, the possibility to solve the base model efficiently may be used as a tool to elaborate wildfire risk maps.
The verification and calibration of the model taking into account data from real wildfires in four of the authors' geographic environment (Chile) is still an open issue, while scenarios related to other geographic regions are successfully simulated in [4] (Germany) and [8,9] (varios regions of Spain). Clearly, realistic scenarios can be considered if geographical information and meteorological data are integrated, which can be achieved through a Geographical Information System (GIS). We refer to [8] for details and mention that we are currently making the effort to access similar information for Chile. On the other hand, and as we mentioned in [1], another omission is the neglect of moisture and pyrolysis effects that can be handled via a multivalued enthalpy function [5,6].