Next Article in Journal
Numerical Solutions of Fractional Differential Equations by Using Fractional Explicit Adams Method
Next Article in Special Issue
A Spatial-Temporal Model for the Evolution of the COVID-19 Pandemic in Spain Including Mobility
Previous Article in Journal
Rapid Consensus Structure: Continuous Common Knowledge in Asynchronous Distributed Systems
Previous Article in Special Issue
On the Existence of Solutions of a Two-Layer Green Roof Mathematical Model
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
CI2MA and Departamento de Ingeniería Matemática, Facultad de Ciencias Físicas y Matemáticas, Universidad de Concepción, Casilla 160-C, Concepción 4030000, Chile
2
Departamento de Silvicultura, Facultad de Ciencias Forestales, Universidad de Concepción, Casilla 160-C, Concepción 4070374, Chile
3
Departament de Matemàtiques, Universitat de València, Av. Vicent Andrés Estellés, E-46100 Burjassot, Spain
4
GIMNAP-Departamento de Matemática, Facultad de Ciencias, Universidad del Bío-Bío, Casilla 5-C, Concepción 4051381, Chile
5
CI2MA, Universidad de Concepción, Casilla 160-C, Concepció n 4030000, Chile
*
Author to whom correspondence should be addressed.
Mathematics 2020, 8(10), 1674; https://doi.org/10.3390/math8101674
Submission received: 20 August 2020 / Revised: 22 September 2020 / Accepted: 24 September 2020 / Published: 1 October 2020
(This article belongs to the Special Issue Modeling and Numerical Analysis of Energy and Environment)

Abstract

:
The propagation of a forest fire 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 efficient 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) 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. In addition, the effect of spatial heterogeneity as described by a variable topography is studied. The variability of topography influences the local velocity and direction of wind. It is demonstrated how this variability affects the direction and speed of propagation of the wildfire and the location and size of the area of fuel consumed. The possibility to solve the base model efficiently 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 quantified by the rapidity of consumption of a given amount of fuel by a fire 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 influential. 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.

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

2. Summary of the Mathematical Model

The governing model is the system of convection–diffusion–reaction partial differential equations (PDEs)
u t + · w ( x , t ) u = · K ( u ) u + f ( u , v , x ) , v t = g ( u , v ) ,
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
u w K ( u ) u · n = 0 , ( x , t ) Ω × ( 0 , + ) ,
where n is the unit normal vector to Ω , and the initial conditions
u ( x , 0 ) = u 0 ( x ) , v ( x , 0 ) = v 0 ( x ) , x Ω .
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 / ( R U ) ) , 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 ε = R U / 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
w ( x , t ) = w 0 ( x , t ) + T ( x ) ,
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
K ( u ) = κ ( 1 + ε u ) 3 + 1 ,
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
q = κ ( 1 + ε u ) 3 u = K ( u ) 1 u .
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
f ( u , v , x ) = v ζ ( u ) α u , ζ ( u ) : = c ( u ) exp u 1 + ε u .
Here the function c ( u ) indicates the phase change between the endothermic (solid) and the exothermic (gaseous) phase,
c ( u ) = 1 if u u pc , 0 if u < u pc ,
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
g ( u , v ) = ε v q react ζ ( u ) .
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.

3. 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)
d u / d t = C ( u ) + D ( u ) u A u + v ζ ( u ) , d v / d t = ( ε / q react ) v ζ ( u ) .
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 ( u ) = ( C ( u ) i ) i M , D ( u ) u = ( ( D ( u ) u ) i ) i M , where
C ( u ) i = l = 1 2 1 Δ x f ^ i + 1 2 e l f ^ i 1 2 e l , D ( u ) u i = l = 1 2 1 2 Δ x 2 ( K i + K i e l ) u i e l ( K i e l + 2 K i + K i + e l ) u i + ( K i + K i + e l ) u i + e l ,
where f ^ i + 1 2 e l is the numerical flux corresponding to the fifth-order WENO spatial discretization of the convective term · ( w u ) [35]. For the boundary conditions (2), we can apply absorbing or reflective boundary conditions in each dimension by imposing ghost nodes.
To summarize the strategy to convert (9) into a fully discrete method, let us consider first the ODE system that results from (1), (5)–(8) by omitting heat transport and diffusion terms, namely
d u / d t = v ζ ( u ) , d v / d t = ( ε / q react ) v ζ ( u ) .
We discretize (10) implicitly in time for v and explicitly for u . The resulting scheme can be written as
u n + 1 = u n + Δ t ζ ( u n ) v n + 1 , v n + 1 = v n 1 + Δ t ( ε / q react ) ζ ( u n ) .
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
ψ Δ t ( u n , v n ) = ( u n + 1 , v n + 1 ) .
The remaining part of the system (9), namely the system of ordinary differential equations
d u / d t = C ( u ) + D ( u ) u A u ,
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
φ Δ t ( u n , v n ) = ( u n + 1 , v n ) ,
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
( u n + 1 , v n + 1 ) = ψ Δ t / 2 φ Δ t ψ Δ t / 2 ( u n , v n ) .
Once again we refer to [1] for details on the numerical method.

4. Numerical Results

4.1. 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 k cal / 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
ρ = 100 kg m 3 , C = 1 kJ kg 1 K 1 , k = 1 W m 1 K 1 .
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
t 0 = 0.03 10 9 s 1 exp ( 1 / 0.03 ) = 8986.8 s , l 0 = ( 0.089868 m 2 ) 1 / 2 = 0.2998 m
(see the formulas in Section 2). In the numerical examples we keep h constant and small, such that α = 10 3 . The inverse of the conductivity coefficient is chosen as κ = 0.1 , and we choose non-dimensional wind velocities in a different way in each example.
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
d u / d t = ( q react / ε ) d v / d t .
On the other hand, the second equation in (10) implies that d v / d t 0 , while v ( t ) 0 if v ( 0 ) 0 . Thus we conclude that
u ( t ) = ( q react / ε ) v ( t ) + ( q react / ε ) v ( 0 ) + u ( 0 ) ( q react / ε ) v ( 0 ) + u ( 0 ) = : u max .
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.9 U = 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.

4.2. 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
u 0 = 21 for ( x , y ) [ 34 , 46 ] 2 , 0 elsewhere
along with
v 0 ( x , y ) = 0.9 for all ( x , y ) Ω .
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
D = 1 u max * 0 u max * K ( s ) d s = 1 + κ 4 ε u max * ( 1 + ε u max * ) 4 1 2.3816 .
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
K ( u ) = κ ( 1 + ε u ) 3 + 1 for u > u crit , 0 for u u crit .
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.

4.3. 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 x- and 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
T ( x , y ) = h 1 s ( x x 1 , y y 1 , γ 1 ) + h 2 s ( x x 2 , y y 2 , γ 2 ) .
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.
For Scenarios 2.1 to 2.4, the initial conditions of temperature and fuel are
u 0 ( x , y ) = 21 χ 9 , 21 2 ( x , y ) , v 0 ( x , y ) = 0.6 χ Ω ( x , y ) .
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 (Figure 5 third and fourth row) or centers the propagation on the diagonal of the domain (Figure 6 first and second row) respectively.

4.4. 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 i j : = ( x , y ) Ω | x ( 55 + 10 i ) | 6 , | y ( 55 + 10 j ) | 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 i j 0 : = ( x , y ) Ω | x ( 55 + 10 i ) | 5 , | y ( 55 + 10 j ) | 5 , i , j = 1 , , 5 ,
we simulate the wildfire model starting from the initial condition
u 0 ( x , y ) = 21 χ D i j ( x , y ) , v 0 ( x , y ) = 0.6 ,
where D i j 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).
For the topography T 1 , all the values obtained for the time in which 5% of fuel consumed is reached are in less time than T risk = 6.45 being the minimum value T risk = 4.92 that is reached on [ 100 , 110 ] 2 and the maximum value T risk = 5.98 that is reached on [ 60 , 70 ] × [ 100 , 110 ] . The evolution of temperature and fuel can be seen in Figure 9 and Figure 10 and the risk map obtained, which is displayed on the top right of Figure 8, corresponds to a slightly higher risk to the case of flat a terrain.
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 Figure 11 and Figure 12 and the risk map obtained which can be seen in the bottom left part of Figure 8 corresponds to similar risk to the flat a terrain.
For topography T 3 , the majority of values are higher than T risk = 6.45 being the minimum value T risk = 5.85 that is attained on [ 90 , 100 ] × [ 60 , 70 ] and the maximum value T risk = 7.90 that is reached on [ 80 , 90 ] × [ 100 , 110 ] . The evolution of temperature and fuel can be seen in Figure 13 and Figure 14 and the risk map obtained which can be seen in the bottom right part of Figure 8 corresponds to a low-risk scenario.

5. 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].

Author Contributions

Conceptualization, R.B. and E.G.; Data curation, D.I., P.M. and L.M.V.; formal analysis, R.B., D.I., P.M. and L.M.V.; investigation, R.B. and E.G.; methodology, R.B., E.G., D.I., P.M. and L.M.V.; software, R.B., E.G., D.I., P.M. and L.M.V.; supervision, R.B.; validation, E.G., D.I., P.M. and L.M.V.; visualization, E.G.; writing—original draft, R.B.; writing—review and editing, E.G. All authors have read and agreed to the published version of the manuscript.

Funding

R.B. is supported by CRHIAM, project ANID/FONDAP/15130015 and Fondecyt project 1170473. D.I. acknowledges CONICYT scholarship CONICYT-PCHA/Doctorado Nacional/2014-21140362. P.M. is supported by Spanish MINECO project MTM2017-83942-P and Conicyt (Chile), project PAI-MEC, folio 80150006. L.M.V. is supported by Fondecyt project 1181511. R.B., D.I. and L.M.V. are also supported by CONICYT/ PIA/Concurso Apoyo a Centros Científicos y Tecnológicos de Excelencia con Financiamiento Basal AFB170001, and by the INRIA Associated Team “Efficient numerical schemes for non-local transport phenomena” (NOLOCO; 2018–2020). We are grateful to one anonymous referee of this work for the suggestion to include Reference [41].

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following symbols and abbreviations are used in this manuscript:
List of symbols:
Latin symbols:
Apre-exponential factor [s−1]
Cspecific heat [kJ/(kg K)]
c ( u ) phase change function [-]
C ( u ) discretization of convective term [-]
D ( u ) discretization of diffusive term [-]
E A reaction energy [ kcal/mol]
f ^ i + 1 2 e l numerical flux in the l-coordinate [-]
f ( u , v , x ) reactive term for u [-]
g ( u , v ) reactive term for v [-]
i = ( i , j ) index vector [-]
kconductivity coefficient [W/(mK)]
K i approximate value K ( u i ) [-]
K ( u ) diffusion coefficient function [-]
Llength of domain x or y coordinate [-]
l 0 length scale [m]
Mnumber points x or y coordinate [-]
q react reaction heat [-]
rArrhenius expression [-]
RUniversal gas constant [cal/(K mol)]
T topography term [-]
ttime variable [-]
t n time discretization [-]
t 0 time scale [s]
uaverage temperature [-]
u pc phase change temperature [-]
u max maximum temperature [-]
Uabsolute temperature [K]
U ref ambient temperature [K]
u n vector of components u i n
u i n approximate value u ( x i , t n ) [-]
vmass fraction of solid fuel [-]
v n ,vector of components v i n [-]
v i n approximate value v ( x i , t n ) [-]
wvector field [-]
w 0 wind speed [-]
x = ( x , y ) spatial variables [-]
x i = ( x i , y j ) nodes in the Cartesian grid [-]
Greek symbols:
α natural convection coefficient [-]
δ length of optical path for radiation [-]
Δ x , Δ t space and time discretizations [-]
ε inverse of the activation energy [-]
κ inverse of conductivity coefficient [-]
ρ density of the fuel [kg/m3]
σ Stefan-Boltzmann constant
φ Δ t ( u , v ) solution operator of (12) [-]
ψ Δ t ( u , v ) solution operator defined in (11) [-]
Ω × ( 0 , ) domain [-]
Abbreviations:
CFLCourant–Friedrichs–Lewy
CPUcentral processing unit
DIRKdiagonally implicit Runge–Kutta
ERKexplicit Runge–Kutta
H-LDIRK3(2,2,2)acronym of particular IMEX-RK scheme
IMEXimplicit–explicit
IMEX-RKimplicit–explicit Runge–Kutta
LI-IMEXlinearly implicit–explicit
NI-IMEXnonlinearly implicit–explicit
ODEordinary differential equation
PDEpartial differential equation
RKRunge–Kutta
S-LIMEXStrang-type splitting scheme
SSPstrong stability-preserving
WENOweighted essentially non-oscillatory

References

  1. Bürger, R.; Gavilán, E.; Inzunza, D.; Mulet, P.; Villada, L.M. Implicit–explicit methods for a convection-diffusion-reaction model of the propagation of forest fires. Mathematics 2020, 8, 1034. [Google Scholar] [CrossRef]
  2. Asensio, M.I.; Ferragut, L. On a wildland fire model with radiation. Int. J. Numer. Methods Eng. 2002, 54, 137–157. [Google Scholar] [CrossRef]
  3. Pastor, E.; Zárate, L.; Planas, E.; Arnaldos, J. Mathematical models and calculation systems for the study of wildland fire behaviour. Prog. Energy Combust. Sci. 2003, 29, 139–153. [Google Scholar] [CrossRef]
  4. Eberle, S.; Freeden, W.; Matthes, U. Forest fire spreading. In Handbook of Geomathematics, 2nd ed.; Freeden, W., Nashed, M.Z., Sonar, T., Eds.; Springer: Berlin, Germany, 2015; pp. 1349–1385. [Google Scholar]
  5. Ferragut, L.; Asensio, M.I.; Monedero, S. Modelling radiation and moisture content in fire spread. Commun. Numer. Methods Eng. 2007, 23, 819–833. [Google Scholar] [CrossRef]
  6. Ferragut, L.; Asensio, M.I.; Monedero, S. A numerical method for solving convection-reaction-diffusion multivalued equations in fire spread modelling. Adv. Eng. Softw. 2007, 38, 366–371. [Google Scholar] [CrossRef]
  7. Asensio, M.I.; Ferragut, L.; Simon, J. A convection model for fire spread simulation. Appl. Math. Lett. 2005, 18, 673–677. [Google Scholar] [CrossRef]
  8. Prieto, D.; Asensio, M.I.; Ferragut, L.; Cascón, J.M.; Morillo, A. A GIS-based fire spread simulator integrating a simplified physical wildland fire model and a wind field model. Int. J. Geogr. Inf. Sci. 2017, 31, 2142–2163. [Google Scholar] [CrossRef]
  9. Ferragut, L.; Asensio, M.I.; Cascón, J.M.; Prieto, D. A simplified wildland fire model applied to a real case. In Advances in Differential Equations and Applications; Casas, F., Martínez, V., Eds.; SEMA Springer Series 4; Springer International Publishing: Cham, Switzerland, 2014; pp. 155–167. [Google Scholar]
  10. San Martín, D.; Torres, C.E. Ngen-Kütral: Toward an open source framework for Chilean wildfire spreading. In Proceedings of the 37th International Conference of the Chilean Computer Science Society (SCCC), Santiago, Chile, 5–9 November 2018; pp. 1–8. [Google Scholar]
  11. San Martín, D.; Torres, C.E. Exploring a spectral numerical algorithm for solving a wildfire mathematical model. In Proceedings of the 38th International Conference of the Chilean Computer Science Society (SCCC), Concepción, Chile, 4–9 November 2019; pp. 1–7. [Google Scholar]
  12. Barovik, D.; Taranchuk, V. Mathematical modelling of running crown forest fires. Math. Model. Anal. 2010, 15, 161–174. [Google Scholar] [CrossRef] [Green Version]
  13. Taranchuk, V.; Barovik, D. Numerical modeling of surface forest fire spread in nonuniform woodland. In Computer Algebra Systems in Teaching and Research Vol. VIII; Prokopenya, A.N., Gil-Świderska, A., Siłuszyk, M., Eds.; Siedlce University of Natural Sciences and Humanities, Institute of Mathematics and Physics: Siedlce, Poland, 2019; pp. 159–168. [Google Scholar]
  14. Bebernes, J.; Eberly, D. Mathematical Problems from Combustion Theory; Springer: New York, NY, USA, 1989. [Google Scholar]
  15. Glassman, I.; Yetter, R.A.; Glumac, N.G. Combustion, 5th ed.; Academic Press: Boston, MA, USA, 2015. [Google Scholar]
  16. Séro-Guillaume, O.; Margerit, J. Modelling forest fires. Part I: A complete set of equations derived by extended irreversible thermodynamics. Int. J. Heat Mass Transf. 2002, 45, 1705–1722. [Google Scholar] [CrossRef]
  17. Margerit, J.; Séro-Guillaume, O. Modelling forest fires. Part II: Reduction to two-dimensional models and simulation of propagation. Int. J. Heat Mass Transf. 2002, 45, 1723–1737. [Google Scholar] [CrossRef]
  18. Linn, R.R.; Canfield, J.M.; Cunningham, P.; Edminster, C.; Dupuy, J.-L.; Pimont, F. Using periodic line fires to gain a new perspective on multi-dimensional aspects of forward fire spread. Agric. For. Meteorol. 2012, 157, 60–76. [Google Scholar] [CrossRef]
  19. Moinuddin, K.A.M.; Sutherland, D. Modeling of tree fires and fires transitioning from the forest floor to the canopy with a physics-based model. Math. Comput. Simul. 2020, 175, 81–95. [Google Scholar] [CrossRef]
  20. Hundsdorfer, W.; Verwer, J.G. Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations; Springer: Berlin, Germany, 2003. [Google Scholar]
  21. Butcher, J.C. Numerical Methods for Ordinary Differential Equations, 2nd ed.; Wiley: Chichester, UK, 2008. [Google Scholar]
  22. Hairer, E.; Nørsett, S.P.; Wanner, G. Solving Ordinary Differential Equations. I. Non-Stiff Problems; Springer Series in Computational Mathematics; Springer: Berlin, Germany, 1993; Volume 8. [Google Scholar]
  23. Ascher, U.; Ruuth, S.; Spiteri, J. Implicit-explicit Runge-Kutta methods for time dependent partial differential equations. Appl. Numer. Math. 1997, 25, 151–167. [Google Scholar] [CrossRef]
  24. Pareschi, L.; Russo, G. Implicit-Explicit Runge-Kutta schemes and applications to hyperbolic systems with relaxation. J. Sci. Comput. 2005, 25, 129–155. [Google Scholar]
  25. Boscarino, S.; Bürger, R.; Mulet, P.; Russo, G.; Villada, L.M. Linearly implicit IMEX Runge-Kutta methods for a class of degenerate convection-diffusion problems. SIAM J. Sci. Comput. 2015, 37, B305–B331. [Google Scholar] [CrossRef]
  26. Boscarino, S.; Bürger, R.; Mulet, P.; Russo, G.; Villada, L.M. On linearly implicit IMEX Runge-Kutta Methods for degenerate convection-diffusion problems modelling polydisperse sedimentation. Bull. Braz. Math. Soc. (New Ser.) 2016, 47, 171–185. [Google Scholar] [CrossRef]
  27. Boscarino, S.; Filbet, F.; Russo, G. High order semi-implicit schemes for time dependent partial differential equations. J. Sci. Comput. 2016, 68, 975–1001. [Google Scholar] [CrossRef] [Green Version]
  28. Boscarino, S.; LeFloch, P.G.; Russo, G. High order asymptotic-preserving methods for fully nonlinear relaxation problems. SIAM J. Sci. Comput. 2014, 36, A377–A395. [Google Scholar] [CrossRef]
  29. Donat, R.; Guerrero, F.; Mulet, P. Implicit-explicit methods for models for vertical equilibrium multiphase flow. Comput. Math. Appl. 2014, 68, 363–383. [Google Scholar] [CrossRef]
  30. Bürger, R.; Chowell, G.; Gavilán, E.; Mulet, P.; Villada, L.M. Numerical solution of a spatio-temporal gender-structured model for hantavirus infection in rodents. Math. Biosci. Eng. 2018, 15, 95–123. [Google Scholar] [CrossRef] [Green Version]
  31. Bürger, R.; Chowell, G.; Gavilán, E.; Mulet, P.; Villada, L.M. Numerical solution of a spatio-temporal predator-prey model with infected prey. Math. Biosci. Eng. 2019, 16, 438–473. [Google Scholar] [CrossRef]
  32. Bürger, R.; Goatin, P.; Inzunza, D.; Villada, L.M. A non-local pedestrian flow model accounting for anisotropic interactions and domain boundaries. Math. Biosci. Eng. 2020, 17, 5883–5906. [Google Scholar]
  33. Bürger, R.; Inzunza, D.; Mulet, P.; Villada, L.M. Implicit-explicit schemes for nonlinear nonlocal equations with a gradient flow structure in one space dimension. Numer. Methods Partial Differ. Equ. 2019, 35, 1008–1034. [Google Scholar] [CrossRef]
  34. Bürger, R.; Inzunza, D.; Mulet, P.; Villada, L.M. Implicit-explicit methods for a class of nonlinear nonlocal gradient flow equations modelling collective behaviour. Appl. Numer. Math. 2019, 144, 234–252. [Google Scholar] [CrossRef]
  35. Bürger, R.; Mulet, P.; Villada, L.M. Regularized nonlinear solvers for IMEX methods applied to diffusively corrected multi-species kinematic flow models. SIAM J. Sci. Comput. 2013, 35, B751–B777. [Google Scholar] [CrossRef]
  36. Abgrall, R.; Shu, C.-W. Handbook of Numerical Analysis for Hyperbolic Problems. Basic and Fundamental Issues; North-Holland: Amsterdam, The Netherlands, 2016. [Google Scholar]
  37. Shu, C.-W. Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws. In Advanced Numerical Approximation of Nonlinear Hyperbolic Equations; Quarteroni, A., Ed.; Lecture Notes in Mathematics; Springer: Berlin, Germany, 1998; Volume 1697, pp. 325–432. [Google Scholar]
  38. Shu, C.-W. High order weighted essentially nonoscillatory schemes for convection dominated problems. SIAM Rev. 2009, 51, 82–126. [Google Scholar] [CrossRef] [Green Version]
  39. LeVeque, R.J. Finite Volume Methods for Hyperbolic Problems; Cambridge University Press: Cambridge, UK, 2002. [Google Scholar]
  40. Hesthaven, J.S. Numerical Methods for Conservation Laws; Society for Industrial and Applied Mathematics: Philadelphia, PA, USA, 2018. [Google Scholar]
  41. Marin, M.; Vlase, S.; Paun, M. Considerations on double porosity structure for micropolar bodies. AIP Adv. 2015, 5, 037113. [Google Scholar] [CrossRef] [Green Version]
  42. Holden, H.; Karlsen, K.H.; Lie, K.-A.; Risebro, N.H. Splitting Methods for Partial Differential Equations with Rough Solutions; European Mathematical Society: Zurich, Switzerland, 2010. [Google Scholar]
  43. Strang, G. On the construction and comparison of difference schemes. SIAM J. Numer. Anal. 1968, 5, 506–517. [Google Scholar] [CrossRef]
  44. Bennett, M.; Fitzgerald, S.A.; Parker, B.; Main, M.L.; Perleberg, A.; Schnepf, C.; Mahoney, R.L. Reduce the Risk of Fire on Your Forest Property; Oregon State University: Corvallis, OR, USA, 2010. [Google Scholar]
Figure 1. Scenarios 1.1 to 1.4: diffusion coefficients K ( u ) .
Figure 1. Scenarios 1.1 to 1.4: diffusion coefficients K ( u ) .
Mathematics 08 01674 g001
Figure 2. Scenarios 1.1 to 1.4: three-dimensional plots of simulated temperature and two-dimensional plots of simulated burnt fuel, at simulated time T = 12 . The red and blue parts correspond to areas with v 0 = 0.9 and v 0 = 0 , respectively. The small white square indicates the location of the initial fire focus. The percentages represent the portion of burnt fuel of the total fuel initially available in Ω .
Figure 2. Scenarios 1.1 to 1.4: three-dimensional plots of simulated temperature and two-dimensional plots of simulated burnt fuel, at simulated time T = 12 . The red and blue parts correspond to areas with v 0 = 0.9 and v 0 = 0 , respectively. The small white square indicates the location of the initial fire focus. The percentages represent the portion of burnt fuel of the total fuel initially available in Ω .
Mathematics 08 01674 g002
Figure 3. Topographies for T ( x , y ) given by (15) with h 1 = h 2 = 60 and (left) x 1 = 0 , y 1 = 150 , γ 1 = 10 , 000 , x 2 = 150 , y 2 = 150 , γ 2 = 300 (topography T 1 for Scenario 2.2), (middle) x 1 = 150 , y 1 = 50 , γ 1 = 1000 , x 2 = 50 , y 2 = 150 , γ 2 = 1000 (topography T 2 for Scenario 2.3), (right) x 1 = 150 , y 1 = 75 , γ 1 = 3000 , x 2 = 50 , y 2 = 150 , γ 2 = 2000 (topography T 3 for Scenario 2.4), shown in each case as three-dimensional plots (top) and as contour maps (bottom).
Figure 3. Topographies for T ( x , y ) given by (15) with h 1 = h 2 = 60 and (left) x 1 = 0 , y 1 = 150 , γ 1 = 10 , 000 , x 2 = 150 , y 2 = 150 , γ 2 = 300 (topography T 1 for Scenario 2.2), (middle) x 1 = 150 , y 1 = 50 , γ 1 = 1000 , x 2 = 50 , y 2 = 150 , γ 2 = 1000 (topography T 2 for Scenario 2.3), (right) x 1 = 150 , y 1 = 75 , γ 1 = 3000 , x 2 = 50 , y 2 = 150 , γ 2 = 2000 (topography T 3 for Scenario 2.4), shown in each case as three-dimensional plots (top) and as contour maps (bottom).
Mathematics 08 01674 g003
Figure 4. Scenario 2.0: accuracy test. Numerical approximations at simulated time T = 0.2 and T = 0.4 obtained with the S-LIMEX scheme with Δ x = 80 / M and M = 160 and M = 640 (which correspond to the reference solution for computation of approximate errors in Table 1).
Figure 4. Scenario 2.0: accuracy test. Numerical approximations at simulated time T = 0.2 and T = 0.4 obtained with the S-LIMEX scheme with Δ x = 80 / M and M = 160 and M = 640 (which correspond to the reference solution for computation of approximate errors in Table 1).
Mathematics 08 01674 g004
Figure 5. Scenarios 2.1 and 2.2: simulation of the propagation of a forest fire on a flat domain (first and second row) and on a domain with topography T 1 (third and fourth row), at indicated simulated times. 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.
Figure 5. Scenarios 2.1 and 2.2: simulation of the propagation of a forest fire on a flat domain (first and second row) and on a domain with topography T 1 (third and fourth row), at indicated simulated times. 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.
Mathematics 08 01674 g005
Figure 6. Scenarios 2.3 and 2.4: simulation of the propagation of a forest fire on domains with topographies T 2 (first and second row) and T 3 (third and fourth row).
Figure 6. Scenarios 2.3 and 2.4: simulation of the propagation of a forest fire on domains with topographies T 2 (first and second row) and T 3 (third and fourth row).
Mathematics 08 01674 g006
Figure 7. Scenario 3.1: principle of construction of a risk map for a flat domain. Here and in Scenarios 3.2 to 3.4, the propagation of a wildfire is simulated starting from an array of 5 × 5 = 25 initial fire foci, marked by small white squares. For each initial position the simulation is stopped when 5 % of the initially available fuel is burnt. The corresponding times T risk are mapped in Figure 8. (For the present flat case, the shapes of the temperature and fuel distributions are the always the same in relation to the initial fire; instead of displaying 2 × 25 identical results, we only show the ‘corner’ cases of the 5 × 5 initial positions).
Figure 7. Scenario 3.1: principle of construction of a risk map for a flat domain. Here and in Scenarios 3.2 to 3.4, the propagation of a wildfire is simulated starting from an array of 5 × 5 = 25 initial fire foci, marked by small white squares. For each initial position the simulation is stopped when 5 % of the initially available fuel is burnt. The corresponding times T risk are mapped in Figure 8. (For the present flat case, the shapes of the temperature and fuel distributions are the always the same in relation to the initial fire; instead of displaying 2 × 25 identical results, we only show the ‘corner’ cases of the 5 × 5 initial positions).
Mathematics 08 01674 g007
Figure 8. Scenarios 3.1 to 3.4: risk maps for the subdomain [ 60 , 110 ] × [ 60 , 110 ] under the assumption of a flat domain (Scenario 3.1, top left) and domains with topography T 1 , T 2 , and T 3 (Scenarios 3.2 to 3.4, top right, bottom left, and bottom right), visualizing for each patch of size 10 × 10 the final time T risk . A low value of T risk implies a high risk.
Figure 8. Scenarios 3.1 to 3.4: risk maps for the subdomain [ 60 , 110 ] × [ 60 , 110 ] under the assumption of a flat domain (Scenario 3.1, top left) and domains with topography T 1 , T 2 , and T 3 (Scenarios 3.2 to 3.4, top right, bottom left, and bottom right), visualizing for each patch of size 10 × 10 the final time T risk . A low value of T risk implies a high risk.
Mathematics 08 01674 g008
Figure 9. Scenario 3.2: simulation of temperature of a wildfire from 5 × 5 = 25 initial fire foci, marked by small white squares, on a domain with topography T 1 . Here and in Figure 10, Figure 11, Figure 12, Figure 13 and Figure 14, for each initial position where the simulation is stopped, and the numerical solutions for u and v are portrayed, when 5 % of the initially available fuel is burnt. The corresponding simulated time determines T risk for each initial position.
Figure 9. Scenario 3.2: simulation of temperature of a wildfire from 5 × 5 = 25 initial fire foci, marked by small white squares, on a domain with topography T 1 . Here and in Figure 10, Figure 11, Figure 12, Figure 13 and Figure 14, for each initial position where the simulation is stopped, and the numerical solutions for u and v are portrayed, when 5 % of the initially available fuel is burnt. The corresponding simulated time determines T risk for each initial position.
Mathematics 08 01674 g009
Figure 10. Scenario 3.2: simulation of fuel of a wildfire from 5 × 5 = 25 initial fire foci, marked by small white squares, on a domain with topography T 1 .
Figure 10. Scenario 3.2: simulation of fuel of a wildfire from 5 × 5 = 25 initial fire foci, marked by small white squares, on a domain with topography T 1 .
Mathematics 08 01674 g010
Figure 11. Scenario 3.3: simulation of temperature of a wildfire from 5 × 5 = 25 initial fire foci, marked by small white squares, on a domain with topography T 2 .
Figure 11. Scenario 3.3: simulation of temperature of a wildfire from 5 × 5 = 25 initial fire foci, marked by small white squares, on a domain with topography T 2 .
Mathematics 08 01674 g011
Figure 12. Scenario 3.3: simulation of fuel of a wildfire from 5 × 5 = 25 initial fire foci, marked by small white squares, on a domain with topography T 2 .
Figure 12. Scenario 3.3: simulation of fuel of a wildfire from 5 × 5 = 25 initial fire foci, marked by small white squares, on a domain with topography T 2 .
Mathematics 08 01674 g012
Figure 13. Scenario 3.4: simulation of temperature of a wildfire from 5 × 5 = 25 initial fire foci, marked by small white squares, on a domain with topography T 3 .
Figure 13. Scenario 3.4: simulation of temperature of a wildfire from 5 × 5 = 25 initial fire foci, marked by small white squares, on a domain with topography T 3 .
Mathematics 08 01674 g013
Figure 14. Scenario 3.4: simulation of fuel of a wildfire from 5 × 5 = 25 initial fire foci, marked by small white squares, on a domain with topography T 3 .
Figure 14. Scenario 3.4: simulation of fuel of a wildfire from 5 × 5 = 25 initial fire foci, marked by small white squares, on a domain with topography T 3 .
Mathematics 08 01674 g014
Table 1. Scenario 2.0: Approximate L 1 errors e M ( u ) and e M ( v ) for Strang LIMEX method with Δ x = 80 / M . Reference solution computed with M = 640 .
Table 1. Scenario 2.0: Approximate L 1 errors e M ( u ) and e M ( v ) for Strang LIMEX method with Δ x = 80 / M . Reference solution computed with M = 640 .
T = 0.2 T = 0.3 T = 0.4
M e M ( u ) e M ( v ) e M ( u ) e M ( v ) e M ( u ) e M ( v )
402.23.7 ×   10 2 3.24.5 ×   10 2 4.96.4 ×   10 2
801.21.6 ×   10 2 2.02.3 ×   10 2 3.33.7 ×   10 2
1600.71.0 ×   10 2 1.21.5 ×   10 2 2.02.4 ×   10 2
3200.45.7 ×   10 3 0.68.9 ×   10 3 1.11.4 ×   10 2

Share and Cite

MDPI and ACS Style

Bürger, R.; Gavilán, E.; Inzunza, D.; Mulet, P.; Villada, L.M. Exploring a Convection–Diffusion–Reaction Model of the Propagation of Forest Fires: Computation of Risk Maps for Heterogeneous Environments. Mathematics 2020, 8, 1674. https://doi.org/10.3390/math8101674

AMA Style

Bürger R, Gavilán E, Inzunza D, Mulet P, Villada LM. Exploring a Convection–Diffusion–Reaction Model of the Propagation of Forest Fires: Computation of Risk Maps for Heterogeneous Environments. Mathematics. 2020; 8(10):1674. https://doi.org/10.3390/math8101674

Chicago/Turabian Style

Bürger, Raimund, Elvis Gavilán, Daniel Inzunza, Pep Mulet, and Luis Miguel Villada. 2020. "Exploring a Convection–Diffusion–Reaction Model of the Propagation of Forest Fires: Computation of Risk Maps for Heterogeneous Environments" Mathematics 8, no. 10: 1674. https://doi.org/10.3390/math8101674

APA Style

Bürger, R., Gavilán, E., Inzunza, D., Mulet, P., & Villada, L. M. (2020). Exploring a Convection–Diffusion–Reaction Model of the Propagation of Forest Fires: Computation of Risk Maps for Heterogeneous Environments. Mathematics, 8(10), 1674. https://doi.org/10.3390/math8101674

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop