Implicit-Explicit Methods for a Convection-Diffusion-Reaction Model of the Propagation of Forest Fires

: Numerical techniques for approximate solution of a system of reaction-diffusion-convection partial differential equations modeling the evolution of temperature and fuel density in a wildﬁre are proposed. These schemes combine linearly implicit-explicit Runge–Kutta (IMEX-RK) methods and Strang-type splitting technique to adequately handle the non-linear parabolic term and the stiffness in the reactive part. Weighted essentially non-oscillatory (WENO) reconstructions are applied to the discretization of the nonlinear convection term. Examples are focused on the applicative problem of determining the width of a ﬁrebreak to prevent the propagation of forest ﬁres. Results illustrate that the model and numerical scheme provide an effective tool for deﬁning that width and the parameters for control strategies of wildland ﬁres.


Scope
It is the purpose of this contribution to provide an efficient numerical method for the solution of a wildland fire model [1]. This model can be expressed as a convection-diffusion-reaction system of the type where t is time, x ∈ Ω is the spatial variable where Ω ⊂ R 2 , and u = u(x, t) and v = v(x, t) are the sought scalar functions. Here u represents the non-dimensionalized temperature and v the non-dimensionalized mass fraction of solid fuel. Moreover, K is a given diffusion coefficient that depends only on u, and w is an advection velocity that represents wind speed. It is assumed that div w = 0. 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 evolution of a wildfire is described by the system (1) together with zero-flux boundary conditions uw − K(u)∇u · n = 0, (x, t) ∈ ∂Ω × (0, +∞), (2) where n is the unit normal vector to ∂Ω, and the initial conditions Explicit numerical schemes for (1) on a uniform Cartesian grid of meshwidth ∆x and time step ∆t are easy to implement but are associated with a Courant-Friedrichs-Lewy (CFL) stability condition [2] that imposes the proportionality ∆t ≈ ∆x 2 [3][4][5][6][7][8][9][10][11][12][13]. This restriction makes long-term simulations of (1)-(3) on a uniform grid unacceptably slow. A well-known remedy consists in handling the diffusive term in (1) by an implicit discretization. The resulting semi-implicit or implicit-explicit (IMEX) schemes (e.g., [3,11]) are associated with an acceptable CFL condition ∆t ≈ ∆x. (Here we recall that semi-implicit (or semi-explicit) and IMEX schemes are not the same; the former class usually comprises one-step schemes such as the implicit Euler scheme or the trapezoidal rule while the latter larger class is based on intertwined computations of the stages of explicit and implicit Runge-Kutta schemes.) However, the case that K depends nonlinearly on u needs to be handled by special nonlinear solvers [11] or by solving linear problems in each time step. These linear problems arise from carefully distinguishing between stiff and nonstiff unknowns in the discretized version of ∇ · K(u)∇u [3]. These variants will be addressed as nonlinearly implicit (NI-IMEX) and linearly implicit (LI-IMEX) schemes, respectively. In this work we focus on the implementation of LI-IMEX schemes.
On the other hand, the reaction term for this model imposes a very severe restriction on the global time step ∆t. As a strategy to avoid this restriction, we propose a variant of the LI-IMEX scheme consisting in a Strang-type splitting technique. This technique consists in solving the convection-diffusion system for u and coupling this result with the solution of the ordinary differential equation (ODE) for u and v defined by the reaction terms (in which the convective and diffusive term is absent). This results in a numerical scheme that allows a less restrictive time step ∆t.
The main research question of this paper is whether LI-IMEX schemes, possibly equipped with Strang-type operator splitting, indeed provide a gain in CPU time, and hence in efficiency if applied to the wildland fire model (1) under realistic choices of model functions and parameters. Here a gain is understood in comparison with the performance of explicit schemes. The result is not obvious, since LI-IMEX schemes permit larger time steps but are associated with additional computational effort related to the solution of linear systems in each iteration. A secondary research question consists in assessing the suitability of model (1) for the simulation of wildfire control strategies.
From a practical point of view, we address the importance of being able to simulate the effects of artificial firebreaks as a measure of prevention and control of forest fires. The effectiveness of a firebreak, which is normally created by cutting a lane into a forest, is strongly determined by its width, as is discussed in [14,15]. This motivates the implementation of scenarios, such as those shown here, in order to help solve the problem of determining an adequate width and thus protect a certain area from forest fire.

Related Work
The literature includes many approaches to studying forest fires such as physical models, cellular automata and artificial intelligence, to name a few. A general survey on the principles of mathematical modeling of forest fires is provided by Pastor et al. [16] and more recently by Eberle et al. [17]. Most continuous-in-space approaches to model and study forest fires lead to partial differential equations (PDEs) as local forms of the balance equations that describe fuel consumption and heat transfer. A discussion of recent developments of PDE-based forest fire models, as opposed to those based on cellular automata, is included in [18]. We herein directly use the two-PDE model formulated by Asensio and Ferragut [1]. This model is widely accepted and has laid the foundation for various developments such as the inclusion of vegetation, moisture and endothermic pyrolysis by means of a multivalued function representing the enthalpy [19,20]. Reports on implementation of the model by Asensio and Ferragut [1] within software packages for real scenarios include [17,21,22]. For sake of completeness, we mention that PDE-based models have also been developed for a refined description of the combustion front of a wildfire; see e.g., [23] and the references cited in that work. Solutions to the model (1) may be obtained through numerical methods such as mixed finite elements, finite differences, or finite volumes. Due to the structure of Equation (1), there is ample theoretical support for a spatial semi-discretization of the model following a method-of-lines strategy. The resulting scheme turns out to be explicit due to the choice of ODE integrators. However, explicit schemes applied to reaction-diffusion problems strongly restrict the time step due to stability conditions. Therefore, our goal is to treat reaction and diffusion implicitly and convection explicitly. To avoid the need to solve nonlinear systems that appear in the implicit treatment of nonlinear diffusion, we propose LIMEX Runge-Kutta schemes introduced in [3,4].
To underline the relevance and novelty of implicit-explicit discretizations for the wildfire model (1), let us briefly review state-of-the-art methods used in some of the references. Asensio and Ferragut [1] employ mixed finite element methods for the spatial discretization combined with the Euler method in time, which leads to a semi-implicit scheme. The versions of this model with moisture and pyrolisis effects (involving enthalpy) that include heat either a non-local [19] or or a local (diffusive) model of heat radiation [20] are based on related semi-implicit [19] or implicit [20] time discretizations, based on explicit or implicit Euler time stepping. A broader discussion of numerical methods for wildfire models, some of them based on work not cited in this paper, is provided by Eberle et al. [17]. Their approach to discretize the model (4), (5) is based on a collocational method (using radial basis functions) for space discretization in combination with a Crank-Nicolson scheme for time stepping, including a stabilization of convective terms. This approach does not involve a weak formulation of the problem [17] (p. 1363). On the other hand, San Martin and Torres [22] devise a spectral differentiation method based on a method-of-lines formulation (discrete in space, continuous in time) for a version of the model (4), (5) with linear diffusion. This method is a formalism of expressing the spatial partial derivatives in terms of Fourier transforms of the unknowns (temperature and fuel). The time integration is done by the classical, explicit fourth-order Runge-Kutta method. Finally, we mention that the treatments by Barovik and Taranchuk [23] and Taranchuk and Barovik [24] employ finite difference methods with explicit Euler time stepping. All these approaches refer to different models but involve at least two unknowns, temperature and fuel. The discretization methods have in common, however, that the numerical solution is evolved over one single time step by discretizing each model ingredient either explicitly or implicitly. For all these approaches a more efficient alternative may arise through the intertwined explicit and implicit evaluations of stages of IMEX schemes, as is explained in Section 1.1. That said, improvements are most likely to accrue for base methods that are fully explicit, such as the one utilized in [23,24].

Outline of the Paper
The remainder of the paper is organized as follows. In Section 2 the governing mathematical model is introduced. To this end, we formulate in Section 2.2 a model system of PDEs which that describes in simplified form the spatio-temporal evolution of temperature and fuel density. Section 3 is devoted to the description of the numerical schemes proposed to solve the PDE model. Values of the model parameters and an estimate of the maximum temperature attainable are provided in Section 2.4. Section 3.1 introduces the notation for two numerical schemes and details the discretization of the convective flux by a fifth-order weighted essentially non-oscillatory (WENO) reconstruction.
In Section 3.2 the linearly implicit IMEX-RK method and in Section 3.3 the linearly implicit method IMEX-RK with Strang-type splitting are outlined. Section 4 is devoted to the presentation of numerical results. The examples correspond to four different scenarios. The first shows the numerical solution of the model considering an initial ignition in a domain with a constant fuel density. Results of this example, which also include a study of approximate L 1 errors, are shown in Section 4.1. In the second case we also impose an initial ignition but now with a terrain separated into two parts with different initial values for the fuel density. In the third we assume that the two regions that differ in the initial value for the fuel density are separated by a firebreak. In the fourth and last case we show the evolution of the variables in a terrain with random fuel density where the wind slightly changes its direction. The corresponding numerical examples are presented in Sections 4.2-4.4, respectively. Conclusions are collected in Section 5.

Model Development
The mathematical model follows the presentation by Asensio and Ferragut [1]. In a simplified setting, we consider the amount of solid fuel, but not its composition. The chemical reaction of fuel and oxidants into products is supposed to take place at the rate r defined by the Arrhenius equation where E A denotes the activation energy, U is absolute temperature, R = 1.987207 cal/(K mol) is the universal gas constant, and A is a constant. A typical value is E A ≈ 20 kcal/mol [26]. Furthermore, if Y denotes the mass fraction of fuel, then the rate of fuel consumption S f can be chosen as S f = −rρY, where ρ is the density. The balance equation for the oxidant need not be considered [1]. Presupposing that the fire source is large (as should be typical for a wildfire), we may assume that the dominant mode of heat transfer is radiation. Heat flux transfer by radiation through a fuel bed can be derived from the Stefan-Boltzmann law (see [1,17] for details), which leads to the equation q = −4σδU 3 ∇U for the radiative heat flux, where σ is the Stefan-Boltzmann constant and δ = δ is the length of the optical path for radiation through the substance. Furthermore, we take into account vertical heat loss due to gravity by a natural convection term within the energy conservation equation, namely h(x)(U − U ∞ ). Here h is the natural convection coefficient and U ∞ stands for the ambient temperature. While h = h(x) usually represents information on the local topography, we here consider h to be constant. On the other hand, we describe convection of heat due to wind by a term ρCw(x, t).
Here ρ and C are the density and specific heat. The vector field w = w(x, t) is the wind velocity vector. We assume that temperature and velocity correspond to averaged values in a turbulent regime but that the reactive term is governed by the laminar regime. Under the latter assumption and denoting by H the heat of combustion, we may model the energy released in an exothermic reaction by Finally, we need to describe the change between the endothermic (solid) and the exothermic (gaseous) phase. This is done by the function that describes phase change, which is supposed to occur at a temperature U pc . Utilizing this function allows us to handle both phases within the same system of equations. The basic processes that are described in this way are fuel combustion and heat generation when U exceeds the phase change temperature U pc .

Simplified Model in Final Form
Combining the previous ingredients, we obtain the following model: We study the model (4), (5) on a two-dimensional domain Ω ⊂ R 2 with a piecewise smooth boundary Γ = ∂Ω. It is assumed that Ω is large enough so that U and Y do not change on Γ during the simulation time (0, t max ). We may therefore impose either zero-flux or Dirichlet boundary conditions. Spatial heterogeneity can be represented by the functions h = h(x) and the initial distribution of It is assumed that the thermal conductivity coefficient k is a positive constant. Thus, (4) is a nonlinear but strictly parabolic PDE whose solutions are smooth. This means that solutions of (4), (5) do not develop discontinuities (shocks).

Dimensionless Variables
The combustion model (4), (5) can be non-dimensionalized in a rational manner in order to elucidate the significant parameters. We use the Frank-Kamenetskii change of variables (see [26]). We fix a reference value for the temperature and the fuel in which an equilibrium can be assumed, namely the ambient temperature U ref = U ∞ and the initial fuel. Since Y 0 may be variable, we set likewise, we define Furthermore we assume that t 0 and l 0 are characteristic magnitudes for temporal and spatial variables to be fixed later. The non-dimensional change is defined as follows, where we assume that x = (x 1 , x 2 ): Here ε is the non-dimensional inverse of the activation energy, that is ε = RT ∞ /E A . Furthermore, we define the non-dimensional reaction heat q react = HY ref /(CT ∞ ) and the time and length scales We define the non-dimensional natural convection coefficient α = α(ξ, η), the non-dimensional inverse of the conductivity coefficient κ and the non-dimensional phase change function c(u) by the respective expressions For simplicity, let us recover the notation t, x 1 and x 2 for non-dimensional temporal and spatial variables, Ω × (0, t max ) for the non-dimensional domain, and maintain u for the non-dimensional average temperature and v for the non-dimensional mass fraction of solid fuel. Then the model can be written in the form (1) if we define

Parameters Used for the Wildfire Model
We assume an ambient temperature of U ∞ = 303 K and an activation energy E A = 20 kcal/mol = 83.68 kJ/mol, which yields ε = 0.02980905 ≈ 0.03. Furthermore, we select A = 10 9 s −1 [26]. According to [1], realistic values of density ρ for different kinds of wood vary from ρ = 420 kg m −3 (for fir) to 640 kg/m 3 (for yellow pine). The specific heat varies from C = 2.4 kJ kg −1 K −1 (for maple or oak) to C = 2.8 kJ kg −1 K −1 (for yellow pine), and the thermal conductivity ranges from k = 0.059 W m −1 K −1 (for maple or oak) to k = 1.17 W m −1 K −1 for sawdust. On the other hand, the values for air at atmospheric pressure and ambient temperature are ρ = 1.1774 kg m −3 , C = 1.0057 kJ kg −1 K −1 and k = 0.024 W m −1 K −1 . The mean magnitudes for the numerical examples in [1], and which we adopt for our numerical experiments, are Asensio and Ferragut [1] argue furthermore that the heat of combustion for cellulose is H = 15,900 kcal kg −1 , but that the presence of fuels with lower heats of combustion justifies smaller values of H. They choose H in such a way that If we set Y 0 = 1, then this corresponds to H = 300 kJ kg −1 = 71.702 kcal kg −1 . Utilizing the parameters (7), we 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.
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 (11), this means On the other hand, the second equation in (11) 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.

Notation and Semi-Discrete Formulation
We take Ω = [0, L] × [0, L] and denote by u : Here it is understood that K, f and g are given by (6). We use a uniform Cartesian grid with nodes ( and for simplicity we denote K i = K(u i ) and for the reactive part we define Using this notation, we may approximate (1), (6) in semi-discrete form (that is, in discrete in space but continuous in time form) by the system of ODEs Here C(u) and D(u)u represent the spatial discretizations of the convective and diffusive terms with entries given by C(u) = (C(u) i ) i∈M , D(u)u = ((D(u)u) i ) i∈M given by is the fifth-order WENO spatial discretization of the convective term w · ∇u [11]. The fifth-order WENO finite difference discretization of convective terms on Cartesian grids 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 [2,27,28] for these and many other alternative schemes) stems from our familiarity with it. Although there are no shock waves due to presence of diffusion (see the comment at the end of Section 2.2), there may be sharp gradients that are dealt with in a robust manner by the WENO The approximate solution of (8) can be obtained by application of Runge-Kutta ODE solvers. Strong Stability Preserving (SSP) explicit Runge-Kutta methods are a popular class of time integrators whose use leads to a stronger stability constraint on ∆t, see [3][4][5][6][9][10][11]13]. An alternative to explicit RK schemes are implicit-explicit Runge-Kutta (IMEX-RK) methods (see [9][10][11]), for which only the diffusion term is treated implicitly. In this case the stability condition on ∆t is less severe than for explicit RK schemes, but a large number of nonlinear systems need to be solved. To overcome the excessive numerical work for the solution of nonlinear systems, an essential gain is obtained by the approach proposed in [3,4] that is based on linearly implicit-explicit Runge-Kutta schemes. The approach for (1) is based on distinguishing in (8) between stiff and non-stiff dependence of various terms on the solution vectors u and v. This distinction allows us to choose 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 occurrence of the solution 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.

Linearly Implicit IMEX-RK Scheme
Next, we define the following functions, which are given by the right-hand sides of (8) where we have replaced the argument u by u * in all terms for which the dependence is non-stiff. These terms are therefore discretized in an explicit form: Then we rewrite (8) in the form Observe that the only stiff terms, which are treated implicitly, are the linear terms that multiply D(u) and A in K u (u * , u, v) and ζ(u) in K v (u * , v). Thus, we treat u * explicitly as argument of C, D and h. The pair of Butcher arrays [12] of IMEX-RK methods for the time integration of (9) is given by (see [3,5]). Both the explicit and diagonally implicit parts define third-order schemes [5]. The particular value of γ guarantees that the implicit part is a third-order DIRK scheme with the best dampening properties [5].
We now aim at applying a partitioned Runge-Kutta scheme, consisting in the application of the explicit part to the first block of equations and the implicit part to the second and third blocks. If both Butcher arrays satisfyb = b (as in the case of H-LDIRK3(2,2,2)), then the step from t n to t n+1 = t n + ∆t of the linearly implicit IMEX-RK scheme is given by the following Algorithm 1.
Observe that K v i in (10) can be computed by evaluating the right-hand in direct form. In contrast, the linear equation for K u i is solved by using standard and efficient block tridiagonal solvers.

Algorithm 1 LI-IMEX Runge-Kutta scheme
Input: approximate solution vector u n and v n for t = t n do i = 1, . . . , s compute the stage values: solve for K u i the linear system Output: approximate solution vectors u n+1 and v n+1 for t = t n+1 = t n + ∆t.

Linearly Implicit IMEX-RK Method with Strang-Type Splitting
An alternative numerical scheme to approximate solutions of (1), (6) consists in applying a splitting technique. To this end we first consider the ODE system that arises when heat transport and diffusion terms are switched off, namely We discretize (11) utilizing a time-implicit discretization for v and an explicit one for u, i.e., In order to compute (12) 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 define ψ ∆t := R M 2 × R M 2 → R M 2 × R M 2 as the solution obtained in (12) in the way Next, we solve the system of ordinary differential equations then we define ϕ ∆t : where u n+1 is the approximation of (13) by using Algorithm 1. Finally, the Strang splitting method (cf. [2,29]) to solve (8) is formulated as

Numerical Results
In the following examples, we solve (8) numerically for 0 ≤ t ≤ T and (x, y) ∈ Ω := (0, L) 2 , under various choices of the initial datum u 0 and v 0 . Numerical results are obtained by the linearly implicit IMEX scheme denoted by LIMEX (see Algorithm 1), and the Strang-type splitting scheme (14) called S-LIMEX. We consider a Cartesian uniform mesh grid of M 2 points and ∆x = ∆y = L/M. We denote by ∆t the time step used to advance the numerical solution from t = t n to t n+1 = t n + ∆t. For each iteration ∆t is determined by the following usual formula (CFL condition) In all numerical tests we have used C cfl = 0.12.

Remark 1.
A fully explicit discretization of reactive term in (11) by using forward Euler scheme is associated with the CFL time step restriction ∆t ∆x where ρ(·) is the spectral radius and J ( f , g) is the Jacobian matrix of functions given by (6). In our cases, |ρ(J ( f , g))| ≤ |ζ (u)| ≈ 10 8 for the chosen parameters, which imposes a severe CFL restriction for ∆t.

Remark 2.
An explicit discretization of the diffusive term in (8) leads to the following CFL condition, which imposes a severe restriction on the time step ∆t, see [3,10,11].

Example 1: Propagation of a Wildfire with Constant Initial Distribution of Fuel
In this example we compare the approximations obtained by both schemes, S-LIMEX and LIMEX, for different discretizations and simulated times. We consider first a configuration that allows us to compare our results with those of [1]. This means that we solve (1) on the domain Ω with L = 50 in dimensionless variables, which corresponds to a square of side length 50l 0 = 14.99 m. We choose the dimensionless wind vector w = (w 1 , w 2 ) = (300, 300), which blows in the north-east direction at physical speed v = (l 0 /t 0 ) w ≈ 0.0142 m s −1 (much larger values seem more realistic). The initial conditions for the non-dimensional temperature and the fuel are u 0 (x, y) = 31χ [4,12] 2 (x, y), v 0 (x, y) = 0.6χ Ω (x, y).
We display the numerical approximation obtained with both schemes in Figures 1 and 2, with ∆x = ∆y = 50/400 at simulated times T = 0.03 and T = 0.06, which represent 4.5 and 9 minutes approximately. The dynamics of the propagation of wildfire is comparable with the results of [1]. In Figure 3(left) we display the evolution of the non-dimensional temperature within a short time interval on the diagonal y = x of the domain Ω. The numerical approximation was computed with S-LIMEX scheme with ∆x = ∆y = 50/800, corresponding to time step ∆t = 7.8125 × 10 −6 . The highest temperature reached during the initial steps does not exceed the maximum U max for the given parameters and initial conditions chosen. We also observe the formation of a fire front traveling in a direction parallel to w. In Figure 3(right) we compare the numerical approximation for both numerical schemes with respect to the reference solution for two discretization levels with ∆x = ∆y = 50/M for M = 100 and M = 200, corresponding to ∆t = 6.25 × 10 −5 and ∆t = 3.125 × 10 −5 , respectively. Here we have represented the non-dimensional temperature at simulated time T = 0.06 on the diagonal y = x. We observe that the numerical approximations are qualitatively equivalent. In Table 1 we compute the approximate L 1 -errors e M (u) and e M (v) for both numerical schemes with different discretizations ∆x = ∆y = 50/M and M = 50, 100, 200, 400 and 800, with respect to a reference solution computed with S-LIMEX scheme with M = 3200. For these cases, we employ ∆t = 1.25 × 10 −4 for the coarse mesh and ∆t = 7.8125 × 10 −6 for the finer mesh. We observe that for the coarse mesh, the S-LIMEX scheme is more accurate than the LIMEX scheme, however for a finer mesh, LIMEX is more accurate and faster than the S-LIMEX scheme.  On the other hand, to demonstrate that the proposed S-LIMEX scheme really provides an advantage, we consider an alternative scenario where L = 10, w = (1, 1), and the initial conditions are with the same parameters as in the above scenario. For this configuration we compare two schemes: the S-LIMEX scheme and a scheme where all terms are discretized explicitly (Fully Explicit scheme). Since especially the last scheme entails very large computation times, we present results for fairly short simulation times only, namely T = 0.02 and T = 0.03. The scenario is of interest only in terms of accuracy and CPU times. The corresponding results, namely the approximate L 1 -errors e M (u) and e M (v) and CPU times for these two methods are provided in Table 2. Here we compare the error of approximation with respect to a reference solution with M ref = 1280 cells per dimension. According to Table 2, the errors and CPU times produced by the S-LIMEX scheme are smaller than the Fully Explicit scheme. This means that the S-LIMEX scheme is more efficient. Theoretical considerations support that the S-LIMEX scheme should become more efficient for sufficiently fine discretizations, but the results of Table 2 provide evidence that there is a clear gain in efficiency even for moderately fine discretizations (e.g., M = 40 or M = 80).

Example 2: Propagation of a Wildfire with Distribution of Fuel Determined by Two Areas with Different Fuel Densities
In this example we describe the dynamics of the wildfire when it enters a region with different fuel density. We consider Ω := (0, 200) × (0, 200), ε = 0.035, u pc = 20 and w = (300, 300). First, we employ the initial conditions u 0 (x, y) = 21χ [9,21] where δ is a function that is constant on each cell of the computational grid, and takes random values with |δ(x, y)| ≤ 0.025. Notice that the initial density of fuel is higher for x > 100 than for x < 100. In Figure 4 we display the numerical solution at four simulated times, namely the initial condition at T = 0 and the solution at T = 0.1, T = 0.2 and T = 0.3. We observe that when the wildfire enters a region with higher fuel density, the non-dimensional temperature increases. Moreover, the width of the fire in the main direction of propagation is higher than in the less dense area, and this width increases with time and the fire opens.
Next we employ the same initial condition for temperature. However, to detect where the fire is extinguished when the availability of fuel is low, we impose v 0 (x, y) = 0.45 for x ≤ 50, 5 × 10 −4 for x > 50.
In Figure 5 we display the initial condition and the numerical solution at T = 0.1, T = 0.2 and T = 0.3. We observe that when the wildfire enters a region with a very low fuel density, the fire focus width decreases with time and tends to disappear since the temperature of the focus remains below u pc when it has already traveled more than 100 (unit lengths). This behavior allows one to identify the horizontal length after which the fire is extinguished (at the given wind speed). The numerical results in

Example 3: Effect of a Firebreak in the Propagation of a Wildfire
In the prevention and control of forest fires, the role played by the implementation or construction of firebreaks (see Figure 6) is very important. Their role is to mitigate the effect of a fire or ideally prevent the spread of fire. This is usually done through the creation of a rectangular space void of fuel, that is in this case, of vegetation whose objective is to separate two zones to protect a certain area. Clearly, it is of great importance that one should be able to estimate the necessary width of the firebreak, as is pointed out in [14,15]. The protected area could be a forest with a large amount of fuel from which wood is obtained for commercial use, a native forest, or a set of households. To implement this fire control strategy using firebreaks, it is essential to study the effects that determine the width of said area with a fuel density close or equal to zero. This situation motivates this third example where we study the effect of a firebreak in the propagation of a wildfire separating two areas with different fuel distribution. We consider the same parameters as in Example 2 with u 0 as in (15) and the initial distribution of fuel v 0 (x, y) = 0.45 + δ(x, y) for x ≤ 100 − l, where again the function δ is constant on each cell of the computational grid, and takes random values with |δ(x, y)| ≤ 0.025. Moreover, l > 0 is the width of the firebreak separating the two areas. Here we consider two cases, l = 25 and l = 50.
In the case l = 25, we display in Figure 7 the dynamics of propagation of the wildfire at four simulated times. We observe that when the wildfire front crosses the firebreak region, the non-dimensional temperature decreases, however the wildfire does not extinguish and propagates into the area with high fuel density. Its temperature increases and the fire focus widens. In Table 3 we compute the total amount of fuel in the denser area for different widths of the firebreak. We observe that amount of fuel burnt for l = 25 is smaller than for l = 0, see also Figure 4 corresponding to the first part of Example 2. On the other hand, when l = 50, we can see in Figure 8 that the wildfire is extinguished and does not propagate from the area with lower fuel density to that with higher density. These values of l are in accordance with those suggested by CONAF in its handbook [15].
The results obtained in these two scenarios agree with those of the second part of the previous example since there it can be seen that the temperature is below u pc = 20 and remains so over time for x < 100, which suggests that the width of the firebreak should be approximately 50 units of length. The numerical results of Figure 9 were

Example 4: Propagation of a Wildfire with Randomized Initial Distribution of Fuel and a Slight Change of Direction in the Wind
On a macroscopic scale, the wind direction is generally variable, but on a smaller scale and for a fixed region in the context of forest fires this direction can be considered constant or at most slightly variable only. In this configuration we consider Ω with L = 200 in dimensionless variables, which corresponds to a square of side length 200l 0 = 59.96 m with a randomized initial distribution of fuel in the whole domain. We choose initially the dimensionless initial wind vector w = (300, 300) and we change the wind direction at time T = 0.085, which represents approximately 14 minutes, to w * = (50, 400). We wish to observe the dynamics in the wildfire until simulated time T = 0.15. The initial values for the non-dimensional temperature and the fuel distribution are u 0 (x, y) = 38χ [2,14] where δ is a function that is constant on each cell of the computational grid, and takes random values with |δ(x, y)| ≤ 0.4. In Figure 9 we can observe the effect of the change of the wind. From the initial time until before the simulated time T = 0.085 it can be seen that the strip corresponding to the burnt fuel has a certain constant thickness and when the direction of the wind changes to w * we observe an increase in the width of the strip of the burnt area. The numerical results in Figure 9 were obtained by the S-LIMEX scheme with ∆x = ∆y = 100/M and M = 400, corresponding to ∆t = 6.25 × 10 −5 .

Conclusions
We have implemented two numerical methods that allow studying the model defined by (4) and (5). The first method is a linearly implicit IMEX-RK scheme (denoted by "LIMEX") and is described in Section 3.2. The second is a linearly implicit IMEX-RK method with Strang-type splitting (denoted by "S-LIMEX") and is described in Section 3.3. Both methods allow us to study the evolution of the temperature and fuel variables over time based on certain initial conditions. In the first instance, we worked with both methods in order to study the evolution of numerical simulations such as those of Example 1. With respect to the original main research question on the gain of efficiency we mention that as Table 2 indicates, the S-LIMEX method is associated with significantly smaller CPU times, especially for fine discretizations of the computational domain in compare to the fully Explicit scheme. Concerning the secondary research question whether the model is indeed suitable for the simulation of wildfire control, we comment first that by considering a certain set of parameters we observe that the fire focus tends to move through the domain according to the direction of the wind speed and that given a certain constant fuel density, this focus increases as time passes. In Example 1 both methods are compared according to the projection of the fire front onto the diagonal of the domain.
To further demonstrate the practical applicability of the second numerical method we have designed a strategy, knowing that when the availability of fuel is very low then the wildfire tends to be extinguished by the effect of diffusion and convection. This property allows us to detect the horizontal distance that the fire front travels before it is completely extinguished. One may then build a firebreak in order to mitigate the effect of the fire or totally prevent its passage towards a certain area that has to be protected and thus avoid what is seen in one of the subcases of Example 2. Then considering the ideal width of said firebreak we proceed to compare two situations, one in which the width is less than the ideal and whose effect is partial in the sense that the fire does propagate into the area one is interested in protecting. The partiality of the effect becomes apparent in that only less fuel is consumed than in the case when the firebreak is absent (as is illustrated by Example 3), but the propagation of the fire is not stopped. On the other hand, when considering the ideal width we can see that this new firebreak does prevent the fire front from passing into the area that one wishes to protect. We have also studied what happens when a slight change in the direction of the wind occurs in a certain domain, which has the effect of considerably widening the width of the focus of the fire.
Clearly, our approach still exhibits a number of limitations. For the moment, and to focus on computational issues, we have considered idealized examples and no real-world scenarios, for which data of the terrain, vegetation, climate and wind effects would be required. Another omission is the neglect of moisture and pyrolysis effects. From a computational point of view, we mention that the present version of the Strang splitting method is still based on using a single time step for the whole computational domain. It could probably be made more efficient by local time stepping so that time-consuming small time steps are employed only where it is necessary.
Overcoming the limitations addressed above will be our main goal of future research. We mention first that our interest in the wildfire propagation model (1) is partially motivated by experience with IMEX methods for related applicative nonlinear convection-diffusion problems. These methods achieve a considerable gain of efficiency compared with easier-to-implement explicit schemes, but require the solution of linear or nonlinear systems of algebraic equations. The results presented herein are basically encouraging and reconfirm that the model can in principle be used for the prediction of wildfire propagation and its mitigation through cutting firebreaks of sufficient width. However, future research is still necessary to establish IMEX methods as a really competitive tool for wildfire simulation. One aspect is related to the discretization of the reaction terms. Another future research direction points at incorporating further model ingredients. These include topography of the terrain (as expressed by the function h(x) that describes spatial heterogeneity but which has been chosen constant here; see Section 2.2), as well as radiation, moisture content and pyrolysis via a multivalued enthalpy function.
Finally, we remark that as the discussion of Sections 1.1 and 1.2 shows, a considerable number of variants of the model (4), (5), as originally proposed in [1], have appeared along with a great variety of numerical methods. All these models refer to realities that vary drastically between different geographic regions, types of vegetation and climate. Moreover forecasts are sought on spatial and temporal scales that differ substantially from case to case. Thus, a comprehensive comparison of the models and their discretizations should start with the definition of a unique "benchmark" scenario for numerical testing that could be handled by all of them. This endeavor is, however, outside the scope of this paper.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: